7ba6e1d69a079551fb918a9edf4cbf4d3875d738
[unres.git] / 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