Initial commit from Polytechnique Montreal

1 files changed, 305 insertions, 0 deletions

diff --git a/Trivac/src/FLDTMX.f b/Trivac/src/FLDTMX.f
new file mode 100755
index 0000000..aaaf33d
--- /dev/null
+++ b/Trivac/src/FLDTMX.f
@@ -0,0 +1,305 @@
+*DECK FLDTMX
+      FUNCTION FLDTMX(F,N,IBLSZ,ITER,IPTRK,IPSYS,IPFLUX) RESULT(X)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Multiplication of A^(-1)B times the harmonic flux in TRIVAC.
+*
+*Copyright:
+* Copyright (C) 2020 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): A. Hebert
+*
+*Parameters: input
+* F       harmonic flux vector.
+* N       number of unknowns in one harmonic.
+* IBLSZ   block size of the Arnoldi Hessenberg matrix.
+* ITER    Arnoldi iteration index.
+* IPTRK   L_TRACK pointer to the tracking information.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* IPFLUX  L_FLUX pointer to the solution.
+*
+*Parameters: output
+* X       result of the multiplication.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER, INTENT(IN) :: N,IBLSZ,ITER
+      COMPLEX(KIND=8), DIMENSION(N,IBLSZ), INTENT(IN) :: F
+      COMPLEX(KIND=8), DIMENSION(N,IBLSZ) :: X
+      TYPE(C_PTR) IPTRK,IPSYS,IPFLUX
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(NSTATE=40)
+      INTEGER ISTATE(NSTATE)
+      REAL EPSCON(5),TIME(2)
+      CHARACTER TEXT12*12,HSMG*131
+      LOGICAL LADJ,LUPS
+      REAL(KIND=8) DERTOL
+      INTERFACE
+        FUNCTION FLDONE_TEMPLATE(X,B,N,IPTRK,IPSYS,IPFLUX) RESULT(Y)
+          USE GANLIB
+          INTEGER, INTENT(IN) :: N
+          REAL(KIND=8), DIMENSION(N), INTENT(IN) :: X, B
+          REAL(KIND=8), DIMENSION(N) :: Y
+          TYPE(C_PTR) IPTRK,IPSYS,IPFLUX
+        END FUNCTION FLDONE_TEMPLATE
+      END INTERFACE
+      PROCEDURE(FLDONE_TEMPLATE) :: FLDONE
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, DIMENSION(:), ALLOCATABLE :: WORK
+      REAL, DIMENSION(:,:), ALLOCATABLE :: GAF1,GRAD
+      REAL, DIMENSION(:), POINTER :: AGAR
+      REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: DWORK1,DWORK2
+      TYPE(C_PTR) AGAR_PTR
+*
+*     TIME(1) : CPU TIME FOR THE SOLUTION OF LINEAR SYSTEMS.
+*     TIME(2) : CPU TIME FOR BILINEAR PRODUCT EVALUATIONS.
+      CALL LCMGET(IPFLUX,'CPU-TIME',TIME)
+      CALL KDRCPU(TK1)
+*----
+*  RECOVER INFORMATION FROM IPTRK, IPSYS AND IPFLUX
+*----
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+      NEL=ISTATE(1)
+      NUN=ISTATE(2)
+      NLF=ISTATE(30)
+      CALL LCMGET(IPSYS,'STATE-VECTOR',ISTATE)
+      NGRP=ISTATE(1)
+      LL4=ISTATE(2)
+      ITY=ISTATE(4)
+      NBMIX=ISTATE(7)
+      NAN=ISTATE(8)
+      IF(ITY.EQ.13) LL4=LL4*NLF/2 ! SPN cases
+      CALL LCMGET(IPFLUX,'STATE-VECTOR',ISTATE)
+      LADJ=ISTATE(3).EQ.10
+      ICL1=ISTATE(8)
+      ICL2=ISTATE(9)
+      IREBAL=ISTATE(10)
+      MAXINR=ISTATE(11)
+      NADI=ISTATE(13)
+      NSTARD=ISTATE(15)
+      IMPX=ISTATE(40)
+      CALL LCMGET(IPFLUX,'EPS-CONVERGE',EPSCON)
+      EPSINR=EPSCON(1)
+      EPSMSR=EPSCON(4)
+      IF(LL4*NGRP.NE.N) CALL XABORT('FLDTMX: INCONSISTENT UNKNOWNS.')
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(WORK(NUN),GAF1(NUN,NGRP),GRAD(NUN,NGRP))
+*----
+*  CHECK FOR UP-SCATTERING.
+*----
+      LUPS=.FALSE.
+      DO 20 IGR=1,NGRP-1
+      DO 10 JGR=IGR+1,NGRP
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.GT.0) THEN
+        LUPS=.TRUE.
+        MAXINR=MAX(MAXINR,10)
+        GO TO 30
+      ENDIF
+   10 CONTINUE
+   20 CONTINUE
+*----
+*  MAIN LOOP OVER MODES.
+*----
+   30 DO 240 IMOD=1,IBLSZ
+      IF(LADJ) THEN
+*       ADJOINT SOLUTION
+*----
+*  COMPUTE B TIMES THE FLUX.
+*----
+        DO 70 IGR=1,NGRP
+        DO 40 I=1,LL4
+        GAF1(I,IGR)=0.0
+   40   CONTINUE
+        DO 60 JGR=1,NGRP
+        WRITE(TEXT12,'(1HB,2I3.3)') JGR,IGR
+        CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+        IF(ILONG.EQ.0) GO TO 60
+        CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+        CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /))
+        DO 50 I=1,ILONG
+        IOF=(JGR-1)*LL4+I
+        GAF1(I,IGR)=GAF1(I,IGR)+AGAR(I)*REAL(F(IOF,IMOD),KIND=4)
+        IF(ABS(AIMAG(F(IOF,IMOD))).GT.1.0E-8) THEN
+          WRITE(HSMG,'(13HFLDTMX: FLUX(,2I8,2H)=,1P,2E12.4,
+     1    12H IS COMPLEX.)') IOF,IMOD,F(IOF,IMOD)
+          CALL XABORT(HSMG)
+        ENDIF
+   50   CONTINUE
+   60   CONTINUE
+   70   CONTINUE
+        CALL KDRCPU(TK2)
+        TIME(2)=TIME(2)+(TK2-TK1)
+*----
+*  COMPUTE A^(-1)B WITHOUT DOWN-SCATTERING.
+*----
+        DO 120 IGR=NGRP,1,-1
+        CALL KDRCPU(TK1)
+        DO 80 I=1,LL4
+        GRAD(I,IGR)=GAF1(I,IGR)
+   80   CONTINUE
+        DO 110 JGR=NGRP,IGR+1,-1
+        WRITE(TEXT12,'(1HA,2I3.3)') JGR,IGR
+        CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+        IF(ILONG.EQ.0) GO TO 110
+        IF(ITY.EQ.13) THEN
+          CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,JGR),WORK)
+          DO 90 I=1,LL4
+          GRAD(I,IGR)=GRAD(I,IGR)+WORK(I)
+   90     CONTINUE
+        ELSE
+          CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+          CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /))
+          DO 100 I=1,ILONG
+          GRAD(I,IGR)=GRAD(I,IGR)+AGAR(I)*GRAD(I,JGR)
+  100     CONTINUE
+        ENDIF
+  110   CONTINUE
+        CALL KDRCPU(TK2)
+        TIME(2)=TIME(2)+(TK2-TK1)
+*
+        CALL KDRCPU(TK1)
+        WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+        IF(NSTARD.EQ.0) THEN
+          WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+          CALL FLDADI(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,IGR),NADI)
+          JTER=NADI
+        ELSE
+*         use a GMRES solution of the linear system
+          DERTOL=EPSMSR
+          ISTATE(39)=IGR
+          CALL LCMPUT(IPFLUX,'STATE-VECTOR',NSTATE,1,ISTATE)
+          ALLOCATE(DWORK1(LL4),DWORK2(LL4))
+          DWORK1(:LL4)=GRAD(:LL4,IGR) ! source
+          DWORK2(:LL4)=0.0            ! estimate of the flux
+          CALL FLDMRA(DWORK1,FLDONE,LL4,DERTOL,NSTARD,NADI,IMPX,IPTRK,
+     1    IPSYS,IPFLUX,DWORK2,JTER)
+          GRAD(:LL4,IGR)=REAL(DWORK2(:LL4))
+          DEALLOCATE(DWORK2,DWORK1)
+        ENDIF
+        CALL KDRCPU(TK2)
+        TIME(1)=TIME(1)+(TK2-TK1)
+  120   CONTINUE
+      ELSE
+*       DIRECT SOLUTION
+*----
+*  COMPUTE B TIMES THE FLUX.
+*----
+        DO 160 IGR=1,NGRP
+        DO 130 I=1,LL4
+        GAF1(I,IGR)=0.0
+  130   CONTINUE
+        DO 150 JGR=1,NGRP
+        WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+        CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+        IF(ILONG.EQ.0) GO TO 150
+        CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+        CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /))
+        DO 140 I=1,ILONG
+        IOF=(JGR-1)*LL4+I
+        GAF1(I,IGR)=GAF1(I,IGR)+AGAR(I)*REAL(F(IOF,IMOD),KIND=4)
+        IF(ABS(AIMAG(F(IOF,IMOD))).GT.1.0E-8) THEN
+          WRITE(HSMG,'(13HFLDTMX: FLUX(,2I8,2H)=,1P,2E12.4,
+     1    12H IS COMPLEX.)') IOF,IMOD,F(IOF,IMOD)
+          CALL XABORT(HSMG)
+        ENDIF
+  140   CONTINUE
+  150   CONTINUE
+  160   CONTINUE
+        CALL KDRCPU(TK2)
+        TIME(2)=TIME(2)+(TK2-TK1)
+*----
+*  COMPUTE A^(-1)B WITHOUT UP-SCATTERING.
+*----
+        DO 210 IGR=1,NGRP
+        CALL KDRCPU(TK1)
+        DO 170 I=1,LL4
+        GRAD(I,IGR)=GAF1(I,IGR)
+  170   CONTINUE
+        DO 200 JGR=1,IGR-1
+        WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+        CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+        IF(ILONG.EQ.0) GO TO 200
+        IF(ITY.EQ.13) THEN
+          CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,JGR),WORK)
+          DO 180 I=1,LL4
+          GRAD(I,IGR)=GRAD(I,IGR)+WORK(I)
+  180     CONTINUE
+        ELSE
+          CALL LCMGPD(IPSYS,TEXT12,AGAR_PTR)
+          CALL C_F_POINTER(AGAR_PTR,AGAR,(/ ILONG /))
+          DO 190 I=1,ILONG
+          GRAD(I,IGR)=GRAD(I,IGR)+AGAR(I)*GRAD(I,JGR)
+  190     CONTINUE
+        ENDIF
+  200   CONTINUE
+        CALL KDRCPU(TK2)
+        TIME(2)=TIME(2)+(TK2-TK1)
+*
+        CALL KDRCPU(TK1)
+        IF(NSTARD.EQ.0) THEN
+          WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+          CALL FLDADI(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,IGR),NADI)
+          JTER=-NADI
+        ELSE
+*         use a GMRES solution of the linear system
+          DERTOL=EPSMSR
+          ISTATE(39)=IGR
+          CALL LCMPUT(IPFLUX,'STATE-VECTOR',NSTATE,1,ISTATE)
+          ALLOCATE(DWORK1(LL4),DWORK2(LL4))
+          DWORK1(:LL4)=GRAD(:LL4,IGR) ! source
+          DWORK2(:LL4)=0.0            ! estimate of the flux
+          CALL FLDMRA(DWORK1,FLDONE,LL4,DERTOL,NSTARD,NADI,IMPX,IPTRK,
+     1    IPSYS,IPFLUX,DWORK2,JTER)
+          GRAD(:LL4,IGR)=REAL(DWORK2(:LL4))
+          DEALLOCATE(DWORK2,DWORK1)
+        ENDIF
+        CALL KDRCPU(TK2)
+        TIME(1)=TIME(1)+(TK2-TK1)
+  210   CONTINUE
+      ENDIF
+*----
+*  PERFORM THERMAL (UP/DOWN-SCATTERING) ITERATIONS.
+*----
+      KTER=0
+      IF((IREBAL.EQ.1).OR.LUPS) THEN
+        CALL FLDTHR(IPTRK,IPSYS,IPFLUX,LADJ,LL4,ITY,NUN,NGRP,ICL1,ICL2,
+     1  IMPX,NADI,NSTARD,MAXINR,EPSINR,KTER,TIME(1),TIME(2),GRAD)
+      ENDIF
+      DO 230 IGR=1,NGRP
+      DO 220 I=1,LL4
+      IOF=(IGR-1)*LL4+I
+      X(IOF,IMOD)=GRAD(I,IGR)
+  220 CONTINUE
+  230 CONTINUE
+*----
+*  END OF LOOP OVER MODES.
+*----
+  240 CONTINUE
+      CALL LCMPUT(IPFLUX,'CPU-TIME',2,2,TIME)
+      IF(IMPX.GT.10) WRITE(6,250) ITER,JTER,KTER
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(GRAD,GAF1,WORK)
+      RETURN
+  250 FORMAT(49H FLDTMX: MATRIX MULTIPLICATION AT IRAM ITERATION=,I5,
+     1 18H INNER ITERATIONS=,I5,20H THERMAL ITERATIONS=,I5)
+      END FUNCTION FLDTMX