Initial commit from Polytechnique Montreal

1 files changed, 192 insertions, 0 deletions

diff --git a/Trivac/src/FLDBMX.f b/Trivac/src/FLDBMX.f
new file mode 100755
index 0000000..8aaa075
--- /dev/null
+++ b/Trivac/src/FLDBMX.f
@@ -0,0 +1,192 @@
+*DECK FLDBMX
+      FUNCTION FLDBMX(F,N,IBLSZ,ITER,IPTRK,IPSYS,IPFLUX) RESULT(X)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Multiplication of A^(-1)B times the harmonic flux in BIVAC.
+*
+*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 LUPS
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, DIMENSION(:), ALLOCATABLE :: WORK1,WORK2
+      REAL, DIMENSION(:,:), ALLOCATABLE :: GAF1,GRAD
+*
+*     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(14)
+      CALL LCMGET(IPSYS,'STATE-VECTOR',ISTATE)
+      NGRP=ISTATE(1)
+      LL4=ISTATE(2)
+      ITY=ISTATE(4)
+      NBMIX=ISTATE(7)
+      NAN=ISTATE(8)
+      IF(ITY.EQ.11) LL4=LL4*NLF/2 ! SPN cases
+      CALL LCMGET(IPFLUX,'STATE-VECTOR',ISTATE)
+      ICL1=ISTATE(8)
+      ICL2=ISTATE(9)
+      IREBAL=ISTATE(10)
+      MAXINR=ISTATE(11)
+      NADI=ISTATE(13)
+      IMPX=ISTATE(40)
+      CALL LCMGET(IPFLUX,'EPS-CONVERGE',EPSCON)
+      EPSINR=EPSCON(1)
+      IF(LL4*NGRP.NE.N) CALL XABORT('FLDBMX: INCONSISTENT UNKNOWNS.')
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(WORK1(NUN),WORK2(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 150 IMOD=1,IBLSZ
+*----
+*  COMPUTE B TIMES THE FLUX.
+*----
+      DO 80 IGR=1,NGRP
+      DO 40 I=1,LL4
+      GAF1(I,IGR)=0.0
+   40 CONTINUE
+      DO 70 JGR=1,NGRP
+      WRITE(TEXT12,'(1HB,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 70
+      DO 50 I=1,LL4
+      IOF=(JGR-1)*LL4+I
+      WORK1(I)=REAL(F(IOF,IMOD),KIND=4)
+      IF(ABS(AIMAG(F(IOF,IMOD))).GT.1.0E-8) THEN
+        WRITE(HSMG,'(13HFLDBMX: FLUX(,2I8,2H)=,1P,2E12.4,
+     1  12H IS COMPLEX.)') IOF,IMOD,F(IOF,IMOD)
+        CALL XABORT(HSMG)
+      ENDIF
+   50 CONTINUE
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,WORK1(1),WORK2(1))
+      DO 60 I=1,LL4
+      GAF1(I,IGR)=GAF1(I,IGR)+WORK2(I)
+   60 CONTINUE
+   70 CONTINUE
+   80 CONTINUE
+      CALL KDRCPU(TK2)
+      TIME(2)=TIME(2)+(TK2-TK1)
+*----
+*  COMPUTE A^(-1)B WITHOUT UP-SCATTERING.
+*----
+      DO 120 IGR=1,NGRP
+      CALL KDRCPU(TK1)
+      DO 90 I=1,LL4
+      GRAD(I,IGR)=GAF1(I,IGR)
+   90 CONTINUE
+      DO 110 JGR=1,IGR-1
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,JGR
+      CALL LCMLEN(IPSYS,TEXT12,ILONG,ITYLCM)
+      IF(ILONG.EQ.0) GO TO 110
+      CALL MTLDLM(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,JGR),WORK2(1))
+      DO 100 I=1,LL4
+      GRAD(I,IGR)=GRAD(I,IGR)+WORK2(I)
+  100 CONTINUE
+  110 CONTINUE
+      CALL KDRCPU(TK2)
+      TIME(2)=TIME(2)+(TK2-TK1)
+*
+      CALL KDRCPU(TK1)
+      WRITE(TEXT12,'(1HA,2I3.3)') IGR,IGR
+      IF(ITY.EQ.11) THEN
+*       SIMPLIFIED PN BIVAC TRACKING.
+        IF(NAN.EQ.0) CALL XABORT('FLDBMX: SPN-ONLY ALGORITHM.')
+        CALL FLDBSS(TEXT12,IPTRK,IPSYS,LL4,NBMIX,NAN,GRAD(1,IGR),NADI)
+      ELSE
+        CALL MTLDLS(TEXT12,IPTRK,IPSYS,LL4,ITY,GRAD(1,IGR))
+      ENDIF
+      CALL KDRCPU(TK2)
+      TIME(1)=TIME(1)+(TK2-TK1)
+  120 CONTINUE
+*----
+*  PERFORM THERMAL (UP-SCATTERING) ITERATIONS.
+*----
+      KTER=0
+      IF((IREBAL.EQ.1).OR.LUPS) THEN
+        CALL FLDBHR(IPTRK,IPSYS,.FALSE.,LL4,ITY,NUN,NGRP,ICL1,ICL2,IMPX,
+     1  MAXINR,EPSINR,KTER,TIME(1),TIME(2),GRAD)
+      ENDIF
+      DO 140 IGR=1,NGRP
+      DO 130 I=1,LL4
+      IOF=(IGR-1)*LL4+I
+      X(IOF,IMOD)=GRAD(I,IGR)
+  130 CONTINUE
+  140 CONTINUE
+*----
+*  END OF LOOP OVER MODES.
+*----
+  150 CONTINUE
+      CALL LCMPUT(IPFLUX,'CPU-TIME',2,2,TIME)
+      IF(IMPX.GT.10) WRITE(6,200) ITER,KTER
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(GRAD,GAF1,WORK2,WORK1)
+      RETURN
+  200 FORMAT(49H FLDBMX: MATRIX MULTIPLICATION AT IRAM ITERATION=,I5,
+     1 20H THERMAL ITERATIONS=,I5)
+      END FUNCTION FLDBMX