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