FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_precond_BILU_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
6!C
7!C***
8!C*** module hecmw_precond_BILU_33
9!C***
10!C
12 use hecmw_util
14
15 private
16
20
21 integer(kind=kint) :: N
22 real(kind=kreal), pointer :: dlu0(:) => null()
23 real(kind=kreal), pointer :: allu0(:) => null()
24 real(kind=kreal), pointer :: aulu0(:) => null()
25 integer(kind=kint), pointer :: inumFI1L(:) => null()
26 integer(kind=kint), pointer :: inumFI1U(:) => null()
27 integer(kind=kint), pointer :: FI1L(:) => null()
28 integer(kind=kint), pointer :: FI1U(:) => null()
29
30 logical, save :: INITIALIZED = .false.
31
32contains
33
34 subroutine hecmw_precond_bilu_33_setup(hecMAT)
35 implicit none
36 type(hecmwst_matrix), intent(inout) :: hecmat
37 integer(kind=kint ) :: np, npu, npl
38 integer(kind=kint ) :: precond
39 real (kind=kreal) :: sigma, sigma_diag
40
41 real(kind=kreal), pointer :: d(:)
42 real(kind=kreal), pointer :: al(:)
43 real(kind=kreal), pointer :: au(:)
44
45 integer(kind=kint ), pointer :: inl(:), inu(:)
46 integer(kind=kint ), pointer :: ial(:)
47 integer(kind=kint ), pointer :: iau(:)
48
49 if (initialized) then
50 if (hecmat%Iarray(98) == 1) then ! need symbolic and numerical setup
52 else if (hecmat%Iarray(97) == 1) then ! need numerical setup only
53 call hecmw_precond_bilu_33_clear() ! TEMPORARY
54 else
55 return
56 endif
57 endif
58
59 n = hecmat%N
60 np = hecmat%NP
61 npl = hecmat%NPL
62 npu = hecmat%NPU
63 d => hecmat%D
64 al => hecmat%AL
65 au => hecmat%AU
66 inl => hecmat%indexL
67 inu => hecmat%indexU
68 ial => hecmat%itemL
69 iau => hecmat%itemU
70 precond = hecmw_mat_get_precond(hecmat)
71 sigma = hecmw_mat_get_sigma(hecmat)
72 sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
73
74 if (precond.eq.10) call form_ilu0_33 &
75 & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
76 & sigma, sigma_diag)
77 if (precond.eq.11) call form_ilu1_33 &
78 & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
79 & sigma, sigma_diag)
80 if (precond.eq.12) call form_ilu2_33 &
81 & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
82 & sigma, sigma_diag)
83
84 initialized = .true.
85 hecmat%Iarray(98) = 0 ! symbolic setup done
86 hecmat%Iarray(97) = 0 ! numerical setup done
87
88 end subroutine hecmw_precond_bilu_33_setup
89
91 implicit none
92 real(kind=kreal), intent(inout) :: ww(:)
93 integer(kind=kint) :: i, j, isl, iel, isu, ieu, k
94 real(kind=kreal) :: sw1, sw2, sw3, x1, x2, x3
95 !C
96 !C-- FORWARD
97
98 do i= 1, n
99 sw1= ww(3*i-2)
100 sw2= ww(3*i-1)
101 sw3= ww(3*i )
102 isl= inumfi1l(i-1)+1
103 iel= inumfi1l(i)
104 do j= isl, iel
105 k= fi1l(j)
106 x1= ww(3*k-2)
107 x2= ww(3*k-1)
108 x3= ww(3*k )
109 sw1= sw1 - allu0(9*j-8)*x1-allu0(9*j-7)*x2-allu0(9*j-6)*x3
110 sw2= sw2 - allu0(9*j-5)*x1-allu0(9*j-4)*x2-allu0(9*j-3)*x3
111 sw3= sw3 - allu0(9*j-2)*x1-allu0(9*j-1)*x2-allu0(9*j )*x3
112 enddo
113
114 x1= sw1
115 x2= sw2
116 x3= sw3
117 x2= x2 - dlu0(9*i-5)*x1
118 x3= x3 - dlu0(9*i-2)*x1 - dlu0(9*i-1)*x2
119 x3= dlu0(9*i )* x3
120 x2= dlu0(9*i-4)*( x2 - dlu0(9*i-3)*x3 )
121 x1= dlu0(9*i-8)*( x1 - dlu0(9*i-6)*x3 - dlu0(9*i-7)*x2)
122 ww(3*i-2)= x1
123 ww(3*i-1)= x2
124 ww(3*i )= x3
125 enddo
126
127 !C
128 !C-- BACKWARD
129
130 do i= n, 1, -1
131 isu= inumfi1u(i-1) + 1
132 ieu= inumfi1u(i)
133 sw1= 0.d0
134 sw2= 0.d0
135 sw3= 0.d0
136 do j= ieu, isu, -1
137 k= fi1u(j)
138 x1= ww(3*k-2)
139 x2= ww(3*k-1)
140 x3= ww(3*k )
141 sw1= sw1 + aulu0(9*j-8)*x1+aulu0(9*j-7)*x2+aulu0(9*j-6)*x3
142 sw2= sw2 + aulu0(9*j-5)*x1+aulu0(9*j-4)*x2+aulu0(9*j-3)*x3
143 sw3= sw3 + aulu0(9*j-2)*x1+aulu0(9*j-1)*x2+aulu0(9*j )*x3
144 enddo
145 x1= sw1
146 x2= sw2
147 x3= sw3
148 x2= x2 - dlu0(9*i-5)*x1
149 x3= x3 - dlu0(9*i-2)*x1 - dlu0(9*i-1)*x2
150 x3= dlu0(9*i )* x3
151 x2= dlu0(9*i-4)*( x2 - dlu0(9*i-3)*x3 )
152 x1= dlu0(9*i-8)*( x1 - dlu0(9*i-6)*x3 - dlu0(9*i-7)*x2)
153 ww(3*i-2)= ww(3*i-2) - x1
154 ww(3*i-1)= ww(3*i-1) - x2
155 ww(3*i )= ww(3*i ) - x3
156 enddo
157 end subroutine hecmw_precond_bilu_33_apply
158
160 implicit none
161 if (associated(dlu0)) deallocate(dlu0)
162 if (associated(allu0)) deallocate(allu0)
163 if (associated(aulu0)) deallocate(aulu0)
164 if (associated(inumfi1l)) deallocate(inumfi1l)
165 if (associated(inumfi1u)) deallocate(inumfi1u)
166 if (associated(fi1l)) deallocate(fi1l)
167 if (associated(fi1u)) deallocate(fi1u)
168 nullify(dlu0)
169 nullify(allu0)
170 nullify(aulu0)
171 nullify(inumfi1l)
172 nullify(inumfi1u)
173 nullify(fi1l)
174 nullify(fi1u)
175 initialized = .false.
176 end subroutine hecmw_precond_bilu_33_clear
177
178 !C
179 !C***
180 !C*** FORM_ILU0_33
181 !C***
182 !C
183 !C form ILU(0) matrix
184 !C
185 subroutine form_ilu0_33 &
186 & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
187 & sigma, sigma_diag)
188 implicit none
189 integer(kind=kint ), intent(in):: n, np, npu, npl
190 real (kind=kreal), intent(in):: sigma, sigma_diag
191
192 real(kind=kreal), dimension(9*NPL), intent(in):: al
193 real(kind=kreal), dimension(9*NPU), intent(in):: au
194 real(kind=kreal), dimension(9*NP ), intent(in):: d
195
196 integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
197 integer(kind=kint ), dimension( NPL),intent(in) :: ial
198 integer(kind=kint ), dimension( NPU),intent(in) :: iau
199
200 integer(kind=kint), dimension(:), allocatable :: iw1, iw2
201 real (kind=kreal), dimension(3,3) :: rhs_aij, dkinv, aik, akj
202 integer(kind=kint) :: i,jj,ij0,kk
203 integer(kind=kint) :: j,k
204 allocate (iw1(np) , iw2(np))
205 allocate(dlu0(9*np), allu0(9*npl), aulu0(9*npu))
206 allocate(inumfi1l(0:np), inumfi1u(0:np), fi1l(npl), fi1u(npu))
207
208 do i=1,9*np
209 dlu0(i) = d(i)
210 end do
211 do i=1,9*npl
212 allu0(i) = al(i)
213 end do
214 do i=1,9*npu
215 aulu0(i) = au(i)
216 end do
217 do i=0,np
218 inumfi1l(i) = inl(i)
219 inumfi1u(i) = inu(i)
220 end do
221 do i=1,npl
222 fi1l(i) = ial(i)
223 end do
224 do i=1,npu
225 fi1u(i) = iau(i)
226 end do
227
228 !C
229 !C +----------------------+
230 !C | ILU(0) factorization |
231 !C +----------------------+
232 !C===
233 do i=1,np
234 dlu0(9*i-8)=dlu0(9*i-8)*sigma_diag
235 dlu0(9*i-4)=dlu0(9*i-4)*sigma_diag
236 dlu0(9*i )=dlu0(9*i )*sigma_diag
237 enddo
238
239 i = 1
240 call ilu1a33 (dkinv, &
241 dlu0(9*i-8), dlu0(9*i-7), dlu0(9*i-6), &
242 dlu0(9*i-5), dlu0(9*i-4), dlu0(9*i-3), &
243 dlu0(9*i-2), dlu0(9*i-1), dlu0(9*i ))
244 dlu0(9*i-8)= dkinv(1,1)
245 dlu0(9*i-7)= dkinv(1,2)
246 dlu0(9*i-6)= dkinv(1,3)
247 dlu0(9*i-5)= dkinv(2,1)
248 dlu0(9*i-4)= dkinv(2,2)
249 dlu0(9*i-3)= dkinv(2,3)
250 dlu0(9*i-2)= dkinv(3,1)
251 dlu0(9*i-1)= dkinv(3,2)
252 dlu0(9*i )= dkinv(3,3)
253
254 do i= 2, np
255 iw1= 0
256 iw2= 0
257
258 do k= inumfi1l(i-1)+1, inumfi1l(i)
259 iw1(fi1l(k))= k
260 enddo
261
262 do k= inumfi1u(i-1)+1, inumfi1u(i)
263 iw2(fi1u(k))= k
264 enddo
265
266 do kk= inl(i-1)+1, inl(i)
267 k= ial(kk)
268
269 dkinv(1,1)= dlu0(9*k-8)
270 dkinv(1,2)= dlu0(9*k-7)
271 dkinv(1,3)= dlu0(9*k-6)
272 dkinv(2,1)= dlu0(9*k-5)
273 dkinv(2,2)= dlu0(9*k-4)
274 dkinv(2,3)= dlu0(9*k-3)
275 dkinv(3,1)= dlu0(9*k-2)
276 dkinv(3,2)= dlu0(9*k-1)
277 dkinv(3,3)= dlu0(9*k )
278
279 aik(1,1)= allu0(9*kk-8)
280 aik(1,2)= allu0(9*kk-7)
281 aik(1,3)= allu0(9*kk-6)
282 aik(2,1)= allu0(9*kk-5)
283 aik(2,2)= allu0(9*kk-4)
284 aik(2,3)= allu0(9*kk-3)
285 aik(3,1)= allu0(9*kk-2)
286 aik(3,2)= allu0(9*kk-1)
287 aik(3,3)= allu0(9*kk )
288
289 do jj= inu(k-1)+1, inu(k)
290 j= iau(jj)
291 if (iw1(j).eq.0.and.iw2(j).eq.0) cycle
292
293 akj(1,1)= aulu0(9*jj-8)
294 akj(1,2)= aulu0(9*jj-7)
295 akj(1,3)= aulu0(9*jj-6)
296 akj(2,1)= aulu0(9*jj-5)
297 akj(2,2)= aulu0(9*jj-4)
298 akj(2,3)= aulu0(9*jj-3)
299 akj(3,1)= aulu0(9*jj-2)
300 akj(3,2)= aulu0(9*jj-1)
301 akj(3,3)= aulu0(9*jj )
302
303 call ilu1b33 (rhs_aij, dkinv, aik, akj)
304
305 if (j.eq.i) then
306 dlu0(9*i-8)= dlu0(9*i-8) - rhs_aij(1,1)
307 dlu0(9*i-7)= dlu0(9*i-7) - rhs_aij(1,2)
308 dlu0(9*i-6)= dlu0(9*i-6) - rhs_aij(1,3)
309 dlu0(9*i-5)= dlu0(9*i-5) - rhs_aij(2,1)
310 dlu0(9*i-4)= dlu0(9*i-4) - rhs_aij(2,2)
311 dlu0(9*i-3)= dlu0(9*i-3) - rhs_aij(2,3)
312 dlu0(9*i-2)= dlu0(9*i-2) - rhs_aij(3,1)
313 dlu0(9*i-1)= dlu0(9*i-1) - rhs_aij(3,2)
314 dlu0(9*i )= dlu0(9*i ) - rhs_aij(3,3)
315 endif
316
317 if (j.lt.i) then
318 ij0= iw1(j)
319 allu0(9*ij0-8)= allu0(9*ij0-8) - rhs_aij(1,1)
320 allu0(9*ij0-7)= allu0(9*ij0-7) - rhs_aij(1,2)
321 allu0(9*ij0-6)= allu0(9*ij0-6) - rhs_aij(1,3)
322 allu0(9*ij0-5)= allu0(9*ij0-5) - rhs_aij(2,1)
323 allu0(9*ij0-4)= allu0(9*ij0-4) - rhs_aij(2,2)
324 allu0(9*ij0-3)= allu0(9*ij0-3) - rhs_aij(2,3)
325 allu0(9*ij0-2)= allu0(9*ij0-2) - rhs_aij(3,1)
326 allu0(9*ij0-1)= allu0(9*ij0-1) - rhs_aij(3,2)
327 allu0(9*ij0 )= allu0(9*ij0 ) - rhs_aij(3,3)
328 endif
329
330 if (j.gt.i) then
331 ij0= iw2(j)
332 aulu0(9*ij0-8)= aulu0(9*ij0-8) - rhs_aij(1,1)
333 aulu0(9*ij0-7)= aulu0(9*ij0-7) - rhs_aij(1,2)
334 aulu0(9*ij0-6)= aulu0(9*ij0-6) - rhs_aij(1,3)
335 aulu0(9*ij0-5)= aulu0(9*ij0-5) - rhs_aij(2,1)
336 aulu0(9*ij0-4)= aulu0(9*ij0-4) - rhs_aij(2,2)
337 aulu0(9*ij0-3)= aulu0(9*ij0-3) - rhs_aij(2,3)
338 aulu0(9*ij0-2)= aulu0(9*ij0-2) - rhs_aij(3,1)
339 aulu0(9*ij0-1)= aulu0(9*ij0-1) - rhs_aij(3,2)
340 aulu0(9*ij0 )= aulu0(9*ij0 ) - rhs_aij(3,3)
341 endif
342
343 enddo
344 enddo
345
346 call ilu1a33 (dkinv, &
347 dlu0(9*i-8), dlu0(9*i-7), dlu0(9*i-6), &
348 dlu0(9*i-5), dlu0(9*i-4), dlu0(9*i-3), &
349 dlu0(9*i-2), dlu0(9*i-1), dlu0(9*i ))
350 dlu0(9*i-8)= dkinv(1,1)
351 dlu0(9*i-7)= dkinv(1,2)
352 dlu0(9*i-6)= dkinv(1,3)
353 dlu0(9*i-5)= dkinv(2,1)
354 dlu0(9*i-4)= dkinv(2,2)
355 dlu0(9*i-3)= dkinv(2,3)
356 dlu0(9*i-2)= dkinv(3,1)
357 dlu0(9*i-1)= dkinv(3,2)
358 dlu0(9*i )= dkinv(3,3)
359 enddo
360
361 deallocate (iw1, iw2)
362 end subroutine form_ilu0_33
363
364 !C
365 !C***
366 !C*** FORM_ILU1_33
367 !C***
368 !C
369 !C form ILU(1) matrix
370 !C
371 subroutine form_ilu1_33 &
372 & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
373 & sigma, sigma_diag)
374 implicit none
375 integer(kind=kint ), intent(in):: n, np, npu, npl
376 real (kind=kreal), intent(in):: sigma, sigma_diag
377
378 real(kind=kreal), dimension(9*NPL), intent(in):: al
379 real(kind=kreal), dimension(9*NPU), intent(in):: au
380 real(kind=kreal), dimension(9*NP ), intent(in):: d
381
382 integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
383 integer(kind=kint ), dimension( NPL),intent(in) :: ial
384 integer(kind=kint ), dimension( NPU),intent(in) :: iau
385
386 integer(kind=kint), dimension(:), allocatable :: iw1, iw2
387 integer(kind=kint), dimension(:), allocatable :: iwsl, iwsu
388 real (kind=kreal), dimension(3,3) :: rhs_aij, dkinv, aik, akj
389 integer(kind=kint) :: nplf1,npuf1
390 integer(kind=kint) :: i,jj,jj1,ij0,kk,ik,kk1,kk2,l,isk,iek,isj,iej
391 integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
392 integer(kind=kint) :: j,k,isl,isu
393 !C
394 !C +--------------+
395 !C | find fill-in |
396 !C +--------------+
397 !C===
398
399 !C
400 !C-- count fill-in
401 allocate (iw1(np) , iw2(np))
402 allocate (inumfi1l(0:np), inumfi1u(0:np))
403
404 inumfi1l= 0
405 inumfi1u= 0
406
407 nplf1= 0
408 npuf1= 0
409 do i= 2, np
410 icou= 0
411 iw1= 0
412 iw1(i)= 1
413 do l= inl(i-1)+1, inl(i)
414 iw1(ial(l))= 1
415 enddo
416 do l= inu(i-1)+1, inu(i)
417 iw1(iau(l))= 1
418 enddo
419
420 isk= inl(i-1) + 1
421 iek= inl(i)
422 do k= isk, iek
423 kk= ial(k)
424 isj= inu(kk-1) + 1
425 iej= inu(kk )
426 do j= isj, iej
427 jj= iau(j)
428 if (iw1(jj).eq.0 .and. jj.lt.i) then
429 inumfi1l(i)= inumfi1l(i)+1
430 iw1(jj)= 1
431 endif
432 if (iw1(jj).eq.0 .and. jj.gt.i) then
433 inumfi1u(i)= inumfi1u(i)+1
434 iw1(jj)= 1
435 endif
436 enddo
437 enddo
438 nplf1= nplf1 + inumfi1l(i)
439 npuf1= npuf1 + inumfi1u(i)
440 enddo
441
442 !C
443 !C-- specify fill-in
444 allocate (iwsl(0:np), iwsu(0:np))
445 allocate (fi1l(npl+nplf1), fi1u(npu+npuf1))
446 allocate (allu0(9*(npl+nplf1)), aulu0(9*(npu+npuf1)))
447
448 fi1l= 0
449 fi1u= 0
450
451 iwsl= 0
452 iwsu= 0
453 do i= 1, np
454 iwsl(i)= inl(i)-inl(i-1) + inumfi1l(i) + iwsl(i-1)
455 iwsu(i)= inu(i)-inu(i-1) + inumfi1u(i) + iwsu(i-1)
456 enddo
457
458 do i= 2, np
459 icoul= 0
460 icouu= 0
461 inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
462 inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
463 icou= 0
464 iw1= 0
465 iw1(i)= 1
466 do l= inl(i-1)+1, inl(i)
467 iw1(ial(l))= 1
468 enddo
469 do l= inu(i-1)+1, inu(i)
470 iw1(iau(l))= 1
471 enddo
472
473 isk= inl(i-1) + 1
474 iek= inl(i)
475 do k= isk, iek
476 kk= ial(k)
477 isj= inu(kk-1) + 1
478 iej= inu(kk )
479 do j= isj, iej
480 jj= iau(j)
481 if (iw1(jj).eq.0 .and. jj.lt.i) then
482 icoul = icoul + 1
483 fi1l(icoul+iwsl(i-1)+inl(i)-inl(i-1))= jj
484 iw1(jj) = 1
485 endif
486 if (iw1(jj).eq.0 .and. jj.gt.i) then
487 icouu = icouu + 1
488 fi1u(icouu+iwsu(i-1)+inu(i)-inu(i-1))= jj
489 iw1(jj) = 1
490 endif
491 enddo
492 enddo
493 enddo
494 !C===
495
496 !C
497 !C +-------------------------------------------------+
498 !C | SORT and RECONSTRUCT matrix considering fill-in |
499 !C +-------------------------------------------------+
500 !C===
501 allu0= 0.d0
502 aulu0= 0.d0
503 isl = 0
504 isu = 0
505 do i= 1, np
506 icoul1= inl(i) - inl(i-1)
507 icoul2= inumfi1l(i) - inumfi1l(i-1)
508 icoul3= icoul1 + icoul2
509 icouu1= inu(i) - inu(i-1)
510 icouu2= inumfi1u(i) - inumfi1u(i-1)
511 icouu3= icouu1 + icouu2
512 !C
513 !C-- LOWER part
514 icou0= 0
515 do k= inl(i-1)+1, inl(i)
516 icou0 = icou0 + 1
517 iw1(icou0)= ial(k)
518 enddo
519
520 do k= inumfi1l(i-1)+1, inumfi1l(i)
521 icou0 = icou0 + 1
522 iw1(icou0)= fi1l(icou0+iwsl(i-1))
523 enddo
524
525 do k= 1, icoul3
526 iw2(k)= k
527 enddo
528 call fill_in_s33_sort (iw1, iw2, icoul3, np)
529
530 do k= 1, icoul3
531 fi1l(k+isl)= iw1(k)
532 ik= iw2(k)
533 if (ik.le.inl(i)-inl(i-1)) then
534 kk1= 9*( k+isl)
535 kk2= 9*(ik+inl(i-1))
536 allu0(kk1-8)= al(kk2-8)
537 allu0(kk1-7)= al(kk2-7)
538 allu0(kk1-6)= al(kk2-6)
539 allu0(kk1-5)= al(kk2-5)
540 allu0(kk1-4)= al(kk2-4)
541 allu0(kk1-3)= al(kk2-3)
542 allu0(kk1-2)= al(kk2-2)
543 allu0(kk1-1)= al(kk2-1)
544 allu0(kk1 )= al(kk2 )
545 endif
546 enddo
547 !C
548 !C-- UPPER part
549 icou0= 0
550 do k= inu(i-1)+1, inu(i)
551 icou0 = icou0 + 1
552 iw1(icou0)= iau(k)
553 enddo
554
555 do k= inumfi1u(i-1)+1, inumfi1u(i)
556 icou0 = icou0 + 1
557 iw1(icou0)= fi1u(icou0+iwsu(i-1))
558 enddo
559
560 do k= 1, icouu3
561 iw2(k)= k
562 enddo
563 call fill_in_s33_sort (iw1, iw2, icouu3, np)
564
565 do k= 1, icouu3
566 fi1u(k+isu)= iw1(k)
567 ik= iw2(k)
568 if (ik.le.inu(i)-inu(i-1)) then
569 kk1= 9*( k+isu)
570 kk2= 9*(ik+inu(i-1))
571 aulu0(kk1-8)= au(kk2-8)
572 aulu0(kk1-7)= au(kk2-7)
573 aulu0(kk1-6)= au(kk2-6)
574 aulu0(kk1-5)= au(kk2-5)
575 aulu0(kk1-4)= au(kk2-4)
576 aulu0(kk1-3)= au(kk2-3)
577 aulu0(kk1-2)= au(kk2-2)
578 aulu0(kk1-1)= au(kk2-1)
579 aulu0(kk1 )= au(kk2 )
580 endif
581 enddo
582
583 isl= isl + icoul3
584 isu= isu + icouu3
585 enddo
586
587 !C===
588 do i= 1, np
589 inumfi1l(i)= iwsl(i)
590 inumfi1u(i)= iwsu(i)
591 enddo
592 deallocate (iwsl, iwsu)
593
594 !C
595 !C +----------------------+
596 !C | ILU(1) factorization |
597 !C +----------------------+
598 !C===
599 allocate (dlu0(9*np))
600 dlu0= d
601 do i=1,np
602 dlu0(9*i-8)=dlu0(9*i-8)*sigma_diag
603 dlu0(9*i-4)=dlu0(9*i-4)*sigma_diag
604 dlu0(9*i )=dlu0(9*i )*sigma_diag
605 enddo
606
607 i = 1
608 call ilu1a33 (dkinv, &
609 dlu0(9*i-8), dlu0(9*i-7), dlu0(9*i-6), &
610 dlu0(9*i-5), dlu0(9*i-4), dlu0(9*i-3), &
611 dlu0(9*i-2), dlu0(9*i-1), dlu0(9*i ))
612 dlu0(9*i-8)= dkinv(1,1)
613 dlu0(9*i-7)= dkinv(1,2)
614 dlu0(9*i-6)= dkinv(1,3)
615 dlu0(9*i-5)= dkinv(2,1)
616 dlu0(9*i-4)= dkinv(2,2)
617 dlu0(9*i-3)= dkinv(2,3)
618 dlu0(9*i-2)= dkinv(3,1)
619 dlu0(9*i-1)= dkinv(3,2)
620 dlu0(9*i )= dkinv(3,3)
621
622 do i= 2, np
623 iw1= 0
624 iw2= 0
625
626 do k= inumfi1l(i-1)+1, inumfi1l(i)
627 iw1(fi1l(k))= k
628 enddo
629
630 do k= inumfi1u(i-1)+1, inumfi1u(i)
631 iw2(fi1u(k))= k
632 enddo
633
634 do kk= inl(i-1)+1, inl(i)
635 k= ial(kk)
636
637 dkinv(1,1)= dlu0(9*k-8)
638 dkinv(1,2)= dlu0(9*k-7)
639 dkinv(1,3)= dlu0(9*k-6)
640 dkinv(2,1)= dlu0(9*k-5)
641 dkinv(2,2)= dlu0(9*k-4)
642 dkinv(2,3)= dlu0(9*k-3)
643 dkinv(3,1)= dlu0(9*k-2)
644 dkinv(3,2)= dlu0(9*k-1)
645 dkinv(3,3)= dlu0(9*k )
646
647 do kk1= inumfi1l(i-1)+1, inumfi1l(i)
648 if (k.eq.fi1l(kk1)) then
649 aik(1,1)= allu0(9*kk1-8)
650 aik(1,2)= allu0(9*kk1-7)
651 aik(1,3)= allu0(9*kk1-6)
652 aik(2,1)= allu0(9*kk1-5)
653 aik(2,2)= allu0(9*kk1-4)
654 aik(2,3)= allu0(9*kk1-3)
655 aik(3,1)= allu0(9*kk1-2)
656 aik(3,2)= allu0(9*kk1-1)
657 aik(3,3)= allu0(9*kk1 )
658 exit
659 endif
660 enddo
661
662 do jj= inu(k-1)+1, inu(k)
663 j= iau(jj)
664 do jj1= inumfi1u(k-1)+1, inumfi1u(k)
665 if (j.eq.fi1u(jj1)) then
666 akj(1,1)= aulu0(9*jj1-8)
667 akj(1,2)= aulu0(9*jj1-7)
668 akj(1,3)= aulu0(9*jj1-6)
669 akj(2,1)= aulu0(9*jj1-5)
670 akj(2,2)= aulu0(9*jj1-4)
671 akj(2,3)= aulu0(9*jj1-3)
672 akj(3,1)= aulu0(9*jj1-2)
673 akj(3,2)= aulu0(9*jj1-1)
674 akj(3,3)= aulu0(9*jj1 )
675 exit
676 endif
677 enddo
678
679 call ilu1b33 (rhs_aij, dkinv, aik, akj)
680
681 if (j.eq.i) then
682 dlu0(9*i-8)= dlu0(9*i-8) - rhs_aij(1,1)
683 dlu0(9*i-7)= dlu0(9*i-7) - rhs_aij(1,2)
684 dlu0(9*i-6)= dlu0(9*i-6) - rhs_aij(1,3)
685 dlu0(9*i-5)= dlu0(9*i-5) - rhs_aij(2,1)
686 dlu0(9*i-4)= dlu0(9*i-4) - rhs_aij(2,2)
687 dlu0(9*i-3)= dlu0(9*i-3) - rhs_aij(2,3)
688 dlu0(9*i-2)= dlu0(9*i-2) - rhs_aij(3,1)
689 dlu0(9*i-1)= dlu0(9*i-1) - rhs_aij(3,2)
690 dlu0(9*i )= dlu0(9*i ) - rhs_aij(3,3)
691 endif
692
693 if (j.lt.i) then
694 ij0= iw1(j)
695 allu0(9*ij0-8)= allu0(9*ij0-8) - rhs_aij(1,1)
696 allu0(9*ij0-7)= allu0(9*ij0-7) - rhs_aij(1,2)
697 allu0(9*ij0-6)= allu0(9*ij0-6) - rhs_aij(1,3)
698 allu0(9*ij0-5)= allu0(9*ij0-5) - rhs_aij(2,1)
699 allu0(9*ij0-4)= allu0(9*ij0-4) - rhs_aij(2,2)
700 allu0(9*ij0-3)= allu0(9*ij0-3) - rhs_aij(2,3)
701 allu0(9*ij0-2)= allu0(9*ij0-2) - rhs_aij(3,1)
702 allu0(9*ij0-1)= allu0(9*ij0-1) - rhs_aij(3,2)
703 allu0(9*ij0 )= allu0(9*ij0 ) - rhs_aij(3,3)
704 endif
705
706 if (j.gt.i) then
707 ij0= iw2(j)
708 aulu0(9*ij0-8)= aulu0(9*ij0-8) - rhs_aij(1,1)
709 aulu0(9*ij0-7)= aulu0(9*ij0-7) - rhs_aij(1,2)
710 aulu0(9*ij0-6)= aulu0(9*ij0-6) - rhs_aij(1,3)
711 aulu0(9*ij0-5)= aulu0(9*ij0-5) - rhs_aij(2,1)
712 aulu0(9*ij0-4)= aulu0(9*ij0-4) - rhs_aij(2,2)
713 aulu0(9*ij0-3)= aulu0(9*ij0-3) - rhs_aij(2,3)
714 aulu0(9*ij0-2)= aulu0(9*ij0-2) - rhs_aij(3,1)
715 aulu0(9*ij0-1)= aulu0(9*ij0-1) - rhs_aij(3,2)
716 aulu0(9*ij0 )= aulu0(9*ij0 ) - rhs_aij(3,3)
717 endif
718
719 enddo
720 enddo
721
722 call ilu1a33 (dkinv, &
723 dlu0(9*i-8), dlu0(9*i-7), dlu0(9*i-6), &
724 dlu0(9*i-5), dlu0(9*i-4), dlu0(9*i-3), &
725 dlu0(9*i-2), dlu0(9*i-1), dlu0(9*i ))
726 dlu0(9*i-8)= dkinv(1,1)
727 dlu0(9*i-7)= dkinv(1,2)
728 dlu0(9*i-6)= dkinv(1,3)
729 dlu0(9*i-5)= dkinv(2,1)
730 dlu0(9*i-4)= dkinv(2,2)
731 dlu0(9*i-3)= dkinv(2,3)
732 dlu0(9*i-2)= dkinv(3,1)
733 dlu0(9*i-1)= dkinv(3,2)
734 dlu0(9*i )= dkinv(3,3)
735 enddo
736
737 deallocate (iw1, iw2)
738 !C===
739 end subroutine form_ilu1_33
740
741 !C
742 !C***
743 !C*** FORM_ILU2_33
744 !C***
745 !C
746 !C form ILU(2) matrix
747 !C
748 subroutine form_ilu2_33 &
749 & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
750 & sigma, sigma_diag)
751 implicit none
752 integer(kind=kint ), intent(in):: n, np, npu, npl
753 real (kind=kreal), intent(in):: sigma, sigma_diag
754
755 real(kind=kreal), dimension(9*NPL), intent(in):: al
756 real(kind=kreal), dimension(9*NPU), intent(in):: au
757 real(kind=kreal), dimension(9*NP ), intent(in):: d
758
759 integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
760 integer(kind=kint ), dimension( NPL),intent(in) :: ial
761 integer(kind=kint ), dimension( NPU),intent(in) :: iau
762
763 integer(kind=kint), dimension(:), allocatable:: iw1 , iw2
764 integer(kind=kint), dimension(:), allocatable:: iwsl, iwsu
765 integer(kind=kint), dimension(:), allocatable:: iconfi1l, iconfi1u
766 integer(kind=kint), dimension(:), allocatable:: inumfi2l, inumfi2u
767 integer(kind=kint), dimension(:), allocatable:: fi2l, fi2u
768 real (kind=kreal), dimension(3,3) :: rhs_aij, dkinv, aik, akj
769 integer(kind=kint) :: nplf1,nplf2,npuf1,npuf2,ias,iconik,iconkj
770 integer(kind=kint) :: i,jj,ij0,kk,ik,kk1,kk2,l,isk,iek,isj,iej
771 integer(kind=kint) :: icou,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
772 integer(kind=kint) :: j,k,isl,isu
773
774 !C
775 !C +------------------+
776 !C | find fill-in (1) |
777 !C +------------------+
778 !C===
779
780 !C
781 !C-- count fill-in
782 allocate (iw1(np) , iw2(np))
783 allocate (inumfi2l(0:np), inumfi2u(0:np))
784
785 inumfi2l= 0
786 inumfi2u= 0
787
788 nplf1= 0
789 npuf1= 0
790 do i= 2, np
791 icou= 0
792 iw1= 0
793 iw1(i)= 1
794 do l= inl(i-1)+1, inl(i)
795 iw1(ial(l))= 1
796 enddo
797 do l= inu(i-1)+1, inu(i)
798 iw1(iau(l))= 1
799 enddo
800
801 isk= inl(i-1) + 1
802 iek= inl(i)
803 do k= isk, iek
804 kk= ial(k)
805 isj= inu(kk-1) + 1
806 iej= inu(kk )
807 do j= isj, iej
808 jj= iau(j)
809 if (iw1(jj).eq.0 .and. jj.lt.i) then
810 inumfi2l(i)= inumfi2l(i)+1
811 iw1(jj)= 1
812 endif
813 if (iw1(jj).eq.0 .and. jj.gt.i) then
814 inumfi2u(i)= inumfi2u(i)+1
815 iw1(jj)= 1
816 endif
817 enddo
818 enddo
819 nplf1= nplf1 + inumfi2l(i)
820 npuf1= npuf1 + inumfi2u(i)
821 enddo
822
823 !C
824 !C-- specify fill-in
825 allocate (iwsl(0:np), iwsu(0:np))
826 allocate (fi2l(nplf1), fi2u(npuf1))
827
828 fi2l= 0
829 fi2u= 0
830
831 do i= 2, np
832 icoul= 0
833 icouu= 0
834 inumfi2l(i)= inumfi2l(i-1) + inumfi2l(i)
835 inumfi2u(i)= inumfi2u(i-1) + inumfi2u(i)
836 icou= 0
837 iw1= 0
838 iw1(i)= 1
839 do l= inl(i-1)+1, inl(i)
840 iw1(ial(l))= 1
841 enddo
842 do l= inu(i-1)+1, inu(i)
843 iw1(iau(l))= 1
844 enddo
845
846 isk= inl(i-1) + 1
847 iek= inl(i)
848 do k= isk, iek
849 kk= ial(k)
850 isj= inu(kk-1) + 1
851 iej= inu(kk )
852 do j= isj, iej
853 jj= iau(j)
854 if (iw1(jj).eq.0 .and. jj.lt.i) then
855 icoul = icoul + 1
856 fi2l(icoul+inumfi2l(i-1))= jj
857 iw1(jj)= 1
858 endif
859 if (iw1(jj).eq.0 .and. jj.gt.i) then
860 icouu = icouu + 1
861 fi2u(icouu+inumfi2u(i-1))= jj
862 iw1(jj)= 1
863 endif
864 enddo
865 enddo
866 enddo
867 !C===
868
869 !C
870 !C +------------------+
871 !C | find fill-in (2) |
872 !C +------------------+
873 !C===
874 allocate (inumfi1l(0:np), inumfi1u(0:np))
875
876 nplf2= 0
877 npuf2= 0
878 inumfi1l= 0
879 inumfi1u= 0
880 !C
881 !C-- count fill-in
882 do i= 2, np
883 iw1= 0
884 iw1(i)= 1
885 do l= inl(i-1)+1, inl(i)
886 iw1(ial(l))= 2
887 enddo
888 do l= inu(i-1)+1, inu(i)
889 iw1(iau(l))= 2
890 enddo
891
892 do l= inumfi2l(i-1)+1, inumfi2l(i)
893 iw1(fi2l(l))= 1
894 enddo
895
896 do l= inumfi2u(i-1)+1, inumfi2u(i)
897 iw1(fi2u(l))= 1
898 enddo
899
900 isk= inl(i-1) + 1
901 iek= inl(i)
902 do k= isk, iek
903 kk= ial(k)
904 isj= inumfi2u(kk-1) + 1
905 iej= inumfi2u(kk)
906 do j= isj, iej
907 jj= fi2u(j)
908 if (iw1(jj).eq.0 .and. jj.lt.i) then
909 inumfi1l(i)= inumfi1l(i) + 1
910 iw1(jj)= 1
911 endif
912 if (iw1(jj).eq.0 .and. jj.gt.i) then
913 inumfi1u(i)= inumfi1u(i) + 1
914 iw1(jj)= 1
915 endif
916 enddo
917 enddo
918
919 isk= inumfi2l(i-1)+1
920 iek= inumfi2l(i)
921 do k= isk, iek
922 kk= fi2l(k)
923 isj= inu(kk-1) + 1
924 iej= inu(kk )
925 do j= isj, iej
926 jj= iau(j)
927 if (iw1(jj).eq.0 .and. jj.lt.i) then
928 inumfi1l(i)= inumfi1l(i) + 1
929 iw1(jj)= 1
930 endif
931 if (iw1(jj).eq.0 .and. jj.gt.i) then
932 inumfi1u(i)= inumfi1u(i) + 1
933 iw1(jj)= 1
934 endif
935 enddo
936 enddo
937 nplf2= nplf2 + inumfi1l(i)
938 npuf2= npuf2 + inumfi1u(i)
939 enddo
940
941 !C
942 !C-- specify fill-in
943 allocate (fi1l(npl+nplf1+nplf2))
944 allocate (fi1u(npu+npuf1+npuf2))
945
946 allocate (iconfi1l(npl+nplf1+nplf2))
947 allocate (iconfi1u(npu+npuf1+npuf2))
948
949 iwsl= 0
950 iwsu= 0
951 do i= 1, np
952 iwsl(i)= inl(i)-inl(i-1) + inumfi2l(i)-inumfi2l(i-1) + &
953 & inumfi1l(i) + iwsl(i-1)
954 iwsu(i)= inu(i)-inu(i-1) + inumfi2u(i)-inumfi2u(i-1) + &
955 & inumfi1u(i) + iwsu(i-1)
956 enddo
957
958 do i= 2, np
959 icoul= 0
960 icouu= 0
961 inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
962 inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
963 icou= 0
964 iw1= 0
965 iw1(i)= 1
966 do l= inl(i-1)+1, inl(i)
967 iw1(ial(l))= 1
968 enddo
969 do l= inu(i-1)+1, inu(i)
970 iw1(iau(l))= 1
971 enddo
972
973 do l= inumfi2l(i-1)+1, inumfi2l(i)
974 iw1(fi2l(l))= 1
975 enddo
976
977 do l= inumfi2u(i-1)+1, inumfi2u(i)
978 iw1(fi2u(l))= 1
979 enddo
980
981 isk= inl(i-1) + 1
982 iek= inl(i)
983 do k= isk, iek
984 kk= ial(k)
985 isj= inumfi2u(kk-1) + 1
986 iej= inumfi2u(kk )
987 do j= isj, iej
988 jj= fi2u(j)
989 if (iw1(jj).eq.0 .and. jj.lt.i) then
990 ias= inl(i)-inl(i-1)+inumfi2l(i)-inumfi2l(i-1)+iwsl(i-1)
991 icoul = icoul + 1
992 fi1l(icoul+ias)= jj
993 iw1(jj) = 1
994 endif
995 if (iw1(jj).eq.0 .and. jj.gt.i) then
996 ias= inu(i)-inu(i-1)+inumfi2u(i)-inumfi2u(i-1)+iwsu(i-1)
997 icouu = icouu + 1
998 fi1u(icouu+ias)= jj
999 iw1(jj) = 1
1000 endif
1001 enddo
1002 enddo
1003
1004 isk= inumfi2l(i-1) + 1
1005 iek= inumfi2l(i)
1006 do k= isk, iek
1007 kk= fi2l(k)
1008 isj= inu(kk-1) + 1
1009 iej= inu(kk )
1010 do j= isj, iej
1011 jj= iau(j)
1012 if (iw1(jj).eq.0 .and. jj.lt.i) then
1013 ias= inl(i)-inl(i-1)+inumfi2l(i)-inumfi2l(i-1)+iwsl(i-1)
1014 icoul = icoul + 1
1015 fi1l(icoul+ias)= jj
1016 iw1(jj) = 1
1017 endif
1018 if (iw1(jj).eq.0 .and. jj.gt.i) then
1019 ias= inu(i)-inu(i-1)+inumfi2u(i)-inumfi2u(i-1)+iwsu(i-1)
1020 icouu = icouu + 1
1021 fi1u(icouu+ias)= jj
1022 iw1(jj) = 1
1023 endif
1024 enddo
1025 enddo
1026 enddo
1027 !C===
1028
1029 !C
1030 !C +-------------------------------------------------+
1031 !C | SORT and RECONSTRUCT matrix considering fill-in |
1032 !C +-------------------------------------------------+
1033 !C===
1034 allocate (allu0(9*(npl+nplf1+nplf2)))
1035 allocate (aulu0(9*(npu+npuf1+npuf2)))
1036
1037 allu0= 0.d0
1038 aulu0= 0.d0
1039 isl = 0
1040 isu = 0
1041
1042 iconfi1l= 0
1043 iconfi1u= 0
1044
1045 do i= 1, np
1046 icoul1= inl(i) - inl(i-1)
1047 icoul2= inumfi2l(i) - inumfi2l(i-1) + icoul1
1048 icoul3= inumfi1l(i) - inumfi1l(i-1) + icoul2
1049
1050 icouu1= inu(i) - inu(i-1)
1051 icouu2= inumfi2u(i) - inumfi2u(i-1) + icouu1
1052 icouu3= inumfi1u(i) - inumfi1u(i-1) + icouu2
1053
1054 !C
1055 !C-- LOWER part
1056 icou= 0
1057 do k= inl(i-1)+1, inl(i)
1058 icou = icou + 1
1059 iw1(icou)= ial(k)
1060 enddo
1061
1062 icou= 0
1063 do k= inumfi2l(i-1)+1, inumfi2l(i)
1064 icou = icou + 1
1065 iw1(icou+icoul1)= fi2l(k)
1066 enddo
1067
1068 icou= 0
1069 do k= inumfi1l(i-1)+1, inumfi1l(i)
1070 icou = icou + 1
1071 iw1(icou+icoul2)= fi1l(icou+icoul2+isl)
1072 enddo
1073
1074 do k= 1, icoul3
1075 iw2(k)= k
1076 enddo
1077
1078 call fill_in_s33_sort (iw1, iw2, icoul3, np)
1079
1080 do k= 1, icoul3
1081 fi1l(k+isl)= iw1(k)
1082 ik= iw2(k)
1083 if (ik.le.inl(i)-inl(i-1)) then
1084 kk1= 9*( k+isl)
1085 kk2= 9*(ik+inl(i-1))
1086 allu0(kk1-8)= al(kk2-8)
1087 allu0(kk1-7)= al(kk2-7)
1088 allu0(kk1-6)= al(kk2-6)
1089 allu0(kk1-5)= al(kk2-5)
1090 allu0(kk1-4)= al(kk2-4)
1091 allu0(kk1-3)= al(kk2-3)
1092 allu0(kk1-2)= al(kk2-2)
1093 allu0(kk1-1)= al(kk2-1)
1094 allu0(kk1 )= al(kk2 )
1095 endif
1096 enddo
1097
1098 icou= 0
1099 do k= inl(i-1)+1, inl(i)
1100 icou = icou + 1
1101 iw1(icou)= 0
1102 enddo
1103
1104 icou= 0
1105 do k= inumfi2l(i-1)+1, inumfi2l(i)
1106 icou = icou + 1
1107 iw1(icou+icoul1)= 1
1108 enddo
1109
1110 icou= 0
1111 do k= inumfi1l(i-1)+1, inumfi1l(i)
1112 icou = icou + 1
1113 iw1(icou+icoul2)= 2
1114 enddo
1115
1116 do k= 1, icoul3
1117 iconfi1l(k+isl)= iw1(iw2(k))
1118 enddo
1119 !C
1120 !C-- UPPER part
1121 icou= 0
1122 do k= inu(i-1)+1, inu(i)
1123 icou = icou + 1
1124 iw1(icou)= iau(k)
1125 enddo
1126
1127 icou= 0
1128 do k= inumfi2u(i-1)+1, inumfi2u(i)
1129 icou = icou + 1
1130 iw1(icou+icouu1)= fi2u(k)
1131 enddo
1132
1133 icou= 0
1134 do k= inumfi1u(i-1)+1, inumfi1u(i)
1135 icou = icou + 1
1136 iw1(icou+icouu2)= fi1u(icou+icouu2+isu)
1137 enddo
1138
1139 do k= 1, icouu3
1140 iw2(k)= k
1141 enddo
1142 call fill_in_s33_sort (iw1, iw2, icouu3, np)
1143
1144 do k= 1, icouu3
1145 fi1u(k+isu)= iw1(k)
1146 ik= iw2(k)
1147 if (ik.le.inu(i)-inu(i-1)) then
1148 kk1= 9*( k+isu)
1149 kk2= 9*(ik+inu(i-1))
1150 aulu0(kk1-8)= au(kk2-8)
1151 aulu0(kk1-7)= au(kk2-7)
1152 aulu0(kk1-6)= au(kk2-6)
1153 aulu0(kk1-5)= au(kk2-5)
1154 aulu0(kk1-4)= au(kk2-4)
1155 aulu0(kk1-3)= au(kk2-3)
1156 aulu0(kk1-2)= au(kk2-2)
1157 aulu0(kk1-1)= au(kk2-1)
1158 aulu0(kk1 )= au(kk2 )
1159 endif
1160 enddo
1161
1162 icou= 0
1163 do k= inu(i-1)+1, inu(i)
1164 icou = icou + 1
1165 iw1(icou)= 0
1166 enddo
1167
1168 icou= 0
1169 do k= inumfi2u(i-1)+1, inumfi2u(i)
1170 icou = icou + 1
1171 iw1(icou+icouu1)= 1
1172 enddo
1173
1174 icou= 0
1175 do k= inumfi1u(i-1)+1, inumfi1u(i)
1176 icou = icou + 1
1177 iw1(icou+icouu2)= 2
1178 enddo
1179
1180 do k= 1, icouu3
1181 iconfi1u(k+isu)= iw1(iw2(k))
1182 enddo
1183
1184 isl= isl + icoul3
1185 isu= isu + icouu3
1186 enddo
1187 !C===
1188 do i= 1, np
1189 inumfi1l(i)= iwsl(i)
1190 inumfi1u(i)= iwsu(i)
1191 enddo
1192
1193 deallocate (iwsl, iwsu)
1194 deallocate (inumfi2l, inumfi2u)
1195 deallocate ( fi2l, fi2u)
1196
1197 !C
1198 !C +----------------------+
1199 !C | ILU(2) factorization |
1200 !C +----------------------+
1201 !C===
1202 allocate (dlu0(9*np))
1203 dlu0= d
1204 do i=1,np
1205 dlu0(9*i-8)=dlu0(9*i-8)*sigma_diag
1206 dlu0(9*i-4)=dlu0(9*i-4)*sigma_diag
1207 dlu0(9*i )=dlu0(9*i )*sigma_diag
1208 enddo
1209
1210 i = 1
1211 call ilu1a33 (dkinv, &
1212 dlu0(9*i-8), dlu0(9*i-7), dlu0(9*i-6), &
1213 dlu0(9*i-5), dlu0(9*i-4), dlu0(9*i-3), &
1214 dlu0(9*i-2), dlu0(9*i-1), dlu0(9*i ))
1215 dlu0(9*i-8)= dkinv(1,1)
1216 dlu0(9*i-7)= dkinv(1,2)
1217 dlu0(9*i-6)= dkinv(1,3)
1218 dlu0(9*i-5)= dkinv(2,1)
1219 dlu0(9*i-4)= dkinv(2,2)
1220 dlu0(9*i-3)= dkinv(2,3)
1221 dlu0(9*i-2)= dkinv(3,1)
1222 dlu0(9*i-1)= dkinv(3,2)
1223 dlu0(9*i )= dkinv(3,3)
1224
1225 do i= 2, np
1226 iw1= 0
1227 iw2= 0
1228
1229 do k= inumfi1l(i-1)+1, inumfi1l(i)
1230 iw1(fi1l(k))= k
1231 enddo
1232
1233 do k= inumfi1u(i-1)+1, inumfi1u(i)
1234 iw2(fi1u(k))= k
1235 enddo
1236
1237 do kk= inumfi1l(i-1)+1, inumfi1l(i)
1238 k= fi1l(kk)
1239 iconik= iconfi1l(kk)
1240
1241 dkinv(1,1)= dlu0(9*k-8)
1242 dkinv(1,2)= dlu0(9*k-7)
1243 dkinv(1,3)= dlu0(9*k-6)
1244 dkinv(2,1)= dlu0(9*k-5)
1245 dkinv(2,2)= dlu0(9*k-4)
1246 dkinv(2,3)= dlu0(9*k-3)
1247 dkinv(3,1)= dlu0(9*k-2)
1248 dkinv(3,2)= dlu0(9*k-1)
1249 dkinv(3,3)= dlu0(9*k )
1250
1251 aik(1,1)= allu0(9*kk-8)
1252 aik(1,2)= allu0(9*kk-7)
1253 aik(1,3)= allu0(9*kk-6)
1254 aik(2,1)= allu0(9*kk-5)
1255 aik(2,2)= allu0(9*kk-4)
1256 aik(2,3)= allu0(9*kk-3)
1257 aik(3,1)= allu0(9*kk-2)
1258 aik(3,2)= allu0(9*kk-1)
1259 aik(3,3)= allu0(9*kk )
1260
1261 do jj= inumfi1u(k-1)+1, inumfi1u(k)
1262 j= fi1u(jj)
1263 iconkj= iconfi1u(jj)
1264
1265 if ((iconik+iconkj).lt.2) then
1266 akj(1,1)= aulu0(9*jj-8)
1267 akj(1,2)= aulu0(9*jj-7)
1268 akj(1,3)= aulu0(9*jj-6)
1269 akj(2,1)= aulu0(9*jj-5)
1270 akj(2,2)= aulu0(9*jj-4)
1271 akj(2,3)= aulu0(9*jj-3)
1272 akj(3,1)= aulu0(9*jj-2)
1273 akj(3,2)= aulu0(9*jj-1)
1274 akj(3,3)= aulu0(9*jj )
1275
1276 call ilu1b33 (rhs_aij, dkinv, aik, akj)
1277
1278 if (j.eq.i) then
1279 dlu0(9*i-8)= dlu0(9*i-8) - rhs_aij(1,1)
1280 dlu0(9*i-7)= dlu0(9*i-7) - rhs_aij(1,2)
1281 dlu0(9*i-6)= dlu0(9*i-6) - rhs_aij(1,3)
1282 dlu0(9*i-5)= dlu0(9*i-5) - rhs_aij(2,1)
1283 dlu0(9*i-4)= dlu0(9*i-4) - rhs_aij(2,2)
1284 dlu0(9*i-3)= dlu0(9*i-3) - rhs_aij(2,3)
1285 dlu0(9*i-2)= dlu0(9*i-2) - rhs_aij(3,1)
1286 dlu0(9*i-1)= dlu0(9*i-1) - rhs_aij(3,2)
1287 dlu0(9*i )= dlu0(9*i ) - rhs_aij(3,3)
1288 endif
1289
1290 if (j.lt.i) then
1291 ij0= iw1(j)
1292 allu0(9*ij0-8)= allu0(9*ij0-8) - rhs_aij(1,1)
1293 allu0(9*ij0-7)= allu0(9*ij0-7) - rhs_aij(1,2)
1294 allu0(9*ij0-6)= allu0(9*ij0-6) - rhs_aij(1,3)
1295 allu0(9*ij0-5)= allu0(9*ij0-5) - rhs_aij(2,1)
1296 allu0(9*ij0-4)= allu0(9*ij0-4) - rhs_aij(2,2)
1297 allu0(9*ij0-3)= allu0(9*ij0-3) - rhs_aij(2,3)
1298 allu0(9*ij0-2)= allu0(9*ij0-2) - rhs_aij(3,1)
1299 allu0(9*ij0-1)= allu0(9*ij0-1) - rhs_aij(3,2)
1300 allu0(9*ij0 )= allu0(9*ij0 ) - rhs_aij(3,3)
1301 endif
1302
1303 if (j.gt.i) then
1304 ij0= iw2(j)
1305 aulu0(9*ij0-8)= aulu0(9*ij0-8) - rhs_aij(1,1)
1306 aulu0(9*ij0-7)= aulu0(9*ij0-7) - rhs_aij(1,2)
1307 aulu0(9*ij0-6)= aulu0(9*ij0-6) - rhs_aij(1,3)
1308 aulu0(9*ij0-5)= aulu0(9*ij0-5) - rhs_aij(2,1)
1309 aulu0(9*ij0-4)= aulu0(9*ij0-4) - rhs_aij(2,2)
1310 aulu0(9*ij0-3)= aulu0(9*ij0-3) - rhs_aij(2,3)
1311 aulu0(9*ij0-2)= aulu0(9*ij0-2) - rhs_aij(3,1)
1312 aulu0(9*ij0-1)= aulu0(9*ij0-1) - rhs_aij(3,2)
1313 aulu0(9*ij0 )= aulu0(9*ij0 ) - rhs_aij(3,3)
1314 endif
1315 endif
1316 enddo
1317 enddo
1318
1319 call ilu1a33 (dkinv, &
1320 dlu0(9*i-8), dlu0(9*i-7), dlu0(9*i-6), &
1321 dlu0(9*i-5), dlu0(9*i-4), dlu0(9*i-3), &
1322 dlu0(9*i-2), dlu0(9*i-1), dlu0(9*i ))
1323 dlu0(9*i-8)= dkinv(1,1)
1324 dlu0(9*i-7)= dkinv(1,2)
1325 dlu0(9*i-6)= dkinv(1,3)
1326 dlu0(9*i-5)= dkinv(2,1)
1327 dlu0(9*i-4)= dkinv(2,2)
1328 dlu0(9*i-3)= dkinv(2,3)
1329 dlu0(9*i-2)= dkinv(3,1)
1330 dlu0(9*i-1)= dkinv(3,2)
1331 dlu0(9*i )= dkinv(3,3)
1332 enddo
1333
1334 deallocate (iw1, iw2)
1335 deallocate (iconfi1l, iconfi1u)
1336 !C===
1337 end subroutine form_ilu2_33
1338
1339
1340 !C
1341 !C***
1342 !C*** fill_in_S33_SORT
1343 !C***
1344 !C
1345 subroutine fill_in_s33_sort (STEM, INUM, N, NP)
1346 use hecmw_util
1347 implicit none
1348 integer(kind=kint) :: n, np
1349 integer(kind=kint), dimension(NP) :: stem
1350 integer(kind=kint), dimension(NP) :: inum
1351 integer(kind=kint), dimension(:), allocatable :: istack
1352 integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
1353
1354 allocate (istack(-np:+np))
1355
1356 m = 100
1357 nstack= np
1358
1359 jstack= 0
1360 l = 1
1361 ir = n
1362
1363 ip= 0
1364 1 continue
1365 ip= ip + 1
1366
1367 if (ir-l.lt.m) then
1368 do j= l+1, ir
1369 ss= stem(j)
1370 ii= inum(j)
1371
1372 do i= j-1,1,-1
1373 if (stem(i).le.ss) goto 2
1374 stem(i+1)= stem(i)
1375 inum(i+1)= inum(i)
1376 end do
1377 i= 0
1378
1379 2 continue
1380 stem(i+1)= ss
1381 inum(i+1)= ii
1382 end do
1383
1384 if (jstack.eq.0) then
1385 deallocate (istack)
1386 return
1387 endif
1388
1389 ir = istack(jstack)
1390 l = istack(jstack-1)
1391 jstack= jstack - 2
1392 else
1393
1394 k= (l+ir) / 2
1395 temp = stem(k)
1396 stem(k) = stem(l+1)
1397 stem(l+1)= temp
1398
1399 it = inum(k)
1400 inum(k) = inum(l+1)
1401 inum(l+1)= it
1402
1403 if (stem(l+1).gt.stem(ir)) then
1404 temp = stem(l+1)
1405 stem(l+1)= stem(ir)
1406 stem(ir )= temp
1407 it = inum(l+1)
1408 inum(l+1)= inum(ir)
1409 inum(ir )= it
1410 endif
1411
1412 if (stem(l).gt.stem(ir)) then
1413 temp = stem(l)
1414 stem(l )= stem(ir)
1415 stem(ir)= temp
1416 it = inum(l)
1417 inum(l )= inum(ir)
1418 inum(ir)= it
1419 endif
1420
1421 if (stem(l+1).gt.stem(l)) then
1422 temp = stem(l+1)
1423 stem(l+1)= stem(l)
1424 stem(l )= temp
1425 it = inum(l+1)
1426 inum(l+1)= inum(l)
1427 inum(l )= it
1428 endif
1429
1430 i= l + 1
1431 j= ir
1432
1433 ss= stem(l)
1434 ii= inum(l)
1435
1436 3 continue
1437 i= i + 1
1438 if (stem(i).lt.ss) goto 3
1439
1440 4 continue
1441 j= j - 1
1442 if (stem(j).gt.ss) goto 4
1443
1444 if (j.lt.i) goto 5
1445
1446 temp = stem(i)
1447 stem(i)= stem(j)
1448 stem(j)= temp
1449
1450 it = inum(i)
1451 inum(i)= inum(j)
1452 inum(j)= it
1453
1454 goto 3
1455
1456 5 continue
1457
1458 stem(l)= stem(j)
1459 stem(j)= ss
1460 inum(l)= inum(j)
1461 inum(j)= ii
1462
1463 jstack= jstack + 2
1464
1465 if (jstack.gt.nstack) then
1466 write (*,*) 'NSTACK overflow'
1467 stop
1468 endif
1469
1470 if (ir-i+1.ge.j-1) then
1471 istack(jstack )= ir
1472 istack(jstack-1)= i
1473 ir= j-1
1474 else
1475 istack(jstack )= j-1
1476 istack(jstack-1)= l
1477 l= i
1478 endif
1479
1480 endif
1481
1482 goto 1
1483
1484 end subroutine fill_in_s33_sort
1485
1486 !C
1487 !C***
1488 !C*** ILU1a33
1489 !C***
1490 !C
1491 !C computes LU factorization of 3*3 Diagonal Block
1492 !C
1493 subroutine ilu1a33 (ALU, D11,D12,D13,D21,D22,D23,D31,D32,D33)
1494 use hecmw_util
1495 implicit none
1496 real(kind=kreal) :: alu(3,3), pw(3)
1497 real(kind=kreal) :: d11,d12,d13,d21,d22,d23,d31,d32,d33
1498 integer(kind=kint) :: i,j,k
1499
1500 alu(1,1)= d11
1501 alu(1,2)= d12
1502 alu(1,3)= d13
1503 alu(2,1)= d21
1504 alu(2,2)= d22
1505 alu(2,3)= d23
1506 alu(3,1)= d31
1507 alu(3,2)= d32
1508 alu(3,3)= d33
1509
1510 do k= 1, 3
1511 if (alu(k,k) == 0.d0) then
1512 !write(*,*) ALU(1:3,1:3)
1513 stop 'ERROR: Divide by zero in ILU setup'
1514 endif
1515 alu(k,k)= 1.d0/alu(k,k)
1516 do i= k+1, 3
1517 alu(i,k)= alu(i,k) * alu(k,k)
1518 do j= k+1, 3
1519 pw(j)= alu(i,j) - alu(i,k)*alu(k,j)
1520 enddo
1521 do j= k+1, 3
1522 alu(i,j)= pw(j)
1523 enddo
1524 enddo
1525 enddo
1526
1527 return
1528 end subroutine ilu1a33
1529
1530 !C
1531 !C***
1532 !C*** ILU1b33
1533 !C***
1534 !C
1535 !C computes L_ik * D_k_INV * U_kj at ILU factorization
1536 !C for 3*3 Block Type Matrix
1537 !C
1538 subroutine ilu1b33 (RHS_Aij, DkINV, Aik, Akj)
1539 use hecmw_util
1540 implicit none
1541 real(kind=kreal) :: rhs_aij(3,3), dkinv(3,3), aik(3,3), akj(3,3)
1542 real(kind=kreal) :: x1,x2,x3
1543
1544 !C
1545 !C-- 1st Col.
1546 x1= akj(1,1)
1547 x2= akj(2,1)
1548 x3= akj(3,1)
1549
1550 x2= x2 - dkinv(2,1)*x1
1551 x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
1552
1553 x3= dkinv(3,3)* x3
1554 x2= dkinv(2,2)*( x2 - dkinv(2,3)*x3 )
1555 x1= dkinv(1,1)*( x1 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
1556
1557 rhs_aij(1,1)= aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3
1558 rhs_aij(2,1)= aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3
1559 rhs_aij(3,1)= aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3
1560
1561 !C
1562 !C-- 2nd Col.
1563 x1= akj(1,2)
1564 x2= akj(2,2)
1565 x3= akj(3,2)
1566
1567 x2= x2 - dkinv(2,1)*x1
1568 x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
1569
1570 x3= dkinv(3,3)* x3
1571 x2= dkinv(2,2)*( x2 - dkinv(2,3)*x3 )
1572 x1= dkinv(1,1)*( x1 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
1573
1574 rhs_aij(1,2)= aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3
1575 rhs_aij(2,2)= aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3
1576 rhs_aij(3,2)= aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3
1577
1578 !C
1579 !C-- 3rd Col.
1580 x1= akj(1,3)
1581 x2= akj(2,3)
1582 x3= akj(3,3)
1583
1584 x2= x2 - dkinv(2,1)*x1
1585 x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
1586
1587 x3= dkinv(3,3)* x3
1588 x2= dkinv(2,2)*( x2 - dkinv(2,3)*x3 )
1589 x1= dkinv(1,1)*( x1 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
1590
1591 rhs_aij(1,3)= aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3
1592 rhs_aij(2,3)= aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3
1593 rhs_aij(3,3)= aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3
1594
1595 return
1596 end subroutine ilu1b33
1597
1598end module hecmw_precond_bilu_33
real(kind=kreal) function, public hecmw_mat_get_sigma(hecmat)
integer(kind=kint) function, public hecmw_mat_get_precond(hecmat)
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecmat)
subroutine, public hecmw_precond_bilu_33_clear()
subroutine, public hecmw_precond_bilu_33_apply(ww)
subroutine, public hecmw_precond_bilu_33_setup(hecmat)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal