C     MATXTST1 - CDC - L.P.GERALDO - AP314
C     TESTS A SYMMETRIC MATRIX FOR POSITIVE DEFINITENESS.
C     THIS PROGRAM ALSO COMPUTATES THE DETERMINANT, INERTIA AND
C     INVERSE MATRIX
C
C    
C     OPTIONS FOR MATRIX INPUT (IOP)
C
C     OPTION 1 - INPUT ERRORS AND CORRELATION MATRIX
C     OPTION 2 - INPUT COVARIANCE MATRIX
C     OPTION 3 - INPUT CROSS SECTIONS AND RELATIVE COVARIANCE MATRIX
C
C     OPTIONS FOR MATRIX INFORMATION OUTPUT (IND)
C     TOGETHER WITH POSITIVE DEFINITENESS TEST.
C
C     OPTION - 00 = ONLY INERTIA
C     OPTION - 10 = DETERMINANT AND INERTIA
C     OPTION - 01 = INERTIA AND INVERSE MATRIX
C     OPTION - 11 = DETERMINANT, INERTIA AND INVERSE MATRIX
C
C
      REAL A(50,50),C(50,50),F(50,50),E(50),X(50)
      REAL WORK(50),DET(3),T,D,RCOND,Z(1)
      INTEGER KPVT(50),INERT(3)
      DATA LDA/50/
      DATA IN/'N '/
      WRITE(*,5)
    5 FORMAT(1X,'DATA INPUT:   ')
      OPEN(3,FILE=' ',STATUS='OLD ')
      READ(3,10)N
   10 FORMAT(I2)
      READ(3,10)IOP
      READ(3,10)IND
      INP=100+IND
   15 FORMAT(6E12.6)
      GO TO (20,40,50) IOP
   20 CONTINUE
      READ(3,15) (E(I), I=1,N)
      DO 25 I=1,N
   25    READ(3,15) (C(I,J), J=1,I)
      DO 30 I=1,N
         DO 30 J=1,I
   30    C(J,I)=C(I,J)
      DO 35 I=1,N
         DO 35 J=1,N
   35    A(I,J) = E(I)*E(J)*C(I,J)
      GO TO 70
   40 CONTINUE
      DO 45 I=1,N
   45    READ(3,15) (A(I,J), J=1,I)
      DO 46 I=1,N
         DO 46 J=1,N
   46    A(I,J)=A(J,I)
      DO 47 I=1,N
         DO 47 J=1,I
   47 C(I,J)=A(I,J)/(A(I,I)*A(J,J))**0.5
      GO TO 70
   50 CONTINUE
      READ(3,15) (X(I), I=1,N)
      DO 55 I=1,N
   55    READ(3,15) (F(I,J), J=1,I)
      DO 60 I=1,N
         DO 60 J=1,N
   60 F(I,J)=F(J,I)
      DO 65 I=1,N
         DO 65 J=1,N
   65 A(I,J)=X(I)*X(J)*F(I,J)
      DO 67 I=1,N
         DO 67 J=1,I
   67 C(I,J)=A(I,J)/(A(I,I)*A(J,J))**0.5
   70 CONTINUE
C     
C     CORRELATION MATRIX TEST
C
      DO 73 I=1,N
         DO 73 J=1,I
         IF(I .EQ. J) GO TO 73
         IF(ABS(C(I,J)) .LE. 1) GO TO 73
         WRITE(*,71)
   71    FORMAT(1X,'... I ....... J ....... C(I,J) ... ')
         WRITE(*,72)I,J,C(I,J)
   72    FORMAT(1X,I5,I10,F18.6)
   73 CONTINUE
      DO 74 I=1,N
         DO 74 J=1,I
         IF(I .EQ. J) GO TO 74
         IF(ABS(C(I,J)) .GT. 1) GO TO 180
   74 CONTINUE
C
C     POSITIVE DEFINITENESS TEST
C
      CALL SSICO(A,LDA,N,KPVT,RCOND,Z)
      WRITE(*,75)RCOND
   75 FORMAT(1X,'RCOND= ',E12.6)
      T=1.0+RCOND
      IF(T .NE. 1.0) GO TO 110
