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
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
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
|
