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/MCGABG.f | 210 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 210 insertions(+) create mode 100644 Dragon/src/MCGABG.f (limited to 'Dragon/src/MCGABG.f') diff --git a/Dragon/src/MCGABG.f b/Dragon/src/MCGABG.f new file mode 100644 index 0000000..c724d04 --- /dev/null +++ b/Dragon/src/MCGABG.f @@ -0,0 +1,210 @@ +*DECK MCGABG + SUBROUTINE MCGABG(IPRINT,LFORW,PACA,N,LC,EPSM,MAXM,IM,MCU,JU, + 1 DIAGF,CF,ILUDF,ILUCF,RHS,F,FAC,LC0,IM0,MCU0) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solve the ACA corrective system using BICGSTAB. +* +*Reference: (p382) +* MEURANT, G. 1999. "Computer Solution of Large Linear Systems". +* Studies in Mathematics and its Applications vol.28. North Holland. +* 776p. +* +*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 +* IPRINT print parameter. +* LFORW flag set to .false. to transpose the coefficient matrix. +* PACA type of preconditioner to solve the ACA corrective system. +* N dimension of the corrective system. +* LC dimension of profiled matrices MCU and CQ. +* IM connection matrix. +* MCU connection matrix. +* DIAGF diagonal elements of the matrix to inverse. +* CF non-diagonal elements of the matrix to inverse. +* RHS right hand-side of the corrective system (already +* preconditioned). +* FAC scaling factor for precision. +* LC0 used in ILU0-ACA acceleration. +* IM0 used in ILU0-ACA acceleration. +* MCU0 used in ILU0-ACA acceleration. +* EPSM stopping criterion. +* MAXM maximum number of iterations allowed. +* +*Parameters: output +* F corrective fluxes and currents. +* +*Parameters: scratch +* JU undefined. +* ILUDF undefined. +* ILUCF undefined. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER IPRINT,PACA,N,LC,IM(N+1),MCU(LC),JU(N),MAXM,LC0,IM0(*), + 1 MCU0(*) + REAL DIAGF(N),CF(LC),ILUDF(N),ILUCF(LC),EPSM,FAC + DOUBLE PRECISION RHS(N),F(N) + LOGICAL LFORW +*---- +* LOCAL VARIABLE +*---- + REAL EPSMAX,EPSINF,EPS2 + PARAMETER (EPSMAX=1E-7) + INTEGER I,J,ITER + DOUBLE PRECISION R,BI,WI,RT1 + DOUBLE PRECISION DDOT,AUX(2),EPS,FNORM,RHSN,ASIN,ASIN2,SIN,CN,SQ2 + LOGICAL DEBUG + INTRINSIC SQRT,ABS,SIGN + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: PI,RI,SI,ROT,API, + 1 ASI +*---- +* SCRATCH STORAGE ALLOCATION +*---- + ALLOCATE(PI(N),RI(N),SI(N),ROT(N),API(N),ASI(N)) +* + SQ2=1.D0/SQRT(2.D0) +*---- + DEBUG=.FALSE. + EPSINF=EPSMAX*FAC + ITER=0 +* + RHSN=0.0 + DO I=1,N + RHSN=MAX(RHSN,ABS(RHS(I))) + ENDDO + IF(RHSN.LT.EPSINF) THEN + DO I=1,N + F(I)=0.0 + ENDDO + IF(DEBUG) WRITE(6,200) RHSN,EPSINF + GO TO 40 + ENDIF + EPS2=EPSMAX*REAL(RHSN) + EPS2=EPS2*EPS2 +*--- +* initial corrective flux is set to rhs +* calculate (P times (D times RHS)) -> RI + CALL MCGPRA(LFORW,3,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGF,CF,ILUDF, + 1 ILUCF,DIAGF,RHS(1),RI,LC0,IM0,MCU0,CF) + DO I=1,N + F(I)=RHS(I) + RI(I)=RHS(I)-RI(I) + PI(I)=RI(I) + ROT(I)=RI(I) + ENDDO + R=DDOT(N,RI,1,RI,1) + FNORM=DDOT(N,F,1,F,1) + EPS=SQRT(R/FNORM) + IF(DEBUG) WRITE(6,100) ITER,EPS,EPSM + IF(EPS.LE.EPSM) THEN + IF(IPRINT.GT.2) WRITE(6,100) ITER,EPS,EPSM + GO TO 40 + ENDIF + AUX(1)=R +* + DO WHILE (ITER.LT.MAXM) +* BiCGSTAB iterations + ITER=ITER+1 +* calculate (P times (D times PI)) -> API + CALL MCGPRA(LFORW,3,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGF,CF, + 1 ILUDF,ILUCF,DIAGF,PI(1),API,LC0,IM0,MCU0,CF) +* + AUX(2)=AUX(1)/DDOT(N,API,1,ROT,1) + DO J=1,N + SI(J)=RI(J)-AUX(2)*API(J) + ENDDO + ITER=ITER+1 +* calculate (P times (D times SI)) -> ASI + CALL MCGPRA(LFORW,3,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGF,CF, + 1 ILUDF,ILUCF,DIAGF,SI(1),ASI,LC0,IM0,MCU0,CF) +* +!!!! ASIN=DDOT(N,ASI,1,ASI,1) +!!!! ASIN2=DDOT(N,ASI,1,SI,1) +!!!! IF(ASIN.GT.EPSMAX*ASIN2) THEN +!!!! WI=ASIN2/ASIN +!!!! ELSE +!!!!* assuming lucky breakdown +!!!! WI=1.0 +!!!! ENDIF +* Modification proposed by Sleijpen and Van der Vorst (Numerical Algorithms, 10:203-223, 1995) + ASIN2=DDOT(N,ASI,1,SI,1) + ASIN=SQRT(DDOT(N,ASI,1,ASI,1)) + SIN=SQRT(DDOT(N,SI,1,SI,1)) + CN=ASIN*SIN + IF(CN.GT.EPSMAX*ASIN2) THEN + CN=ASIN2/CN + WI=MAX(ABS(CN),SQ2)*SIN/ASIN + WI=SIGN(WI,CN) + ELSE +* assuming lucky breakdown + WI=1.0 + ENDIF +* calculate new iterate + DO J=1,N + F(J)=F(J)+AUX(2)*PI(J)+WI*SI(J) + RI(J)=SI(J)-WI*ASI(J) + ENDDO + R=DDOT(N,RI,1,RI,1) + FNORM=DDOT(N,F,1,F,1) + IF(FNORM.LT.EPS2) GOTO 30 + EPS=SQRT(R/FNORM) + IF(DEBUG) WRITE(6,100) ITER,EPS,EPSM + IF(EPS.LE.EPSM) GO TO 20 + RT1=AUX(1) + AUX(1)=DDOT(N,RI,1,ROT,1) + BI=AUX(1)/RT1*AUX(2)/WI + DO J=1,N + PI(J)=RI(J)+BI*(PI(J)-WI*API(J)) + ENDDO + ENDDO + 20 CONTINUE +* determine final residual norm + ITER=ITER+1 +* calculate (P times (D times F)) -> RI + CALL MCGPRA(LFORW,3,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGF,CF,ILUDF, + 1 ILUCF,DIAGF,F(1),RI,LC0,IM0,MCU0,CF) + DO I=1,N + RI(I)=RHS(I)-RI(I) + ENDDO +* + R=DDOT(N,RI,1,RI,1) + FNORM=DDOT(N,F,1,F,1) + IF(FNORM.LT.EPS2) GOTO 30 + EPS=SQRT(R/FNORM) + IF(IPRINT.GT.2) WRITE(6,100) ITER,EPS,EPSM +!!!! IF(EPS.GT.EPSM) THEN +!!!! DO I=1,N +!!!! PI(I)=RI(I) +!!!! ROT(I)=RI(I) +!!!! ENDDO +!!!! GO TO 10 +!!!! ENDIF + GO TO 40 +* + 30 IF(DEBUG) WRITE(6,300) ITER,FNORM,EPS2 + F(:N)=0.0 +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + 40 DEALLOCATE(ASI,API,ROT,SI,RI,PI) + RETURN +* + 100 FORMAT(12X,14H MCGABG: ITER=,I3,5H EPS=,E9.2,5H TAR=,E9.2) + 200 FORMAT(12X,27H MCGABG: RHS INFINITE NORM=,E9.2,5H LIM=,E9.2/ + 1 12X,33H -> ACA CORRECTION IS SET TO ZERO) + 300 FORMAT(12X,14H MCGABG: ITER=,I3,7H FNORM=,E9.2,5H LIM=,E9.2) + END -- cgit v1.2.3