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