Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid

[unres.git] / source / unres / src_MD-M-newcorr / banach.f

diff --git a/source/unres/src_MD-M-newcorr/banach.f b/source/unres/src_MD-M-newcorr/banach.f
new file mode 100644 (file)
index 0000000..7c43d77
--- /dev/null
+++ b/source/unres/src_MD-M-newcorr/banach.f
@@ -0,0 +1,99 @@
+C
+C**********************
+      SUBROUTINE BANACH(N,NMAX,A,X,osob)
+C**********************
+C     Banachiewicz
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      DIMENSION A(NMAX,NMAX),X(NMAX),D(MAXRES6)
+      COMMON /BANII/ D
+      logical osob
+      osob=.false.
+      if (dabs(a(1,1)).lt.1.0d-15) then
+        osob=.true.
+        return
+      endif
+      D(1)=1./A(1,1)
+      DO 80 I=2,N
+      A(I,1)=A(1,I)
+      DO 81 J=2,I-1
+      XX=A(J,I)
+      DO 82 K=1,J-1
+      XX=XX-A(I,K)*A(J,K)
+   82 CONTINUE
+      A(I,J)=XX
+   81 CONTINUE
+      XX=A(I,I)
+      JJJJ=I-1
+      DO 83 J=1,JJJJ
+      AIJ=A(I,J)
+      AIJD=AIJ*D(J)
+      A(I,J)=AIJD
+      XX=XX-AIJ*AIJD
+   83 CONTINUE 
+      if (dabs(xx).lt.1.0d-15) then
+        osob=.true.
+        return
+      endif
+      D(I)=1./XX
+   80 CONTINUE
+C
+      CALL BANAII(N,NMAX,A,X)
+      RETURN
+      END
+C************************
+      SUBROUTINE BANAII(N,NMAX,A,X)
+C************************
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      DIMENSION A(NMAX,NMAX),X(NMAX),D(MAXRES6)
+      COMMON /BANII/ D
+      DO 90 I=1,N
+      Z=X(I)
+      JJJJ=I-1
+      DO 91 J=JJJJ,1,-1
+      Z=Z-A(I,J)*X(J)
+   91 CONTINUE
+      X(I)=Z
+   90 CONTINUE
+      DO 92 I=N,1,-1
+      Z=X(I)*D(I)
+      JJJJ=I+1
+      DO 93 J=JJJJ,N
+      Z=Z-A(J,I)*X(J)
+   93 CONTINUE
+      X(I)=Z
+   92 CONTINUE
+      RETURN
+      END
+C
+      SUBROUTINE MATINVERT(N,NMAX,A,A1,osob)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS' 
+      DIMENSION A(NMAX,NMAX),A1(NMAX,NMAX),D(MAXRES6)
+      COMMON /BANII/ D
+      DIMENSION X(NMAX)
+      logical osob
+      DO I=1,N
+        X(I)=0.0
+      ENDDO
+      X(1)=1.0
+      CALL BANACH(N,NMAX,A,X,osob)
+      if (osob) return
+      DO I=1,N
+        A1(I,1)=X(I)
+      ENDDO
+      DO I=2,N
+        DO J=1,N
+          X(J)=0.0
+        ENDDO
+        X(I)=1.0
+        CALL BANAII(N,NMAX,A,X)
+        DO J=1,N
+          A1(J,I)=X(J)
+        ENDDO
+      ENDDO
+      RETURN
+      END
+
+