
      SUBROUTINE SMSINV(V,B,N,MFLAG,NRANK)
*     Inverse symmetric matrix and solution of system of linear  equations
*     V(*) = symmetric matrix of N cols and rows, replaced by inverse
*     B(N) = right hand side, for MFLAG .ne. 0 replaced by solution
*     B is a dummy argument for MFLAG .eq. 0 (inversion only)
*     NRANK is returned rank of matrix; only part of matrix is inverted 
*     for NRANK < N (remaining elements set to 0)    
*     Limitation to N <= NDIMLD!
      INTEGER NDIMLD
      REAL EPS
      PARAMETER (NDIMLD=100,EPS=1.E-6)
      INTEGER I,J,K,L,M,N,II,JJ,KK,JK,LK,MFLAF,NRANK 
      REAL V(*),B(*),VKK,VJK,PIVOT,DIAG(NDIMLD)
      LOGICAL DONE(NDIMLD)
      IF(N.GT.NDIMLD) STOP ' SMSINV: argument N too large '
*     reset flags
      II=0
      DO I=1,N              
       DONE(I)=.FALSE.
       II=II+I
       DIAG(I)=ABS(V(II))
      END DO 
*     loop 
      NRANK=0
      DO I=1,N               
*      search for pivot and test for linearity and zero matrix
       K  =0
       VKK=0.0
       JJ=0 
       DO J=1,N              
        JJ=JJ+J
        IF(.NOT.DONE(J)) THEN 
           IF(K.EQ.0) K=J
           IF(ABS(V(JJ)).GE.MAX(VKK,EPS*DIAG(J))) THEN 
              VKK=ABS(V(JJ))
              K  =J
           END IF 
        END IF     
       END DO                
       DONE(K)=.TRUE.
*      KK is previous diagonal element
       KK=(K*K-K)/2 
       IF(VKK.NE.0.0) THEN
*         preparation for exchange step 
          NRANK=NRANK+1
          PIVOT=1.0/V(KK+K)
          V(KK+K)=-PIVOT
          IF(MFLAG.NE.0) B(K)=PIVOT*B(K)
          JJ=0
*         exchange step
          JK=KK
          DO J=1,N          
           IF(J.LE.K) JK=JK+1
           IF(J.GT.K) JK=JK+J-1
           IF(J.NE.K) THEN
              VJK  =V(JK)
              V(JK)=PIVOT*VJK
              IF(MFLAG.NE.0) B(J) =B(J)-B(K)*VJK
*             index of pivot row/col 
              LK=KK 
              DO L=1,J      
               IF(L.LE.K) LK=LK+1
               IF(L.GT.K) LK=LK+L-1
               IF(L.NE.K) THEN
                  V(JJ+L)=V(JJ+L)-V(LK)*VJK
               END IF
              END DO
           END IF
           JJ=JJ+J 
          END DO  
       ELSE
*         no pivot found, clear 
          JK=KK
          DO J=1,N          
           IF(J.LE.K) JK=JK+1
           IF(J.GT.K) JK=JK+J-1
           V(JK)=0.0
          END DO 
       END IF 
      END DO 
*     change sign
      DO I=1,(N*N+N)/2
       V(I)=-V(I)
      END DO 
      END

      SUBROUTINE GNLRD(A,N,NDIM1,IDX)
*     LR decomposition of general N*N quadratic matrix in A(NDIM1,N)
*     with result overwriting matrix A. 
*     IDX(N) is an array, returned with information about the order,   
*     with IDX(1)=0 in case of singular matrix.
*     Limitation to N <= NDIMP! 
      INTEGER N,IDX(N),NDIM1,NDIMP,I,J,K,IMAX 
      PARAMETER (NDIMP=100) 
      REAL A(NDIM1,N),SCALE(NDIMP),ELMAX,TMP,SUM
