Initial commit from Polytechnique Montreal

1 files changed, 210 insertions, 0 deletions

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