1 files changed, 98 insertions, 0 deletions
diff --git a/Utilib/src/MSRLUS.f b/Utilib/src/MSRLUS.f
new file mode 100644
index 0000000..3b25713
--- /dev/null
+++ b/Utilib/src/MSRLUS.f
@@ -0,0 +1,98 @@
+*DECK MSRLUS
+      SUBROUTINE MSRLUS(LFORW,N,LC,IM,MCU,JU,ILUDF,ILUCF,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 IM.
+* IM      elements of row i in MCU(IM(i):IM(i+1)-1)
+* MCU
+* 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 L and U.
+* XIN     input vector x. 
+*    
+*Parameters: output
+* XOUT    output vector y.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*---
+* SUBROUTINE ARGUMENTS
+*---
+      INTEGER N,LC,IM(N+1),MCU(LC),JU(N)
+      REAL ILUDF(N),ILUCF(LC)
+      DOUBLE PRECISION XIN(N),XOUT(N)
+      LOGICAL LFORW
+*---
+* LOCAL VARIABLES
+*---
+      INTEGER I,J,IJ
+      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)) XOUT(I)=XOUT(I)
+     1        -ILUCF(IJ)*XOUT(J)
+           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)) XOUT(I)=XOUT(I)
+     1        -ILUCF(IJ)*XOUT(J)
+           ENDDO
+           XOUT(I)=ILUDF(I)*XOUT(I)
+        ENDDO
+      ELSE
+*---
+* FORWARD SOLVE (TRANSPOSED LINEAR SYSTEM)
+*---
+        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)) XOUT(J)=XOUT(J)
+     1        -ILUDF(I)*ILUCF(IJ)*XOUT(I)
+           ENDDO
+        ENDDO
+*---
+* BACKWARD SOLVE (TRANSPOSED LINEAR SYSTEM)
+*---
+        DO I=N,1,-1
+           DO IJ=IM(I)+1,JU(I)-1
+              J=MCU(IJ)
+              IF ((J.GT.0).AND.(J.LE.N)) XOUT(J)=XOUT(J)
+     1        -ILUCF(IJ)*XOUT(I)
+           ENDDO
+        ENDDO
+      ENDIF
+*
+      RETURN
+      END