FrontISTR  5.7.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, 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 = 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 
170  subroutine hecmw_precond_bilu_nn_clear()
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 
1568 end module hecmw_precond_bilu_nn
hecmw_precond_bilu_nn::ilu1ann
subroutine ilu1ann(ALU, D, NDOF)
Definition: hecmw_precond_BILU_nn.f90:1501
hecmw_matrix_misc::hecmw_mat_get_precond
integer(kind=kint) function, public hecmw_mat_get_precond(hecMAT)
Definition: hecmw_matrix_misc.f90:321
hecmw_precond_bilu_nn::ilu1bnn
subroutine ilu1bnn(RHS_Aij, DkINV, Aik, Akj, NDOF)
Definition: hecmw_precond_BILU_nn.f90:1539
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_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_matrix_misc::hecmw_mat_get_sigma_diag
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
Definition: hecmw_matrix_misc.f90:694
hecmw_matrix_misc
Definition: hecmw_matrix_misc.f90:6
hecmw_precond_bilu_nn::ilu1a33
subroutine ilu1a33(ALU, D11, D12, D13, D21, D22, D23, D31, D32, D33)
Definition: hecmw_precond_BILU_nn.f90:1465
hecmw_precond_bilu_nn::hecmw_precond_bilu_nn_apply
subroutine, public hecmw_precond_bilu_nn_apply(WW, NDOF)
Definition: hecmw_precond_BILU_nn.f90:91
hecmw_precond_bilu_nn::hecmw_precond_bilu_nn_clear
subroutine, public hecmw_precond_bilu_nn_clear()
Definition: hecmw_precond_BILU_nn.f90:171
hecmw_precond_bilu_nn
Definition: hecmw_precond_BILU_nn.f90:11
hecmw_matrix_misc::hecmw_mat_get_sigma
real(kind=kreal) function, public hecmw_mat_get_sigma(hecMAT)
Definition: hecmw_matrix_misc.f90:714
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444