C     IF(N.GT.NDIMLD) STOP ' GNLRD: argument N too large ' ! changed!
      IF(N.GT.NDIMP) STOP ' GNLRD: argument N too large '
      DO I=1,N
       ELMAX=0.0
       DO J=1,N 
C       IF(ABS(A(I,J)).GT.ELMAX) ELMAX=ABS(A(I,I))  ! changed!
        IF(ABS(A(I,J)).GT.ELMAX) ELMAX=ABS(A(I,J))
       END DO 
       IF(ELMAX.EQ.0.0) THEN
          IDX(1)=0 ! matrix singular
          RETURN   
       END IF 
       SCALE(I)=1.0/ELMAX 
      END DO 
      DO J=1,N
       DO I=1,J-1
        SUM=A(I,J)
        DO K=1,I-1
         SUM=SUM-A(I,K)*A(K,J)
        END DO 
        A(I,J)=SUM 
       END DO 
       ELMAX=0.0
       DO I=J,N
        SUM=A(I,J)
        DO K=1,J-1
         SUM=SUM-A(I,K)*A(K,J)
        END DO  
        A(I,J)=SUM
        TMP=SCALE(I)*ABS(SUM)
        IF(TMP.GE.ELMAX) THEN 
           IMAX=I
           ELMAX=TMP 
        END IF 
       END DO  
       IF(J.NE.IMAX) THEN 
          DO K=1,N
           TMP=A(IMAX,K)
           A(IMAX,K)=A(J,K)
           A(J,K)=TMP 
          END DO  
          SCALE(IMAX)=SCALE(J)  
       END IF
       IDX(J)=IMAX
       IF(J.NE.N) THEN 
          TMP=1.0/A(J,J)
          DO I=J+1,N
           A(I,J)=A(I,J)*TMP 
          END DO  
       END IF 
      END DO  
      write(*,*) 'before return' 
      END

      SUBROUTINE GNSOL(A,N,NDIM1,IDX,B)
*     Solution of matrix equation using LR decomposition, with   
*     arrays A and IDX from previous call to GNLRD.
*     B(N) is array of right-hand side, replaced by result.
      REAL A(NDIM1,N),B(N),SUM
      INTEGER IDX(N),NDIM1,N,I,J,K 
      DO I=1,N
       K=IDX(I)
       SUM=B(K)
       B(K)=B(I)
       DO J=1,I-1
        SUM=SUM-A(I,J)*B(J)
       END DO   
       B(I)=SUM
      END DO 
      DO I=N,1,-1
       SUM=B(I)
       DO J=I+1,N
        SUM=SUM-A(I,J)*B(J)
       END DO   
       B(I)=SUM/A(I,I)
      END DO  
      END

      SUBROUTINE GNINV(A,AINV,N,NDIM1,IDX) 
*     Inverse of matrix AINV returned, using LR decomposition with   
*     arrays A and IDX from previous call to GNLRD.
      REAL A(NDIM1,N),AINV(NDIM1,N)
      INTEGER IDX(N),I,J,N,NDIM1
      DO J=1,N
       DO I=1,N
        AINV(I,J)=0.0 
       END DO 
       AINV(J,J)=1.0
      END DO
      DO J=1,N
       CALL GNSOL(A,N,NDIM1,IDX,AINV(1,J))
      END DO
      END

      SUBROUTINE BMSOL(W,B,NBAND,NROW)      
*     Solution of band matrix equation W * X = B
*     with symmetric band matrix in special storage mode in array W.
*     Array B(NROW) is right-hand side, replaced by result.
      INTEGER I,J,N,NN,IMAX,JMAX,NDIM,NBAND,NROW
      PARAMETER (NDIM=100)
      REAL W(NBAND,NROW),B(NROW),DIAG(NDIM),RATIO
      IF(NROW.GT.NDIM) STOP ' BMSOL: argument NROW too large ' ! added
      DO N=1,NROW
       DIAG(N)=W(1,N)
      END DO
