FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
12 use hecmw_util
14
15 private
16
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
31contains
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
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
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
2377end module hecmw_precond_bilu_66
real(kind=kreal) function, public hecmw_mat_get_sigma(hecmat)
integer(kind=kint) function, public hecmw_mat_get_precond(hecmat)
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecmat)
subroutine, public hecmw_precond_bilu_66_setup(hecmat)
subroutine, public hecmw_precond_bilu_66_clear()
subroutine, public hecmw_precond_bilu_66_apply(ww)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal