summaryrefslogtreecommitdiff
path: root/Trivac/src/FLDTMX.f
diff options
context:
space:
mode:
Diffstat (limited to 'Trivac/src/FLDTMX.f')
-rwxr-xr-xTrivac/src/FLDTMX.f305
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
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 22:55:06 +0000