source/unres/src_MD/src/gauss.f

   1       subroutine gauss(RO,AP,MT,M,N,*)
   2 c
   3 c CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
   4 c RO IS A SQUARE MATRIX
   5 c THE CALCULATED PRODUCT IS STORED IN AP
   6 c ABNORMAL EXIT IF RO IS SINGULAR
   7 c
   8       integer MT, M, N, M1,I,J,IM,
   9      & I1,MI,MI1
  10       double precision RO(MT,M),AP(MT,N),X,RM,PR,
  11      &  Y
  12       if(M.ne.1)goto 10
  13       X=RO(1,1)
  14       if(dabs(X).le.1.0D-13) return 1
  15       X=1.0/X
  16       do 16 I=1,N
  17 16     AP(1,I)=AP(1,I)*X
  18        return
  19 10     continue
  20         M1=M-1
  21         DO1 I=1,M1
  22         IM=I
  23         RM=DABS(RO(I,I))
  24         I1=I+1
  25         do 2 J=I1,M
  26         if(DABS(RO(J,I)).LE.RM) goto 2
  27         RM=DABS(RO(J,I))
  28         IM=J
  29 2       continue
  30         If(IM.eq.I)goto 17
  31         do 3 J=1,N
  32         PR=AP(I,J)
  33         AP(I,J)=AP(IM,J)
  34 3       AP(IM,J)=PR
  35         do 4 J=I,M
  36         PR=RO(I,J)
  37         RO(I,J)=RO(IM,J)
  38 4       RO(IM,J)=PR
  39 17      X=RO(I,I)
  40         if(dabs(X).le.1.0E-13) return 1
  41         X=1.0/X
  42         do 5 J=1,N
  43 5       AP(I,J)=X*AP(I,J)
  44         do 6 J=I1,M
  45 6       RO(I,J)=X*RO(I,J)
  46         do 7 J=I1,M
  47         Y=RO(J,I)
  48         do 8 K=1,N
  49 8       AP(J,K)=AP(J,K)-Y*AP(I,K)
  50         do 9 K=I1,M
  51 9       RO(J,K)=RO(J,K)-Y*RO(I,K)
  52 7       continue
  53 1       continue
  54         X=RO(M,M)
  55         if(dabs(X).le.1.0E-13) return 1
  56         X=1.0/X
  57         do 11 J=1,N
  58 11      AP(M,J)=X*AP(M,J)
  59         do 12 I=1,M1
  60         MI=M-I
  61         MI1=MI+1
  62         do 14 J=1,N
  63         X=AP(MI,J)
  64         do 15 K=MI1,M
  65 15      X=X-AP(K,J)*RO(MI,K)
  66 14      AP(MI,J)=X
  67 12      continue
  68         return
  69         end