      SUBROUTINE SMSINV(V,B,N,MFLAG,NRANK)
      IMPLICIT NONE
*     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,MFLAG,LK,NRANK 
      REAL V(*),B(*),VKK,VJK,PIVOT,DIAG(NDIMLD)
      LOGICAL DONE(NDIMLD)
*     ...
      IF(N.GT.NDIMLD) STOP ' SMSINV: argument N too large '
      II=0
      DO I=1,N              ! reset flags
       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
