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