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