C
C
C            OUTPUT DATA
C
      WRITE(*,76)
   76 FORMAT(1X,'This matrix is singular to working precision')
      STOP
  110 CALL SSIDI(A,LDA,N,KPVT,DET,INERT,WORK,INP)
      IF (INP .EQ. 101) GO TO 112
      IF (INP .EQ. 110) GO TO 130
      IF (INP .EQ. 100) GO TO 154
  112 WRITE(*,114)
  114 FORMAT(1X,'Do you want to save the inverse matrix?Y/N ')
      READ(*,115,ERR=112)IOP
  115 FORMAT(A1)
      IF (IOP .EQ. IN) GO TO 120
      WRITE(*,116)
  116 FORMAT(1X,'Please type the output file name')
      OPEN(5,FILE=' ',STATUS='NEW ')
      WRITE(5,117)
  117 FORMAT(1X,'xxxxxxxxxx INVERSE MATRIX xxxxxxxxxx ')
      DO 118 I=1,N
         DO 118 J=1,N
  118 A(J,I)=A(I,J)
      DO 119 J=1,N
  119    WRITE(5,128) (A(J,I), I=1,J)
      GO TO 130
  120 WRITE(*,122)
  122 FORMAT(1X,'xxxxxxxxxx INVERSE MATRIX xxxxxxxxxx ')
      DO 124 I=1,N
         DO 124 J=1,N
  124 A(J,I)=A(I,J)
      DO 126 J=1,N
  126    WRITE(*,128) (A(I,J),I=1,J)
  128 FORMAT(1X,50E12.6)
      IF (INP .NE. 111) GO TO 154
  130 CONTINUE
  140 DT=DET(1)*10.0**DET(2)
  149 WRITE(*,150) DT
  150 FORMAT(1X,'DET= ',E12.6)
  154 WRITE(*,155)
  155 FORMAT(1X,'xxxxx INERTIA=Number of Eigenvalues xxxxx ')
      WRITE(*,160) (INERT(L),L=1,3)
  160 FORMAT(1X,'Positive= ',I2,'   Negative= ',I2,'   Zero= 'I2)
      IF (INERT(1) .EQ. N) GO TO 162
      IF (INERT(2) .EQ. N) GO TO 164
      IF (INERT(1)+INERT(3) .EQ. N) GO TO 166
      IF (INERT(2)+INERT(3) .EQ. N) GO TO 168
      IF (INERT(1)+INERT(2)+INERT(3) .EQ. N) GO TO 170
  162 WRITE(*,163)
  163 FORMAT(1X,'This matrix is Positive Definite')
      STOP
  164 WRITE(*,165)
  165 FORMAT(1X,'This matrix is Negative Definite')
      STOP
  166 WRITE(*,167)
  167 FORMAT(1X,'This matrix is Positive ')
      STOP
  168 WRITE(*,169)
  169 FORMAT(1X,'This matrix is Negative ')
      STOP
  170 WRITE(*,171)
  171 FORMAT(1X,'This matrix is Indefinite ')
      STOP
  180 WRITE(*,185)
  185 FORMAT(1X,'This matrix is unreal because the correlation coefficie
     1nts above are higher than 1')
      STOP
      END
      SUBROUTINE SSICO(A,LDA,N,KPVT,RCOND,Z)
      INTEGER LDA,N,KPVT(1)
      REAL A(LDA,1),Z(1)
      REAL RCOND
C
C     SSICO FACTORS A REAL SYMMETRIC MATRIX BY ELIMINATION WITH
C     SYMMETRIC PIVOTING AND ESTIMATES THE CONDITION OF THE MATRIX
C
C     INTERNAL VARIABLES
C
      REAL AK,AKM1,BK,BKM1,SDOT,DENOM,EK,T
      REAL ANORM,S,SASUM,YNORM
      INTEGER I,INFO,J,JM1,K,KP,KPS,KS
C
      DO 30 J=1,N
         Z(J)=SASUM(J,A(1,J),1)
         JM1=J-1
         IF(JM1 .LT. 1) GO TO 20
         DO 10 I=1,JM1
            Z(I)=Z(I)+ABS(A(I,J))
   10    CONTINUE
   20    CONTINUE
   30 CONTINUE
      ANORM=0.0E0
      DO 40 J=1,N
         ANORM=AMAX1(ANORM,Z(J))
   40 CONTINUE
      CALL SSIFA(A,LDA,N,KPVT,INFO)
      EK=1.0E0
      DO 50 J=1,N
         Z(J)=0.0E0
   50 CONTINUE
      K=N
   60 IF(K .EQ. 0) GO TO 120
      KS=1
      IF(KPVT(K) .LT. 0) KS=2
      KP=IABS(KPVT(K))
      KPS=K+1-KS
      IF(KP .EQ. KPS) GO TO 70
      T=Z(KPS)
      Z(KPS)=Z(KP)
      Z(KP)=T
   70 CONTINUE
      IF(Z(K) .NE. 0.0E0) EK=SIGN(EK,Z(K))
      Z(K)=Z(K)+EK
      CALL SAXPY(K-KS,Z(K),A(1,K),1,Z(1),1)
      IF(KS .EQ. 1) GO TO 80
      IF(Z(K-1) .NE. 0.0E0) EK=SIGN(EK,Z(K-1))
      Z(K-1)=Z(K-1)+EK
      CALL SAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1)
   80 CONTINUE
      IF(KS .EQ. 2) GO TO 100
      IF(ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 90
      S=ABS(A(K,K))/ABS(Z(K))
      CALL SSCAL(N,S,Z,1)
      EK=S*EK
   90 CONTINUE
      IF(A(K,K) .NE. 0.0E0) Z(K)=Z(K)/A(K,K)
      IF(A(K,K) .EQ. 0.0E0) Z(K)=1.0E0
      GO TO 110
  100 CONTINUE
      AK=A(K,K)/A(K-1,K)
      AKM1=A(K-1,K-1)/A(K-1,K)
      BK=Z(K)/A(K-1,K)
      BKM1=Z(K-1)/A(K-1,K)
      DENOM=AK*AKM1-1.0E0
      Z(K)=(AKM1*BK-BKM1)/DENOM
  110 CONTINUE
      K=K-KS
      GO TO 60
  120 CONTINUE
      S=1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
      K=1
  130 IF(K .GT. N) GO TO 160
      KS=1
      IF(KPVT(K) .LT. 0) KS=2
      IF(K .EQ. 1) GO TO 150
      Z(K)=Z(K)+SDOT(K-1,A(1,K),1,Z(1),1)
      IF(KS .EQ. 2)
     *  Z(K+1)=Z(K+1)+SDOT(K-1,A(1,K+1),1,Z(1),1)
      KP=IABS(KPVT(K))
      IF(KP .EQ. K) GO TO 140
      T=Z(K)
      Z(K)=Z(KP)
      Z(KP)=T
  140 CONTINUE
  150 CONTINUE
      K=K+KS
      GO TO 130
  160 CONTINUE
      S=1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
      YNORM=1.0E0
      K=N
  170 IF(K .EQ. 0) GO TO 230
      KS=1
      IF(KPVT(K) .LT. 0) KS=2
      IF(K .EQ. KS) GO TO 190
      KP=IABS(KPVT(K))
      KPS=K+1-KS
      IF(KP .EQ. KPS) GO TO 180
      T=Z(KPS)
      Z(KPS)=Z(KP)
      Z(KP)=T
  180 CONTINUE
      CALL SAXPY(K-KS,Z(K),A(1,K),1,Z(1),1)
      IF(KS .EQ. 2) CALL SAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1)
  190 CONTINUE
      IF(KS .EQ. 2) GO TO 210
      IF(ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 200
        S=ABS(A(K,K))/ABS(Z(K))
      CALL SSCAL(N,S,Z,1)
      YNORM=S*YNORM
  200 CONTINUE
      IF(A(K,K) .NE. 0.0E0) Z(K)=Z(K)/A(K,K)
      IF(A(K,K) .EQ. 0.0E0) Z(K)=1.0E0
      GO TO 220
  210 CONTINUE
      AK=A(K,K)/A(K-1,K)
      AKM1=A(K-1,K-1)/A(K-1,K)
      BK=Z(K)/A(K-1,K)
      BKM1=Z(K-1)/A(K-1,K)
      DENOM=AK*AKM1-1.0E0
      Z(K)=(AKM1*BK-BKM1)/DENOM
      Z(K-1)=(AK*BKM1-BK)/DENOM
  220 CONTINUE
      K=K-KS
      GO TO 170
  230 CONTINUE
      S=1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
      YNORM=S*YNORM
      K=1
  240 IF(K .GT. N) GO TO 270
      KS=1
      IF(KPVT(K) .LT. 0) KS=2
      IF(K .EQ. 1) GO TO 260
      Z(K)=Z(K)+SDOT(K-1,A(1,K),1,Z(1),1)
      IF(KS .EQ. 2)
     *  Z(K+1)=Z(K+1)+SDOT(K-1,A(1,K+1),1,Z(1),1)
      KP=IABS(KPVT(K))
      IF(KP .EQ. K) GO TO 250
      T=Z(K)
      Z(K)=Z(KP)
      Z(KP)=T
  250 CONTINUE
  260 CONTINUE
      K=K+KS
      GO TO 240
  270 CONTINUE
      S=1.0E0/SASUM(N,Z,1)
      CALL SSCAL(N,S,Z,1)
      YNORM=S*YNORM
      IF(ANORM .NE. 0.0E0) RCOND=YNORM/ANORM
      IF(ANORM .EQ. 0.0E0) RCOND=0.0E0
      RETURN
      END
C
C
C
      SUBROUTINE SSIFA(A,LDA,N,KPVT,INFO)
C     SSIFA FACTORS A REAL SYMMETRIC MATRIX BY ELIMINATION
C     WITH SYMMETRIC PIVOTING
C
      INTEGER LDA,N,KPVT(1),INFO
      REAL A(LDA,1)
      REAL AK,AKM1,BK,BKM1,DENOM,MULK,MULKM1,T
      REAL ABSAKK,ALPHA,COLMAX,ROWMAX
      INTEGER IMAX,IMAXP1,J,JJ,JMAX,K,KM1,KM2,KSTEP,ISAMAX
      LOGICAL SWAP
C
      ALPHA = (1.0E0 + SQRT(17.0E0))/8.0E0
      INFO = 0
      K=N
   10 CONTINUE
      IF(K .EQ. 0) GO TO 200
      IF(K .GT. 1) GO TO 20
      KPVT(1) = 1
      IF (A(1,1) .EQ. 0.0E0) INFO =1
      GO TO 200
   20 CONTINUE
      KM1 = K - 1
      ABSAKK = ABS(A(K,K))
      IMAX = ISAMAX(K-1,A(1,K),1)
      COLMAX = ABS(A(IMAX,K))
      IF (ABSAKK .LT. ALPHA*COLMAX) GO TO 30
         KSTEP = 1
         SWAP = .FALSE.
      GO TO 90
   30 CONTINUE
      ROWMAX = 0.0E0
      IMAXP1 = IMAX +1
      DO 40 J = IMAXP1, K
         ROWMAX = AMAX1(ROWMAX,ABS(A(IMAX,J)))
   40 CONTINUE
      IF (IMAX .EQ. 1) GO TO 50
      JMAX = ISAMAX(IMAX-1,A(1,IMAX),1)
      ROWMAX = AMAX1(ROWMAX,ABS(A(JMAX,IMAX)))
   50 CONTINUE
      IF (ABS(A(IMAX,IMAX)) .LT. ALPHA*ROWMAX) GO TO 60
         KSTEP = 1
         SWAP = .TRUE.
      GO TO 80
   60 CONTINUE
      IF (ABSAKK .LT. ALPHA*COLMAX*(COLMAX/ROWMAX)) GO TO 70
         KSTEP = 1
         SWAP = .FALSE.
      GO TO 80
   70 CONTINUE
      KSTEP = 2
      SWAP = IMAX .NE. KM1
   80 CONTINUE
   90 CONTINUE
      IF (AMAX1(ABSAKK,COLMAX) .NE. 0.0E0) GO TO 100
      KPVT(K) = K
      INFO = K
      GO TO 190
  100 CONTINUE
      IF (KSTEP .EQ. 2) GO TO 140
      IF (.NOT.SWAP) GO TO 120
      CALL SSWAP(IMAX,A(1,IMAX),1,A(1,K),1)
      DO 110 JJ = IMAX, K
         J = K + IMAX - JJ
         T = A(J,K)
         A(J,K) = A(IMAX,J)
         A(IMAX,J) = T
  110 CONTINUE
  120 CONTINUE
      DO 130 JJ = 1, KM1
         J = K - JJ
         MULK = -A(J,K)/A(K,K)
         T = MULK
         CALL SAXPY(J,T,A(1,K),1,A(1,J),1)
         A(J,K) = MULK
  130 CONTINUE
      KPVT(K) = K
      IF (SWAP) KPVT(K) = IMAX
      GO TO 190
  140 CONTINUE
      IF (.NOT.SWAP) GO TO 160
      CALL SSWAP(IMAX,A(1,IMAX),1,A(1,K-1),1)
      DO 150 JJ = IMAX, KM1
         J = KM1 + IMAX - JJ
         T = A(J,K-1)
         A(J,K-1) = A(IMAX,J)
         A(IMAX,J) = T
  150 CONTINUE
      T = A(K-1,K)
      A(K-1,K) = A(IMAX,K)
      A(IMAX,K) = T
  160 CONTINUE
      KM2 = K - 2
      IF (KM2 .EQ. 0) GO TO 180
      AK = A(K,K)/A(K-1,K)
      AKM1 = A(K-1,K-1)/A(K-1,K)
      DENOM = 1.0E0 - AK*AKM1
      DO 170 JJ = 1, KM2
         J = KM1 - JJ
         BK = A(J,K)/A(K-1,K)
         BKM1 = A(J,K-1)/A(K-1,K)
         MULK = (AKM1*BK - BKM1)/DENON
         MULKM1 = (AK*BKM1 - BK)/DENOM
         T = MULK
         CALL SAXPY(J,T,A(1,K),1,A(1,J),1)
         T = MULKM1
         CALL SAXPY(J,T,A(1,K-1),1,A(1,J),1)
         A(J,K) = MULK
         A(J,K-1) = MULKM1
  170 CONTINUE
  180 CONTINUE
      KPVT(K) = 1 - K
      IF (SWAP) KPVT(K) = -IMAX
      KPVT(K-1) = KPVT(K)
  190 CONTINUE
      K = K - KSTEP
      GO TO 10
  200 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE SSIDI(A,LDA,N,KPVT,DET,INERT,WORK,JOB)
C    
C     SSIDI COMPUTES THE DETERMINANT, INERTIA AND INVERSE OF A REAL
C     SYMMETRIC MATRIX USING THE FACTORS FROM SSIFA.
C
      INTEGER LDA,N,JOB
      REAL A(LDA,1),WORK(1),DET(2)
      INTEGER KPVT(1),INERT(3)
      REAL AKKP1,SDOT,TEMP
      REAL TEN,D,T,AK,AKP1
      INTEGER J,JB,K,KM1,KS,KSTEP
      LOGICAL NOINV,NODET,NOERT
C
      NOINV = MOD(JOB,10) .EQ. 0
      NODET = MOD(JOB,100)/10 .EQ. 0
      NOERT = MOD(JOB,1000)/100 .EQ. 0
      IF (NODET .AND. NOERT) GO TO 140
      IF (NOERT) GO TO 10
      INERT(1) = 0
      INERT(2) = 0
      INERT(3) = 0
   10 CONTINUE
      IF(NODET) GO TO 20
      DET(1) = 1.0E0
      DET(2) = 0.0E0
      TEN = 10.0E0
   20 CONTINUE
      T = 0.0E0
      DO 130 K = 1, N
         D = A(K,K)
         IF (KPVT(K) .GT. 0) GO TO 50
         IF (T .NE. 0.0E0) GO TO 30
         T = ABS(A(K,K+1))
         D = (D/T)*A(K+1,K+1) - T
         GO TO 40
   30    CONTINUE
         D = T
         T = 0.0E0
   40    CONTINUE
   50    CONTINUE
         IF (NOERT) GO TO 60
         IF (D .GT. 0.0E0) INERT(1) = INERT(1) +1
         IF (D .LT. 0.0E0) INERT(2) = INERT(2) +1
         IF (D .EQ. 0.0E0) INERT(3) = INERT(3) +1
   60    CONTINUE
         IF (NODET) GO TO 120
         DET(1) = D*DET(1)
         IF (DET(1) .EQ. 0.0E0) GO TO 110
   70    IF (ABS(DET(1)) .GE. 1.0E0) GO TO 80
         DET(1) = TEN*DET(1)
         DET(2) = DET(2) - 1.0E0
         GO TO 70
   80    CONTINUE
   90    IF(ABS(DET(1)) .LT. TEN) GO TO 100
         DET(1) = DET(1)/TEN
         DET(2) = DET(2) + 1.0E0
         GO TO 90
  100    CONTINUE
  110    CONTINUE
  120    CONTINUE
  130 CONTINUE
  140 CONTINUE
      IF (NOINV) GO TO 270
      K = 1
  150 IF (K .GT. N) GO TO 260
      KM1 = K - 1
      IF (KPVT(K) .LT. 0) GO TO 180
      A(K,K) = 1.0E0/A(K,K)
      IF (KM1 .LT. 1) GO TO 170
      CALL SCOPY(KM1,A(1,K),1,WORK,1)
      DO 160 J=1, KM1
         A(J,K) = SDOT(J,A(1,J),1,WORK,1)
         CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K),1)
  160 CONTINUE
      A(K,K) = A(K,K) + SDOT(KM1,WORK,1,A(1,K),1)
  170 CONTINUE
      KSTEP = 1
      GO TO 220
  180 CONTINUE
      T = ABS(A(K,K+1))
      AK = A(K,K)/T
      AKP1 = A(K+1,K+1)/T
      AKKP1 = A(K,K+1)/T
      D = T*(AK*AKP1 - 1.0E0)
      A(K,K) = AKP1/D
      A(K+1,K+1) = AK/D
      A(K,K+1) = -AKKP1/D
      IF (KM1 .LT. 1) GO TO 210
      CALL SCOPY(KM1,A(1,K+1),1,WORK,1)
      DO 190 J =1, KM1
         A(J,K+1) = SDOT(J,A(1,J),1,WORK,1)
         CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K+1),1)
  190 CONTINUE
      A(K+1,K+1) = A(K+1,K+1) + SDOT(KM1,WORK,1,A(1,K+1),1)
      A(K,K+1) = A(K,K+1) + SDOT(KM1,A(1,K),1,A(1,K+1),1)
      CALL SCOPY(KM1,A(1,K),1,WORK,1)
      DO 200 J = 1, KM1
         A(J,K) = SDOT(J,A(1,J),1,WORK,1)
         CALL SAXPY(J-1,WORK(J),A(1,J),1,A(1,K),1)
  200 CONTINUE
      A(K,K) = A(K,K) + SDOT(KM1,WORK,1,A(1,K),1)
  210 CONTINUE
      KSTEP = 2
  220 CONTINUE
      KS = IABS(KPVT(K))
      IF (KS .EQ. K) GO TO 250
      CALL SSWAP(KS,A(1,KS),1,A(1,K),1)
      DO 230 JB = KS, K
         J = K + KS - JB
         TEMP = A(J,K)
         A(J,K) = A(KS,J)
         A(KS,J) = TEMP
  230 CONTINUE
      IF (KSTEP .EQ. 1) GO TO 240
      TEMP = A(KS,K+1)
      A(KS,K+1) = A(K,K+1)
      A(K,K+1) = TEMP
  240 CONTINUE
  250 CONTINUE
      K = K + KSTEP
      GO TO 150
  260 CONTINUE
  270 CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE SSWAP(N,SX,INCX,SY,INCY)
