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