Initial commit from Polytechnique Montreal

1 files changed, 184 insertions, 0 deletions

diff --git a/Trivac/src/FLDARN.f b/Trivac/src/FLDARN.f
new file mode 100755
index 0000000..471e250
--- /dev/null
+++ b/Trivac/src/FLDARN.f
@@ -0,0 +1,184 @@
+*DECK FLDARN
+      SUBROUTINE FLDARN (FLDATV,IPTRK,IPSYS,IPFLUX,LL4,NUN,NGRP,LMOD,
+     1 IBLSZ,ADJ,IMPX,EPSOUT,MAXOUT,EVECT,FKEFFV)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of a multigroup eigenvalue system for the calculation of the
+* LMOD first orthogonal harmonics of the diffusion or SPN equation.
+* Use the implicit restarted Arnoldi method (IRAM).
+*
+*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
+* FLDATV  function pointer for the multiplication of A^(-1)B times the
+*         harmonic flux
+* IPTRK   L_TRACK pointer to the BIVAC tracking information.
+* IPSYS   L_SYSTEM pointer to system matrices.
+* IPFLUX  L_FLUX pointer to the solution.
+* LL4     order of the system matrices.
+* NUN     number of unknowns in each energy group.
+* NGRP    number of energy groups.
+* LMOD    number of orthogonal harmonics to compute.
+* IBLSZ   block size of the Arnoldi Hessenberg matrix.
+* ADJ     adjoint calculation flag.
+* IMPX    print parameter: =0: no print; =1: minimum printing.
+* EPSOUT  convergence criteria for the flux.
+* MAXOUT  maximum number of outer iterations.
+* EVECT   initial estimate of the unknown vector.
+*
+*Parameters: output
+* EVECT   converged unknown vector.
+* FKEFFV  effective multiplication factor of each harmonic.
+*
+*Reference:
+* J. BAGLAMA, "Augmented Block Householder Arnoldi Method,"
+* Linear Algebra Appl., 429, Issue 10, 2315-2334 (2008).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPSYS,IPFLUX
+      INTEGER LL4,NUN,NGRP,LMOD,IBLSZ,IMPX,MAXOUT
+      LOGICAL ADJ
+      REAL EPSOUT
+      COMPLEX EVECT(NUN,NGRP,LMOD),FKEFFV(LMOD)
+*----
+*  LOCAL VARIABLES
+*----
+      INTERFACE
+        FUNCTION FLDATV(F,N,IBLSZ,ITER,IPTRK,IPSYS,IPFLUX) RESULT(X)
+          USE GANLIB
+          INTEGER, INTENT(IN) :: N,IBLSZ,ITER
+          REAL(KIND=8), DIMENSION(N,IBLSZ), INTENT(IN) :: F
+          REAL(KIND=8), DIMENSION(N,IBLSZ) :: X
+          TYPE(C_PTR) IPTRK,IPSYS,IPFLUX
+        END FUNCTION FLDATV
+      END INTERFACE
+      REAL TIME(2)
+      REAL(KIND=8) DEPSOUT
+      CHARACTER(LEN=8) TEXT8
+      TYPE(C_PTR) JPFLUX,KPFLUX,MPFLUX
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:) :: GAR
+      COMPLEX(KIND=8), ALLOCATABLE, DIMENSION(:,:) :: V, D
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      N=LL4*NGRP
+      ALLOCATE(V(N,LMOD),D(LMOD,LMOD),GAR(NUN))
+*----
+*  SET TIMER
+*----
+*     TIME(1) : CPU TIME FOR THE SOLUTION OF LINEAR SYSTEMS.
+*     TIME(2) : CPU TIME FOR BILINEAR PRODUCT EVALUATIONS.
+      TIME(1)=0.0
+      TIME(2)=0.0
+      CALL LCMPUT(IPFLUX,'CPU-TIME',2,2,TIME)
+*----
+*  FLUX INITIALIZATION
+*----
+      DO IMOD=1,LMOD
+        V(:N,IMOD)=1.0D0
+        V(1:MIN(IBLSZ,IMOD)-1,IMOD)=0.0D0
+      ENDDO
+      CALL LCMLEN(IPFLUX,'MODE',ILONG,ITYLCM)
+      IF(ILONG.GT.0) THEN
+        DO IMOD=1,LMOD
+          JPFLUX=LCMGID(IPFLUX,'MODE')
+          CALL LCMLEL(JPFLUX,IMOD,ILONG,ITYLCM)
+          IF(ILONG.EQ.0) CYCLE
+          KPFLUX=LCMGIL(JPFLUX,IMOD)
+          IF(ADJ) THEN
+            CALL LCMLEN(KPFLUX,'AFLUX',LENA,ITYLCM)
+            IF(LENA.EQ.0) CYCLE
+            MPFLUX=LCMGID(KPFLUX,'AFLUX')
+            DO IGR=1,NGRP
+              IF(ITYLCM.EQ.2) THEN
+                CALL LCMGDL(MPFLUX,IGR,GAR)
+                EVECT(:NUN,IGR,IMOD)=GAR(:NUN)
+              ELSE IF(ITYLCM.EQ.6) THEN
+                CALL LCMGDL(MPFLUX,IGR,EVECT(1,IGR,IMOD))
+              ENDIF
+            ENDDO
+          ELSE
+            CALL LCMLEN(KPFLUX,'FLUX',LEND,ITYLCM)
+            IF(LEND.EQ.0) CYCLE
+            MPFLUX=LCMGID(KPFLUX,'FLUX')
+            DO IGR=1,NGRP
+              IF(ITYLCM.EQ.2) THEN
+                CALL LCMGDL(MPFLUX,IGR,GAR)
+                EVECT(:NUN,IGR,IMOD)=GAR(:NUN)
+              ELSE IF(ITYLCM.EQ.6) THEN
+                CALL LCMGDL(MPFLUX,IGR,EVECT(1,IGR,IMOD))
+              ENDIF
+            ENDDO
+          ENDIF
+          DO IGR=1,NGRP
+            DO IUN=1,LL4
+              IOF=(IGR-1)*LL4+IUN
+              V(IOF,IMOD)=EVECT(IUN,IGR,IMOD)
+            ENDDO
+          ENDDO
+        ENDDO
+      ENDIF
+*----
+*  CALL IRAM SOLVER
+*----
+      DEPSOUT=EPSOUT
+      CALL ALBEIGS(FLDATV,N,IBLSZ,LMOD,MAXOUT,DEPSOUT,IMPX,ITER,V,D,
+     1 IPTRK,IPSYS,IPFLUX)
+      DO IMOD=1,LMOD
+        FKEFFV(IMOD)=CMPLX(D(IMOD,IMOD),KIND=4)
+        DO IGR=1,NGRP
+          DO IUN=1,LL4
+            IOF=(IGR-1)*LL4+IUN
+            EVECT(IUN,IGR,IMOD)=CMPLX(V(IOF,IMOD),KIND=4)
+          ENDDO
+        ENDDO
+      ENDDO
+*----
+*  PRINTOUTS
+*----
+      IF(IMPX.GE.1) THEN
+        CALL LCMGET(IPFLUX,'CPU-TIME',TIME)
+        WRITE (6,650) ITER,TIME(1),TIME(2),TIME(1)+TIME(2)
+        WRITE (6,670) (FKEFFV(IMOD),IMOD=1,LMOD)
+      ENDIF
+      IF(IMPX.GE.3) THEN
+        TEXT8=' DIRECT'
+        IF(ADJ) TEXT8=' ADJOINT'
+        DO IMOD=1,LMOD
+          WRITE (6,'(/A8,13H HARMONIC NB.,I3/)') TEXT8,IMOD
+          DO IGR=1,NGRP
+            WRITE (6,680) IGR,(REAL(EVECT(I,IGR,IMOD)),I=1,LL4)
+          ENDDO
+        ENDDO
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(GAR,D,V)
+      RETURN
+*
+  650 FORMAT(/31H FLDARN: CONVERGENCE OF IRAM IN,I5,11H ITERATIONS/
+     1 9X,54HCPU TIME USED TO SOLVE THE TRIANGULAR LINEAR SYSTEMS =,
+     2 F10.3/23X,34HTO COMPUTE THE BILINEAR PRODUCTS =,F10.3,20X,
+     3 16HTOTAL CPU TIME =,F10.3)
+  670 FORMAT(//21H FLDARN: EIGENVALUES:/(5X,1P,E17.10,3H + ,E17.10,1Hi))
+  680 FORMAT(43H FLDARN: EIGENVECTOR CORRESPONDING TO GROUP,I4//
+     1 (5X,1P,8E14.5))
+      END