FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_precond_RIF_33.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
7 use hecmw_util
11
12 private
13
17
18 integer(4),parameter :: krealp = 8
19
20 integer(kind=kint) :: NPFIU, NPFIL
21 integer(kind=kint) :: N
22 integer(kind=kint), pointer :: inumFI1L(:) => null()
23 integer(kind=kint), pointer :: inumFI1U(:) => null()
24 integer(kind=kint), pointer :: FI1L(:) => null()
25 integer(kind=kint), pointer :: FI1U(:) => null()
26
27 integer(kind=kint), pointer :: indexL(:) => null()
28 integer(kind=kint), pointer :: indexU(:) => null()
29 integer(kind=kint), pointer :: itemL(:) => null()
30 integer(kind=kint), pointer :: itemU(:) => null()
31 real(kind=kreal), pointer :: d(:) => null()
32 real(kind=kreal), pointer :: al(:) => null()
33 real(kind=kreal), pointer :: au(:) => null()
34
35 real(kind=krealp), pointer :: sainvl(:) => null()
36 real(kind=krealp), pointer :: sainvd(:) => null()
37 real(kind=krealp), pointer :: rifu(:) => null()
38 real(kind=krealp), pointer :: rifl(:) => null()
39 real(kind=krealp), pointer :: rifd(:) => null()
40
41contains
42
43 !C***
44 !C*** hecmw_precond_33_sainv_setup
45 !C***
46 subroutine hecmw_precond_rif_33_setup(hecMAT)
47 implicit none
48 type(hecmwst_matrix) :: hecmat
49
50 integer(kind=kint ) :: precond
51
52 real(kind=krealp) :: filter
53
54 n = hecmat%N
55 precond = hecmw_mat_get_precond(hecmat)
56
57 d => hecmat%D
58 au=> hecmat%AU
59 al=> hecmat%AL
60 indexl => hecmat%indexL
61 indexu => hecmat%indexU
62 iteml => hecmat%itemL
63 itemu => hecmat%itemU
64
65 if (precond.eq.21) call form_ilu1_rif_33(hecmat)
66
67 allocate (sainvd(9*hecmat%NP))
68 allocate (sainvl(9*npfiu))
69 allocate (rifd(9*hecmat%NP))
70 allocate (rifu(9*npfiu))
71 sainvd = 0.0d0
72 sainvl = 0.0d0
73 rifd = 0.0d0
74 rifu = 0.0d0
75
76 filter= hecmat%Rarray(5)
77
78 write(*,"(a,F15.8)")"### RIF FILTER :",filter
79
80 call hecmw_rif_33(hecmat)
81
82 allocate (rifl(9*npfiu))
83 rifl = 0.0d0
84
85 call hecmw_rif_make_u_33(hecmat)
86
87 end subroutine hecmw_precond_rif_33_setup
88
90 implicit none
91 real(kind=kreal), intent(inout) :: zp(:)
92 integer(kind=kint) :: in, i, j, isl, iel
93 real(kind=kreal) :: sw1, sw2, sw3, x1, x2, x3
94
95 !C-- FORWARD
96 do i= 1, n
97 sw1= zp(3*i-2)
98 sw2= zp(3*i-1)
99 sw3= zp(3*i )
100
101 isl= inumfi1l(i-1)+1
102 iel= inumfi1l(i)
103 do j= isl, iel
104 in= fi1l(j)
105 x1= zp(3*in-2)
106 x2= zp(3*in-1)
107 x3= zp(3*in )
108 sw1= sw1 - rifl(9*j-8)*x1 - rifl(9*j-7)*x2 - rifl(9*j-6)*x3
109 sw2= sw2 - rifl(9*j-5)*x1 - rifl(9*j-4)*x2 - rifl(9*j-3)*x3
110 sw3= sw3 - rifl(9*j-2)*x1 - rifl(9*j-1)*x2 - rifl(9*j )*x3
111 enddo
112 x1= sw1
113 x2= sw2
114 x3= sw3
115 x2= x2 - rifd(9*i-5)*x1
116 x3= x3 - rifd(9*i-2)*x1 - rifd(9*i-1)*x2
117 zp(3*i-2)= x1
118 zp(3*i-1)= x2
119 zp(3*i )= x3
120 enddo
121
122 do i= 1, n
123 zp(3*i-2)= zp(3*i-2)*rifd(9*i-8)
124 zp(3*i-1)= zp(3*i-1)*rifd(9*i-4)
125 zp(3*i )= zp(3*i )*rifd(9*i )
126 enddo
127
128 !C-- BACKWARD
129 do i= n, 1, -1
130 sw1= zp(3*i-2)
131 sw2= zp(3*i-1)
132 sw3= zp(3*i )
133
134 isl= inumfi1u(i-1) + 1
135 iel= inumfi1u(i)
136 do j= iel, isl, -1
137 in= fi1u(j)
138 x1= zp(3*in-2)
139 x2= zp(3*in-1)
140 x3= zp(3*in )
141 sw1= sw1 - rifu(9*j-8)*x1 - rifu(9*j-7)*x2 - rifu(9*j-6)*x3
142 sw2= sw2 - rifu(9*j-5)*x1 - rifu(9*j-4)*x2 - rifu(9*j-3)*x3
143 sw3= sw3 - rifu(9*j-2)*x1 - rifu(9*j-1)*x2 - rifu(9*j )*x3
144 enddo
145
146 x1= sw1
147 x2= sw2
148 x3= sw3
149 x2= x2 - rifd(9*i-1)*x3
150 x1= x1 - rifd(9*i-2)*x3 - rifd(9*i-5)*x2
151 zp(3*i-2)= x1
152 zp(3*i-1)= x2
153 zp(3*i )= x3
154 enddo
155 end subroutine hecmw_precond_rif_33_apply
156
157
158 !C***
159 !C*** hecmw_rif_33
160 !C***
161 subroutine hecmw_rif_33(hecMAT)
162 implicit none
163 type (hecmwst_matrix) :: hecmat
164 integer(kind=kint) :: i, j, js, je, in, itr, np
165 real(kind=krealp) :: x1, x2, x3, dd, dd1, dd2, dd3, dtemp(3)
166 real(kind=krealp) :: filter
167 real(kind=krealp), allocatable :: zz(:), vv(:)
168
169 filter= hecmat%Rarray(5)
170
171 np = hecmat%NP
172 allocate(vv(3*np))
173 allocate(zz(3*np))
174
175 do itr=1,np
176
177 !------------------------------ iitr = 1 ----------------------------------------
178
179 zz(:) = 0.0d0
180 vv(:) = 0.0d0
181
182 !{v}=[A]{zi}
183
184 zz(3*itr-2)= sainvd(9*itr-8)
185 zz(3*itr-1)= sainvd(9*itr-5)
186 zz(3*itr )= sainvd(9*itr-2)
187
188 zz(3*itr-2)= 1.0d0
189
190 js= inumfi1l(itr-1) + 1
191 je= inumfi1l(itr )
192 do j= js, je
193 in = fi1l(j)
194 zz(3*in-2)= sainvl(9*j-8)
195 zz(3*in-1)= sainvl(9*j-7)
196 zz(3*in )= sainvl(9*j-6)
197 enddo
198
199 do i= 1, itr
200 x1= zz(3*i-2)
201 x2= zz(3*i-1)
202 x3= zz(3*i )
203 vv(3*i-2) = vv(3*i-2) + d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
204 vv(3*i-1) = vv(3*i-1) + d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
205 vv(3*i ) = vv(3*i ) + d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i )*x3
206
207 js= indexl(i-1) + 1
208 je= indexl(i )
209 do j=js,je
210 in = iteml(j)
211 vv(3*in-2)= vv(3*in-2) + al(9*j-8)*x1 + al(9*j-5)*x2 + al(9*j-2)*x3
212 vv(3*in-1)= vv(3*in-1) + al(9*j-7)*x1 + al(9*j-4)*x2 + al(9*j-1)*x3
213 vv(3*in )= vv(3*in ) + al(9*j-6)*x1 + al(9*j-3)*x2 + al(9*j )*x3
214 enddo
215
216 js= indexu(i-1) + 1
217 je= indexu(i )
218 do j= js, je
219 in = itemu(j)
220 vv(3*in-2)= vv(3*in-2) + au(9*j-8)*x1 + au(9*j-5)*x2 + au(9*j-2)*x3
221 vv(3*in-1)= vv(3*in-1) + au(9*j-7)*x1 + au(9*j-4)*x2 + au(9*j-1)*x3
222 vv(3*in )= vv(3*in ) + au(9*j-6)*x1 + au(9*j-3)*x2 + au(9*j )*x3
223 enddo
224 enddo
225
226 !{d}={v^t}{z_j}
227 do i= itr,np
228 sainvd(9*i-8) = vv(3*i-2)
229 sainvd(9*i-4) = vv(3*i-2)*sainvd(9*i-7) + vv(3*i-1)
230 sainvd(9*i ) = vv(3*i-2)*sainvd(9*i-6) + vv(3*i-1)*sainvd(9*i-3) + vv(3*i)
231 js= inumfi1l(i-1) + 1
232 je= inumfi1l(i )
233 do j= js, je
234 in = fi1l(j)
235 x1= vv(3*in-2)
236 x2= vv(3*in-1)
237 x3= vv(3*in )
238 sainvd(9*i-8)= sainvd(9*i-8) + x1*sainvl(9*j-8) + x2*sainvl(9*j-7) + x3*sainvl(9*j-6)
239 sainvd(9*i-4)= sainvd(9*i-4) + x1*sainvl(9*j-5) + x2*sainvl(9*j-4) + x3*sainvl(9*j-3)
240 sainvd(9*i )= sainvd(9*i ) + x1*sainvl(9*j-2) + x2*sainvl(9*j-1) + x3*sainvl(9*j )
241 enddo
242 enddo
243
244 !Update D
245 dd = 1.0d0/sainvd(9*itr-8)
246
247 sainvd(9*itr-4) =sainvd(9*itr-4)*dd
248 sainvd(9*itr ) =sainvd(9*itr )*dd
249
250 do i =itr+1,np
251 sainvd(9*i-8) = sainvd(9*i-8)*dd
252 sainvd(9*i-4) = sainvd(9*i-4)*dd
253 sainvd(9*i ) = sainvd(9*i )*dd
254 enddo
255
256 rifd(9*itr-8) = dd !RIF UPDATE
257 rifd(9*itr-5) = sainvd(9*itr-4) !RIF UPDATE
258 rifd(9*itr-2) = sainvd(9*itr ) !RIF UPDATE
259
260 js= inumfi1u(itr-1) + 1
261 je= inumfi1u(itr )
262 do j= js, je
263 rifu(9*j-8) = sainvd(9*fi1u(j)-8)
264 rifu(9*j-7) = sainvd(9*fi1u(j)-4)
265 rifu(9*j-6) = sainvd(9*fi1u(j) )
266 enddo
267
268 !Update Z
269
270 dd2=sainvd(9*itr-4)
271 if(abs(dd2) > filter)then
272 sainvd(9*itr-7)= sainvd(9*itr-7) - dd2*zz(3*itr-2)
273 js= inumfi1l(itr-1) + 1
274 je= inumfi1l(itr )
275 do j= js, je
276 in = fi1l(j)
277 sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
278 sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
279 sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
280 enddo
281 endif
282
283 dd3=sainvd(9*itr )
284 if(abs(dd3) > filter)then
285 sainvd(9*itr-6)= sainvd(9*itr-6) - dd3*zz(3*itr-2)
286 js= inumfi1l(itr-1) + 1
287 je= inumfi1l(itr )
288 do j= js, je
289 in = fi1l(j)
290 sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
291 sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
292 sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
293 enddo
294 endif
295
296 do i= itr +1,np
297 js= inumfi1l(i-1) + 1
298 je= inumfi1l(i )
299 dd1=sainvd(9*i-8)
300 if(abs(dd1) > filter)then
301 do j= js, je
302 in = fi1l(j)
303 if (in > itr) exit
304 sainvl(9*j-8) = sainvl(9*j-8)-dd1*zz(3*in-2)
305 sainvl(9*j-7) = sainvl(9*j-7)-dd1*zz(3*in-1)
306 sainvl(9*j-6) = sainvl(9*j-6)-dd1*zz(3*in )
307 enddo
308 endif
309 dd2=sainvd(9*i-4)
310 if(abs(dd2) > filter)then
311 do j= js, je
312 in = fi1l(j)
313 if (in > itr) exit
314 sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
315 sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
316 sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
317 enddo
318 endif
319 dd3=sainvd(9*i )
320 if(abs(dd3) > filter)then
321 do j= js, je
322 in = fi1l(j)
323 if (in > itr) exit
324 sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
325 sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
326 sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
327 enddo
328 endif
329 enddo
330
331 !------------------------------ iitr = 1 ----------------------------------------
332
333 zz(:) = 0.0d0
334 vv(:) = 0.0d0
335
336 !{v}=[A]{zi}
337
338 zz(3*itr-2)= sainvd(9*itr-7)
339 zz(3*itr-1)= sainvd(9*itr-4)
340 zz(3*itr )= sainvd(9*itr-1)
341
342 zz(3*itr-1)= 1.0d0
343
344 js= inumfi1l(itr-1) + 1
345 je= inumfi1l(itr )
346 do j= js, je
347 in = fi1l(j)
348 zz(3*in-2)= sainvl(9*j-5)
349 zz(3*in-1)= sainvl(9*j-4)
350 zz(3*in )= sainvl(9*j-3)
351 enddo
352
353 do i= 1, itr
354 x1= zz(3*i-2)
355 x2= zz(3*i-1)
356 x3= zz(3*i )
357 vv(3*i-2) = vv(3*i-2) + d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
358 vv(3*i-1) = vv(3*i-1) + d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
359 vv(3*i ) = vv(3*i ) + d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i )*x3
360
361 js= indexl(i-1) + 1
362 je= indexl(i )
363 do j=js,je
364 in = iteml(j)
365 vv(3*in-2)= vv(3*in-2) + al(9*j-8)*x1 + al(9*j-5)*x2 + al(9*j-2)*x3
366 vv(3*in-1)= vv(3*in-1) + al(9*j-7)*x1 + al(9*j-4)*x2 + al(9*j-1)*x3
367 vv(3*in )= vv(3*in ) + al(9*j-6)*x1 + al(9*j-3)*x2 + al(9*j )*x3
368 enddo
369
370 js= indexu(i-1) + 1
371 je= indexu(i )
372 do j= js, je
373 in = itemu(j)
374 vv(3*in-2)= vv(3*in-2) + au(9*j-8)*x1 + au(9*j-5)*x2 + au(9*j-2)*x3
375 vv(3*in-1)= vv(3*in-1) + au(9*j-7)*x1 + au(9*j-4)*x2 + au(9*j-1)*x3
376 vv(3*in )= vv(3*in ) + au(9*j-6)*x1 + au(9*j-3)*x2 + au(9*j )*x3
377 enddo
378 enddo
379
380 !{d}={v^t}{z_j}
381 dtemp(1) = sainvd(9*itr-8)
382
383 do i=itr,np
384 sainvd(9*i-8) = vv(3*i-2)
385 sainvd(9*i-4) = vv(3*i-2)*sainvd(9*i-7) + vv(3*i-1)
386 sainvd(9*i ) = vv(3*i-2)*sainvd(9*i-6) + vv(3*i-1)*sainvd(9*i-3) + vv(3*i)
387 js= inumfi1l(i-1) + 1
388 je= inumfi1l(i )
389 do j= js, je
390 in = fi1l(j)
391 x1= vv(3*in-2)
392 x2= vv(3*in-1)
393 x3= vv(3*in )
394 sainvd(9*i-8)= sainvd(9*i-8) + x1*sainvl(9*j-8) + x2*sainvl(9*j-7) + x3*sainvl(9*j-6)
395 sainvd(9*i-4)= sainvd(9*i-4) + x1*sainvl(9*j-5) + x2*sainvl(9*j-4) + x3*sainvl(9*j-3)
396 sainvd(9*i )= sainvd(9*i ) + x1*sainvl(9*j-2) + x2*sainvl(9*j-1) + x3*sainvl(9*j )
397 enddo
398 enddo
399
400 !Update D
401 dd = 1.0d0/sainvd(9*itr-4)
402
403 sainvd(9*itr-8) = dtemp(1)
404 sainvd(9*itr ) =sainvd(9*itr )*dd
405
406 do i =itr+1,np
407 sainvd(9*i-8) = sainvd(9*i-8)*dd
408 sainvd(9*i-4) = sainvd(9*i-4)*dd
409 sainvd(9*i ) = sainvd(9*i )*dd
410 enddo
411
412 rifd(9*itr-4) = dd !RIF UPDATE
413 rifd(9*itr-1) = sainvd(9*itr ) !RIF UPDATE
414
415 js= inumfi1u(itr-1) + 1
416 je= inumfi1u(itr )
417 do j= js, je
418 rifu(9*j-5) = sainvd(9*fi1u(j)-8)
419 rifu(9*j-4) = sainvd(9*fi1u(j)-4)
420 rifu(9*j-3) = sainvd(9*fi1u(j) )
421 enddo
422
423 !Update Z
424
425 dd3=sainvd(9*itr )
426 if(abs(dd3) > filter)then
427 sainvd(9*itr-6)= sainvd(9*itr-6) - dd3*zz(3*itr-2)
428 sainvd(9*itr-3)= sainvd(9*itr-3) - dd3*zz(3*itr-1)
429
430 js= inumfi1l(itr-1) + 1
431 je= inumfi1l(itr )
432 do j= js, je
433 in = fi1l(j)
434 sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
435 sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
436 sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
437 enddo
438 endif
439
440 do i= itr +1,np
441 js= inumfi1l(i-1) + 1
442 je= inumfi1l(i )
443 dd1=sainvd(9*i-8)
444 if(abs(dd1) > filter)then
445 do j= js, je
446 in = fi1l(j)
447 if (in > itr) exit
448 sainvl(9*j-8) = sainvl(9*j-8)-dd1*zz(3*in-2)
449 sainvl(9*j-7) = sainvl(9*j-7)-dd1*zz(3*in-1)
450 sainvl(9*j-6) = sainvl(9*j-6)-dd1*zz(3*in )
451 enddo
452 endif
453 dd2=sainvd(9*i-4)
454 if(abs(dd2) > filter)then
455 do j= js, je
456 in = fi1l(j)
457 if (in > itr) exit
458 sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
459 sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
460 sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
461 enddo
462 endif
463 dd3=sainvd(9*i )
464 if(abs(dd3) > filter)then
465 do j= js, je
466 in = fi1l(j)
467 if (in > itr) exit
468 sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
469 sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
470 sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
471 enddo
472 endif
473 enddo
474
475
476 !------------------------------ iitr = 1 ----------------------------------------
477
478 zz(:) = 0.0d0
479 vv(:) = 0.0d0
480
481 !{v}=[A]{zi}
482
483 zz(3*itr-2)= sainvd(9*itr-6)
484 zz(3*itr-1)= sainvd(9*itr-3)
485 zz(3*itr )= sainvd(9*itr )
486
487 zz(3*itr )= 1.0d0
488
489 js= inumfi1l(itr-1) + 1
490 je= inumfi1l(itr )
491 do j= js, je
492 in = fi1l(j)
493 zz(3*in-2)= sainvl(9*j-2)
494 zz(3*in-1)= sainvl(9*j-1)
495 zz(3*in )= sainvl(9*j )
496 enddo
497
498 do i= 1, itr
499 x1= zz(3*i-2)
500 x2= zz(3*i-1)
501 x3= zz(3*i )
502 vv(3*i-2) = vv(3*i-2) + d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
503 vv(3*i-1) = vv(3*i-1) + d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
504 vv(3*i ) = vv(3*i ) + d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i )*x3
505
506 js= indexl(i-1) + 1
507 je= indexl(i )
508 do j=js,je
509 in = iteml(j)
510 vv(3*in-2)= vv(3*in-2) + al(9*j-8)*x1 + al(9*j-5)*x2 + al(9*j-2)*x3
511 vv(3*in-1)= vv(3*in-1) + al(9*j-7)*x1 + al(9*j-4)*x2 + al(9*j-1)*x3
512 vv(3*in )= vv(3*in ) + al(9*j-6)*x1 + al(9*j-3)*x2 + al(9*j )*x3
513 enddo
514
515 js= indexu(i-1) + 1
516 je= indexu(i )
517 do j= js, je
518 in = itemu(j)
519 vv(3*in-2)= vv(3*in-2) + au(9*j-8)*x1 + au(9*j-5)*x2 + au(9*j-2)*x3
520 vv(3*in-1)= vv(3*in-1) + au(9*j-7)*x1 + au(9*j-4)*x2 + au(9*j-1)*x3
521 vv(3*in )= vv(3*in ) + au(9*j-6)*x1 + au(9*j-3)*x2 + au(9*j )*x3
522 enddo
523 enddo
524
525 !{d}={v^t}{z_j}
526
527 dtemp(1) = sainvd(9*itr-8)
528 dtemp(2) = sainvd(9*itr-4)
529
530 do i=itr,np
531 sainvd(9*i-8) = vv(3*i-2)
532 sainvd(9*i-4) = vv(3*i-2)*sainvd(9*i-7) + vv(3*i-1)
533 sainvd(9*i ) = vv(3*i-2)*sainvd(9*i-6) + vv(3*i-1)*sainvd(9*i-3) + vv(3*i)
534 js= inumfi1l(i-1) + 1
535 je= inumfi1l(i )
536 do j= js, je
537 in = fi1l(j)
538 x1= vv(3*in-2)
539 x2= vv(3*in-1)
540 x3= vv(3*in )
541 sainvd(9*i-8)= sainvd(9*i-8) + x1*sainvl(9*j-8) + x2*sainvl(9*j-7) + x3*sainvl(9*j-6)
542 sainvd(9*i-4)= sainvd(9*i-4) + x1*sainvl(9*j-5) + x2*sainvl(9*j-4) + x3*sainvl(9*j-3)
543 sainvd(9*i )= sainvd(9*i ) + x1*sainvl(9*j-2) + x2*sainvl(9*j-1) + x3*sainvl(9*j )
544 enddo
545 enddo
546
547 !Update D
548 dd = 1.0d0/sainvd(9*itr )
549
550 sainvd(9*itr-8) = dtemp(1)
551 sainvd(9*itr-4) = dtemp(2)
552
553 do i =itr+1,np
554 sainvd(9*i-8) = sainvd(9*i-8)*dd
555 sainvd(9*i-4) = sainvd(9*i-4)*dd
556 sainvd(9*i ) = sainvd(9*i )*dd
557 enddo
558
559 rifd(9*itr) = dd !RIF UPDATE
560
561 js= inumfi1u(itr-1) + 1
562 je= inumfi1u(itr )
563 do j= js, je
564 rifu(9*j-2) = sainvd(9*fi1u(j)-8)
565 rifu(9*j-1) = sainvd(9*fi1u(j)-4)
566 rifu(9*j ) = sainvd(9*fi1u(j) )
567 enddo
568
569 !Update Z
570
571 do i= itr +1,np
572 js= inumfi1l(i-1) + 1
573 je= inumfi1l(i )
574 dd1=sainvd(9*i-8)
575 if(abs(dd1) > filter)then
576 do j= js, je
577 in = fi1l(j)
578 if (in > itr) exit
579 sainvl(9*j-8) = sainvl(9*j-8)-dd1*zz(3*in-2)
580 sainvl(9*j-7) = sainvl(9*j-7)-dd1*zz(3*in-1)
581 sainvl(9*j-6) = sainvl(9*j-6)-dd1*zz(3*in )
582 enddo
583 endif
584 dd2=sainvd(9*i-4)
585 if(abs(dd2) > filter)then
586 do j= js, je
587 in = fi1l(j)
588 if (in > itr) exit
589 sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
590 sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
591 sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
592 enddo
593 endif
594 dd3=sainvd(9*i )
595 if(abs(dd3) > filter)then
596 do j= js, je
597 in = fi1l(j)
598 if (in > itr) exit
599 sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
600 sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
601 sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
602 enddo
603 endif
604 enddo
605 enddo
606 deallocate(vv)
607 deallocate(zz)
608
609 end subroutine hecmw_rif_33
610
611 subroutine hecmw_rif_make_u_33(hecMAT)
612 implicit none
613 type (hecmwst_matrix) :: hecmat
614 integer(kind=kint) i,j,k,n,m,o
615 integer(kind=kint) is,ie,js,je
616
617 n = 1
618 do i= 1, hecmat%NP
619 is=inumfi1l(i-1) + 1
620 ie=inumfi1l(i )
621 flag1:do k= is, ie
622 m = fi1l(k)
623 js=inumfi1u(m-1) + 1
624 je=inumfi1u(m )
625 do j= js,je
626 o = fi1u(j)
627 if (o == i)then
628 rifl(9*n-8)=rifu(9*j-8)
629 rifl(9*n-7)=rifu(9*j-5)
630 rifl(9*n-6)=rifu(9*j-2)
631 rifl(9*n-5)=rifu(9*j-7)
632 rifl(9*n-4)=rifu(9*j-4)
633 rifl(9*n-3)=rifu(9*j-1)
634 rifl(9*n-2)=rifu(9*j-6)
635 rifl(9*n-1)=rifu(9*j-3)
636 rifl(9*n )=rifu(9*j )
637 n = n + 1
638 cycle flag1
639 endif
640 enddo
641 enddo flag1
642 enddo
643 end subroutine hecmw_rif_make_u_33
644
645 !C***
646 !C*** FORM_ILU1_33
647 !C*** form ILU(1) matrix
648 subroutine form_ilu0_rif_33(hecMAT)
649 implicit none
650 type(hecmwst_matrix) :: hecmat
651
652 allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
653 allocate (fi1l(hecmat%NPL), fi1u(hecmat%NPU))
654
655 inumfi1l = 0
656 inumfi1u = 0
657 fi1l = 0
658 fi1u = 0
659
660 inumfi1l(:) = hecmat%indexL(:)
661 inumfi1u(:) = hecmat%indexU(:)
662 fi1l(:) = hecmat%itemL(:)
663 fi1u(:) = hecmat%itemU(:)
664
665 npfiu = hecmat%NPU
666 npfil = hecmat%NPL
667
668 end subroutine form_ilu0_rif_33
669
670 !C***
671 !C*** FORM_ILU1_33
672 !C*** form ILU(1) matrix
673 subroutine form_ilu1_rif_33(hecMAT)
674 implicit none
675 type(hecmwst_matrix) :: hecmat
676
677 integer(kind=kint),allocatable :: iwsl(:), iwsu(:), iw1(:), iw2(:)
678 integer(kind=kint) :: nplf1,npuf1
679 integer(kind=kint) :: i,jj,kk,l,isk,iek,isj,iej
680 integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
681 integer(kind=kint) :: j,k,isl,isu
682 !C
683 !C +--------------+
684 !C | find fill-in |
685 !C +--------------+
686 !C===
687
688 !C
689 !C-- count fill-in
690 allocate (iw1(hecmat%NP) , iw2(hecmat%NP))
691 allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
692
693 inumfi1l= 0
694 inumfi1u= 0
695
696 nplf1= 0
697 npuf1= 0
698 do i= 2, hecmat%NP
699 icou= 0
700 iw1= 0
701 iw1(i)= 1
702 do l= indexl(i-1)+1, indexl(i)
703 iw1(iteml(l))= 1
704 enddo
705 do l= indexu(i-1)+1, indexu(i)
706 iw1(itemu(l))= 1
707 enddo
708
709 isk= indexl(i-1) + 1
710 iek= indexl(i)
711 do k= isk, iek
712 kk= iteml(k)
713 isj= indexu(kk-1) + 1
714 iej= indexu(kk )
715 do j= isj, iej
716 jj= itemu(j)
717 if (iw1(jj).eq.0 .and. jj.lt.i) then
718 inumfi1l(i)= inumfi1l(i)+1
719 iw1(jj)= 1
720 endif
721 if (iw1(jj).eq.0 .and. jj.gt.i) then
722 inumfi1u(i)= inumfi1u(i)+1
723 iw1(jj)= 1
724 endif
725 enddo
726 enddo
727 nplf1= nplf1 + inumfi1l(i)
728 npuf1= npuf1 + inumfi1u(i)
729 enddo
730
731 !C
732 !C-- specify fill-in
733 allocate (iwsl(0:hecmat%NP), iwsu(0:hecmat%NP))
734 allocate (fi1l(hecmat%NPL+nplf1), fi1u(hecmat%NPU+npuf1))
735
736 npfiu = hecmat%NPU+npuf1
737 npfil = hecmat%NPL+nplf1
738
739 fi1l= 0
740 fi1u= 0
741
742 iwsl= 0
743 iwsu= 0
744 do i= 1, hecmat%NP
745 iwsl(i)= indexl(i)-indexl(i-1) + inumfi1l(i) + iwsl(i-1)
746 iwsu(i)= indexu(i)-indexu(i-1) + inumfi1u(i) + iwsu(i-1)
747 enddo
748
749 do i= 2, hecmat%NP
750 icoul= 0
751 icouu= 0
752 inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
753 inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
754 icou= 0
755 iw1= 0
756 iw1(i)= 1
757 do l= indexl(i-1)+1, indexl(i)
758 iw1(iteml(l))= 1
759 enddo
760 do l= indexu(i-1)+1, indexu(i)
761 iw1(itemu(l))= 1
762 enddo
763
764 isk= indexl(i-1) + 1
765 iek= indexl(i)
766 do k= isk, iek
767 kk= iteml(k)
768 isj= indexu(kk-1) + 1
769 iej= indexu(kk )
770 do j= isj, iej
771 jj= itemu(j)
772 if (iw1(jj).eq.0 .and. jj.lt.i) then
773 icoul = icoul + 1
774 fi1l(icoul+iwsl(i-1)+indexl(i)-indexl(i-1))= jj
775 iw1(jj) = 1
776 endif
777 if (iw1(jj).eq.0 .and. jj.gt.i) then
778 icouu = icouu + 1
779 fi1u(icouu+iwsu(i-1)+indexu(i)-indexu(i-1))= jj
780 iw1(jj) = 1
781 endif
782 enddo
783 enddo
784 enddo
785
786 isl = 0
787 isu = 0
788 do i= 1, hecmat%NP
789 icoul1= indexl(i) - indexl(i-1)
790 icoul2= inumfi1l(i) - inumfi1l(i-1)
791 icoul3= icoul1 + icoul2
792 icouu1= indexu(i) - indexu(i-1)
793 icouu2= inumfi1u(i) - inumfi1u(i-1)
794 icouu3= icouu1 + icouu2
795 !C
796 !C-- LOWER part
797 icou0= 0
798 do k= indexl(i-1)+1, indexl(i)
799 icou0 = icou0 + 1
800 iw1(icou0)= iteml(k)
801 enddo
802
803 do k= inumfi1l(i-1)+1, inumfi1l(i)
804 icou0 = icou0 + 1
805 iw1(icou0)= fi1l(icou0+iwsl(i-1))
806 enddo
807
808 do k= 1, icoul3
809 iw2(k)= k
810 enddo
811 call rif_sort_33 (iw1, iw2, icoul3, hecmat%NP)
812
813 do k= 1, icoul3
814 fi1l(k+isl)= iw1(k)
815 enddo
816 !C
817 !C-- UPPER part
818 icou0= 0
819 do k= indexu(i-1)+1, indexu(i)
820 icou0 = icou0 + 1
821 iw1(icou0)= itemu(k)
822 enddo
823
824 do k= inumfi1u(i-1)+1, inumfi1u(i)
825 icou0 = icou0 + 1
826 iw1(icou0)= fi1u(icou0+iwsu(i-1))
827 enddo
828
829 do k= 1, icouu3
830 iw2(k)= k
831 enddo
832 call rif_sort_33 (iw1, iw2, icouu3, hecmat%NP)
833
834 do k= 1, icouu3
835 fi1u(k+isu)= iw1(k)
836 enddo
837
838 isl= isl + icoul3
839 isu= isu + icouu3
840 enddo
841
842 !C===
843 do i= 1, hecmat%NP
844 inumfi1l(i)= iwsl(i)
845 inumfi1u(i)= iwsu(i)
846 enddo
847
848 deallocate (iw1, iw2)
849 deallocate (iwsl, iwsu)
850 !C===
851 end subroutine form_ilu1_rif_33
852
853 !C
854 !C***
855 !C*** fill_in_S33_SORT
856 !C***
857 !C
858 subroutine rif_sort_33(STEM, INUM, N, NP)
859 use hecmw_util
860 implicit none
861 integer(kind=kint) :: n, np
862 integer(kind=kint), dimension(NP) :: stem
863 integer(kind=kint), dimension(NP) :: inum
864 integer(kind=kint), dimension(:), allocatable :: istack
865 integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
866
867 allocate (istack(-np:+np))
868
869 m = 100
870 nstack= np
871
872 jstack= 0
873 l = 1
874 ir = n
875
876 ip= 0
877 1 continue
878 ip= ip + 1
879
880 if (ir-l.lt.m) then
881 do j= l+1, ir
882 ss= stem(j)
883 ii= inum(j)
884
885 do i= j-1,1,-1
886 if (stem(i).le.ss) goto 2
887 stem(i+1)= stem(i)
888 inum(i+1)= inum(i)
889 end do
890 i= 0
891
892 2 continue
893 stem(i+1)= ss
894 inum(i+1)= ii
895 end do
896
897 if (jstack.eq.0) then
898 deallocate (istack)
899 return
900 endif
901
902 ir = istack(jstack)
903 l = istack(jstack-1)
904 jstack= jstack - 2
905 else
906
907 k= (l+ir) / 2
908 temp = stem(k)
909 stem(k) = stem(l+1)
910 stem(l+1)= temp
911
912 it = inum(k)
913 inum(k) = inum(l+1)
914 inum(l+1)= it
915
916 if (stem(l+1).gt.stem(ir)) then
917 temp = stem(l+1)
918 stem(l+1)= stem(ir)
919 stem(ir )= temp
920 it = inum(l+1)
921 inum(l+1)= inum(ir)
922 inum(ir )= it
923 endif
924
925 if (stem(l).gt.stem(ir)) then
926 temp = stem(l)
927 stem(l )= stem(ir)
928 stem(ir)= temp
929 it = inum(l)
930 inum(l )= inum(ir)
931 inum(ir)= it
932 endif
933
934 if (stem(l+1).gt.stem(l)) then
935 temp = stem(l+1)
936 stem(l+1)= stem(l)
937 stem(l )= temp
938 it = inum(l+1)
939 inum(l+1)= inum(l)
940 inum(l )= it
941 endif
942
943 i= l + 1
944 j= ir
945
946 ss= stem(l)
947 ii= inum(l)
948
949 3 continue
950 i= i + 1
951 if (stem(i).lt.ss) goto 3
952
953 4 continue
954 j= j - 1
955 if (stem(j).gt.ss) goto 4
956
957 if (j.lt.i) goto 5
958
959 temp = stem(i)
960 stem(i)= stem(j)
961 stem(j)= temp
962
963 it = inum(i)
964 inum(i)= inum(j)
965 inum(j)= it
966
967 goto 3
968
969 5 continue
970
971 stem(l)= stem(j)
972 stem(j)= ss
973 inum(l)= inum(j)
974 inum(j)= ii
975
976 jstack= jstack + 2
977
978 if (jstack.gt.nstack) then
979 write (*,*) 'NSTACK overflow'
980 stop
981 endif
982
983 if (ir-i+1.ge.j-1) then
984 istack(jstack )= ir
985 istack(jstack-1)= i
986 ir= j-1
987 else
988 istack(jstack )= j-1
989 istack(jstack-1)= l
990 l= i
991 endif
992
993 endif
994
995 goto 1
996
997 end subroutine rif_sort_33
998
1000 implicit none
1001
1002 if (associated(sainvd)) deallocate(sainvd)
1003 if (associated(sainvl)) deallocate(sainvl)
1004 if (associated(rifd)) deallocate(rifd)
1005 if (associated(rifu)) deallocate(rifu)
1006 if (associated(rifl)) deallocate(rifl)
1007 if (associated(inumfi1l)) deallocate(inumfi1l)
1008 if (associated(inumfi1u)) deallocate(inumfi1u)
1009 if (associated(fi1l)) deallocate(fi1l)
1010 if (associated(fi1u)) deallocate(fi1u)
1011 nullify(inumfi1l)
1012 nullify(inumfi1u)
1013 nullify(fi1l)
1014 nullify(fi1u)
1015 nullify(d)
1016 nullify(al)
1017 nullify(au)
1018 nullify(indexl)
1019 nullify(indexu)
1020 nullify(iteml)
1021 nullify(itemu)
1022
1023 end subroutine hecmw_precond_rif_33_clear
1024end module hecmw_precond_rif_33
integer(kind=kint) function, public hecmw_mat_get_precond(hecmat)
subroutine, public hecmw_precond_rif_33_setup(hecmat)
subroutine, public hecmw_precond_rif_33_apply(zp)
subroutine, public hecmw_precond_rif_33_clear()
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal