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