*     decomposition
      DO N=1,NROW
       IF(W(1,N)+DIAG(N).LE.DIAG(N)) THEN
*         singular case
          DO J=1,NBAND
           W(J,N)=0.0
          END DO
       ELSE
          W(1,N)=1.0/W(1,N)
          IMAX=MIN(NBAND-1,NROW-N)
          JMAX=IMAX
          DO I=1,IMAX
           RATIO=W(I+1,N)*W(1,N)
           DO J=1,JMAX
            W(J,N+I)=W(J,N+I)-W(J+I,N)*RATIO
           END DO
           JMAX=JMAX-1
           W(I+1,N)=RATIO
          END DO
       END IF
      END DO
*     solution by forward ...
      DO N=1,NROW
       DO J=1,MIN(NBAND-1,NROW-N)
        B(J+N)=B(J+N)-W(J+1,N)*B(N)
       END DO
      END DO
*     ... and backward substitution
      DO NN=1,NROW
       N=NROW+1-NN
       B(N)=B(N)*W(1,N)
       DO J=1,MIN(NBAND-1,NROW-N)
        B(N)=B(N)-W(J+1,N)*B(J+N)
       END DO
      END DO
      END

      SUBROUTINE CHLRD(W,N)
*     Cholesky decomposition of symmetric N*N matrix W, 
*     overwritten by result of decomposition.
*     W(1)=0 returned in case of non-positive diagonal element  
      REAL W(*)
      II=0
      DO I=1,N
       II=II+I
       IF(W(II).LE.0.0) THEN
*         non-positive diagonal element          
          W(1)=0.0
          RETURN           
       END IF
       W(II)=SQRT(W(II))      ! diagonal element of R 
       DO J=I+1,N 
        IJ=I+(J*J-J)/2
        W(IJ)=W(IJ)/W(II)     ! other elements of R
       END DO  
       DO J=I+1,N
        IJ=I+(J*J-J)/2
        DO K=J,N   
         JK=J+(K*K-K)/2
         IK=I+(K*K-K)/2
         W(JK)=W(JK)-W(IJ)*W(IK)
        END DO        
       END DO
      END DO
      END

      SUBROUTINE CHSOL(W,N,B)
*     Solution of matrix equation using LR decomposition, with   
*     arrays W from previous call to CHLRD.
*     B(N) is array of right-hand side, replaced by result.
      REAL W(*),B(N)
      DO I=1,N      ! forward solution
       SUM=B(I)
       DO J=1,I-1  
        SUM=SUM+W(J+(I*I-I)/2)*B(J)
       END DO        
       B(I)=-SUM/W((I*I+I)/2)
      END DO
      DO I=1,N,-1   ! backward solution
       SUM=B(I)
       DO J=I+1,N
        SUM=SUM-W(I+(J*J-J)/2)*B(J)
       END DO
       B(I)=SUM/W((I*I+I)/2)
      END DO        
      END  

      SUBROUTINE SMJROT(V,A,N,NDIM1,DIAG,NROT)
*     Eigenvalues and eigenvectors of symmetric N*N matrix V(*), which
*     is destroyed during computation.
*     Returned: the N*N matrix of eigenvectors in array A(NDIM1,N)
*     and the diagonal elements (eigenvalues) in array DIAG(N).  
*     NROT is the performed number of rotations 
*     Limitation to N <= NDD! 
      INTEGER NDD
      PARAMETER (NDD=100)
      REAL V(*),A(NDIM1,N),DIAG(N),DDAG(NDD) 
      REAL AIP,AIQ,VIP,VIQ,AM,SM,DIFF,EPSFAC,TEST
      REAL COSPHI,SINPHI,TANPHI,TANPH2,THETA
      INTEGER I,J,N,II,IJ,P,Q,PQ,NDIM1,NROT,MSY,MDI,ITER
      PARAMETER (EPSFAC=100.0)
