X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fgauss.f;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fgauss.f;h=7ba6e1d69a079551fb918a9edf4cbf4d3875d738;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git

diff --git a/source/unres/src_MD-M-newcorr/gauss.f b/source/unres/src_MD-M-newcorr/gauss.f
new file mode 100644
index 0000000..7ba6e1d
--- /dev/null
+++ b/source/unres/src_MD-M-newcorr/gauss.f
@@ -0,0 +1,69 @@
+      subroutine gauss(RO,AP,MT,M,N,*)
+c
+c CALCULATES (RO**(-1))*AP BY GAUSS ELIMINATION
+c RO IS A SQUARE MATRIX
+c THE CALCULATED PRODUCT IS STORED IN AP
+c ABNORMAL EXIT IF RO IS SINGULAR
+c       
+      integer MT, M, N, M1,I,J,IM,
+     & I1,MI,MI1    
+      double precision RO(MT,M),AP(MT,N),X,RM,PR,
+     &  Y  
+      if(M.ne.1)goto 10
+      X=RO(1,1)
+      if(dabs(X).le.1.0D-13) return 1
+      X=1.0/X
+      do 16 I=1,N
+16     AP(1,I)=AP(1,I)*X
+       return
+10     continue
+        M1=M-1
+        DO1 I=1,M1
+        IM=I
+        RM=DABS(RO(I,I))
+        I1=I+1
+        do 2 J=I1,M
+        if(DABS(RO(J,I)).LE.RM) goto 2
+        RM=DABS(RO(J,I))
+        IM=J
+2       continue
+        If(IM.eq.I)goto 17
+        do 3 J=1,N
+        PR=AP(I,J)
+        AP(I,J)=AP(IM,J)
+3       AP(IM,J)=PR
+        do 4 J=I,M
+        PR=RO(I,J)
+        RO(I,J)=RO(IM,J)
+4       RO(IM,J)=PR
+17      X=RO(I,I)
+        if(dabs(X).le.1.0E-13) return 1
+        X=1.0/X
+        do 5 J=1,N
+5       AP(I,J)=X*AP(I,J)
+        do 6 J=I1,M
+6       RO(I,J)=X*RO(I,J)
+        do 7 J=I1,M
+        Y=RO(J,I)
+        do 8 K=1,N
+8       AP(J,K)=AP(J,K)-Y*AP(I,K)
+        do 9 K=I1,M
+9       RO(J,K)=RO(J,K)-Y*RO(I,K)
+7       continue
+1       continue
+        X=RO(M,M)
+        if(dabs(X).le.1.0E-13) return 1
+        X=1.0/X
+        do 11 J=1,N
+11      AP(M,J)=X*AP(M,J)
+        do 12 I=1,M1
+        MI=M-I
+        MI1=MI+1
+        do 14 J=1,N
+        X=AP(MI,J)
+        do 15 K=MI1,M
+15      X=X-AP(K,J)*RO(MI,K)
+14      AP(MI,J)=X
+12      continue
+        return
+        end