C
C     SSWAP INTERCHANGES TWO VECTORS
C
      REAL SX(1),SY(1),STEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
      IF(N .LE. 0 ) RETURN
      M = MOD(N,3)
      IF (M .EQ. 0) GO TO 40
      DO 30 I = 1,M
         STEMP = SX(I)
         SX(I) = SY(I)
         SY(I) = STEMP
   30 CONTINUE
      IF (N .LT. 3) RETURN
   40 MP1 = M+1
      DO 50 I = MP1,N,3
         STEMP = SX(I)
         SX(I) = SY(I)
         SY(I) = STEMP
         STEMP = SX(I + 1)
         SX(I + 1) = SY(I + 1)
         SY(I + 1) = STEMP
         STEMP = SX(I + 2)
         SX(I + 2) = SY(I + 2)
         SY(I + 2) = STEMP
   50 CONTINUE
      RETURN
      END
C
C
C
C
      SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
C
C     SAXPY COMPUTES THE OPERATION: CONSTANT TIMES A VECTOR
C     PLUS A VECTOR.
C
      REAL SX(1),SY(1),SA
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF(N.LE.0)RETURN
      IF(SA .EQ. 0.0) RETURN
      M=MOD(N,4)
      IF(M .EQ. 0) GO TO 40
      DO 30 I=1,M
         SY(I)=SY(I)+SA*SX(I)
   30 CONTINUE
      IF(N .LT. 4) RETURN
   40 MP1=M+1
      DO 50 I=MP1,N,4
         SY(I)=SY(I)+SA*SX(I)
         SY(I+1)=SY(I+1)+SA*SX(I+1)
         SY(I+2)=SY(I+2)+SA*SX(I+2)
         SY(I+3)=SY(I+3)+SA*SX(I+3)
   50 CONTINUE
      RETURN
      END
