From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Dragon/src/MCGPRA.f | 159 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) create mode 100644 Dragon/src/MCGPRA.f (limited to 'Dragon/src/MCGPRA.f') 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 -- cgit v1.2.3