FrontISTR  6.0.0
Large-scale structural analysis program with finit element method
hecmw_precond_BILU_66.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_66
9 !C***
10 !C
11 module hecmw_precond_bilu_66
12  use hecmw_util
13  use hecmw_matrix_misc
14 
15  private
16 
17  public:: hecmw_precond_bilu_66_setup
18  public:: hecmw_precond_bilu_66_apply
19  public:: hecmw_precond_bilu_66_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  real(kind=kreal), pointer :: alu(:) => null()
30 
31 contains
32 
33  subroutine hecmw_precond_bilu_66_setup(hecMAT)
34  implicit none
35  type(hecmwst_matrix), intent(in) :: hecmat
36  integer(kind=kint ) :: np, npu, npl
37  integer(kind=kint ) :: precond
38  real (kind=kreal) :: sigma_diag
39 
40  real(kind=kreal), pointer :: d(:)
41  real(kind=kreal), pointer :: al(:)
42  real(kind=kreal), pointer :: au(:)
43 
44  integer(kind=kint ), pointer :: inl(:), inu(:)
45  integer(kind=kint ), pointer :: ial(:)
46  integer(kind=kint ), pointer :: iau(:)
47 
48  real (kind=kreal):: alutmp(6,6)
49  integer(kind=kint ):: ip
50 
51  n = hecmat%N
52  np = hecmat%NP
53  npl = hecmat%NPL
54  npu = hecmat%NPU
55  d => hecmat%D
56  al => hecmat%AL
57  au => hecmat%AU
58  inl => hecmat%indexL
59  inu => hecmat%indexU
60  ial => hecmat%itemL
61  iau => hecmat%itemU
62  precond = hecmw_mat_get_precond(hecmat)
63  sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
64 
65  if (precond.eq.10) call form_ilu0_66 &
66  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
67  & sigma_diag)
68  if (precond.eq.11) call form_ilu1_66 &
69  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
70  & sigma_diag)
71  if (precond.eq.12) call form_ilu2_66 &
72  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
73  & sigma_diag)
74 
75  allocate(alu(36*np))
76 
77  do ip= 1, n
78  call ilu1a66 (alutmp, &
79  & dlu0(36*ip-35), dlu0(36*ip-34), dlu0(36*ip-33), dlu0(36*ip-32),&
80  & dlu0(36*ip-31), dlu0(36*ip-30), dlu0(36*ip-29), dlu0(36*ip-28),&
81  & dlu0(36*ip-27), dlu0(36*ip-26), dlu0(36*ip-25), dlu0(36*ip-24),&
82  & dlu0(36*ip-23), dlu0(36*ip-22), dlu0(36*ip-21), dlu0(36*ip-20),&
83  & dlu0(36*ip-19), dlu0(36*ip-18), dlu0(36*ip-17), dlu0(36*ip-16),&
84  & dlu0(36*ip-15), dlu0(36*ip-14), dlu0(36*ip-13), dlu0(36*ip-12),&
85  & dlu0(36*ip-11), dlu0(36*ip-10), dlu0(36*ip-9 ), dlu0(36*ip-8 ),&
86  & dlu0(36*ip-7), dlu0(36*ip-6), dlu0(36*ip-5), dlu0(36*ip-4),&
87  & dlu0(36*ip-3), dlu0(36*ip-2), dlu0(36*ip-1), dlu0(36*ip ))
88  alu(36*ip-35)= alutmp(1,1)
89  alu(36*ip-34)= alutmp(1,2)
90  alu(36*ip-33)= alutmp(1,3)
91  alu(36*ip-32)= alutmp(1,4)
92  alu(36*ip-31)= alutmp(1,5)
93  alu(36*ip-30)= alutmp(1,6)
94  alu(36*ip-29)= alutmp(2,1)
95  alu(36*ip-28)= alutmp(2,2)
96  alu(36*ip-27)= alutmp(2,3)
97  alu(36*ip-26)= alutmp(2,4)
98  alu(36*ip-25)= alutmp(2,5)
99  alu(36*ip-24)= alutmp(2,6)
100  alu(36*ip-23)= alutmp(3,1)
101  alu(36*ip-22)= alutmp(3,2)
102  alu(36*ip-21)= alutmp(3,3)
103  alu(36*ip-20)= alutmp(3,4)
104  alu(36*ip-19)= alutmp(3,5)
105  alu(36*ip-18)= alutmp(3,6)
106  alu(36*ip-17)= alutmp(4,1)
107  alu(36*ip-16)= alutmp(4,2)
108  alu(36*ip-15)= alutmp(4,3)
109  alu(36*ip-14)= alutmp(4,4)
110  alu(36*ip-13)= alutmp(4,5)
111  alu(36*ip-12)= alutmp(4,6)
112  alu(36*ip-11)= alutmp(5,1)
113  alu(36*ip-10)= alutmp(5,2)
114  alu(36*ip-9 )= alutmp(5,3)
115  alu(36*ip-8 )= alutmp(5,4)
116  alu(36*ip-7 )= alutmp(5,5)
117  alu(36*ip-6 )= alutmp(5,6)
118  alu(36*ip-5 )= alutmp(6,1)
119  alu(36*ip-4 )= alutmp(6,2)
120  alu(36*ip-3 )= alutmp(6,3)
121  alu(36*ip-2 )= alutmp(6,4)
122  alu(36*ip-1 )= alutmp(6,5)
123  alu(36*ip )= alutmp(6,6)
124  enddo
125  end subroutine hecmw_precond_bilu_66_setup
126 
127  subroutine hecmw_precond_bilu_66_apply(WW)
128  implicit none
129  real(kind=kreal), intent(inout) :: ww(:)
130  integer(kind=kint) :: i, j, isl, iel, isu, ieu, k
131  real(kind=kreal) :: sw1, sw2, sw3, x1, x2, x3
132  real(kind=kreal) :: sw4, sw5, sw6, x4, x5, x6
133  !C
134  !C-- FORWARD
135 
136  do i= 1, n
137  sw1= ww(6*i-5)
138  sw2= ww(6*i-4)
139  sw3= ww(6*i-3)
140  sw4= ww(6*i-2)
141  sw5= ww(6*i-1)
142  sw6= ww(6*i )
143  isl= inumfi1l(i-1)+1
144  iel= inumfi1l(i)
145  do j= isl, iel
146  k= fi1l(j)
147  x1= ww(6*k-5)
148  x2= ww(6*k-4)
149  x3= ww(6*k-3)
150  x4= ww(6*k-2)
151  x5= ww(6*k-1)
152  x6= ww(6*k )
153  sw1= sw1 -allu0(36*j-35)*x1-allu0(36*j-34)*x2-allu0(36*j-33)*x3-allu0(36*j-32)*x4-allu0(36*j-31)*x5-allu0(36*j-30)*x6
154  sw2= sw2 -allu0(36*j-29)*x1-allu0(36*j-28)*x2-allu0(36*j-27)*x3-allu0(36*j-26)*x4-allu0(36*j-25)*x5-allu0(36*j-24)*x6
155  sw3= sw3 -allu0(36*j-23)*x1-allu0(36*j-22)*x2-allu0(36*j-21)*x3-allu0(36*j-20)*x4-allu0(36*j-19)*x5-allu0(36*j-18)*x6
156  sw4= sw4 -allu0(36*j-17)*x1-allu0(36*j-16)*x2-allu0(36*j-15)*x3-allu0(36*j-14)*x4-allu0(36*j-13)*x5-allu0(36*j-12)*x6
157  sw5= sw5 -allu0(36*j-11)*x1-allu0(36*j-10)*x2-allu0(36*j-9 )*x3-allu0(36*j-8 )*x4-allu0(36*j-7 )*x5-allu0(36*j-6 )*x6
158  sw6= sw6 -allu0(36*j-5 )*x1-allu0(36*j-4 )*x2-allu0(36*j-3 )*x3-allu0(36*j-2 )*x4-allu0(36*j-1 )*x5-allu0(36*j )*x6
159  enddo
160 
161  x1= sw1
162  x2= sw2
163  x3= sw3
164  x4= sw4
165  x5= sw5
166  x6= sw6
167  x2= x2 -alu(36*i-29)*x1
168  x3= x3 -alu(36*i-23)*x1 -alu(36*i-22)*x2
169  x4= x4 -alu(36*i-17)*x1 -alu(36*i-16)*x2 -alu(36*i-15)*x3
170  x5= x5 -alu(36*i-11)*x1 -alu(36*i-10)*x2 -alu(36*i-9 )*x3 -alu(36*i-8)*x4
171  x6= x6 -alu(36*i-5 )*x1 -alu(36*i-4 )*x2 -alu(36*i-3 )*x3 -alu(36*i-2)*x4 -alu(36*i-1)*x5
172  x6= alu(36*i )* x6
173  x5= alu(36*i-7 )*( x5 -alu(36*i-6 )*x6 )
174  x4= alu(36*i-14)*( x4 -alu(36*i-12)*x6 -alu(36*i-13)*x5)
175  x3= alu(36*i-21)*( x3 -alu(36*i-18)*x6 -alu(36*i-19)*x5 -alu(36*i-20)*x4)
176  x2= alu(36*i-28)*( x2 -alu(36*i-24)*x6 -alu(36*i-25)*x5 -alu(36*i-26)*x4 -alu(36*i-27)*x3)
177  x1= alu(36*i-35)*( x1 -alu(36*i-30)*x6 -alu(36*i-31)*x5 -alu(36*i-32)*x4 -alu(36*i-33)*x3 -alu(36*i-34)*x2)
178  ww(6*i-5)= x1
179  ww(6*i-4)= x2
180  ww(6*i-3)= x3
181  ww(6*i-2)= x4
182  ww(6*i-1)= x5
183  ww(6*i )= x6
184  enddo
185 
186  !C
187  !C-- BACKWARD
188 
189  do i= n, 1, -1
190  isu= inumfi1u(i-1) + 1
191  ieu= inumfi1u(i)
192  sw1= 0.d0
193  sw2= 0.d0
194  sw3= 0.d0
195  sw4= 0.d0
196  sw5= 0.d0
197  sw6= 0.d0
198  do j= ieu, isu, -1
199  k= fi1u(j)
200  x1= ww(6*k-5)
201  x2= ww(6*k-4)
202  x3= ww(6*k-3)
203  x4= ww(6*k-2)
204  x5= ww(6*k-1)
205  x6= ww(6*k )
206  sw1= sw1 +aulu0(36*j-35)*x1+aulu0(36*j-34)*x2+aulu0(36*j-33)*x3+aulu0(36*j-32)*x4+aulu0(36*j-31)*x5+aulu0(36*j-30)*x6
207  sw2= sw2 +aulu0(36*j-29)*x1+aulu0(36*j-28)*x2+aulu0(36*j-27)*x3+aulu0(36*j-26)*x4+aulu0(36*j-25)*x5+aulu0(36*j-24)*x6
208  sw3= sw3 +aulu0(36*j-23)*x1+aulu0(36*j-22)*x2+aulu0(36*j-21)*x3+aulu0(36*j-20)*x4+aulu0(36*j-19)*x5+aulu0(36*j-18)*x6
209  sw4= sw4 +aulu0(36*j-17)*x1+aulu0(36*j-16)*x2+aulu0(36*j-15)*x3+aulu0(36*j-14)*x4+aulu0(36*j-13)*x5+aulu0(36*j-12)*x6
210  sw5= sw5 +aulu0(36*j-11)*x1+aulu0(36*j-10)*x2+aulu0(36*j-9 )*x3+aulu0(36*j-8 )*x4+aulu0(36*j-7 )*x5+aulu0(36*j-6 )*x6
211  sw6= sw6 +aulu0(36*j-5 )*x1+aulu0(36*j-4 )*x2+aulu0(36*j-3 )*x3+aulu0(36*j-2 )*x4+aulu0(36*j-1 )*x5+aulu0(36*j )*x6
212  enddo
213  x1= sw1
214  x2= sw2
215  x3= sw3
216  x4= sw4
217  x5= sw5
218  x6= sw6
219  x2= x2 -alu(36*i-29)*x1
220  x3= x3 -alu(36*i-23)*x1 -alu(36*i-22)*x2
221  x4= x4 -alu(36*i-17)*x1 -alu(36*i-16)*x2 -alu(36*i-15)*x3
222  x5= x5 -alu(36*i-11)*x1 -alu(36*i-10)*x2 -alu(36*i-9 )*x3 -alu(36*i-8)*x4
223  x6= x6 -alu(36*i-5 )*x1 -alu(36*i-4 )*x2 -alu(36*i-3 )*x3 -alu(36*i-2)*x4 -alu(36*i-1)*x5
224  x6= alu(36*i )* x6
225  x5= alu(36*i-7 )*( x5 -alu(36*i-6 )*x6 )
226  x4= alu(36*i-14)*( x4 -alu(36*i-12)*x6 -alu(36*i-13)*x5)
227  x3= alu(36*i-21)*( x3 -alu(36*i-18)*x6 -alu(36*i-19)*x5 -alu(36*i-20)*x4)
228  x2= alu(36*i-28)*( x2 -alu(36*i-24)*x6 -alu(36*i-25)*x5 -alu(36*i-26)*x4 -alu(36*i-27)*x3)
229  x1= alu(36*i-35)*( x1 -alu(36*i-30)*x6 -alu(36*i-31)*x5 -alu(36*i-32)*x4 -alu(36*i-33)*x3 -alu(36*i-34)*x2)
230  ww(6*i-5)= ww(6*i-5) -x1
231  ww(6*i-4)= ww(6*i-4) -x2
232  ww(6*i-3)= ww(6*i-3) -x3
233  ww(6*i-2)= ww(6*i-2) -x4
234  ww(6*i-1)= ww(6*i-1) -x5
235  ww(6*i )= ww(6*i ) -x6
236  enddo
237  end subroutine hecmw_precond_bilu_66_apply
238 
239  subroutine hecmw_precond_bilu_66_clear()
240  implicit none
241  if (associated(dlu0)) deallocate(dlu0)
242  if (associated(allu0)) deallocate(allu0)
243  if (associated(aulu0)) deallocate(aulu0)
244  if (associated(inumfi1l)) deallocate(inumfi1l)
245  if (associated(inumfi1u)) deallocate(inumfi1u)
246  if (associated(fi1l)) deallocate(fi1l)
247  if (associated(fi1u)) deallocate(fi1u)
248  if (associated(alu)) deallocate(alu)
249  nullify(dlu0)
250  nullify(allu0)
251  nullify(aulu0)
252  nullify(inumfi1l)
253  nullify(inumfi1u)
254  nullify(fi1l)
255  nullify(fi1u)
256  nullify(alu)
257  end subroutine hecmw_precond_bilu_66_clear
258 
259  !C
260  !C***
261  !C*** FORM_ILU1_66
262  !C***
263  !C
264  !C form ILU(1) matrix
265  !C
266  subroutine form_ilu0_66 &
267  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
268  & sigma_diag)
269  implicit none
270  integer(kind=kint ), intent(in):: n, np, npu, npl
271  real (kind=kreal), intent(in):: sigma_diag
272 
273  real(kind=kreal), dimension(36*NPL), intent(in):: al
274  real(kind=kreal), dimension(36*NPU), intent(in):: au
275  real(kind=kreal), dimension(36*NP ), intent(in):: d
276 
277  integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
278  integer(kind=kint ), dimension( NPL),intent(in) :: ial
279  integer(kind=kint ), dimension( NPU),intent(in) :: iau
280 
281  integer(kind=kint), dimension(:), allocatable :: iw1, iw2
282  real (kind=kreal), dimension(6,6) :: rhs_aij, dkinv, aik, akj
283  real(kind=kreal) :: d11,d12,d13,d14,d15,d16
284  real(kind=kreal) :: d21,d22,d23,d24,d25,d26
285  real(kind=kreal) :: d31,d32,d33,d34,d35,d36
286  real(kind=kreal) :: d41,d42,d43,d44,d45,d46
287  real(kind=kreal) :: d51,d52,d53,d54,d55,d56
288  real(kind=kreal) :: d61,d62,d63,d64,d65,d66
289  integer(kind=kint) :: i,jj,ij0,kk
290  integer(kind=kint) :: j,k
291  allocate (iw1(np) , iw2(np))
292  allocate(dlu0(36*np), allu0(36*npl), aulu0(36*npu))
293  allocate(inumfi1l(0:np), inumfi1u(0:np), fi1l(npl), fi1u(npu))
294 
295  do i=1,36*np
296  dlu0(i) = d(i)
297  end do
298  do i=1,36*npl
299  allu0(i) = al(i)
300  end do
301  do i=1,36*npu
302  aulu0(i) = au(i)
303  end do
304  do i=0,np
305  inumfi1l(i) = inl(i)
306  inumfi1u(i) = inu(i)
307  end do
308  do i=1,npl
309  fi1l(i) = ial(i)
310  end do
311  do i=1,npu
312  fi1u(i) = iau(i)
313  end do
314 
315  !C
316  !C +----------------------+
317  !C | ILU(0) factorization |
318  !C +----------------------+
319  !C===
320  do i=1,np
321  dlu0(36*i-35)=dlu0(36*i-35)*sigma_diag
322  dlu0(36*i-28)=dlu0(36*i-28)*sigma_diag
323  dlu0(36*i-21)=dlu0(36*i-21)*sigma_diag
324  dlu0(36*i-14)=dlu0(36*i-14)*sigma_diag
325  dlu0(36*i-7 )=dlu0(36*i-7 )*sigma_diag
326  dlu0(36*i )=dlu0(36*i )*sigma_diag
327  enddo
328 
329  do i= 2, np
330  iw1= 0
331  iw2= 0
332 
333  do k= inumfi1l(i-1)+1, inumfi1l(i)
334  iw1(fi1l(k))= k
335  enddo
336 
337  do k= inumfi1u(i-1)+1, inumfi1u(i)
338  iw2(fi1u(k))= k
339  enddo
340 
341  do kk= inl(i-1)+1, inl(i)
342  k= ial(kk)
343  d11= dlu0(36*k-35)
344  d12= dlu0(36*k-34)
345  d13= dlu0(36*k-33)
346  d14= dlu0(36*k-32)
347  d15= dlu0(36*k-31)
348  d16= dlu0(36*k-30)
349  d21= dlu0(36*k-29)
350  d22= dlu0(36*k-28)
351  d23= dlu0(36*k-27)
352  d24= dlu0(36*k-26)
353  d25= dlu0(36*k-25)
354  d26= dlu0(36*k-24)
355  d31= dlu0(36*k-23)
356  d32= dlu0(36*k-22)
357  d33= dlu0(36*k-21)
358  d34= dlu0(36*k-20)
359  d35= dlu0(36*k-19)
360  d36= dlu0(36*k-18)
361  d41= dlu0(36*k-17)
362  d42= dlu0(36*k-16)
363  d43= dlu0(36*k-15)
364  d44= dlu0(36*k-14)
365  d45= dlu0(36*k-13)
366  d46= dlu0(36*k-12)
367  d51= dlu0(36*k-11)
368  d52= dlu0(36*k-10)
369  d53= dlu0(36*k-9 )
370  d54= dlu0(36*k-8 )
371  d55= dlu0(36*k-7 )
372  d56= dlu0(36*k-6 )
373  d61= dlu0(36*k-5 )
374  d62= dlu0(36*k-4 )
375  d63= dlu0(36*k-3 )
376  d64= dlu0(36*k-2 )
377  d65= dlu0(36*k-1 )
378  d66= dlu0(36*k )
379 
380  call ilu1a66 (dkinv,d11,d12,d13,d14,d15,d16,d21,d22,d23,d24,d25,d26, &
381  & d31,d32,d33,d34,d35,d36,d41,d42,d43,d44,d45,d46,d51,d52,d53,d54,d55,d56, &
382  & d61,d62,d63,d64,d65,d66)
383 
384  aik(1,1)= allu0(36*kk-35)
385  aik(1,2)= allu0(36*kk-34)
386  aik(1,3)= allu0(36*kk-33)
387  aik(1,4)= allu0(36*kk-32)
388  aik(1,5)= allu0(36*kk-31)
389  aik(1,6)= allu0(36*kk-30)
390  aik(2,1)= allu0(36*kk-29)
391  aik(2,2)= allu0(36*kk-28)
392  aik(2,3)= allu0(36*kk-27)
393  aik(2,4)= allu0(36*kk-26)
394  aik(2,5)= allu0(36*kk-25)
395  aik(2,6)= allu0(36*kk-24)
396  aik(3,1)= allu0(36*kk-23)
397  aik(3,2)= allu0(36*kk-22)
398  aik(3,3)= allu0(36*kk-21)
399  aik(3,4)= allu0(36*kk-20)
400  aik(3,5)= allu0(36*kk-19)
401  aik(3,6)= allu0(36*kk-18)
402  aik(4,1)= allu0(36*kk-17)
403  aik(4,2)= allu0(36*kk-16)
404  aik(4,3)= allu0(36*kk-15)
405  aik(4,4)= allu0(36*kk-14)
406  aik(4,5)= allu0(36*kk-13)
407  aik(4,6)= allu0(36*kk-12)
408  aik(5,1)= allu0(36*kk-11)
409  aik(5,2)= allu0(36*kk-10)
410  aik(5,3)= allu0(36*kk-9)
411  aik(5,4)= allu0(36*kk-8)
412  aik(5,5)= allu0(36*kk-7)
413  aik(5,6)= allu0(36*kk-6)
414  aik(6,1)= allu0(36*kk-5)
415  aik(6,2)= allu0(36*kk-4)
416  aik(6,3)= allu0(36*kk-3)
417  aik(6,4)= allu0(36*kk-2)
418  aik(6,5)= allu0(36*kk-1)
419  aik(6,6)= allu0(36*kk )
420 
421  do jj= inu(k-1)+1, inu(k)
422  j= iau(jj)
423  if (iw1(j).eq.0.and.iw2(j).eq.0) cycle
424 
425  akj(1,1)= aulu0(36*jj-35)
426  akj(1,2)= aulu0(36*jj-34)
427  akj(1,3)= aulu0(36*jj-33)
428  akj(1,4)= aulu0(36*jj-32)
429  akj(1,5)= aulu0(36*jj-31)
430  akj(1,6)= aulu0(36*jj-30)
431  akj(2,1)= aulu0(36*jj-29)
432  akj(2,2)= aulu0(36*jj-28)
433  akj(2,3)= aulu0(36*jj-27)
434  akj(2,4)= aulu0(36*jj-26)
435  akj(2,5)= aulu0(36*jj-25)
436  akj(2,6)= aulu0(36*jj-24)
437  akj(3,1)= aulu0(36*jj-23)
438  akj(3,2)= aulu0(36*jj-22)
439  akj(3,3)= aulu0(36*jj-21)
440  akj(3,4)= aulu0(36*jj-20)
441  akj(3,5)= aulu0(36*jj-19)
442  akj(3,6)= aulu0(36*jj-18)
443  akj(4,1)= aulu0(36*jj-17)
444  akj(4,2)= aulu0(36*jj-16)
445  akj(4,3)= aulu0(36*jj-15)
446  akj(4,4)= aulu0(36*jj-14)
447  akj(4,5)= aulu0(36*jj-13)
448  akj(4,6)= aulu0(36*jj-12)
449  akj(5,1)= aulu0(36*jj-11)
450  akj(5,2)= aulu0(36*jj-10)
451  akj(5,3)= aulu0(36*jj-9)
452  akj(5,4)= aulu0(36*jj-8)
453  akj(5,5)= aulu0(36*jj-7)
454  akj(5,6)= aulu0(36*jj-6)
455  akj(6,1)= aulu0(36*jj-5)
456  akj(6,2)= aulu0(36*jj-4)
457  akj(6,3)= aulu0(36*jj-3)
458  akj(6,4)= aulu0(36*jj-2)
459  akj(6,5)= aulu0(36*jj-1)
460  akj(6,6)= aulu0(36*jj )
461 
462  call ilu1b66 (rhs_aij, dkinv, aik, akj)
463 
464  if (j.eq.i) then
465  dlu0(36*i-35)= dlu0(36*i-35) - rhs_aij(1,1)
466  dlu0(36*i-34)= dlu0(36*i-34) - rhs_aij(1,2)
467  dlu0(36*i-33)= dlu0(36*i-33) - rhs_aij(1,3)
468  dlu0(36*i-32)= dlu0(36*i-32) - rhs_aij(1,4)
469  dlu0(36*i-31)= dlu0(36*i-31) - rhs_aij(1,5)
470  dlu0(36*i-30)= dlu0(36*i-30) - rhs_aij(1,6)
471  dlu0(36*i-29)= dlu0(36*i-29) - rhs_aij(2,1)
472  dlu0(36*i-28)= dlu0(36*i-28) - rhs_aij(2,2)
473  dlu0(36*i-27)= dlu0(36*i-27) - rhs_aij(2,3)
474  dlu0(36*i-26)= dlu0(36*i-26) - rhs_aij(2,4)
475  dlu0(36*i-25)= dlu0(36*i-25) - rhs_aij(2,5)
476  dlu0(36*i-24)= dlu0(36*i-24) - rhs_aij(2,6)
477  dlu0(36*i-23)= dlu0(36*i-23) - rhs_aij(3,1)
478  dlu0(36*i-22)= dlu0(36*i-22) - rhs_aij(3,2)
479  dlu0(36*i-21)= dlu0(36*i-21) - rhs_aij(3,3)
480  dlu0(36*i-20)= dlu0(36*i-20) - rhs_aij(3,4)
481  dlu0(36*i-19)= dlu0(36*i-19) - rhs_aij(3,5)
482  dlu0(36*i-18)= dlu0(36*i-18) - rhs_aij(3,6)
483  dlu0(36*i-17)= dlu0(36*i-17) - rhs_aij(4,1)
484  dlu0(36*i-16)= dlu0(36*i-16) - rhs_aij(4,2)
485  dlu0(36*i-15)= dlu0(36*i-15) - rhs_aij(4,3)
486  dlu0(36*i-14)= dlu0(36*i-14) - rhs_aij(4,4)
487  dlu0(36*i-13)= dlu0(36*i-13) - rhs_aij(4,5)
488  dlu0(36*i-12)= dlu0(36*i-12) - rhs_aij(4,6)
489  dlu0(36*i-11)= dlu0(36*i-11) - rhs_aij(5,1)
490  dlu0(36*i-10)= dlu0(36*i-10) - rhs_aij(5,2)
491  dlu0(36*i-9 )= dlu0(36*i-9 ) - rhs_aij(5,3)
492  dlu0(36*i-8 )= dlu0(36*i-8 ) - rhs_aij(5,4)
493  dlu0(36*i-7 )= dlu0(36*i-7 ) - rhs_aij(5,5)
494  dlu0(36*i-6 )= dlu0(36*i-6 ) - rhs_aij(5,6)
495  dlu0(36*i-5 )= dlu0(36*i-5 ) - rhs_aij(6,1)
496  dlu0(36*i-4 )= dlu0(36*i-4 ) - rhs_aij(6,2)
497  dlu0(36*i-3 )= dlu0(36*i-3 ) - rhs_aij(6,3)
498  dlu0(36*i-2 )= dlu0(36*i-2 ) - rhs_aij(6,4)
499  dlu0(36*i-1 )= dlu0(36*i-1 ) - rhs_aij(6,5)
500  dlu0(36*i )= dlu0(36*i ) - rhs_aij(6,6)
501  endif
502 
503  if (j.lt.i) then
504  ij0= iw1(j)
505  allu0(36*ij0-35)= allu0(36*ij0-35) - rhs_aij(1,1)
506  allu0(36*ij0-34)= allu0(36*ij0-34) - rhs_aij(1,2)
507  allu0(36*ij0-33)= allu0(36*ij0-33) - rhs_aij(1,3)
508  allu0(36*ij0-32)= allu0(36*ij0-32) - rhs_aij(1,4)
509  allu0(36*ij0-31)= allu0(36*ij0-31) - rhs_aij(1,5)
510  allu0(36*ij0-30)= allu0(36*ij0-30) - rhs_aij(1,6)
511  allu0(36*ij0-29)= allu0(36*ij0-29) - rhs_aij(2,1)
512  allu0(36*ij0-28)= allu0(36*ij0-28) - rhs_aij(2,2)
513  allu0(36*ij0-27)= allu0(36*ij0-27) - rhs_aij(2,3)
514  allu0(36*ij0-26)= allu0(36*ij0-26) - rhs_aij(2,4)
515  allu0(36*ij0-25)= allu0(36*ij0-25) - rhs_aij(2,5)
516  allu0(36*ij0-24)= allu0(36*ij0-24) - rhs_aij(2,6)
517  allu0(36*ij0-23)= allu0(36*ij0-23) - rhs_aij(3,1)
518  allu0(36*ij0-22)= allu0(36*ij0-22) - rhs_aij(3,2)
519  allu0(36*ij0-21)= allu0(36*ij0-21) - rhs_aij(3,3)
520  allu0(36*ij0-20)= allu0(36*ij0-20) - rhs_aij(3,4)
521  allu0(36*ij0-19)= allu0(36*ij0-19) - rhs_aij(3,5)
522  allu0(36*ij0-18)= allu0(36*ij0-18) - rhs_aij(3,6)
523  allu0(36*ij0-17)= allu0(36*ij0-17) - rhs_aij(4,1)
524  allu0(36*ij0-16)= allu0(36*ij0-16) - rhs_aij(4,2)
525  allu0(36*ij0-15)= allu0(36*ij0-15) - rhs_aij(4,3)
526  allu0(36*ij0-14)= allu0(36*ij0-14) - rhs_aij(4,4)
527  allu0(36*ij0-13)= allu0(36*ij0-13) - rhs_aij(4,5)
528  allu0(36*ij0-12)= allu0(36*ij0-12) - rhs_aij(4,6)
529  allu0(36*ij0-11)= allu0(36*ij0-11) - rhs_aij(5,1)
530  allu0(36*ij0-10)= allu0(36*ij0-10) - rhs_aij(5,2)
531  allu0(36*ij0-9 )= allu0(36*ij0-9 ) - rhs_aij(5,3)
532  allu0(36*ij0-8 )= allu0(36*ij0-8 ) - rhs_aij(5,4)
533  allu0(36*ij0-7 )= allu0(36*ij0-7 ) - rhs_aij(5,5)
534  allu0(36*ij0-6 )= allu0(36*ij0-6 ) - rhs_aij(5,6)
535  allu0(36*ij0-5 )= allu0(36*ij0-5 ) - rhs_aij(6,1)
536  allu0(36*ij0-4 )= allu0(36*ij0-4 ) - rhs_aij(6,2)
537  allu0(36*ij0-3 )= allu0(36*ij0-3 ) - rhs_aij(6,3)
538  allu0(36*ij0-2 )= allu0(36*ij0-2 ) - rhs_aij(6,4)
539  allu0(36*ij0-1 )= allu0(36*ij0-1 ) - rhs_aij(6,5)
540  allu0(36*ij0 )= allu0(36*ij0 ) - rhs_aij(6,6)
541  endif
542 
543  if (j.gt.i) then
544  ij0= iw2(j)
545  aulu0(36*ij0-35)= aulu0(36*ij0-35) - rhs_aij(1,1)
546  aulu0(36*ij0-34)= aulu0(36*ij0-34) - rhs_aij(1,2)
547  aulu0(36*ij0-33)= aulu0(36*ij0-33) - rhs_aij(1,3)
548  aulu0(36*ij0-32)= aulu0(36*ij0-32) - rhs_aij(1,4)
549  aulu0(36*ij0-31)= aulu0(36*ij0-31) - rhs_aij(1,5)
550  aulu0(36*ij0-30)= aulu0(36*ij0-30) - rhs_aij(1,6)
551  aulu0(36*ij0-29)= aulu0(36*ij0-29) - rhs_aij(2,1)
552  aulu0(36*ij0-28)= aulu0(36*ij0-28) - rhs_aij(2,2)
553  aulu0(36*ij0-27)= aulu0(36*ij0-27) - rhs_aij(2,3)
554  aulu0(36*ij0-26)= aulu0(36*ij0-26) - rhs_aij(2,4)
555  aulu0(36*ij0-25)= aulu0(36*ij0-25) - rhs_aij(2,5)
556  aulu0(36*ij0-24)= aulu0(36*ij0-24) - rhs_aij(2,6)
557  aulu0(36*ij0-23)= aulu0(36*ij0-23) - rhs_aij(3,1)
558  aulu0(36*ij0-22)= aulu0(36*ij0-22) - rhs_aij(3,2)
559  aulu0(36*ij0-21)= aulu0(36*ij0-21) - rhs_aij(3,3)
560  aulu0(36*ij0-20)= aulu0(36*ij0-20) - rhs_aij(3,4)
561  aulu0(36*ij0-19)= aulu0(36*ij0-19) - rhs_aij(3,5)
562  aulu0(36*ij0-18)= aulu0(36*ij0-18) - rhs_aij(3,6)
563  aulu0(36*ij0-17)= aulu0(36*ij0-17) - rhs_aij(4,1)
564  aulu0(36*ij0-16)= aulu0(36*ij0-16) - rhs_aij(4,2)
565  aulu0(36*ij0-15)= aulu0(36*ij0-15) - rhs_aij(4,3)
566  aulu0(36*ij0-14)= aulu0(36*ij0-14) - rhs_aij(4,4)
567  aulu0(36*ij0-13)= aulu0(36*ij0-13) - rhs_aij(4,5)
568  aulu0(36*ij0-12)= aulu0(36*ij0-12) - rhs_aij(4,6)
569  aulu0(36*ij0-11)= aulu0(36*ij0-11) - rhs_aij(5,1)
570  aulu0(36*ij0-10)= aulu0(36*ij0-10) - rhs_aij(5,2)
571  aulu0(36*ij0-9 )= aulu0(36*ij0-9 ) - rhs_aij(5,3)
572  aulu0(36*ij0-8 )= aulu0(36*ij0-8 ) - rhs_aij(5,4)
573  aulu0(36*ij0-7 )= aulu0(36*ij0-7 ) - rhs_aij(5,5)
574  aulu0(36*ij0-6 )= aulu0(36*ij0-6 ) - rhs_aij(5,6)
575  aulu0(36*ij0-5 )= aulu0(36*ij0-5 ) - rhs_aij(6,1)
576  aulu0(36*ij0-4 )= aulu0(36*ij0-4 ) - rhs_aij(6,2)
577  aulu0(36*ij0-3 )= aulu0(36*ij0-3 ) - rhs_aij(6,3)
578  aulu0(36*ij0-2 )= aulu0(36*ij0-2 ) - rhs_aij(6,4)
579  aulu0(36*ij0-1 )= aulu0(36*ij0-1 ) - rhs_aij(6,5)
580  aulu0(36*ij0 )= aulu0(36*ij0 ) - rhs_aij(6,6)
581  endif
582 
583  enddo
584  enddo
585  enddo
586 
587  deallocate (iw1, iw2)
588  end subroutine form_ilu0_66
589 
590  !C
591  !C***
592  !C*** FORM_ILU1_66
593  !C***
594  !C
595  !C form ILU(1) matrix
596  !C
597  subroutine form_ilu1_66 &
598  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
599  & sigma_diag)
600  implicit none
601  integer(kind=kint ), intent(in):: n, np, npu, npl
602  real (kind=kreal), intent(in):: sigma_diag
603 
604  real(kind=kreal), dimension(36*NPL), intent(in):: al
605  real(kind=kreal), dimension(36*NPU), intent(in):: au
606  real(kind=kreal), dimension(36*NP ), intent(in):: d
607 
608  integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
609  integer(kind=kint ), dimension( NPL),intent(in) :: ial
610  integer(kind=kint ), dimension( NPU),intent(in) :: iau
611 
612  integer(kind=kint), dimension(:), allocatable :: iw1, iw2
613  integer(kind=kint), dimension(:), allocatable :: iwsl, iwsu
614  real (kind=kreal), dimension(6,6) :: rhs_aij, dkinv, aik, akj
615  real(kind=kreal) :: d11,d12,d13,d14,d15,d16
616  real(kind=kreal) :: d21,d22,d23,d24,d25,d26
617  real(kind=kreal) :: d31,d32,d33,d34,d35,d36
618  real(kind=kreal) :: d41,d42,d43,d44,d45,d46
619  real(kind=kreal) :: d51,d52,d53,d54,d55,d56
620  real(kind=kreal) :: d61,d62,d63,d64,d65,d66
621  integer(kind=kint) :: nplf1,npuf1
622  integer(kind=kint) :: i,jj,jj1,ij0,kk,ik,kk1,kk2,l,isk,iek,isj,iej
623  integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
624  integer(kind=kint) :: j,k,isl,isu
625  !C
626  !C +--------------+
627  !C | find fill-in |
628  !C +--------------+
629  !C===
630 
631  !C
632  !C-- count fill-in
633  allocate (iw1(np) , iw2(np))
634  allocate (inumfi1l(0:np), inumfi1u(0:np))
635 
636  inumfi1l= 0
637  inumfi1u= 0
638 
639  nplf1= 0
640  npuf1= 0
641  do i= 2, np
642  icou= 0
643  iw1= 0
644  iw1(i)= 1
645  do l= inl(i-1)+1, inl(i)
646  iw1(ial(l))= 1
647  enddo
648  do l= inu(i-1)+1, inu(i)
649  iw1(iau(l))= 1
650  enddo
651 
652  isk= inl(i-1) + 1
653  iek= inl(i)
654  do k= isk, iek
655  kk= ial(k)
656  isj= inu(kk-1) + 1
657  iej= inu(kk )
658  do j= isj, iej
659  jj= iau(j)
660  if (iw1(jj).eq.0 .and. jj.lt.i) then
661  inumfi1l(i)= inumfi1l(i)+1
662  iw1(jj)= 1
663  endif
664  if (iw1(jj).eq.0 .and. jj.gt.i) then
665  inumfi1u(i)= inumfi1u(i)+1
666  iw1(jj)= 1
667  endif
668  enddo
669  enddo
670  nplf1= nplf1 + inumfi1l(i)
671  npuf1= npuf1 + inumfi1u(i)
672  enddo
673 
674  !C
675  !C-- specify fill-in
676  allocate (iwsl(0:np), iwsu(0:np))
677  allocate (fi1l(npl+nplf1), fi1u(npu+npuf1))
678  allocate (allu0(36*(npl+nplf1)), aulu0(36*(npu+npuf1)))
679 
680  fi1l= 0
681  fi1u= 0
682 
683  iwsl= 0
684  iwsu= 0
685  do i= 1, np
686  iwsl(i)= inl(i)-inl(i-1) + inumfi1l(i) + iwsl(i-1)
687  iwsu(i)= inu(i)-inu(i-1) + inumfi1u(i) + iwsu(i-1)
688  enddo
689 
690  do i= 2, np
691  icoul= 0
692  icouu= 0
693  inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
694  inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
695  icou= 0
696  iw1= 0
697  iw1(i)= 1
698  do l= inl(i-1)+1, inl(i)
699  iw1(ial(l))= 1
700  enddo
701  do l= inu(i-1)+1, inu(i)
702  iw1(iau(l))= 1
703  enddo
704 
705  isk= inl(i-1) + 1
706  iek= inl(i)
707  do k= isk, iek
708  kk= ial(k)
709  isj= inu(kk-1) + 1
710  iej= inu(kk )
711  do j= isj, iej
712  jj= iau(j)
713  if (iw1(jj).eq.0 .and. jj.lt.i) then
714  icoul = icoul + 1
715  fi1l(icoul+iwsl(i-1)+inl(i)-inl(i-1))= jj
716  iw1(jj) = 1
717  endif
718  if (iw1(jj).eq.0 .and. jj.gt.i) then
719  icouu = icouu + 1
720  fi1u(icouu+iwsu(i-1)+inu(i)-inu(i-1))= jj
721  iw1(jj) = 1
722  endif
723  enddo
724  enddo
725  enddo
726  !C===
727 
728  !C
729  !C +-------------------------------------------------+
730  !C | SORT and RECONSTRUCT matrix considering fill-in |
731  !C +-------------------------------------------------+
732  !C===
733  allu0= 0.d0
734  aulu0= 0.d0
735  isl = 0
736  isu = 0
737  do i= 1, np
738  icoul1= inl(i) - inl(i-1)
739  icoul2= inumfi1l(i) - inumfi1l(i-1)
740  icoul3= icoul1 + icoul2
741  icouu1= inu(i) - inu(i-1)
742  icouu2= inumfi1u(i) - inumfi1u(i-1)
743  icouu3= icouu1 + icouu2
744  !C
745  !C-- LOWER part
746  icou0= 0
747  do k= inl(i-1)+1, inl(i)
748  icou0 = icou0 + 1
749  iw1(icou0)= ial(k)
750  enddo
751 
752  do k= inumfi1l(i-1)+1, inumfi1l(i)
753  icou0 = icou0 + 1
754  iw1(icou0)= fi1l(icou0+iwsl(i-1))
755  enddo
756 
757  do k= 1, icoul3
758  iw2(k)= k
759  enddo
760  call fill_in_s66_sort (iw1, iw2, icoul3, np)
761 
762  do k= 1, icoul3
763  fi1l(k+isl)= iw1(k)
764  ik= iw2(k)
765  if (ik.le.inl(i)-inl(i-1)) then
766  kk1= 36*( k+isl)
767  kk2= 36*(ik+inl(i-1))
768  allu0(kk1-35)= al(kk2-35)
769  allu0(kk1-34)= al(kk2-34)
770  allu0(kk1-33)= al(kk2-33)
771  allu0(kk1-32)= al(kk2-32)
772  allu0(kk1-31)= al(kk2-31)
773  allu0(kk1-30)= al(kk2-30)
774  allu0(kk1-29)= al(kk2-29)
775  allu0(kk1-28)= al(kk2-28)
776  allu0(kk1-27)= al(kk2-27)
777  allu0(kk1-26)= al(kk2-26)
778  allu0(kk1-25)= al(kk2-25)
779  allu0(kk1-24)= al(kk2-24)
780  allu0(kk1-23)= al(kk2-23)
781  allu0(kk1-22)= al(kk2-22)
782  allu0(kk1-21)= al(kk2-21)
783  allu0(kk1-20)= al(kk2-20)
784  allu0(kk1-19)= al(kk2-19)
785  allu0(kk1-18)= al(kk2-18)
786  allu0(kk1-17)= al(kk2-17)
787  allu0(kk1-16)= al(kk2-16)
788  allu0(kk1-15)= al(kk2-15)
789  allu0(kk1-14)= al(kk2-14)
790  allu0(kk1-13)= al(kk2-13)
791  allu0(kk1-12)= al(kk2-12)
792  allu0(kk1-11)= al(kk2-11)
793  allu0(kk1-10)= al(kk2-10)
794  allu0(kk1-9 )= al(kk2-9 )
795  allu0(kk1-8 )= al(kk2-8 )
796  allu0(kk1-7 )= al(kk2-7 )
797  allu0(kk1-6 )= al(kk2-6 )
798  allu0(kk1-5 )= al(kk2-5 )
799  allu0(kk1-4 )= al(kk2-4 )
800  allu0(kk1-3 )= al(kk2-3 )
801  allu0(kk1-2 )= al(kk2-2 )
802  allu0(kk1-1 )= al(kk2-1 )
803  allu0(kk1 )= al(kk2 )
804  endif
805  enddo
806  !C
807  !C-- UPPER part
808  icou0= 0
809  do k= inu(i-1)+1, inu(i)
810  icou0 = icou0 + 1
811  iw1(icou0)= iau(k)
812  enddo
813 
814  do k= inumfi1u(i-1)+1, inumfi1u(i)
815  icou0 = icou0 + 1
816  iw1(icou0)= fi1u(icou0+iwsu(i-1))
817  enddo
818 
819  do k= 1, icouu3
820  iw2(k)= k
821  enddo
822  call fill_in_s66_sort (iw1, iw2, icouu3, np)
823 
824  do k= 1, icouu3
825  fi1u(k+isu)= iw1(k)
826  ik= iw2(k)
827  if (ik.le.inu(i)-inu(i-1)) then
828  kk1= 36*( k+isu)
829  kk2= 36*(ik+inu(i-1))
830  aulu0(kk1-35)= au(kk2-35)
831  aulu0(kk1-34)= au(kk2-34)
832  aulu0(kk1-33)= au(kk2-33)
833  aulu0(kk1-32)= au(kk2-32)
834  aulu0(kk1-31)= au(kk2-31)
835  aulu0(kk1-30)= au(kk2-30)
836  aulu0(kk1-29)= au(kk2-29)
837  aulu0(kk1-28)= au(kk2-28)
838  aulu0(kk1-27)= au(kk2-27)
839  aulu0(kk1-26)= au(kk2-26)
840  aulu0(kk1-25)= au(kk2-25)
841  aulu0(kk1-24)= au(kk2-24)
842  aulu0(kk1-23)= au(kk2-23)
843  aulu0(kk1-22)= au(kk2-22)
844  aulu0(kk1-21)= au(kk2-21)
845  aulu0(kk1-20)= au(kk2-20)
846  aulu0(kk1-19)= au(kk2-19)
847  aulu0(kk1-18)= au(kk2-18)
848  aulu0(kk1-17)= au(kk2-17)
849  aulu0(kk1-16)= au(kk2-16)
850  aulu0(kk1-15)= au(kk2-15)
851  aulu0(kk1-14)= au(kk2-14)
852  aulu0(kk1-13)= au(kk2-13)
853  aulu0(kk1-12)= au(kk2-12)
854  aulu0(kk1-11)= au(kk2-11)
855  aulu0(kk1-10)= au(kk2-10)
856  aulu0(kk1-9 )= au(kk2-9 )
857  aulu0(kk1-8 )= au(kk2-8 )
858  aulu0(kk1-7 )= au(kk2-7 )
859  aulu0(kk1-6 )= au(kk2-6 )
860  aulu0(kk1-5 )= au(kk2-5 )
861  aulu0(kk1-4 )= au(kk2-4 )
862  aulu0(kk1-3 )= au(kk2-3 )
863  aulu0(kk1-2 )= au(kk2-2 )
864  aulu0(kk1-1 )= au(kk2-1 )
865  aulu0(kk1 )= au(kk2 )
866  endif
867  enddo
868 
869  isl= isl + icoul3
870  isu= isu + icouu3
871  enddo
872 
873  !C===
874  do i= 1, np
875  inumfi1l(i)= iwsl(i)
876  inumfi1u(i)= iwsu(i)
877  enddo
878  deallocate (iwsl, iwsu)
879 
880  !C
881  !C +----------------------+
882  !C | ILU(1) factorization |
883  !C +----------------------+
884  !C===
885  allocate (dlu0(36*np))
886  dlu0= d
887  do i=1,np
888  dlu0(36*i-35)=dlu0(36*i-35)*sigma_diag
889  dlu0(36*i-28)=dlu0(36*i-28)*sigma_diag
890  dlu0(36*i-21)=dlu0(36*i-21)*sigma_diag
891  dlu0(36*i-14)=dlu0(36*i-14)*sigma_diag
892  dlu0(36*i-7 )=dlu0(36*i-7 )*sigma_diag
893  dlu0(36*i )=dlu0(36*i )*sigma_diag
894  enddo
895 
896  do i= 2, np
897  iw1= 0
898  iw2= 0
899 
900  do k= inumfi1l(i-1)+1, inumfi1l(i)
901  iw1(fi1l(k))= k
902  enddo
903 
904  do k= inumfi1u(i-1)+1, inumfi1u(i)
905  iw2(fi1u(k))= k
906  enddo
907 
908  do kk= inl(i-1)+1, inl(i)
909  k= ial(kk)
910  d11= dlu0(36*k-35)
911  d12= dlu0(36*k-34)
912  d13= dlu0(36*k-33)
913  d14= dlu0(36*k-32)
914  d15= dlu0(36*k-31)
915  d16= dlu0(36*k-30)
916  d21= dlu0(36*k-29)
917  d22= dlu0(36*k-28)
918  d23= dlu0(36*k-27)
919  d24= dlu0(36*k-26)
920  d25= dlu0(36*k-25)
921  d26= dlu0(36*k-24)
922  d31= dlu0(36*k-23)
923  d32= dlu0(36*k-22)
924  d33= dlu0(36*k-21)
925  d34= dlu0(36*k-20)
926  d35= dlu0(36*k-19)
927  d36= dlu0(36*k-18)
928  d41= dlu0(36*k-17)
929  d42= dlu0(36*k-16)
930  d43= dlu0(36*k-15)
931  d44= dlu0(36*k-14)
932  d45= dlu0(36*k-13)
933  d46= dlu0(36*k-12)
934  d51= dlu0(36*k-11)
935  d52= dlu0(36*k-10)
936  d53= dlu0(36*k-9 )
937  d54= dlu0(36*k-8 )
938  d55= dlu0(36*k-7 )
939  d56= dlu0(36*k-6 )
940  d61= dlu0(36*k-5 )
941  d62= dlu0(36*k-4 )
942  d63= dlu0(36*k-3 )
943  d64= dlu0(36*k-2 )
944  d65= dlu0(36*k-1 )
945  d66= dlu0(36*k )
946 
947  call ilu1a66 (dkinv,d11,d12,d13,d14,d15,d16,d21,d22,d23,d24,d25,d26, &
948  & d31,d32,d33,d34,d35,d36,d41,d42,d43,d44,d45,d46,d51,d52,d53,d54,d55,d56, &
949  & d61,d62,d63,d64,d65,d66)
950 
951  do kk1= inumfi1l(i-1)+1, inumfi1l(i)
952  if (k.eq.fi1l(kk1)) then
953  aik(1,1)= allu0(36*kk1-35)
954  aik(1,2)= allu0(36*kk1-34)
955  aik(1,3)= allu0(36*kk1-33)
956  aik(1,4)= allu0(36*kk1-32)
957  aik(1,5)= allu0(36*kk1-31)
958  aik(1,6)= allu0(36*kk1-30)
959  aik(2,1)= allu0(36*kk1-29)
960  aik(2,2)= allu0(36*kk1-28)
961  aik(2,3)= allu0(36*kk1-27)
962  aik(2,4)= allu0(36*kk1-26)
963  aik(2,5)= allu0(36*kk1-25)
964  aik(2,6)= allu0(36*kk1-24)
965  aik(3,1)= allu0(36*kk1-23)
966  aik(3,2)= allu0(36*kk1-22)
967  aik(3,3)= allu0(36*kk1-21)
968  aik(3,4)= allu0(36*kk1-20)
969  aik(3,5)= allu0(36*kk1-19)
970  aik(3,6)= allu0(36*kk1-18)
971  aik(4,1)= allu0(36*kk1-17)
972  aik(4,2)= allu0(36*kk1-16)
973  aik(4,3)= allu0(36*kk1-15)
974  aik(4,4)= allu0(36*kk1-14)
975  aik(4,5)= allu0(36*kk1-13)
976  aik(4,6)= allu0(36*kk1-12)
977  aik(5,1)= allu0(36*kk1-11)
978  aik(5,2)= allu0(36*kk1-10)
979  aik(5,3)= allu0(36*kk1-9)
980  aik(5,4)= allu0(36*kk1-8)
981  aik(5,5)= allu0(36*kk1-7)
982  aik(5,6)= allu0(36*kk1-6)
983  aik(6,1)= allu0(36*kk1-5)
984  aik(6,2)= allu0(36*kk1-4)
985  aik(6,3)= allu0(36*kk1-3)
986  aik(6,4)= allu0(36*kk1-2)
987  aik(6,5)= allu0(36*kk1-1)
988  aik(6,6)= allu0(36*kk1 )
989  exit
990  endif
991  enddo
992 
993  do jj= inu(k-1)+1, inu(k)
994  j= iau(jj)
995  do jj1= inumfi1u(k-1)+1, inumfi1u(k)
996  if (j.eq.fi1u(jj1)) then
997  akj(1,1)= aulu0(36*jj1-35)
998  akj(1,2)= aulu0(36*jj1-34)
999  akj(1,3)= aulu0(36*jj1-33)
1000  akj(1,4)= aulu0(36*jj1-32)
1001  akj(1,5)= aulu0(36*jj1-31)
1002  akj(1,6)= aulu0(36*jj1-30)
1003  akj(2,1)= aulu0(36*jj1-29)
1004  akj(2,2)= aulu0(36*jj1-28)
1005  akj(2,3)= aulu0(36*jj1-27)
1006  akj(2,4)= aulu0(36*jj1-26)
1007  akj(2,5)= aulu0(36*jj1-25)
1008  akj(2,6)= aulu0(36*jj1-24)
1009  akj(3,1)= aulu0(36*jj1-23)
1010  akj(3,2)= aulu0(36*jj1-22)
1011  akj(3,3)= aulu0(36*jj1-21)
1012  akj(3,4)= aulu0(36*jj1-20)
1013  akj(3,5)= aulu0(36*jj1-19)
1014  akj(3,6)= aulu0(36*jj1-18)
1015  akj(4,1)= aulu0(36*jj1-17)
1016  akj(4,2)= aulu0(36*jj1-16)
1017  akj(4,3)= aulu0(36*jj1-15)
1018  akj(4,4)= aulu0(36*jj1-14)
1019  akj(4,5)= aulu0(36*jj1-13)
1020  akj(4,6)= aulu0(36*jj1-12)
1021  akj(5,1)= aulu0(36*jj1-11)
1022  akj(5,2)= aulu0(36*jj1-10)
1023  akj(5,3)= aulu0(36*jj1-9)
1024  akj(5,4)= aulu0(36*jj1-8)
1025  akj(5,5)= aulu0(36*jj1-7)
1026  akj(5,6)= aulu0(36*jj1-6)
1027  akj(6,1)= aulu0(36*jj1-5)
1028  akj(6,2)= aulu0(36*jj1-4)
1029  akj(6,3)= aulu0(36*jj1-3)
1030  akj(6,4)= aulu0(36*jj1-2)
1031  akj(6,5)= aulu0(36*jj1-1)
1032  akj(6,6)= aulu0(36*jj1 )
1033  exit
1034  endif
1035  enddo
1036 
1037  call ilu1b66 (rhs_aij, dkinv, aik, akj)
1038 
1039  if (j.eq.i) then
1040  dlu0(36*i-35)= dlu0(36*i-35) - rhs_aij(1,1)
1041  dlu0(36*i-34)= dlu0(36*i-34) - rhs_aij(1,2)
1042  dlu0(36*i-33)= dlu0(36*i-33) - rhs_aij(1,3)
1043  dlu0(36*i-32)= dlu0(36*i-32) - rhs_aij(1,4)
1044  dlu0(36*i-31)= dlu0(36*i-31) - rhs_aij(1,5)
1045  dlu0(36*i-30)= dlu0(36*i-30) - rhs_aij(1,6)
1046  dlu0(36*i-29)= dlu0(36*i-29) - rhs_aij(2,1)
1047  dlu0(36*i-28)= dlu0(36*i-28) - rhs_aij(2,2)
1048  dlu0(36*i-27)= dlu0(36*i-27) - rhs_aij(2,3)
1049  dlu0(36*i-26)= dlu0(36*i-26) - rhs_aij(2,4)
1050  dlu0(36*i-25)= dlu0(36*i-25) - rhs_aij(2,5)
1051  dlu0(36*i-24)= dlu0(36*i-24) - rhs_aij(2,6)
1052  dlu0(36*i-23)= dlu0(36*i-23) - rhs_aij(3,1)
1053  dlu0(36*i-22)= dlu0(36*i-22) - rhs_aij(3,2)
1054  dlu0(36*i-21)= dlu0(36*i-21) - rhs_aij(3,3)
1055  dlu0(36*i-20)= dlu0(36*i-20) - rhs_aij(3,4)
1056  dlu0(36*i-19)= dlu0(36*i-19) - rhs_aij(3,5)
1057  dlu0(36*i-18)= dlu0(36*i-18) - rhs_aij(3,6)
1058  dlu0(36*i-17)= dlu0(36*i-17) - rhs_aij(4,1)
1059  dlu0(36*i-16)= dlu0(36*i-16) - rhs_aij(4,2)
1060  dlu0(36*i-15)= dlu0(36*i-15) - rhs_aij(4,3)
1061  dlu0(36*i-14)= dlu0(36*i-14) - rhs_aij(4,4)
1062  dlu0(36*i-13)= dlu0(36*i-13) - rhs_aij(4,5)
1063  dlu0(36*i-12)= dlu0(36*i-12) - rhs_aij(4,6)
1064  dlu0(36*i-11)= dlu0(36*i-11) - rhs_aij(5,1)
1065  dlu0(36*i-10)= dlu0(36*i-10) - rhs_aij(5,2)
1066  dlu0(36*i-9 )= dlu0(36*i-9 ) - rhs_aij(5,3)
1067  dlu0(36*i-8 )= dlu0(36*i-8 ) - rhs_aij(5,4)
1068  dlu0(36*i-7 )= dlu0(36*i-7 ) - rhs_aij(5,5)
1069  dlu0(36*i-6 )= dlu0(36*i-6 ) - rhs_aij(5,6)
1070  dlu0(36*i-5 )= dlu0(36*i-5 ) - rhs_aij(6,1)
1071  dlu0(36*i-4 )= dlu0(36*i-4 ) - rhs_aij(6,2)
1072  dlu0(36*i-3 )= dlu0(36*i-3 ) - rhs_aij(6,3)
1073  dlu0(36*i-2 )= dlu0(36*i-2 ) - rhs_aij(6,4)
1074  dlu0(36*i-1 )= dlu0(36*i-1 ) - rhs_aij(6,5)
1075  dlu0(36*i )= dlu0(36*i ) - rhs_aij(6,6)
1076  endif
1077 
1078  if (j.lt.i) then
1079  ij0= iw1(j)
1080  allu0(36*ij0-35)= allu0(36*ij0-35) - rhs_aij(1,1)
1081  allu0(36*ij0-34)= allu0(36*ij0-34) - rhs_aij(1,2)
1082  allu0(36*ij0-33)= allu0(36*ij0-33) - rhs_aij(1,3)
1083  allu0(36*ij0-32)= allu0(36*ij0-32) - rhs_aij(1,4)
1084  allu0(36*ij0-31)= allu0(36*ij0-31) - rhs_aij(1,5)
1085  allu0(36*ij0-30)= allu0(36*ij0-30) - rhs_aij(1,6)
1086  allu0(36*ij0-29)= allu0(36*ij0-29) - rhs_aij(2,1)
1087  allu0(36*ij0-28)= allu0(36*ij0-28) - rhs_aij(2,2)
1088  allu0(36*ij0-27)= allu0(36*ij0-27) - rhs_aij(2,3)
1089  allu0(36*ij0-26)= allu0(36*ij0-26) - rhs_aij(2,4)
1090  allu0(36*ij0-25)= allu0(36*ij0-25) - rhs_aij(2,5)
1091  allu0(36*ij0-24)= allu0(36*ij0-24) - rhs_aij(2,6)
1092  allu0(36*ij0-23)= allu0(36*ij0-23) - rhs_aij(3,1)
1093  allu0(36*ij0-22)= allu0(36*ij0-22) - rhs_aij(3,2)
1094  allu0(36*ij0-21)= allu0(36*ij0-21) - rhs_aij(3,3)
1095  allu0(36*ij0-20)= allu0(36*ij0-20) - rhs_aij(3,4)
1096  allu0(36*ij0-19)= allu0(36*ij0-19) - rhs_aij(3,5)
1097  allu0(36*ij0-18)= allu0(36*ij0-18) - rhs_aij(3,6)
1098  allu0(36*ij0-17)= allu0(36*ij0-17) - rhs_aij(4,1)
1099  allu0(36*ij0-16)= allu0(36*ij0-16) - rhs_aij(4,2)
1100  allu0(36*ij0-15)= allu0(36*ij0-15) - rhs_aij(4,3)
1101  allu0(36*ij0-14)= allu0(36*ij0-14) - rhs_aij(4,4)
1102  allu0(36*ij0-13)= allu0(36*ij0-13) - rhs_aij(4,5)
1103  allu0(36*ij0-12)= allu0(36*ij0-12) - rhs_aij(4,6)
1104  allu0(36*ij0-11)= allu0(36*ij0-11) - rhs_aij(5,1)
1105  allu0(36*ij0-10)= allu0(36*ij0-10) - rhs_aij(5,2)
1106  allu0(36*ij0-9 )= allu0(36*ij0-9 ) - rhs_aij(5,3)
1107  allu0(36*ij0-8 )= allu0(36*ij0-8 ) - rhs_aij(5,4)
1108  allu0(36*ij0-7 )= allu0(36*ij0-7 ) - rhs_aij(5,5)
1109  allu0(36*ij0-6 )= allu0(36*ij0-6 ) - rhs_aij(5,6)
1110  allu0(36*ij0-5 )= allu0(36*ij0-5 ) - rhs_aij(6,1)
1111  allu0(36*ij0-4 )= allu0(36*ij0-4 ) - rhs_aij(6,2)
1112  allu0(36*ij0-3 )= allu0(36*ij0-3 ) - rhs_aij(6,3)
1113  allu0(36*ij0-2 )= allu0(36*ij0-2 ) - rhs_aij(6,4)
1114  allu0(36*ij0-1 )= allu0(36*ij0-1 ) - rhs_aij(6,5)
1115  allu0(36*ij0 )= allu0(36*ij0 ) - rhs_aij(6,6)
1116  endif
1117 
1118  if (j.gt.i) then
1119  ij0= iw2(j)
1120  aulu0(36*ij0-35)= aulu0(36*ij0-35) - rhs_aij(1,1)
1121  aulu0(36*ij0-34)= aulu0(36*ij0-34) - rhs_aij(1,2)
1122  aulu0(36*ij0-33)= aulu0(36*ij0-33) - rhs_aij(1,3)
1123  aulu0(36*ij0-32)= aulu0(36*ij0-32) - rhs_aij(1,4)
1124  aulu0(36*ij0-31)= aulu0(36*ij0-31) - rhs_aij(1,5)
1125  aulu0(36*ij0-30)= aulu0(36*ij0-30) - rhs_aij(1,6)
1126  aulu0(36*ij0-29)= aulu0(36*ij0-29) - rhs_aij(2,1)
1127  aulu0(36*ij0-28)= aulu0(36*ij0-28) - rhs_aij(2,2)
1128  aulu0(36*ij0-27)= aulu0(36*ij0-27) - rhs_aij(2,3)
1129  aulu0(36*ij0-26)= aulu0(36*ij0-26) - rhs_aij(2,4)
1130  aulu0(36*ij0-25)= aulu0(36*ij0-25) - rhs_aij(2,5)
1131  aulu0(36*ij0-24)= aulu0(36*ij0-24) - rhs_aij(2,6)
1132  aulu0(36*ij0-23)= aulu0(36*ij0-23) - rhs_aij(3,1)
1133  aulu0(36*ij0-22)= aulu0(36*ij0-22) - rhs_aij(3,2)
1134  aulu0(36*ij0-21)= aulu0(36*ij0-21) - rhs_aij(3,3)
1135  aulu0(36*ij0-20)= aulu0(36*ij0-20) - rhs_aij(3,4)
1136  aulu0(36*ij0-19)= aulu0(36*ij0-19) - rhs_aij(3,5)
1137  aulu0(36*ij0-18)= aulu0(36*ij0-18) - rhs_aij(3,6)
1138  aulu0(36*ij0-17)= aulu0(36*ij0-17) - rhs_aij(4,1)
1139  aulu0(36*ij0-16)= aulu0(36*ij0-16) - rhs_aij(4,2)
1140  aulu0(36*ij0-15)= aulu0(36*ij0-15) - rhs_aij(4,3)
1141  aulu0(36*ij0-14)= aulu0(36*ij0-14) - rhs_aij(4,4)
1142  aulu0(36*ij0-13)= aulu0(36*ij0-13) - rhs_aij(4,5)
1143  aulu0(36*ij0-12)= aulu0(36*ij0-12) - rhs_aij(4,6)
1144  aulu0(36*ij0-11)= aulu0(36*ij0-11) - rhs_aij(5,1)
1145  aulu0(36*ij0-10)= aulu0(36*ij0-10) - rhs_aij(5,2)
1146  aulu0(36*ij0-9 )= aulu0(36*ij0-9 ) - rhs_aij(5,3)
1147  aulu0(36*ij0-8 )= aulu0(36*ij0-8 ) - rhs_aij(5,4)
1148  aulu0(36*ij0-7 )= aulu0(36*ij0-7 ) - rhs_aij(5,5)
1149  aulu0(36*ij0-6 )= aulu0(36*ij0-6 ) - rhs_aij(5,6)
1150  aulu0(36*ij0-5 )= aulu0(36*ij0-5 ) - rhs_aij(6,1)
1151  aulu0(36*ij0-4 )= aulu0(36*ij0-4 ) - rhs_aij(6,2)
1152  aulu0(36*ij0-3 )= aulu0(36*ij0-3 ) - rhs_aij(6,3)
1153  aulu0(36*ij0-2 )= aulu0(36*ij0-2 ) - rhs_aij(6,4)
1154  aulu0(36*ij0-1 )= aulu0(36*ij0-1 ) - rhs_aij(6,5)
1155  aulu0(36*ij0 )= aulu0(36*ij0 ) - rhs_aij(6,6)
1156  endif
1157 
1158  enddo
1159  enddo
1160  enddo
1161 
1162  deallocate (iw1, iw2)
1163  !C===
1164  end subroutine form_ilu1_66
1165 
1166  !C
1167  !C***
1168  !C*** FORM_ILU2_66
1169  !C***
1170  !C
1171  !C form ILU(2) matrix
1172  !C
1173  subroutine form_ilu2_66 &
1174  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
1175  & sigma_diag)
1176  implicit none
1177  integer(kind=kint ), intent(in):: n, np, npu, npl
1178  real (kind=kreal), intent(in):: sigma_diag
1179 
1180  real(kind=kreal), dimension(36*NPL), intent(in):: al
1181  real(kind=kreal), dimension(36*NPU), intent(in):: au
1182  real(kind=kreal), dimension(36*NP ), intent(in):: d
1183 
1184  integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
1185  integer(kind=kint ), dimension( NPL),intent(in) :: ial
1186  integer(kind=kint ), dimension( NPU),intent(in) :: iau
1187 
1188  integer(kind=kint), dimension(:), allocatable:: iw1 , iw2
1189  integer(kind=kint), dimension(:), allocatable:: iwsl, iwsu
1190  integer(kind=kint), dimension(:), allocatable:: iconfi1l, iconfi1u
1191  integer(kind=kint), dimension(:), allocatable:: inumfi2l, inumfi2u
1192  integer(kind=kint), dimension(:), allocatable:: fi2l, fi2u
1193  real (kind=kreal), dimension(6,6) :: rhs_aij, dkinv, aik, akj
1194  real(kind=kreal) :: d11,d12,d13,d14,d15,d16
1195  real(kind=kreal) :: d21,d22,d23,d24,d25,d26
1196  real(kind=kreal) :: d31,d32,d33,d34,d35,d36
1197  real(kind=kreal) :: d41,d42,d43,d44,d45,d46
1198  real(kind=kreal) :: d51,d52,d53,d54,d55,d56
1199  real(kind=kreal) :: d61,d62,d63,d64,d65,d66
1200  integer(kind=kint) :: nplf1,nplf2,npuf1,npuf2,ias,iconik,iconkj
1201  integer(kind=kint) :: i,jj,ij0,kk,ik,kk1,kk2,l,isk,iek,isj,iej
1202  integer(kind=kint) :: icou,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
1203  integer(kind=kint) :: j,k,isl,isu
1204 
1205  !C
1206  !C +------------------+
1207  !C | find fill-in (1) |
1208  !C +------------------+
1209  !C===
1210 
1211  !C
1212  !C-- count fill-in
1213  allocate (iw1(np) , iw2(np))
1214  allocate (inumfi2l(0:np), inumfi2u(0:np))
1215 
1216  inumfi2l= 0
1217  inumfi2u= 0
1218 
1219  nplf1= 0
1220  npuf1= 0
1221  do i= 2, np
1222  icou= 0
1223  iw1= 0
1224  iw1(i)= 1
1225  do l= inl(i-1)+1, inl(i)
1226  iw1(ial(l))= 1
1227  enddo
1228  do l= inu(i-1)+1, inu(i)
1229  iw1(iau(l))= 1
1230  enddo
1231 
1232  isk= inl(i-1) + 1
1233  iek= inl(i)
1234  do k= isk, iek
1235  kk= ial(k)
1236  isj= inu(kk-1) + 1
1237  iej= inu(kk )
1238  do j= isj, iej
1239  jj= iau(j)
1240  if (iw1(jj).eq.0 .and. jj.lt.i) then
1241  inumfi2l(i)= inumfi2l(i)+1
1242  iw1(jj)= 1
1243  endif
1244  if (iw1(jj).eq.0 .and. jj.gt.i) then
1245  inumfi2u(i)= inumfi2u(i)+1
1246  iw1(jj)= 1
1247  endif
1248  enddo
1249  enddo
1250  nplf1= nplf1 + inumfi2l(i)
1251  npuf1= npuf1 + inumfi2u(i)
1252  enddo
1253 
1254  !C
1255  !C-- specify fill-in
1256  allocate (iwsl(0:np), iwsu(0:np))
1257  allocate (fi2l(nplf1), fi2u(npuf1))
1258 
1259  fi2l= 0
1260  fi2u= 0
1261 
1262  do i= 2, np
1263  icoul= 0
1264  icouu= 0
1265  inumfi2l(i)= inumfi2l(i-1) + inumfi2l(i)
1266  inumfi2u(i)= inumfi2u(i-1) + inumfi2u(i)
1267  icou= 0
1268  iw1= 0
1269  iw1(i)= 1
1270  do l= inl(i-1)+1, inl(i)
1271  iw1(ial(l))= 1
1272  enddo
1273  do l= inu(i-1)+1, inu(i)
1274  iw1(iau(l))= 1
1275  enddo
1276 
1277  isk= inl(i-1) + 1
1278  iek= inl(i)
1279  do k= isk, iek
1280  kk= ial(k)
1281  isj= inu(kk-1) + 1
1282  iej= inu(kk )
1283  do j= isj, iej
1284  jj= iau(j)
1285  if (iw1(jj).eq.0 .and. jj.lt.i) then
1286  icoul = icoul + 1
1287  fi2l(icoul+inumfi2l(i-1))= jj
1288  iw1(jj)= 1
1289  endif
1290  if (iw1(jj).eq.0 .and. jj.gt.i) then
1291  icouu = icouu + 1
1292  fi2u(icouu+inumfi2u(i-1))= jj
1293  iw1(jj)= 1
1294  endif
1295  enddo
1296  enddo
1297  enddo
1298  !C===
1299 
1300  !C
1301  !C +------------------+
1302  !C | find fill-in (2) |
1303  !C +------------------+
1304  !C===
1305  allocate (inumfi1l(0:np), inumfi1u(0:np))
1306 
1307  nplf2= 0
1308  npuf2= 0
1309  inumfi1l= 0
1310  inumfi1u= 0
1311  !C
1312  !C-- count fill-in
1313  do i= 2, np
1314  iw1= 0
1315  iw1(i)= 1
1316  do l= inl(i-1)+1, inl(i)
1317  iw1(ial(l))= 2
1318  enddo
1319  do l= inu(i-1)+1, inu(i)
1320  iw1(iau(l))= 2
1321  enddo
1322 
1323  do l= inumfi2l(i-1)+1, inumfi2l(i)
1324  iw1(fi2l(l))= 1
1325  enddo
1326 
1327  do l= inumfi2u(i-1)+1, inumfi2u(i)
1328  iw1(fi2u(l))= 1
1329  enddo
1330 
1331  isk= inl(i-1) + 1
1332  iek= inl(i)
1333  do k= isk, iek
1334  kk= ial(k)
1335  isj= inumfi2u(kk-1) + 1
1336  iej= inumfi2u(kk)
1337  do j= isj, iej
1338  jj= fi2u(j)
1339  if (iw1(jj).eq.0 .and. jj.lt.i) then
1340  inumfi1l(i)= inumfi1l(i) + 1
1341  iw1(jj)= 1
1342  endif
1343  if (iw1(jj).eq.0 .and. jj.gt.i) then
1344  inumfi1u(i)= inumfi1u(i) + 1
1345  iw1(jj)= 1
1346  endif
1347  enddo
1348  enddo
1349 
1350  isk= inumfi2l(i-1)+1
1351  iek= inumfi2l(i)
1352  do k= isk, iek
1353  kk= fi2l(k)
1354  isj= inu(kk-1) + 1
1355  iej= inu(kk )
1356  do j= isj, iej
1357  jj= iau(j)
1358  if (iw1(jj).eq.0 .and. jj.lt.i) then
1359  inumfi1l(i)= inumfi1l(i) + 1
1360  iw1(jj)= 1
1361  endif
1362  if (iw1(jj).eq.0 .and. jj.gt.i) then
1363  inumfi1u(i)= inumfi1u(i) + 1
1364  iw1(jj)= 1
1365  endif
1366  enddo
1367  enddo
1368  nplf2= nplf2 + inumfi1l(i)
1369  npuf2= npuf2 + inumfi1u(i)
1370  enddo
1371 
1372  !C
1373  !C-- specify fill-in
1374  allocate (fi1l(npl+nplf1+nplf2))
1375  allocate (fi1u(npu+npuf1+npuf2))
1376 
1377  allocate (iconfi1l(npl+nplf1+nplf2))
1378  allocate (iconfi1u(npu+npuf1+npuf2))
1379 
1380  iwsl= 0
1381  iwsu= 0
1382  do i= 1, np
1383  iwsl(i)= inl(i)-inl(i-1) + inumfi2l(i)-inumfi2l(i-1) + &
1384  & inumfi1l(i) + iwsl(i-1)
1385  iwsu(i)= inu(i)-inu(i-1) + inumfi2u(i)-inumfi2u(i-1) + &
1386  & inumfi1u(i) + iwsu(i-1)
1387  enddo
1388 
1389  do i= 2, np
1390  icoul= 0
1391  icouu= 0
1392  inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
1393  inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
1394  icou= 0
1395  iw1= 0
1396  iw1(i)= 1
1397  do l= inl(i-1)+1, inl(i)
1398  iw1(ial(l))= 1
1399  enddo
1400  do l= inu(i-1)+1, inu(i)
1401  iw1(iau(l))= 1
1402  enddo
1403 
1404  do l= inumfi2l(i-1)+1, inumfi2l(i)
1405  iw1(fi2l(l))= 1
1406  enddo
1407 
1408  do l= inumfi2u(i-1)+1, inumfi2u(i)
1409  iw1(fi2u(l))= 1
1410  enddo
1411 
1412  isk= inl(i-1) + 1
1413  iek= inl(i)
1414  do k= isk, iek
1415  kk= ial(k)
1416  isj= inumfi2u(kk-1) + 1
1417  iej= inumfi2u(kk )
1418  do j= isj, iej
1419  jj= fi2u(j)
1420  if (iw1(jj).eq.0 .and. jj.lt.i) then
1421  ias= inl(i)-inl(i-1)+inumfi2l(i)-inumfi2l(i-1)+iwsl(i-1)
1422  icoul = icoul + 1
1423  fi1l(icoul+ias)= jj
1424  iw1(jj) = 1
1425  endif
1426  if (iw1(jj).eq.0 .and. jj.gt.i) then
1427  ias= inu(i)-inu(i-1)+inumfi2u(i)-inumfi2u(i-1)+iwsu(i-1)
1428  icouu = icouu + 1
1429  fi1u(icouu+ias)= jj
1430  iw1(jj) = 1
1431  endif
1432  enddo
1433  enddo
1434 
1435  isk= inumfi2l(i-1) + 1
1436  iek= inumfi2l(i)
1437  do k= isk, iek
1438  kk= fi2l(k)
1439  isj= inu(kk-1) + 1
1440  iej= inu(kk )
1441  do j= isj, iej
1442  jj= iau(j)
1443  if (iw1(jj).eq.0 .and. jj.lt.i) then
1444  ias= inl(i)-inl(i-1)+inumfi2l(i)-inumfi2l(i-1)+iwsl(i-1)
1445  icoul = icoul + 1
1446  fi1l(icoul+ias)= jj
1447  iw1(jj) = 1
1448  endif
1449  if (iw1(jj).eq.0 .and. jj.gt.i) then
1450  ias= inu(i)-inu(i-1)+inumfi2u(i)-inumfi2u(i-1)+iwsu(i-1)
1451  icouu = icouu + 1
1452  fi1u(icouu+ias)= jj
1453  iw1(jj) = 1
1454  endif
1455  enddo
1456  enddo
1457  enddo
1458  !C===
1459 
1460  !C
1461  !C +-------------------------------------------------+
1462  !C | SORT and RECONSTRUCT matrix considering fill-in |
1463  !C +-------------------------------------------------+
1464  !C===
1465  allocate (allu0(9*(npl+nplf1+nplf2)))
1466  allocate (aulu0(9*(npu+npuf1+npuf2)))
1467 
1468  allu0= 0.d0
1469  aulu0= 0.d0
1470  isl = 0
1471  isu = 0
1472 
1473  iconfi1l= 0
1474  iconfi1u= 0
1475 
1476  do i= 1, np
1477  icoul1= inl(i) - inl(i-1)
1478  icoul2= inumfi2l(i) - inumfi2l(i-1) + icoul1
1479  icoul3= inumfi1l(i) - inumfi1l(i-1) + icoul2
1480 
1481  icouu1= inu(i) - inu(i-1)
1482  icouu2= inumfi2u(i) - inumfi2u(i-1) + icouu1
1483  icouu3= inumfi1u(i) - inumfi1u(i-1) + icouu2
1484 
1485  !C
1486  !C-- LOWER part
1487  icou= 0
1488  do k= inl(i-1)+1, inl(i)
1489  icou = icou + 1
1490  iw1(icou)= ial(k)
1491  enddo
1492 
1493  icou= 0
1494  do k= inumfi2l(i-1)+1, inumfi2l(i)
1495  icou = icou + 1
1496  iw1(icou+icoul1)= fi2l(k)
1497  enddo
1498 
1499  icou= 0
1500  do k= inumfi1l(i-1)+1, inumfi1l(i)
1501  icou = icou + 1
1502  iw1(icou+icoul2)= fi1l(icou+icoul2+isl)
1503  enddo
1504 
1505  do k= 1, icoul3
1506  iw2(k)= k
1507  enddo
1508 
1509  call fill_in_s66_sort (iw1, iw2, icoul3, np)
1510 
1511  do k= 1, icoul3
1512  fi1l(k+isl)= iw1(k)
1513  ik= iw2(k)
1514  if (ik.le.inl(i)-inl(i-1)) then
1515  kk1= 36*( k+isl)
1516  kk2= 36*(ik+inl(i-1))
1517  allu0(kk1-35)= al(kk2-35)
1518  allu0(kk1-34)= al(kk2-34)
1519  allu0(kk1-33)= al(kk2-33)
1520  allu0(kk1-32)= al(kk2-32)
1521  allu0(kk1-31)= al(kk2-31)
1522  allu0(kk1-30)= al(kk2-30)
1523  allu0(kk1-29)= al(kk2-29)
1524  allu0(kk1-28)= al(kk2-28)
1525  allu0(kk1-27)= al(kk2-27)
1526  allu0(kk1-26)= al(kk2-26)
1527  allu0(kk1-25)= al(kk2-25)
1528  allu0(kk1-24)= al(kk2-24)
1529  allu0(kk1-23)= al(kk2-23)
1530  allu0(kk1-22)= al(kk2-22)
1531  allu0(kk1-21)= al(kk2-21)
1532  allu0(kk1-20)= al(kk2-20)
1533  allu0(kk1-19)= al(kk2-19)
1534  allu0(kk1-18)= al(kk2-18)
1535  allu0(kk1-17)= al(kk2-17)
1536  allu0(kk1-16)= al(kk2-16)
1537  allu0(kk1-15)= al(kk2-15)
1538  allu0(kk1-14)= al(kk2-14)
1539  allu0(kk1-13)= al(kk2-13)
1540  allu0(kk1-12)= al(kk2-12)
1541  allu0(kk1-11)= al(kk2-11)
1542  allu0(kk1-10)= al(kk2-10)
1543  allu0(kk1-9 )= al(kk2-9 )
1544  allu0(kk1-8 )= al(kk2-8 )
1545  allu0(kk1-7 )= al(kk2-7 )
1546  allu0(kk1-6 )= al(kk2-6 )
1547  allu0(kk1-5 )= al(kk2-5 )
1548  allu0(kk1-4 )= al(kk2-4 )
1549  allu0(kk1-3 )= al(kk2-3 )
1550  allu0(kk1-2 )= al(kk2-2 )
1551  allu0(kk1-1 )= al(kk2-1 )
1552  allu0(kk1 )= al(kk2 )
1553  endif
1554  enddo
1555 
1556  icou= 0
1557  do k= inl(i-1)+1, inl(i)
1558  icou = icou + 1
1559  iw1(icou)= 0
1560  enddo
1561 
1562  icou= 0
1563  do k= inumfi2l(i-1)+1, inumfi2l(i)
1564  icou = icou + 1
1565  iw1(icou+icoul1)= 1
1566  enddo
1567 
1568  icou= 0
1569  do k= inumfi1l(i-1)+1, inumfi1l(i)
1570  icou = icou + 1
1571  iw1(icou+icoul2)= 2
1572  enddo
1573 
1574  do k= 1, icoul3
1575  iconfi1l(k+isl)= iw1(iw2(k))
1576  enddo
1577  !C
1578  !C-- UPPER part
1579  icou= 0
1580  do k= inu(i-1)+1, inu(i)
1581  icou = icou + 1
1582  iw1(icou)= iau(k)
1583  enddo
1584 
1585  icou= 0
1586  do k= inumfi2u(i-1)+1, inumfi2u(i)
1587  icou = icou + 1
1588  iw1(icou+icouu1)= fi2u(k)
1589  enddo
1590 
1591  icou= 0
1592  do k= inumfi1u(i-1)+1, inumfi1u(i)
1593  icou = icou + 1
1594  iw1(icou+icouu2)= fi1u(icou+icouu2+isu)
1595  enddo
1596 
1597  do k= 1, icouu3
1598  iw2(k)= k
1599  enddo
1600  call fill_in_s66_sort (iw1, iw2, icouu3, np)
1601 
1602  do k= 1, icouu3
1603  fi1u(k+isu)= iw1(k)
1604  ik= iw2(k)
1605  if (ik.le.inu(i)-inu(i-1)) then
1606  kk1= 36*( k+isu)
1607  kk2= 36*(ik+inu(i-1))
1608  aulu0(kk1-35)= au(kk2-35)
1609  aulu0(kk1-34)= au(kk2-34)
1610  aulu0(kk1-33)= au(kk2-33)
1611  aulu0(kk1-32)= au(kk2-32)
1612  aulu0(kk1-31)= au(kk2-31)
1613  aulu0(kk1-30)= au(kk2-30)
1614  aulu0(kk1-29)= au(kk2-29)
1615  aulu0(kk1-28)= au(kk2-28)
1616  aulu0(kk1-27)= au(kk2-27)
1617  aulu0(kk1-26)= au(kk2-26)
1618  aulu0(kk1-25)= au(kk2-25)
1619  aulu0(kk1-24)= au(kk2-24)
1620  aulu0(kk1-23)= au(kk2-23)
1621  aulu0(kk1-22)= au(kk2-22)
1622  aulu0(kk1-21)= au(kk2-21)
1623  aulu0(kk1-20)= au(kk2-20)
1624  aulu0(kk1-19)= au(kk2-19)
1625  aulu0(kk1-18)= au(kk2-18)
1626  aulu0(kk1-17)= au(kk2-17)
1627  aulu0(kk1-16)= au(kk2-16)
1628  aulu0(kk1-15)= au(kk2-15)
1629  aulu0(kk1-14)= au(kk2-14)
1630  aulu0(kk1-13)= au(kk2-13)
1631  aulu0(kk1-12)= au(kk2-12)
1632  aulu0(kk1-11)= au(kk2-11)
1633  aulu0(kk1-10)= au(kk2-10)
1634  aulu0(kk1-9 )= au(kk2-9 )
1635  aulu0(kk1-8 )= au(kk2-8 )
1636  aulu0(kk1-7 )= au(kk2-7 )
1637  aulu0(kk1-6 )= au(kk2-6 )
1638  aulu0(kk1-5 )= au(kk2-5 )
1639  aulu0(kk1-4 )= au(kk2-4 )
1640  aulu0(kk1-3 )= au(kk2-3 )
1641  aulu0(kk1-2 )= au(kk2-2 )
1642  aulu0(kk1-1 )= au(kk2-1 )
1643  aulu0(kk1 )= au(kk2 )
1644  endif
1645  enddo
1646 
1647  icou= 0
1648  do k= inu(i-1)+1, inu(i)
1649  icou = icou + 1
1650  iw1(icou)= 0
1651  enddo
1652 
1653  icou= 0
1654  do k= inumfi2u(i-1)+1, inumfi2u(i)
1655  icou = icou + 1
1656  iw1(icou+icouu1)= 1
1657  enddo
1658 
1659  icou= 0
1660  do k= inumfi1u(i-1)+1, inumfi1u(i)
1661  icou = icou + 1
1662  iw1(icou+icouu2)= 2
1663  enddo
1664 
1665  do k= 1, icouu3
1666  iconfi1u(k+isu)= iw1(iw2(k))
1667  enddo
1668 
1669  isl= isl + icoul3
1670  isu= isu + icouu3
1671  enddo
1672  !C===
1673  do i= 1, np
1674  inumfi1l(i)= iwsl(i)
1675  inumfi1u(i)= iwsu(i)
1676  enddo
1677 
1678  deallocate (iwsl, iwsu)
1679  deallocate (inumfi2l, inumfi2u)
1680  deallocate ( fi2l, fi2u)
1681 
1682  !C
1683  !C +----------------------+
1684  !C | ILU(2) factorization |
1685  !C +----------------------+
1686  !C===
1687  allocate (dlu0(36*np))
1688  dlu0= d
1689  do i=1,np
1690  dlu0(36*i-35)=dlu0(36*i-35)*sigma_diag
1691  dlu0(36*i-28)=dlu0(36*i-28)*sigma_diag
1692  dlu0(36*i-21)=dlu0(36*i-21)*sigma_diag
1693  dlu0(36*i-14)=dlu0(36*i-14)*sigma_diag
1694  dlu0(36*i-7 )=dlu0(36*i-7 )*sigma_diag
1695  dlu0(36*i )=dlu0(36*i )*sigma_diag
1696  enddo
1697 
1698  do i= 2, np
1699  iw1= 0
1700  iw2= 0
1701 
1702  do k= inumfi1l(i-1)+1, inumfi1l(i)
1703  iw1(fi1l(k))= k
1704  enddo
1705 
1706  do k= inumfi1u(i-1)+1, inumfi1u(i)
1707  iw2(fi1u(k))= k
1708  enddo
1709 
1710  do kk= inumfi1l(i-1)+1, inumfi1l(i)
1711  k= fi1l(kk)
1712  iconik= iconfi1l(kk)
1713 
1714  d11= dlu0(36*k-35)
1715  d12= dlu0(36*k-34)
1716  d13= dlu0(36*k-33)
1717  d14= dlu0(36*k-32)
1718  d15= dlu0(36*k-31)
1719  d16= dlu0(36*k-30)
1720  d21= dlu0(36*k-29)
1721  d22= dlu0(36*k-28)
1722  d23= dlu0(36*k-27)
1723  d24= dlu0(36*k-26)
1724  d25= dlu0(36*k-25)
1725  d26= dlu0(36*k-24)
1726  d31= dlu0(36*k-23)
1727  d32= dlu0(36*k-22)
1728  d33= dlu0(36*k-21)
1729  d34= dlu0(36*k-20)
1730  d35= dlu0(36*k-19)
1731  d36= dlu0(36*k-18)
1732  d41= dlu0(36*k-17)
1733  d42= dlu0(36*k-16)
1734  d43= dlu0(36*k-15)
1735  d44= dlu0(36*k-14)
1736  d45= dlu0(36*k-13)
1737  d46= dlu0(36*k-12)
1738  d51= dlu0(36*k-11)
1739  d52= dlu0(36*k-10)
1740  d53= dlu0(36*k-9 )
1741  d54= dlu0(36*k-8 )
1742  d55= dlu0(36*k-7 )
1743  d56= dlu0(36*k-6 )
1744  d61= dlu0(36*k-5 )
1745  d62= dlu0(36*k-4 )
1746  d63= dlu0(36*k-3 )
1747  d64= dlu0(36*k-2 )
1748  d65= dlu0(36*k-1 )
1749  d66= dlu0(36*k )
1750 
1751  call ilu1a66 (dkinv,d11,d12,d13,d14,d15,d16,d21,d22,d23,d24,d25,d26, &
1752  & d31,d32,d33,d34,d35,d36,d41,d42,d43,d44,d45,d46,d51,d52,d53,d54,d55,d56, &
1753  & d61,d62,d63,d64,d65,d66)
1754 
1755  aik(1,1)= allu0(36*kk-35)
1756  aik(1,2)= allu0(36*kk-34)
1757  aik(1,3)= allu0(36*kk-33)
1758  aik(1,4)= allu0(36*kk-32)
1759  aik(1,5)= allu0(36*kk-31)
1760  aik(1,6)= allu0(36*kk-30)
1761  aik(2,1)= allu0(36*kk-29)
1762  aik(2,2)= allu0(36*kk-28)
1763  aik(2,3)= allu0(36*kk-27)
1764  aik(2,4)= allu0(36*kk-26)
1765  aik(2,5)= allu0(36*kk-25)
1766  aik(2,6)= allu0(36*kk-24)
1767  aik(3,1)= allu0(36*kk-23)
1768  aik(3,2)= allu0(36*kk-22)
1769  aik(3,3)= allu0(36*kk-21)
1770  aik(3,4)= allu0(36*kk-20)
1771  aik(3,5)= allu0(36*kk-19)
1772  aik(3,6)= allu0(36*kk-18)
1773  aik(4,1)= allu0(36*kk-17)
1774  aik(4,2)= allu0(36*kk-16)
1775  aik(4,3)= allu0(36*kk-15)
1776  aik(4,4)= allu0(36*kk-14)
1777  aik(4,5)= allu0(36*kk-13)
1778  aik(4,6)= allu0(36*kk-12)
1779  aik(5,1)= allu0(36*kk-11)
1780  aik(5,2)= allu0(36*kk-10)
1781  aik(5,3)= allu0(36*kk-9)
1782  aik(5,4)= allu0(36*kk-8)
1783  aik(5,5)= allu0(36*kk-7)
1784  aik(5,6)= allu0(36*kk-6)
1785  aik(6,1)= allu0(36*kk-5)
1786  aik(6,2)= allu0(36*kk-4)
1787  aik(6,3)= allu0(36*kk-3)
1788  aik(6,4)= allu0(36*kk-2)
1789  aik(6,5)= allu0(36*kk-1)
1790  aik(6,6)= allu0(36*kk )
1791 
1792  do jj= inumfi1u(k-1)+1, inumfi1u(k)
1793  j= fi1u(jj)
1794  iconkj= iconfi1u(jj)
1795 
1796  if ((iconik+iconkj).lt.2) then
1797  akj(1,1)= aulu0(36*jj-35)
1798  akj(1,2)= aulu0(36*jj-34)
1799  akj(1,3)= aulu0(36*jj-33)
1800  akj(1,4)= aulu0(36*jj-32)
1801  akj(1,5)= aulu0(36*jj-31)
1802  akj(1,6)= aulu0(36*jj-30)
1803  akj(2,1)= aulu0(36*jj-29)
1804  akj(2,2)= aulu0(36*jj-28)
1805  akj(2,3)= aulu0(36*jj-27)
1806  akj(2,4)= aulu0(36*jj-26)
1807  akj(2,5)= aulu0(36*jj-25)
1808  akj(2,6)= aulu0(36*jj-24)
1809  akj(3,1)= aulu0(36*jj-23)
1810  akj(3,2)= aulu0(36*jj-22)
1811  akj(3,3)= aulu0(36*jj-21)
1812  akj(3,4)= aulu0(36*jj-20)
1813  akj(3,5)= aulu0(36*jj-19)
1814  akj(3,6)= aulu0(36*jj-18)
1815  akj(4,1)= aulu0(36*jj-17)
1816  akj(4,2)= aulu0(36*jj-16)
1817  akj(4,3)= aulu0(36*jj-15)
1818  akj(4,4)= aulu0(36*jj-14)
1819  akj(4,5)= aulu0(36*jj-13)
1820  akj(4,6)= aulu0(36*jj-12)
1821  akj(5,1)= aulu0(36*jj-11)
1822  akj(5,2)= aulu0(36*jj-10)
1823  akj(5,3)= aulu0(36*jj-9)
1824  akj(5,4)= aulu0(36*jj-8)
1825  akj(5,5)= aulu0(36*jj-7)
1826  akj(5,6)= aulu0(36*jj-6)
1827  akj(6,1)= aulu0(36*jj-5)
1828  akj(6,2)= aulu0(36*jj-4)
1829  akj(6,3)= aulu0(36*jj-3)
1830  akj(6,4)= aulu0(36*jj-2)
1831  akj(6,5)= aulu0(36*jj-1)
1832  akj(6,6)= aulu0(36*jj )
1833 
1834  call ilu1b66 (rhs_aij, dkinv, aik, akj)
1835 
1836  if (j.eq.i) then
1837  dlu0(36*i-35)= dlu0(36*i-35) - rhs_aij(1,1)
1838  dlu0(36*i-34)= dlu0(36*i-34) - rhs_aij(1,2)
1839  dlu0(36*i-33)= dlu0(36*i-33) - rhs_aij(1,3)
1840  dlu0(36*i-32)= dlu0(36*i-32) - rhs_aij(1,4)
1841  dlu0(36*i-31)= dlu0(36*i-31) - rhs_aij(1,5)
1842  dlu0(36*i-30)= dlu0(36*i-30) - rhs_aij(1,6)
1843  dlu0(36*i-29)= dlu0(36*i-29) - rhs_aij(2,1)
1844  dlu0(36*i-28)= dlu0(36*i-28) - rhs_aij(2,2)
1845  dlu0(36*i-27)= dlu0(36*i-27) - rhs_aij(2,3)
1846  dlu0(36*i-26)= dlu0(36*i-26) - rhs_aij(2,4)
1847  dlu0(36*i-25)= dlu0(36*i-25) - rhs_aij(2,5)
1848  dlu0(36*i-24)= dlu0(36*i-24) - rhs_aij(2,6)
1849  dlu0(36*i-23)= dlu0(36*i-23) - rhs_aij(3,1)
1850  dlu0(36*i-22)= dlu0(36*i-22) - rhs_aij(3,2)
1851  dlu0(36*i-21)= dlu0(36*i-21) - rhs_aij(3,3)
1852  dlu0(36*i-20)= dlu0(36*i-20) - rhs_aij(3,4)
1853  dlu0(36*i-19)= dlu0(36*i-19) - rhs_aij(3,5)
1854  dlu0(36*i-18)= dlu0(36*i-18) - rhs_aij(3,6)
1855  dlu0(36*i-17)= dlu0(36*i-17) - rhs_aij(4,1)
1856  dlu0(36*i-16)= dlu0(36*i-16) - rhs_aij(4,2)
1857  dlu0(36*i-15)= dlu0(36*i-15) - rhs_aij(4,3)
1858  dlu0(36*i-14)= dlu0(36*i-14) - rhs_aij(4,4)
1859  dlu0(36*i-13)= dlu0(36*i-13) - rhs_aij(4,5)
1860  dlu0(36*i-12)= dlu0(36*i-12) - rhs_aij(4,6)
1861  dlu0(36*i-11)= dlu0(36*i-11) - rhs_aij(5,1)
1862  dlu0(36*i-10)= dlu0(36*i-10) - rhs_aij(5,2)
1863  dlu0(36*i-9 )= dlu0(36*i-9 ) - rhs_aij(5,3)
1864  dlu0(36*i-8 )= dlu0(36*i-8 ) - rhs_aij(5,4)
1865  dlu0(36*i-7 )= dlu0(36*i-7 ) - rhs_aij(5,5)
1866  dlu0(36*i-6 )= dlu0(36*i-6 ) - rhs_aij(5,6)
1867  dlu0(36*i-5 )= dlu0(36*i-5 ) - rhs_aij(6,1)
1868  dlu0(36*i-4 )= dlu0(36*i-4 ) - rhs_aij(6,2)
1869  dlu0(36*i-3 )= dlu0(36*i-3 ) - rhs_aij(6,3)
1870  dlu0(36*i-2 )= dlu0(36*i-2 ) - rhs_aij(6,4)
1871  dlu0(36*i-1 )= dlu0(36*i-1 ) - rhs_aij(6,5)
1872  dlu0(36*i )= dlu0(36*i ) - rhs_aij(6,6)
1873  endif
1874 
1875  if (j.lt.i) then
1876  ij0= iw1(j)
1877  allu0(36*ij0-35)= allu0(36*ij0-35) - rhs_aij(1,1)
1878  allu0(36*ij0-34)= allu0(36*ij0-34) - rhs_aij(1,2)
1879  allu0(36*ij0-33)= allu0(36*ij0-33) - rhs_aij(1,3)
1880  allu0(36*ij0-32)= allu0(36*ij0-32) - rhs_aij(1,4)
1881  allu0(36*ij0-31)= allu0(36*ij0-31) - rhs_aij(1,5)
1882  allu0(36*ij0-30)= allu0(36*ij0-30) - rhs_aij(1,6)
1883  allu0(36*ij0-29)= allu0(36*ij0-29) - rhs_aij(2,1)
1884  allu0(36*ij0-28)= allu0(36*ij0-28) - rhs_aij(2,2)
1885  allu0(36*ij0-27)= allu0(36*ij0-27) - rhs_aij(2,3)
1886  allu0(36*ij0-26)= allu0(36*ij0-26) - rhs_aij(2,4)
1887  allu0(36*ij0-25)= allu0(36*ij0-25) - rhs_aij(2,5)
1888  allu0(36*ij0-24)= allu0(36*ij0-24) - rhs_aij(2,6)
1889  allu0(36*ij0-23)= allu0(36*ij0-23) - rhs_aij(3,1)
1890  allu0(36*ij0-22)= allu0(36*ij0-22) - rhs_aij(3,2)
1891  allu0(36*ij0-21)= allu0(36*ij0-21) - rhs_aij(3,3)
1892  allu0(36*ij0-20)= allu0(36*ij0-20) - rhs_aij(3,4)
1893  allu0(36*ij0-19)= allu0(36*ij0-19) - rhs_aij(3,5)
1894  allu0(36*ij0-18)= allu0(36*ij0-18) - rhs_aij(3,6)
1895  allu0(36*ij0-17)= allu0(36*ij0-17) - rhs_aij(4,1)
1896  allu0(36*ij0-16)= allu0(36*ij0-16) - rhs_aij(4,2)
1897  allu0(36*ij0-15)= allu0(36*ij0-15) - rhs_aij(4,3)
1898  allu0(36*ij0-14)= allu0(36*ij0-14) - rhs_aij(4,4)
1899  allu0(36*ij0-13)= allu0(36*ij0-13) - rhs_aij(4,5)
1900  allu0(36*ij0-12)= allu0(36*ij0-12) - rhs_aij(4,6)
1901  allu0(36*ij0-11)= allu0(36*ij0-11) - rhs_aij(5,1)
1902  allu0(36*ij0-10)= allu0(36*ij0-10) - rhs_aij(5,2)
1903  allu0(36*ij0-9 )= allu0(36*ij0-9 ) - rhs_aij(5,3)
1904  allu0(36*ij0-8 )= allu0(36*ij0-8 ) - rhs_aij(5,4)
1905  allu0(36*ij0-7 )= allu0(36*ij0-7 ) - rhs_aij(5,5)
1906  allu0(36*ij0-6 )= allu0(36*ij0-6 ) - rhs_aij(5,6)
1907  allu0(36*ij0-5 )= allu0(36*ij0-5 ) - rhs_aij(6,1)
1908  allu0(36*ij0-4 )= allu0(36*ij0-4 ) - rhs_aij(6,2)
1909  allu0(36*ij0-3 )= allu0(36*ij0-3 ) - rhs_aij(6,3)
1910  allu0(36*ij0-2 )= allu0(36*ij0-2 ) - rhs_aij(6,4)
1911  allu0(36*ij0-1 )= allu0(36*ij0-1 ) - rhs_aij(6,5)
1912  allu0(36*ij0 )= allu0(36*ij0 ) - rhs_aij(6,6)
1913  endif
1914 
1915  if (j.gt.i) then
1916  ij0= iw2(j)
1917  aulu0(36*ij0-35)= aulu0(36*ij0-35) - rhs_aij(1,1)
1918  aulu0(36*ij0-34)= aulu0(36*ij0-34) - rhs_aij(1,2)
1919  aulu0(36*ij0-33)= aulu0(36*ij0-33) - rhs_aij(1,3)
1920  aulu0(36*ij0-32)= aulu0(36*ij0-32) - rhs_aij(1,4)
1921  aulu0(36*ij0-31)= aulu0(36*ij0-31) - rhs_aij(1,5)
1922  aulu0(36*ij0-30)= aulu0(36*ij0-30) - rhs_aij(1,6)
1923  aulu0(36*ij0-29)= aulu0(36*ij0-29) - rhs_aij(2,1)
1924  aulu0(36*ij0-28)= aulu0(36*ij0-28) - rhs_aij(2,2)
1925  aulu0(36*ij0-27)= aulu0(36*ij0-27) - rhs_aij(2,3)
1926  aulu0(36*ij0-26)= aulu0(36*ij0-26) - rhs_aij(2,4)
1927  aulu0(36*ij0-25)= aulu0(36*ij0-25) - rhs_aij(2,5)
1928  aulu0(36*ij0-24)= aulu0(36*ij0-24) - rhs_aij(2,6)
1929  aulu0(36*ij0-23)= aulu0(36*ij0-23) - rhs_aij(3,1)
1930  aulu0(36*ij0-22)= aulu0(36*ij0-22) - rhs_aij(3,2)
1931  aulu0(36*ij0-21)= aulu0(36*ij0-21) - rhs_aij(3,3)
1932  aulu0(36*ij0-20)= aulu0(36*ij0-20) - rhs_aij(3,4)
1933  aulu0(36*ij0-19)= aulu0(36*ij0-19) - rhs_aij(3,5)
1934  aulu0(36*ij0-18)= aulu0(36*ij0-18) - rhs_aij(3,6)
1935  aulu0(36*ij0-17)= aulu0(36*ij0-17) - rhs_aij(4,1)
1936  aulu0(36*ij0-16)= aulu0(36*ij0-16) - rhs_aij(4,2)
1937  aulu0(36*ij0-15)= aulu0(36*ij0-15) - rhs_aij(4,3)
1938  aulu0(36*ij0-14)= aulu0(36*ij0-14) - rhs_aij(4,4)
1939  aulu0(36*ij0-13)= aulu0(36*ij0-13) - rhs_aij(4,5)
1940  aulu0(36*ij0-12)= aulu0(36*ij0-12) - rhs_aij(4,6)
1941  aulu0(36*ij0-11)= aulu0(36*ij0-11) - rhs_aij(5,1)
1942  aulu0(36*ij0-10)= aulu0(36*ij0-10) - rhs_aij(5,2)
1943  aulu0(36*ij0-9 )= aulu0(36*ij0-9 ) - rhs_aij(5,3)
1944  aulu0(36*ij0-8 )= aulu0(36*ij0-8 ) - rhs_aij(5,4)
1945  aulu0(36*ij0-7 )= aulu0(36*ij0-7 ) - rhs_aij(5,5)
1946  aulu0(36*ij0-6 )= aulu0(36*ij0-6 ) - rhs_aij(5,6)
1947  aulu0(36*ij0-5 )= aulu0(36*ij0-5 ) - rhs_aij(6,1)
1948  aulu0(36*ij0-4 )= aulu0(36*ij0-4 ) - rhs_aij(6,2)
1949  aulu0(36*ij0-3 )= aulu0(36*ij0-3 ) - rhs_aij(6,3)
1950  aulu0(36*ij0-2 )= aulu0(36*ij0-2 ) - rhs_aij(6,4)
1951  aulu0(36*ij0-1 )= aulu0(36*ij0-1 ) - rhs_aij(6,5)
1952  aulu0(36*ij0 )= aulu0(36*ij0 ) - rhs_aij(6,6)
1953  endif
1954  endif
1955  enddo
1956  enddo
1957  enddo
1958 
1959  deallocate (iw1, iw2)
1960  deallocate (iconfi1l, iconfi1u)
1961  !C===
1962  end subroutine form_ilu2_66
1963 
1964 
1965  !C
1966  !C***
1967  !C*** fill_in_S66_SORT
1968  !C***
1969  !C
1970  subroutine fill_in_s66_sort (STEM, INUM, N, NP)
1971  use hecmw_util
1972  implicit none
1973  integer(kind=kint) :: n, np
1974  integer(kind=kint), dimension(NP) :: stem
1975  integer(kind=kint), dimension(NP) :: inum
1976  integer(kind=kint), dimension(:), allocatable :: istack
1977  integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
1978 
1979  allocate (istack(-np:+np))
1980 
1981  m = 100
1982  nstack= np
1983 
1984  jstack= 0
1985  l = 1
1986  ir = n
1987 
1988  ip= 0
1989  1 continue
1990  ip= ip + 1
1991 
1992  if (ir-l.lt.m) then
1993  do j= l+1, ir
1994  ss= stem(j)
1995  ii= inum(j)
1996 
1997  do i= j-1,1,-1
1998  if (stem(i).le.ss) goto 2
1999  stem(i+1)= stem(i)
2000  inum(i+1)= inum(i)
2001  end do
2002  i= 0
2003 
2004  2 continue
2005  stem(i+1)= ss
2006  inum(i+1)= ii
2007  end do
2008 
2009  if (jstack.eq.0) then
2010  deallocate (istack)
2011  return
2012  endif
2013 
2014  ir = istack(jstack)
2015  l = istack(jstack-1)
2016  jstack= jstack - 2
2017  else
2018 
2019  k= (l+ir) / 2
2020  temp = stem(k)
2021  stem(k) = stem(l+1)
2022  stem(l+1)= temp
2023 
2024  it = inum(k)
2025  inum(k) = inum(l+1)
2026  inum(l+1)= it
2027 
2028  if (stem(l+1).gt.stem(ir)) then
2029  temp = stem(l+1)
2030  stem(l+1)= stem(ir)
2031  stem(ir )= temp
2032  it = inum(l+1)
2033  inum(l+1)= inum(ir)
2034  inum(ir )= it
2035  endif
2036 
2037  if (stem(l).gt.stem(ir)) then
2038  temp = stem(l)
2039  stem(l )= stem(ir)
2040  stem(ir)= temp
2041  it = inum(l)
2042  inum(l )= inum(ir)
2043  inum(ir)= it
2044  endif
2045 
2046  if (stem(l+1).gt.stem(l)) then
2047  temp = stem(l+1)
2048  stem(l+1)= stem(l)
2049  stem(l )= temp
2050  it = inum(l+1)
2051  inum(l+1)= inum(l)
2052  inum(l )= it
2053  endif
2054 
2055  i= l + 1
2056  j= ir
2057 
2058  ss= stem(l)
2059  ii= inum(l)
2060 
2061  3 continue
2062  i= i + 1
2063  if (stem(i).lt.ss) goto 3
2064 
2065  4 continue
2066  j= j - 1
2067  if (stem(j).gt.ss) goto 4
2068 
2069  if (j.lt.i) goto 5
2070 
2071  temp = stem(i)
2072  stem(i)= stem(j)
2073  stem(j)= temp
2074 
2075  it = inum(i)
2076  inum(i)= inum(j)
2077  inum(j)= it
2078 
2079  goto 3
2080 
2081  5 continue
2082 
2083  stem(l)= stem(j)
2084  stem(j)= ss
2085  inum(l)= inum(j)
2086  inum(j)= ii
2087 
2088  jstack= jstack + 2
2089 
2090  if (jstack.gt.nstack) then
2091  write (*,*) 'NSTACK overflow'
2092  stop
2093  endif
2094 
2095  if (ir-i+1.ge.j-1) then
2096  istack(jstack )= ir
2097  istack(jstack-1)= i
2098  ir= j-1
2099  else
2100  istack(jstack )= j-1
2101  istack(jstack-1)= l
2102  l= i
2103  endif
2104 
2105  endif
2106 
2107  goto 1
2108 
2109  end subroutine fill_in_s66_sort
2110 
2111  !C
2112  !C***
2113  !C*** ILU1a66
2114  !C***
2115  !C
2116  !C computes LU factorization of 3*3 Diagonal Block
2117  !C
2118  subroutine ilu1a66 (ALU, D11,D12,D13,D14,D15,D16,D21,D22,D23,D24,D25,D26, &
2119  & D31,D32,D33,D34,D35,D36,D41,D42,D43,D44,D45,D46,D51,D52,D53,D54,D55,D56, &
2120  & D61,D62,D63,D64,D65,D66)
2121  use hecmw_util
2122  implicit none
2123  real(kind=kreal) :: alu(6,6), pw(6)
2124  real(kind=kreal) :: d11,d12,d13,d14,d15,d16
2125  real(kind=kreal) :: d21,d22,d23,d24,d25,d26
2126  real(kind=kreal) :: d31,d32,d33,d34,d35,d36
2127  real(kind=kreal) :: d41,d42,d43,d44,d45,d46
2128  real(kind=kreal) :: d51,d52,d53,d54,d55,d56
2129  real(kind=kreal) :: d61,d62,d63,d64,d65,d66
2130  integer(kind=kint) :: i,j,k
2131 
2132  alu(1,1)= d11
2133  alu(1,2)= d12
2134  alu(1,3)= d13
2135  alu(1,4)= d14
2136  alu(1,5)= d15
2137  alu(1,6)= d16
2138  alu(2,1)= d21
2139  alu(2,2)= d22
2140  alu(2,3)= d23
2141  alu(2,4)= d24
2142  alu(2,5)= d25
2143  alu(2,6)= d26
2144  alu(3,1)= d31
2145  alu(3,2)= d32
2146  alu(3,3)= d33
2147  alu(3,4)= d34
2148  alu(3,5)= d35
2149  alu(3,6)= d36
2150  alu(4,1)= d41
2151  alu(4,2)= d42
2152  alu(4,3)= d43
2153  alu(4,4)= d44
2154  alu(4,5)= d45
2155  alu(4,6)= d46
2156  alu(5,1)= d51
2157  alu(5,2)= d52
2158  alu(5,3)= d53
2159  alu(5,4)= d54
2160  alu(5,5)= d55
2161  alu(5,6)= d56
2162  alu(6,1)= d61
2163  alu(6,2)= d62
2164  alu(6,3)= d63
2165  alu(6,4)= d64
2166  alu(6,5)= d65
2167  alu(6,6)= d66
2168 
2169  do k= 1, 6
2170  alu(k,k)= 1.d0/alu(k,k)
2171  do i= k+1, 6
2172  alu(i,k)= alu(i,k) * alu(k,k)
2173  do j= k+1, 6
2174  pw(j)= alu(i,j) - alu(i,k)*alu(k,j)
2175  enddo
2176  do j= k+1, 6
2177  alu(i,j)= pw(j)
2178  enddo
2179  enddo
2180  enddo
2181 
2182  return
2183  end subroutine ilu1a66
2184 
2185  !C
2186  !C***
2187  !C*** ILU1b66
2188  !C***
2189  !C
2190  !C computes L_ik * D_k_INV * U_kj at ILU factorization
2191  !C for 3*3 Block Type Matrix
2192  !C
2193  subroutine ilu1b66 (RHS_Aij, DkINV, Aik, Akj)
2194  use hecmw_util
2195  implicit none
2196  real(kind=kreal) :: rhs_aij(6,6), dkinv(6,6), aik(6,6), akj(6,6)
2197  real(kind=kreal) :: x1,x2,x3,x4,x5,x6
2198 
2199  !C
2200  !C-- 1st Col.
2201  x1= akj(1,1)
2202  x2= akj(2,1)
2203  x3= akj(3,1)
2204  x4= akj(4,1)
2205  x5= akj(5,1)
2206  x6= akj(6,1)
2207 
2208  x2= x2 -dkinv(2,1)*x1
2209  x3= x3 -dkinv(3,1)*x1 -dkinv(3,2)*x2
2210  x4= x4 -dkinv(4,1)*x1 -dkinv(4,2)*x2 -dkinv(4,3)*x3
2211  x5= x5 -dkinv(5,1)*x1 -dkinv(5,2)*x2 -dkinv(5,3)*x3 -dkinv(5,4)*x4
2212  x6= x6 -dkinv(6,1)*x1 -dkinv(6,2)*x2 -dkinv(6,3)*x3 -dkinv(6,4)*x4 -dkinv(6,5)*x5
2213 
2214  x6= dkinv(6,6)* x6
2215  x5= dkinv(5,5)*( x5 -dkinv(5,6)*x6 )
2216  x4= dkinv(4,4)*( x4 -dkinv(4,6)*x6 -dkinv(4,5)*x5)
2217  x3= dkinv(3,3)*( x3 -dkinv(3,6)*x6 -dkinv(3,5)*x5 -dkinv(3,4)*x4)
2218  x2= dkinv(2,2)*( x2 -dkinv(2,6)*x6 -dkinv(2,5)*x5 -dkinv(2,4)*x4 -dkinv(2,3)*x3)
2219  x1= dkinv(1,1)*( x1 -dkinv(1,6)*x6 -dkinv(1,5)*x5 -dkinv(1,4)*x4 -dkinv(1,3)*x3 -dkinv(1,2)*x2)
2220 
2221  rhs_aij(1,1)= aik(1,1)*x1 +aik(1,2)*x2 +aik(1,3)*x3 +aik(1,4)*x4 +aik(1,5)*x5 +aik(1,6)*x6
2222  rhs_aij(2,1)= aik(2,1)*x1 +aik(2,2)*x2 +aik(2,3)*x3 +aik(2,4)*x4 +aik(2,5)*x5 +aik(2,6)*x6
2223  rhs_aij(3,1)= aik(3,1)*x1 +aik(3,2)*x2 +aik(3,3)*x3 +aik(3,4)*x4 +aik(3,5)*x5 +aik(3,6)*x6
2224  rhs_aij(4,1)= aik(4,1)*x1 +aik(4,2)*x2 +aik(4,3)*x3 +aik(4,4)*x4 +aik(4,5)*x5 +aik(4,6)*x6
2225  rhs_aij(5,1)= aik(5,1)*x1 +aik(5,2)*x2 +aik(5,3)*x3 +aik(5,4)*x4 +aik(5,5)*x5 +aik(5,6)*x6
2226  rhs_aij(6,1)= aik(6,1)*x1 +aik(6,2)*x2 +aik(6,3)*x3 +aik(6,4)*x4 +aik(6,5)*x5 +aik(6,6)*x6
2227 
2228  !C
2229  !C-- 2nd Col.
2230  x1= akj(1,2)
2231  x2= akj(2,2)
2232  x3= akj(3,2)
2233  x4= akj(4,2)
2234  x5= akj(5,2)
2235  x6= akj(6,2)
2236 
2237  x2= x2 -dkinv(2,1)*x1
2238  x3= x3 -dkinv(3,1)*x1 -dkinv(3,2)*x2
2239  x4= x4 -dkinv(4,1)*x1 -dkinv(4,2)*x2 -dkinv(4,3)*x3
2240  x5= x5 -dkinv(5,1)*x1 -dkinv(5,2)*x2 -dkinv(5,3)*x3 -dkinv(5,4)*x4
2241  x6= x6 -dkinv(6,1)*x1 -dkinv(6,2)*x2 -dkinv(6,3)*x3 -dkinv(6,4)*x4 -dkinv(6,5)*x5
2242 
2243  x6= dkinv(6,6)* x6
2244  x5= dkinv(5,5)*( x5 -dkinv(5,6)*x6 )
2245  x4= dkinv(4,4)*( x4 -dkinv(4,6)*x6 -dkinv(4,5)*x5)
2246  x3= dkinv(3,3)*( x3 -dkinv(3,6)*x6 -dkinv(3,5)*x5 -dkinv(3,4)*x4)
2247  x2= dkinv(2,2)*( x2 -dkinv(2,6)*x6 -dkinv(2,5)*x5 -dkinv(2,4)*x4 -dkinv(2,3)*x3)
2248  x1= dkinv(1,1)*( x1 -dkinv(1,6)*x6 -dkinv(1,5)*x5 -dkinv(1,4)*x4 -dkinv(1,3)*x3 -dkinv(1,2)*x2)
2249 
2250  rhs_aij(1,2)= aik(1,1)*x1 +aik(1,2)*x2 +aik(1,3)*x3 +aik(1,4)*x4 +aik(1,5)*x5 +aik(1,6)*x6
2251  rhs_aij(2,2)= aik(2,1)*x1 +aik(2,2)*x2 +aik(2,3)*x3 +aik(2,4)*x4 +aik(2,5)*x5 +aik(2,6)*x6
2252  rhs_aij(3,2)= aik(3,1)*x1 +aik(3,2)*x2 +aik(3,3)*x3 +aik(3,4)*x4 +aik(3,5)*x5 +aik(3,6)*x6
2253  rhs_aij(4,2)= aik(4,1)*x1 +aik(4,2)*x2 +aik(4,3)*x3 +aik(4,4)*x4 +aik(4,5)*x5 +aik(4,6)*x6
2254  rhs_aij(5,2)= aik(5,1)*x1 +aik(5,2)*x2 +aik(5,3)*x3 +aik(5,4)*x4 +aik(5,5)*x5 +aik(5,6)*x6
2255  rhs_aij(6,2)= aik(6,1)*x1 +aik(6,2)*x2 +aik(6,3)*x3 +aik(6,4)*x4 +aik(6,5)*x5 +aik(6,6)*x6
2256 
2257  !C
2258  !C-- 3rd Col.
2259  x1= akj(1,3)
2260  x2= akj(2,3)
2261  x3= akj(3,3)
2262  x4= akj(4,3)
2263  x5= akj(5,3)
2264  x6= akj(6,3)
2265 
2266  x2= x2 -dkinv(2,1)*x1
2267  x3= x3 -dkinv(3,1)*x1 -dkinv(3,2)*x2
2268  x4= x4 -dkinv(4,1)*x1 -dkinv(4,2)*x2 -dkinv(4,3)*x3
2269  x5= x5 -dkinv(5,1)*x1 -dkinv(5,2)*x2 -dkinv(5,3)*x3 -dkinv(5,4)*x4
2270  x6= x6 -dkinv(6,1)*x1 -dkinv(6,2)*x2 -dkinv(6,3)*x3 -dkinv(6,4)*x4 -dkinv(6,5)*x5
2271 
2272  x6= dkinv(6,6)* x6
2273  x5= dkinv(5,5)*( x5 -dkinv(5,6)*x6 )
2274  x4= dkinv(4,4)*( x4 -dkinv(4,6)*x6 -dkinv(4,5)*x5)
2275  x3= dkinv(3,3)*( x3 -dkinv(3,6)*x6 -dkinv(3,5)*x5 -dkinv(3,4)*x4)
2276  x2= dkinv(2,2)*( x2 -dkinv(2,6)*x6 -dkinv(2,5)*x5 -dkinv(2,4)*x4 -dkinv(2,3)*x3)
2277  x1= dkinv(1,1)*( x1 -dkinv(1,6)*x6 -dkinv(1,5)*x5 -dkinv(1,4)*x4 -dkinv(1,3)*x3 -dkinv(1,2)*x2)
2278 
2279  rhs_aij(1,3)= aik(1,1)*x1 +aik(1,2)*x2 +aik(1,3)*x3 +aik(1,4)*x4 +aik(1,5)*x5 +aik(1,6)*x6
2280  rhs_aij(2,3)= aik(2,1)*x1 +aik(2,2)*x2 +aik(2,3)*x3 +aik(2,4)*x4 +aik(2,5)*x5 +aik(2,6)*x6
2281  rhs_aij(3,3)= aik(3,1)*x1 +aik(3,2)*x2 +aik(3,3)*x3 +aik(3,4)*x4 +aik(3,5)*x5 +aik(3,6)*x6
2282  rhs_aij(4,3)= aik(4,1)*x1 +aik(4,2)*x2 +aik(4,3)*x3 +aik(4,4)*x4 +aik(4,5)*x5 +aik(4,6)*x6
2283  rhs_aij(5,3)= aik(5,1)*x1 +aik(5,2)*x2 +aik(5,3)*x3 +aik(5,4)*x4 +aik(5,5)*x5 +aik(5,6)*x6
2284  rhs_aij(6,3)= aik(6,1)*x1 +aik(6,2)*x2 +aik(6,3)*x3 +aik(6,4)*x4 +aik(6,5)*x5 +aik(6,6)*x6
2285 
2286  !C
2287  !C-- 4th Col.
2288  x1= akj(1,4)
2289  x2= akj(2,4)
2290  x3= akj(3,4)
2291  x4= akj(4,4)
2292  x5= akj(5,4)
2293  x6= akj(6,4)
2294 
2295  x2= x2 -dkinv(2,1)*x1
2296  x3= x3 -dkinv(3,1)*x1 -dkinv(3,2)*x2
2297  x4= x4 -dkinv(4,1)*x1 -dkinv(4,2)*x2 -dkinv(4,3)*x3
2298  x5= x5 -dkinv(5,1)*x1 -dkinv(5,2)*x2 -dkinv(5,3)*x3 -dkinv(5,4)*x4
2299  x6= x6 -dkinv(6,1)*x1 -dkinv(6,2)*x2 -dkinv(6,3)*x3 -dkinv(6,4)*x4 -dkinv(6,5)*x5
2300 
2301  x6= dkinv(6,6)* x6
2302  x5= dkinv(5,5)*( x5 -dkinv(5,6)*x6 )
2303  x4= dkinv(4,4)*( x4 -dkinv(4,6)*x6 -dkinv(4,5)*x5)
2304  x3= dkinv(3,3)*( x3 -dkinv(3,6)*x6 -dkinv(3,5)*x5 -dkinv(3,4)*x4)
2305  x2= dkinv(2,2)*( x2 -dkinv(2,6)*x6 -dkinv(2,5)*x5 -dkinv(2,4)*x4 -dkinv(2,3)*x3)
2306  x1= dkinv(1,1)*( x1 -dkinv(1,6)*x6 -dkinv(1,5)*x5 -dkinv(1,4)*x4 -dkinv(1,3)*x3 -dkinv(1,2)*x2)
2307 
2308  rhs_aij(1,4)= aik(1,1)*x1 +aik(1,2)*x2 +aik(1,3)*x3 +aik(1,4)*x4 +aik(1,5)*x5 +aik(1,6)*x6
2309  rhs_aij(2,4)= aik(2,1)*x1 +aik(2,2)*x2 +aik(2,3)*x3 +aik(2,4)*x4 +aik(2,5)*x5 +aik(2,6)*x6
2310  rhs_aij(3,4)= aik(3,1)*x1 +aik(3,2)*x2 +aik(3,3)*x3 +aik(3,4)*x4 +aik(3,5)*x5 +aik(3,6)*x6
2311  rhs_aij(4,4)= aik(4,1)*x1 +aik(4,2)*x2 +aik(4,3)*x3 +aik(4,4)*x4 +aik(4,5)*x5 +aik(4,6)*x6
2312  rhs_aij(5,4)= aik(5,1)*x1 +aik(5,2)*x2 +aik(5,3)*x3 +aik(5,4)*x4 +aik(5,5)*x5 +aik(5,6)*x6
2313  rhs_aij(6,4)= aik(6,1)*x1 +aik(6,2)*x2 +aik(6,3)*x3 +aik(6,4)*x4 +aik(6,5)*x5 +aik(6,6)*x6
2314 
2315  !C
2316  !C-- 5th Col.
2317  x1= akj(1,5)
2318  x2= akj(2,5)
2319  x3= akj(3,5)
2320  x4= akj(4,5)
2321  x5= akj(5,5)
2322  x6= akj(6,5)
2323 
2324  x2= x2 -dkinv(2,1)*x1
2325  x3= x3 -dkinv(3,1)*x1 -dkinv(3,2)*x2
2326  x4= x4 -dkinv(4,1)*x1 -dkinv(4,2)*x2 -dkinv(4,3)*x3
2327  x5= x5 -dkinv(5,1)*x1 -dkinv(5,2)*x2 -dkinv(5,3)*x3 -dkinv(5,4)*x4
2328  x6= x6 -dkinv(6,1)*x1 -dkinv(6,2)*x2 -dkinv(6,3)*x3 -dkinv(6,4)*x4 -dkinv(6,5)*x5
2329 
2330  x6= dkinv(6,6)* x6
2331  x5= dkinv(5,5)*( x5 -dkinv(5,6)*x6 )
2332  x4= dkinv(4,4)*( x4 -dkinv(4,6)*x6 -dkinv(4,5)*x5)
2333  x3= dkinv(3,3)*( x3 -dkinv(3,6)*x6 -dkinv(3,5)*x5 -dkinv(3,4)*x4)
2334  x2= dkinv(2,2)*( x2 -dkinv(2,6)*x6 -dkinv(2,5)*x5 -dkinv(2,4)*x4 -dkinv(2,3)*x3)
2335  x1= dkinv(1,1)*( x1 -dkinv(1,6)*x6 -dkinv(1,5)*x5 -dkinv(1,4)*x4 -dkinv(1,3)*x3 -dkinv(1,2)*x2)
2336 
2337  rhs_aij(1,5)= aik(1,1)*x1 +aik(1,2)*x2 +aik(1,3)*x3 +aik(1,4)*x4 +aik(1,5)*x5 +aik(1,6)*x6
2338  rhs_aij(2,5)= aik(2,1)*x1 +aik(2,2)*x2 +aik(2,3)*x3 +aik(2,4)*x4 +aik(2,5)*x5 +aik(2,6)*x6
2339  rhs_aij(3,5)= aik(3,1)*x1 +aik(3,2)*x2 +aik(3,3)*x3 +aik(3,4)*x4 +aik(3,5)*x5 +aik(3,6)*x6
2340  rhs_aij(4,5)= aik(4,1)*x1 +aik(4,2)*x2 +aik(4,3)*x3 +aik(4,4)*x4 +aik(4,5)*x5 +aik(4,6)*x6
2341  rhs_aij(5,5)= aik(5,1)*x1 +aik(5,2)*x2 +aik(5,3)*x3 +aik(5,4)*x4 +aik(5,5)*x5 +aik(5,6)*x6
2342  rhs_aij(6,5)= aik(6,1)*x1 +aik(6,2)*x2 +aik(6,3)*x3 +aik(6,4)*x4 +aik(6,5)*x5 +aik(6,6)*x6
2343 
2344  !C
2345  !C-- 6th Col.
2346  x1= akj(1,6)
2347  x2= akj(2,6)
2348  x3= akj(3,6)
2349  x4= akj(4,6)
2350  x5= akj(5,6)
2351  x6= akj(6,6)
2352 
2353  x2= x2 -dkinv(2,1)*x1
2354  x3= x3 -dkinv(3,1)*x1 -dkinv(3,2)*x2
2355  x4= x4 -dkinv(4,1)*x1 -dkinv(4,2)*x2 -dkinv(4,3)*x3
2356  x5= x5 -dkinv(5,1)*x1 -dkinv(5,2)*x2 -dkinv(5,3)*x3 -dkinv(5,4)*x4
2357  x6= x6 -dkinv(6,1)*x1 -dkinv(6,2)*x2 -dkinv(6,3)*x3 -dkinv(6,4)*x4 -dkinv(6,5)*x5
2358 
2359  x6= dkinv(6,6)* x6
2360  x5= dkinv(5,5)*( x5 -dkinv(5,6)*x6 )
2361  x4= dkinv(4,4)*( x4 -dkinv(4,6)*x6 -dkinv(4,5)*x5)
2362  x3= dkinv(3,3)*( x3 -dkinv(3,6)*x6 -dkinv(3,5)*x5 -dkinv(3,4)*x4)
2363  x2= dkinv(2,2)*( x2 -dkinv(2,6)*x6 -dkinv(2,5)*x5 -dkinv(2,4)*x4 -dkinv(2,3)*x3)
2364  x1= dkinv(1,1)*( x1 -dkinv(1,6)*x6 -dkinv(1,5)*x5 -dkinv(1,4)*x4 -dkinv(1,3)*x3 -dkinv(1,2)*x2)
2365 
2366  rhs_aij(1,6)= aik(1,1)*x1 +aik(1,2)*x2 +aik(1,3)*x3 +aik(1,4)*x4 +aik(1,5)*x5 +aik(1,6)*x6
2367  rhs_aij(2,6)= aik(2,1)*x1 +aik(2,2)*x2 +aik(2,3)*x3 +aik(2,4)*x4 +aik(2,5)*x5 +aik(2,6)*x6
2368  rhs_aij(3,6)= aik(3,1)*x1 +aik(3,2)*x2 +aik(3,3)*x3 +aik(3,4)*x4 +aik(3,5)*x5 +aik(3,6)*x6
2369  rhs_aij(4,6)= aik(4,1)*x1 +aik(4,2)*x2 +aik(4,3)*x3 +aik(4,4)*x4 +aik(4,5)*x5 +aik(4,6)*x6
2370  rhs_aij(5,6)= aik(5,1)*x1 +aik(5,2)*x2 +aik(5,3)*x3 +aik(5,4)*x4 +aik(5,5)*x5 +aik(5,6)*x6
2371  rhs_aij(6,6)= aik(6,1)*x1 +aik(6,2)*x2 +aik(6,3)*x3 +aik(6,4)*x4 +aik(6,5)*x5 +aik(6,6)*x6
2372 
2373  return
2374  end subroutine ilu1b66
2375 
2376 end module hecmw_precond_bilu_66
hecmw_matrix_misc
Definition: hecmw_matrix_misc.F90:6
hecmw_matrix_misc::hecmw_mat_get_sigma_diag
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
Definition: hecmw_matrix_misc.F90:811
hecmw_matrix_misc::hecmw_mat_get_precond
integer(kind=kint) function, public hecmw_mat_get_precond(hecMAT)
Definition: hecmw_matrix_misc.F90:381
hecmw_precond_bilu_66
Definition: hecmw_precond_BILU_66.f90:11
hecmw_precond_bilu_66::hecmw_precond_bilu_66_apply
subroutine, public hecmw_precond_bilu_66_apply(WW)
Definition: hecmw_precond_BILU_66.f90:128
hecmw_precond_bilu_66::ilu1a66
subroutine ilu1a66(ALU, D11, D12, D13, D14, D15, D16, D21, D22, D23, D24, D25, D26, D31, D32, D33, D34, D35, D36, D41, D42, D43, D44, D45, D46, D51, D52, D53, D54, D55, D56, D61, D62, D63, D64, D65, D66)
Definition: hecmw_precond_BILU_66.f90:2121
hecmw_precond_bilu_66::hecmw_precond_bilu_66_setup
subroutine, public hecmw_precond_bilu_66_setup(hecMAT)
Definition: hecmw_precond_BILU_66.f90:34
hecmw_precond_bilu_66::hecmw_precond_bilu_66_clear
subroutine, public hecmw_precond_bilu_66_clear()
Definition: hecmw_precond_BILU_66.f90:240
hecmw_precond_bilu_66::ilu1b66
subroutine ilu1b66(RHS_Aij, DkINV, Aik, Akj)
Definition: hecmw_precond_BILU_66.f90:2194
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:465