C
C
C
      REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
C
C     SDOT FORMS THE DOT PRODUCT OF TWO VECTORS
C
      REAL SX(1),SY(1),STEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      STEMP=0.0E0
      SDOT=0.0E0
      IF(N .LE. 0)RETURN
   20 M=MOD(N,5)
      IF(M .EQ. 0) GO TO 40
      DO 30 I=1,M
         STEMP=STEMP+SX(I)*SY(I)
   30 CONTINUE
      IF(N .LT. 5) GO TO 60
   40 MP1=M+1
      DO 50 I=MP1,N,5
         STEMP=STEMP+SX(I)*SY(I)+SX(I+1)*SY(I+1)+SX(I+2)*SY(I+2)+
     *    SX(I+3)*SY(I+3)+SX(I+4)*SY(I+4)
   50 CONTINUE
   60 SDOT=STEMP
      RETURN
      END
C
C
      INTEGER FUNCTION ISAMAX(N,SX,INCX)
C
C     ISAMAX FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE
C
      REAL SX(1),SMAX
      INTEGER I,INCX,IX,N
C
      ISAMAX = 0
      IF (N .LT. 1) RETURN
      ISAMAX = 1
      IF (N .EQ. 1) RETURN
      SMAX = ABS(SX(1))
      DO 30 I = 2,N
         IF(ABS(SX(I)) .LE. SMAX) GO TO 30
         ISAMAX = I
         SMAX = ABS(SX(I))
   30 CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE SCOPY(N,SX,INCX,SY,INCY)
