Initial commit from Polytechnique Montreal

1 files changed, 149 insertions, 0 deletions

diff --git a/Dragon/src/MCGACA.f b/Dragon/src/MCGACA.f
new file mode 100644
index 0000000..a69b169
--- /dev/null
+++ b/Dragon/src/MCGACA.f
@@ -0,0 +1,149 @@
+*DECK MCGACA
+      SUBROUTINE MCGACA(LFORW,PACA,N,NG,NFIRST,NGEFF,M,LC,NGIND,NGINDV,
+     1                  NCONV,KPSYS,JPMACR,NZON,IPERM,IM,MCU,JU,XIN,
+     2                  LC0,IM0,MCU0,XOUT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the product of the left-hand side ACA matrix in its multigroup
+* form with a vector and apply group per group left preconditioner.
+*
+*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.
+* PACA    type of preconditioner to solve the ACA corrective system.
+* N       number of unknowns per group.
+* NG      total number of groups.
+* NFIRST  first group to proceed.
+* NGEFF   number of  unconverged groups.
+* M       number of material mixtures.
+* LC      dimension of profiled matrices MCU and CQ.
+* NGIND   index of the groups to process.
+* NGINDV  index to pass from "NGEFF format" to "NG format".
+* NCONV   logical array of convergence status for each group (.TRUE.
+*         not converged).
+* KPSYS   pointer array for each group properties.
+* JPMACR  pointer to the macrolib LCM object ('GROUP' directory).
+* NZON    index-number of the mixture type assigned to each volume.
+* IPERM   permutation array for ACA.
+* IM      connection matrix.
+* MCU     connection matrix.
+* JU      used for ilu0 preconditioner.
+* XIN     undefined.
+* LC0     used in ILU0-ACA acceleration.
+* IM0     used in ILU0-ACA acceleration.
+* MCU0    used in ILU0-ACA acceleration.
+*
+*Parameters: output
+* XOUT    product.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) KPSYS(NGEFF),JPMACR
+      INTEGER PACA,N,NFIRST,NGEFF,NG,M,LC,NGIND(NGEFF),NGINDV(NG),
+     1 NZON(N),IPERM(N),IM(N+1),MCU(LC),JU(N),LC0,IM0(*),MCU0(*)
+      DOUBLE PRECISION XIN(N,NGEFF),XOUT(N,NGEFF)
+      LOGICAL LFORW,NCONV(NGEFF)
+*----
+* LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPSYS,KPMACR
+      INTEGER I,J,II,IG,JG,JJ,JND,IBM
+      REAL, TARGET, SAVE, DIMENSION(1) :: DUMMY
+*----
+* ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: NJJ,IJJ,IPOS
+      REAL, ALLOCATABLE, DIMENSION(:) :: XSCAT
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: TEMP
+*
+      TYPE(C_PTR) DIAGF_PTR,CF_PTR,LUDF_PTR,LUCF_PTR,DIAGQ_PTR,CQ_PTR
+      REAL, POINTER, DIMENSION(:) :: DIAGF,CF,LUDF,LUCF,DIAGQ,CQ
+*----
+*  INITIALIZE POINTERS
+*----
+      LUDF=>DUMMY
+      LUCF=>DUMMY
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(NJJ(0:M),IJJ(0:M),IPOS(0:M),XSCAT(0:M*NG),TEMP(N))
+*
+      DO II=NFIRST,NGEFF
+         IF(NCONV(II)) THEN
+*           compute temp=sum_{g' ne g} Sigma_s^{g<-g'} XIN^{g'}
+            IG=NGIND(II)
+            JPSYS=KPSYS(II)
+            CALL LCMGPD(JPSYS,'DIAGF$MCCG',DIAGF_PTR)
+            CALL LCMGPD(JPSYS,'CF$MCCG',CF_PTR)
+            CALL C_F_POINTER(DIAGF_PTR,DIAGF,(/ N /))
+            CALL C_F_POINTER(CF_PTR,CF,(/ LC /))
+            IF(PACA.GE.2) THEN
+               CALL LCMGPD(JPSYS,'ILUDF$MCCG',LUDF_PTR)
+               CALL C_F_POINTER(LUDF_PTR,LUDF,(/ N /))
+               IF(PACA.LT.4) THEN
+                  CALL LCMGPD(JPSYS,'ILUCF$MCCG',LUCF_PTR)
+                  CALL C_F_POINTER(LUCF_PTR,LUCF,(/ LC /))
+               ENDIF
+            ENDIF
+            KPMACR=LCMGIL(JPMACR,IG)
+            CALL LCMGPD(JPSYS,'DIAGQ$MCCG',DIAGQ_PTR)
+            CALL LCMGPD(JPSYS,'CQ$MCCG',CQ_PTR)
+            CALL C_F_POINTER(DIAGQ_PTR,DIAGQ,(/ N /))
+            CALL C_F_POINTER(CQ_PTR,CQ,(/ LC /))
+            CALL LCMGET(KPMACR,'NJJS00',NJJ(1))
+            CALL LCMGET(KPMACR,'IJJS00',IJJ(1))
+            CALL LCMGET(KPMACR,'IPOS00',IPOS(1))
+            CALL LCMGET(KPMACR,'SCAT00',XSCAT(1))
+            DO I=1,N
+               TEMP(I)=0.0D0
+               J=IPERM(I)
+               IBM=NZON(J)
+               IF(IBM.GT.0) THEN
+                  JG=IJJ(IBM)
+                  DO 10 JND=1,NJJ(IBM)
+                  IF(JG.NE.IG) THEN
+                     JJ=NGINDV(JG)
+                     IF(JJ.GE.NFIRST) THEN
+                        TEMP(I)=TEMP(I)+XSCAT(IPOS(IBM)+JND-1)*XIN(I,JJ)
+                     ENDIF
+                  ENDIF
+                  JG=JG-1
+ 10               CONTINUE
+               ENDIF
+            ENDDO
+*           compute E^{g}*temp
+            CALL MCGPRA(LFORW,1,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGQ,CQ,
+     1           LUDF,LUCF,DIAGF,TEMP(1),XOUT(1,II),LC0,IM0,MCU0,CF)
+*           compute D^{g}*XIN^{g}
+            CALL MCGPRA(LFORW,1,PACA,.FALSE.,N,LC,IM,MCU,JU,DIAGF,CF,
+     1           LUDF,LUCF,DIAGF,XIN(1,II),TEMP(1),LC0,IM0,MCU0,CF)
+*           temp=D^{g}*XIN^{g}-E^{g}*sum_{g' ne g} Sigma_s^{g<-g'} XIN^{g'}
+            DO I=1,N
+               TEMP(I)=TEMP(I)-XOUT(I,II)
+            ENDDO
+*           apply single-group preconditioner XOUT^{g}=P^{g}*temp
+            CALL MCGPRA(LFORW,2,PACA,.TRUE.,N,LC,IM,MCU,JU,DIAGF,CF,
+     1           LUDF,LUCF,DIAGF,XOUT(1,II),TEMP(1),LC0,IM0,MCU0,CF)
+         ENDIF
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(TEMP,XSCAT,IPOS,IJJ,NJJ)
+      RETURN
+      END