1 files changed, 141 insertions, 0 deletions
diff --git a/Utilib/src/MSRLUS2.f b/Utilib/src/MSRLUS2.f
new file mode 100644
index 0000000..d50c33b
--- /dev/null
+++ b/Utilib/src/MSRLUS2.f
@@ -0,0 +1,141 @@
+*DECK MSRLUS
+      SUBROUTINE MSRLUS2(LFORW,N,LC,LC0,IM,MCU,IM0,MCU0,JU,ILUDF,ILUCF,
+     1                   CF,XIN,XOUT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solve LU y = x where LU is stored in the "L\U-I" form in MSR format.
+* Can be use "in-place" i.e. XOUT=XIN.
+*
+*Copyright:
+* Copyright (C) 2002 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): R. Le Tellier
+*
+*Parameters: input
+* LFORW   flag set to .false. to transpose the coefficient matrix.
+* N       order of the system.
+* LC      dimension of MCU.
+* LC0     dimension of MCU0.
+* IM      elements of row i in MCU(IM(i):IM(i+1)-1) for CF.
+* MCU
+* IM0     elements of row i in MCU0(IM0(i):IM0(i+1)-1) for ILUCF.
+* MCU0
+* JU      MCU(JU(i):IM(i+1)) corresponds to U.
+*         MCU(IM(i)+1:JU(i)-1) correspond to L.
+* ILUDF   diagonal elements of U (inversed diagonal).
+* ILUCF   non-diagonal elements of U which differs from CF.
+* CF      non-diagonal elements of the original matrix. 
+* XIN     input vector x. 
+*    
+*Parameters: output
+* XOUT    output vector y.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*---
+* SUBROUTINE ARGUMENTS
+*---
+      INTEGER N,LC,LC0,IM(N+1),MCU(LC),IM0(N+1),MCU0(LC0),JU(N)
+      REAL ILUDF(N),ILUCF(LC0),CF(LC)
+      DOUBLE PRECISION XIN(N),XOUT(N)
+      LOGICAL LFORW
+*---
+* LOCAL VARIABLES
+*---
+      INTEGER I,J,IJ,IK
+      REAL ICFIJ
+      IF(LFORW) THEN
+*---
+* FORWARD SOLVE
+*---
+        DO I=1,N
+           XOUT(I)=XIN(I)
+           DO IJ=IM(I)+1,JU(I)-1
+              J=MCU(IJ)
+              IF ((J.GT.0).AND.(J.LE.N)) THEN
+                 DO IK=IM0(I)+1,IM0(I+1)
+                    IF (MCU0(IK).EQ.J) THEN
+                       ICFIJ=ILUCF(IK)
+                       GOTO 10
+                    ENDIF
+                 ENDDO
+                 ICFIJ=ILUDF(J)*CF(IJ)
+ 10              CONTINUE
+                 XOUT(I)=XOUT(I)-ICFIJ*XOUT(J)
+              ENDIF
+           ENDDO
+        ENDDO
+*---
+* BACKWARD SOLVE
+*---
+        DO I=N,1,-1
+           DO IJ=JU(I),IM(I+1)
+              J=MCU(IJ)
+              IF ((J.GT.0).AND.(J.LE.N)) THEN
+                 DO IK=IM0(I)+1,IM0(I+1)
+                    IF (MCU0(IK).EQ.J) THEN
+                       ICFIJ=ILUCF(IK)
+                       GOTO 20
+                    ENDIF
+                 ENDDO
+                 ICFIJ=CF(IJ)
+ 20              CONTINUE
+                 XOUT(I)=XOUT(I)-ICFIJ*XOUT(J)
+              ENDIF
+           ENDDO
+           XOUT(I)=ILUDF(I)*XOUT(I)
+        ENDDO
+      ELSE
+*---
+* FORWARD SOLVE
+*---
+        DO I=1,N
+           XOUT(I)=XIN(I)
+        ENDDO
+        DO I=1,N
+           DO IJ=JU(I),IM(I+1)
+              J=MCU(IJ)
+              IF ((J.GT.0).AND.(J.LE.N)) THEN
+                 DO IK=IM0(I)+1,IM0(I+1)
+                    IF (MCU0(IK).EQ.J) THEN
+                       ICFIJ=ILUCF(IK)
+                       GOTO 30
+                    ENDIF
+                 ENDDO
+                 ICFIJ=ILUDF(I)*CF(IJ)
+ 30              CONTINUE
+                 XOUT(J)=XOUT(J)-ICFIJ*XOUT(I)
+              ENDIF
+           ENDDO
+        ENDDO
+*---
+* BACKWARD SOLVE
+*---
+        DO I=N,1,-1
+           XOUT(I)=ILUDF(I)*XOUT(I)
+           DO IJ=IM(I)+1,JU(I)-1
+              J=MCU(IJ)
+              IF ((J.GT.0).AND.(J.LE.N)) THEN
+                 DO IK=IM0(I)+1,IM0(I+1)
+                    IF (MCU0(IK).EQ.J) THEN
+                       ICFIJ=ILUCF(IK)
+                       GOTO 40
+                    ENDIF
+                 ENDDO
+                 ICFIJ=CF(IJ)
+ 40              CONTINUE
+                 XOUT(J)=XOUT(J)-ILUDF(J)*ICFIJ*XOUT(I)
+              ENDIF
+           ENDDO
+        ENDDO
+      ENDIF
+*
+      RETURN
+      END