C
C     SCOPY COPIES A VECTOR X TO A VECTOR Y
C
      REAL SX(1),SY(1)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      IF (N .LE. 0) RETURN
      M = MOD(N,7)
      IF(M .EQ. 0) GO TO 40
      DO 30 I=1,M
         SY(I) = SX(I)
   30 CONTINUE
      IF(N .LT. 7) RETURN
   40 MP1 = M+1
      DO 50 I = MP1,N,7
         SY(I)=SX(I)
         SY(I+1)=SX(I+1)
         SY(I+2)=SX(I+2)
         SY(I+3)=SX(I+3)
         SY(I+4)=SX(I+4)
         SY(I+5)=SX(I+5)
         SY(I+6)=SX(I+6)
   50 CONTINUE
      RETURN
      END
C
C
      REAL FUNCTION SASUM(N,SX,INCX)
C
C     SASUM TAKES THE SUM OF THE ABSOLUTE VALUES.
C
      REAL SX(1),STEMP
      INTEGER I,INCX,M,MP1,N,NINCX
C
      SASUM=0.0E0
      STEMP=0.0E0
      IF(N .LE. 0) RETURN
   20 M=MOD(N,6)
      IF(M .EQ. 0) GO TO 40
      DO 30 I=1,M
         STEMP=STEMP+ABS(SX(I))
   30 CONTINUE
      IF(N .LT. 6) GO TO 60
   40 MP1=M+1
      DO 50 I=MP1,N,6
         STEMP=STEMP+ABS(SX(I))+ABS(SX(I+1))+ABS(SX(I+2))+
     *   ABS(SX(I+3))+ABS(SX(I+4))+ABS(SX(I+5))
   50 CONTINUE
   60 SASUM=STEMP
      RETURN
      END
C
C
C
      SUBROUTINE SSCAL(N,SA,SX,INCX)
C
C     SSCAL SCALES A VECTOR BY A CONSTANT.
C
      REAL SA,SX(1)
      INTEGER I,INCX,M,MP1,N,NINCX
C
      IF(N .LE. 0) RETURN
   20 M=MOD(N,5)
      IF(M .EQ. 0) GO TO 40
      DO 30 I=1,M
         SX(I)=SA*SX(I)
   30 CONTINUE
      IF(N .LT. 5) RETURN
   40 MP1=M+1
      DO 50 I=MP1,N,5
         SX(I)=SA*SX(I)
         SX(I+1)=SA*SX(I+1)
         SX(I+2)=SA*SX(I+2)
         SX(I+3)=SA*SX(I+3)
         SX(I+4)=SA*SX(I+4)
   50 CONTINUE
      RETURN
      END