1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
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
|
