Initial commit from Polytechnique Montreal

1 files changed, 159 insertions, 0 deletions

diff --git a/Dragon/src/MCGPRA.f b/Dragon/src/MCGPRA.f
new file mode 100644
index 0000000..169f0f3
--- /dev/null
+++ b/Dragon/src/MCGPRA.f
@@ -0,0 +1,159 @@
+*DECK MCGPRA
+      SUBROUTINE MCGPRA(LFORW,OPT,PACA,FLOUT,NLONG,LC,IM,MCU,JU,DIAGM,
+     1                  CM,ILUDF,ILUCF,DIAGF,XIN,XOUT,LC0,IM0,MCU0,CF)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Multiply a matrix stored in MSR format with a vector
+* and apply preconditioner (diagonal / ILUO stored in MSR format) 
+* to this vector.
+*
+*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.
+* OPT     1 product matrix-vector only;
+*         2 preconditioner only;
+*         3 both.
+* PACA    type of preconditioner.
+* FLOUT   flag for output:
+*         .false. matrix-vector product to -preconditioned vector
+*         to XOUT;
+*         .true.  matrix-vector product to XOUT and
+*         preconditioned vector to XIN.
+* NLONG   size of the corrective system.
+* LC      dimension of profiled matrices MCU and CM.
+* IM      MSR indexes vector.
+* MCU     MSR indexes vector.
+* JU      used in ACA  acceleration for ilu0.
+* DIAGM   diagonal of matrix to multiply by.
+* CM      non-diagonal elements of matrix to multiply by.
+* ILUDF   diagonal of ilu0 matrix.
+* ILUCF   non-diagonal elements of ilu0 matrix.
+* DIAGF   vector of diagonal preconditioning.
+* LC0     used in ILU0-ACA acceleration.
+* IM0     used in ILU0-ACA acceleration.
+* MCU0    used in ILU0-ACA acceleration.
+* CF      non-diagonal elements of the matrix
+*         corresponding to the ILU0 decomposition ILUCF.
+*
+*Parameters: input/output
+* XIN     input vector to be precondition.
+*
+*Parameters: output
+* XOUT    output vector preconditioned.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+* SUBROUTINE ARGUMENTS
+*----
+      INTEGER OPT,PACA,NLONG,LC,IM(NLONG),MCU(LC),JU(NLONG),LC0,IM0(*),
+     1 MCU0(*)
+      REAL DIAGM(NLONG),CM(LC),ILUDF(NLONG),ILUCF(LC),DIAGF(NLONG),
+     1 CF(LC)
+      DOUBLE PRECISION XIN(NLONG),XOUT(NLONG)
+      LOGICAL LFORW,FLOUT
+*---
+* LOCAL VARIABLES
+*---
+      INTEGER I,J,IJ
+      DOUBLE PRECISION FF
+*
+      IF(MOD(OPT,2).EQ.1) THEN
+*---
+* MATRIX VECTOR PRODUCT
+*---
+        IF(LFORW) THEN
+*         direct calculation
+          DO I=1,NLONG
+             FF=XIN(I)*DIAGM(I)
+             DO IJ=IM(I)+1,IM(I+1)
+                J=MCU(IJ)
+                IF((J.GT.0).AND.(J.LE.NLONG)) FF=FF+CM(IJ)*XIN(J)
+             ENDDO
+             XOUT(I)=FF
+          ENDDO
+        ELSE
+*         adjoint calculation
+          DO I=1,NLONG
+             XOUT(I)=XIN(I)*DIAGM(I)
+          ENDDO
+          DO I=1,NLONG
+             DO IJ=IM(I)+1,IM(I+1)
+                J=MCU(IJ)
+                IF((J.GT.0).AND.(J.LE.NLONG)) XOUT(J)=XOUT(J)
+     1                                               +CM(IJ)*XIN(I)
+             ENDDO
+          ENDDO
+        ENDIF
+      ENDIF
+*
+      IF((OPT/2).EQ.1) THEN
+*---
+* APPLY PRECONDITIONER
+*---
+      IF(PACA.EQ.4) THEN
+*     apply ILU0 preconditioner (optimized storage = no extra-storage)
+         IF(FLOUT) THEN
+            CALL MSRLUS1(LFORW,NLONG,LC,IM,MCU,JU,ILUDF,CF,XOUT,XIN)
+         ELSE
+            CALL MSRLUS1(LFORW,NLONG,LC,IM,MCU,JU,ILUDF,CF,XOUT,XOUT)
+         ENDIF
+      ELSE IF(PACA.EQ.3) THEN
+*     apply ILU0 preconditioner (optimized storage)
+         IF(FLOUT) THEN
+            CALL MSRLUS2(LFORW,NLONG,LC,LC0,IM,MCU,IM0,MCU0,JU,ILUDF,
+     1           ILUCF,CF,XOUT,XIN)
+         ELSE
+            CALL MSRLUS2(LFORW,NLONG,LC,LC0,IM,MCU,IM0,MCU0,JU,ILUDF,
+     1           ILUCF,CF,XOUT,XOUT)
+         ENDIF
+      ELSEIF(PACA.EQ.2) THEN
+*     apply ILU0 preconditioner (complete storage)
+         IF(FLOUT) THEN
+            CALL MSRLUS(LFORW,NLONG,LC,IM,MCU,JU,ILUDF,ILUCF,XOUT,XIN)
+         ELSE
+            CALL MSRLUS(LFORW,NLONG,LC,IM,MCU,JU,ILUDF,ILUCF,XOUT,XOUT)
+         ENDIF
+      ELSEIF(PACA.EQ.1) THEN
+*     apply Diagonal preconditioner         
+         IF(FLOUT) THEN
+            DO I=1,NLONG
+               XIN(I)=XOUT(I)/DIAGF(I)
+            ENDDO
+         ELSE
+            DO I=1,NLONG
+               XOUT(I)=XOUT(I)/DIAGF(I)
+            ENDDO
+         ENDIF
+      ELSEIF(PACA.EQ.0) THEN
+*     no preconditioner
+         IF(FLOUT) THEN
+            DO I=1,NLONG
+               XIN(I)=XOUT(I)
+            ENDDO
+         ENDIF
+      ENDIF
+      ELSE
+*---
+* DO NOT APPLY PRECONDITIONER
+*---
+         IF(FLOUT) THEN
+            DO I=1,NLONG
+               XIN(I)=XOUT(I)
+            ENDDO
+         ENDIF
+      ENDIF
+*
+      RETURN
+      END