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