       SUBROUTINE DKMVHF(A,N,V,IFLAG)                                   
C                                             KARIN MEYER               
C                                             NOVEMBER 1983             
C-----------------------------------------------------------------------
      DOUBLE PRECISION A(1),V(1),XX,DMAX,AMAX,ZERO,DIMAX                
      INTEGER IFLAG(1)                                                  
      zero=1.D-12
      IF(N.EQ.1)THEN                                                    
      XX=A(1)                                                           
      IF(DABS(XX).GT.ZERO)THEN                                          
      A(1)=1.D0/XX                                                      
      ELSE                                                              
      A(1)=0.D0                                                         
      END IF                                                            
      RETURN                                                            
      END IF                                                            
                                                                        
      N1=N+1                                                            
      NN=(N*N1)/2                                                         
      DO 1 I=1,N                                                        
    1 IFLAG(I)=0                                                        
C                                                                       
C      SET MINIMUM ABSOLUTE VALUE OF DIAGONAL ELEMENTS FOR              
C      NON-SINGULARITY (MACHINE SPECIFIC!)                              
       ZERO=1.D-12                                                      
C-----------------------------------------------------------------------
C      START LOOP OVER ROWS/COLS                                        
C-----------------------------------------------------------------------
      DO 8 II=1,N                                                       
C      ... FIND DIAGONAL ELEMENT WITH BIGGEST ABSOLUTE VALUE            
          DMAX=0.D0                                                     
          AMAX=0.D0                                                     
          KK=-N                                                         
          DO 2 I=1,N                                                    
C       ... CHECK THAT THIS ROW/COL HAS NOT BEEN PROCESSED              
              KK=KK+N1-I                                                
              IF(IFLAG(I).NE.0)GO TO 2                                  
              K=KK+I                                                    
              IF(DABS(A(K)).GT.AMAX)THEN                                
                  DMAX=A(K)                                             
                  AMAX=DABS(DMAX)                                       
                  IMAX=I                                                
              END IF                                                    
    2     CONTINUE                                                      
C      ... CHECK FOR SINGULARITY                                        
          IF(AMAX.LE.ZERO)GO TO 11                                      
C      ... ALL ELEMENTS SCANNED,SET FLAG                                
          IFLAG(IMAX)=II                                                
                                                                        
C      ... INVERT DIAGONAL                                              
          DIMAX=1.D0/DMAX                                               
C      ... DEVIDE ELEMENTS IN ROW/COL PERTAINING TO THE BIGGEST         
C      DIAGONAL ELEMENT BY DMAX                                         
          IL=IMAX-N                                                     
          DO 3 I=1,IMAX-1                                               
              IL=IL+N1-I                                                
              XX=A(IL)                                                  
              A(IL)=XX*DIMAX                                            
          IF(DABS(XX).LT.0.1D-17)XX=0.D0                                
    3     V(I)=XX                                                       
C      ... NEW DIAGONAL ELEMENT                                         
          IL=IL+N1-IMAX                                                 
          A(IL)=-DIMAX                                                  
          DO 4 I=IMAX+1,N                                               
              IL=IL+1                                                   
              XX=A(IL)                                                  
              A(IL)=XX*DIMAX                                            
          IF(DABS(XX).LT.0.1D-17)XX=0.D0                                
    4     V(I)=XX                                                       
C      ... ADJUST THE OTHER ROWS/COLS :                                 
C      A(I,J)=A(I,J)-A(I,IMAX)*A(J,IMAX)/A(IMAX,IMAX)                   
          IJ=0                                                          
          DO 7 I=1,N                                                    
              IF(I.EQ.IMAX)THEN                                         
      IJ=IJ+N1-I                                                        
      GO TO 7                                                           
      END IF                                                            
              XX=V(I)                                                   
              IF(XX.NE.0.D0)THEN                                        
                  XX=XX*DIMAX                                           
                  DO 5 J=I,N                                            
                      IJ=IJ+1                                           
                      IF(J.NE.IMAX)A(IJ)=A(IJ)-XX*V(J)                  
    5             CONTINUE                                              
              ELSE                                                      
    6             IJ=IJ+N1-I                                            
              END IF                                                    
    7     CONTINUE                                                      
C      ... REPEAT UNTIL ALL ROWS/COLS ARE PROCESSED                     
    8 CONTINUE                                                          
C-----------------------------------------------------------------------
C      END LOOP OVER ROWS/COLS                                          
C-----------------------------------------------------------------------
C      ... REVERSE SIGN                                                 
      DO 9 I=1,NN                                                       
    9 A(I)=-A(I)                                                        
C      ... AND THAT'S IT !                                              
C     PRINT 10,N                                                        
   10 FORMAT(' FULL RANK MATRIX INVERTED, ORDER =',I5)                  
C      RETURN RANK AS LAST ELEMENT OF FLAG VECTOR                       
       IFLAG(N)=N                                                       
      RETURN                                                            
C-----------------------------------------------------------------------
C      MATRIX NOT OF FULL RANK, RETURN GENERALISED INVERSE              
C-----------------------------------------------------------------------
   11 IRANK=II-1                                                        
      IJ=0                                                              
      DO 14 I=1,N                                                       
          IF(IFLAG(I).EQ.0)THEN                                         
C      ... SET REMAINING N-II ROWS/COLS TO ZERO                         
              DO 12 J=I,N                                               
                  IJ=IJ+1                                               
                  A(IJ)=0.D0                                            
   12         CONTINUE                                                  
          ELSE                                                          
              DO 13 J=I,N                                               
                  IJ=IJ+1                                               
                  IF(IFLAG(J).NE.0)THEN                                 
C      ... REVERSE SIGN FOR II-1 ROWS/COLS PREVIOUSLY PROCESSED         
                      A(IJ)=-A(IJ)                                      
                  ELSE                                                  
                      A(IJ)=0.D0                                        
                  END IF                                                
   13         CONTINUE                                                  
          END IF                                                        
   14 CONTINUE                                                          
C     PRINT 15,N,IRANK                                                  
C  15 FORMAT(' GENERALISED INVERSE OF MATRIX WITH ORDER =',I5,          
C    1  '   AND RANK =',I5)                                             
      IFLAG(N)=IRANK                                                    
      RETURN                                                            
      END                                                               
      FUNCTION IHMSSF(I,J,N)
      IF(I-J)1,1,2
 1    IHMSSF=((N+N-I)*(I-1))/2+J
      RETURN
 2    IHMSSF=((N+N-J)*(J-1))/2+I
      RETURN
      END