*     indices in symmetric matrix
      MSY(I,J)=MIN(I,J)+(MAX(I,J)*MAX(I,J)-MAX(I,J))/2
      MDI(I)  =(I*I+I)/2
      NROT=0 
      DO J=1,N   
       DO I=1,N
        A(I,J)=0.0
       END DO
       A(J,J)=1.0 ! A = unit matrix
      END DO
      II=0
      DO I=1,N
       II=II+I
       DIAG(I)=V(II) ! copy diagonal elements to DIAG
       IF(I.LE.NDD) DDAG(I)=0.0
      END DO 
      DO ITER=-N,50   ! iterations
       SM=0.0
       AM=0.0
       IJ=0
       DO I=1,N
        DO J=1,I-1
         IJ=IJ+1
         SM=SM+ABS(V(IJ))                   ! sum 
         IF(ABS(V(IJ)).GT.AM) AM=ABS(V(IJ)) ! maximum
        END DO 
        IJ=IJ+1 
       END DO  
       IF(SM.EQ.0.0) RETURN  ! test convergence
       PQ=0
       DO Q=1,N
        DO P=1,Q-1       
         PQ=PQ+1
         IF(ITER.GE.0.OR.ABS(V(PQ)).EQ.AM) THEN
            TEST=EPSFAC*ABS(V(PQ))
            IF(ITER.LT.0.OR.
     +         ABS(DIAG(P)).NE.ABS(DIAG(P))+TEST.OR.
     +         ABS(DIAG(Q)).NE.ABS(DIAG(Q))+TEST) THEN
               NROT=NROT+1       ! begin/Jacobi rotation  
               DIFF=DIAG(Q)-DIAG(P) 
               IF(ABS(DIFF).EQ.ABS(DIFF)+TEST) THEN
                  TANPHI=V(PQ)/DIFF
               ELSE 
                  THETA =0.5*DIFF/V(PQ)
                  TANPHI=1.0/(THETA+SIGN(SQRT(1.0+THETA*THETA),THETA))
               END IF
               COSPHI=1.0/SQRT(1.0+TANPHI*TANPHI)
               SINPHI=COSPHI*TANPHI
               TANPH2=SINPHI/(1.0+COSPHI)
               DO I=1,N
                IF(I.NE.Q.AND.I.NE.P) THEN
                   VIP=V(MSY(I,P)) ! transform symmetric
                   VIQ=V(MSY(I,Q)) !  matrix
                   V(MSY(I,P))=VIP-SINPHI*(VIQ+TANPH2*VIP)
                   V(MSY(I,Q))=VIQ+SINPHI*(VIP-TANPH2*VIQ)
                END IF
                AIP=A(I,P) ! update rotation
                AIQ=A(I,Q) !  matrix
                A(I,P)=AIP-SINPHI*(AIQ+TANPH2*AIP)
                A(I,Q)=AIQ+SINPHI*(AIP-TANPH2*AIQ)
               END DO 
               DIAG(P)=DIAG(P)-TANPHI*V(PQ) ! update special elements
               DIAG(Q)=DIAG(Q)+TANPHI*V(PQ) !  of symmetric matrix   
               IF(P.LE.NDD) DDAG(P)=DDAG(P)-TANPHI*V(PQ) 
               IF(Q.LE.NDD) DDAG(Q)=DDAG(Q)+TANPHI*V(PQ) 
            END IF
            V(PQ)=0.0  ! off-diagonal element becomes zero  
         END IF        ! end/Jacobi rotation  
        END DO 
        PQ=PQ+1 
       END DO
       II=0       
       DO I=1,MIN(N,NDD)
        II=II+I
        V(II)=V(II)+DDAG(I) ! improve diagonal elements  
        DIAG(I)=V(II)
        DDAG(I)=0.0
       END DO
      END DO  
      PAUSE 'SMJROT: no convergence' 
      END







