FrontISTR  6.0.0
Large-scale structural analysis program with finit element method
hecmw_precond_BILU_44.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5 
6 !C
7 !C***
8 !C*** module hecmw_precond_BILU_44
9 !C***
10 !C
11 module hecmw_precond_bilu_44
12  use hecmw_util
13  use hecmw_matrix_misc
14 
15  private
16 
17  public:: hecmw_precond_bilu_44_setup
18  public:: hecmw_precond_bilu_44_apply
19  public:: hecmw_precond_bilu_44_clear
20 
21  integer(kind=kint) :: N
22  real(kind=kreal), pointer :: dlu0(:) => null()
23  real(kind=kreal), pointer :: allu0(:) => null()
24  real(kind=kreal), pointer :: aulu0(:) => null()
25  integer(kind=kint), pointer :: inumFI1L(:) => null()
26  integer(kind=kint), pointer :: inumFI1U(:) => null()
27  integer(kind=kint), pointer :: FI1L(:) => null()
28  integer(kind=kint), pointer :: FI1U(:) => null()
29 
30  logical, save :: INITIALIZED = .false.
31 
32 contains
33 
34  subroutine hecmw_precond_bilu_44_setup(hecMAT)
35  implicit none
36  type(hecmwst_matrix), intent(inout) :: hecmat
37  integer(kind=kint ) :: np, npu, npl
38  integer(kind=kint ) :: precond
39  real (kind=kreal) :: sigma_diag
40 
41  real(kind=kreal), pointer :: d(:)
42  real(kind=kreal), pointer :: al(:)
43  real(kind=kreal), pointer :: au(:)
44 
45  integer(kind=kint ), pointer :: inl(:), inu(:)
46  integer(kind=kint ), pointer :: ial(:)
47  integer(kind=kint ), pointer :: iau(:)
48 
49  if (initialized) then
50  if (hecmat%Iarray(98) == 1) then ! need symbolic and numerical setup
51  call hecmw_precond_bilu_44_clear()
52  else if (hecmat%Iarray(97) == 1) then ! need numerical setup only
53  call hecmw_precond_bilu_44_clear() ! TEMPORARY
54  else
55  return
56  endif
57  endif
58 
59  n = hecmat%N
60  np = hecmat%NP
61  npl = hecmat%NPL
62  npu = hecmat%NPU
63  d => hecmat%D
64  al => hecmat%AL
65  au => hecmat%AU
66  inl => hecmat%indexL
67  inu => hecmat%indexU
68  ial => hecmat%itemL
69  iau => hecmat%itemU
70  precond = hecmw_mat_get_precond(hecmat)
71  sigma_diag = hecmw_mat_get_sigma_diag(hecmat)
72 
73  !if (PRECOND.eq.10) call FORM_ILU0_44 &
74  call form_ilu0_44 &
75  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
76  & sigma_diag)
77 
78  initialized = .true.
79  hecmat%Iarray(98) = 0 ! symbolic setup done
80  hecmat%Iarray(97) = 0 ! numerical setup done
81 
82  end subroutine hecmw_precond_bilu_44_setup
83 
84  subroutine hecmw_precond_bilu_44_apply(WW)
85  implicit none
86  real(kind=kreal), intent(inout) :: ww(:)
87  integer(kind=kint) :: i, j, isl, iel, isu, ieu, k
88  real(kind=kreal) :: sw1, sw2, sw3, sw4, x1, x2, x3, x4
89  !C
90  !C-- FORWARD
91 
92  do i= 1, n
93  sw1= ww(4*i-3)
94  sw2= ww(4*i-2)
95  sw3= ww(4*i-1)
96  sw4= ww(4*i )
97  isl= inumfi1l(i-1)+1
98  iel= inumfi1l(i)
99  do j= isl, iel
100  k= fi1l(j)
101  x1= ww(4*k-3)
102  x2= ww(4*k-2)
103  x3= ww(4*k-1)
104  x4= ww(4*k )
105  sw1= sw1 - allu0(16*j-15)*x1-allu0(16*j-14)*x2-allu0(16*j-13)*x3-allu0(16*j-12)*x4
106  sw2= sw2 - allu0(16*j-11)*x1-allu0(16*j-10)*x2-allu0(16*j- 9)*x3-allu0(16*j- 8)*x4
107  sw3= sw3 - allu0(16*j- 7)*x1-allu0(16*j- 6)*x2-allu0(16*j- 5)*x3-allu0(16*j- 4)*x4
108  sw4= sw4 - allu0(16*j- 3)*x1-allu0(16*j- 2)*x2-allu0(16*j- 1)*x3-allu0(16*j )*x4
109  enddo
110 
111  x1= sw1
112  x2= sw2
113  x3= sw3
114  x4= sw4
115  x2= x2 - dlu0(16*i-11)*x1
116  x3= x3 - dlu0(16*i- 7)*x1 - dlu0(16*i-6)*x2
117  x4= x4 - dlu0(16*i- 3)*x1 - dlu0(16*i-2)*x2 - dlu0(16*i-1)*x3
118  x4= dlu0(16*i )* x4
119  x3= dlu0(16*i- 5)*(x3 - dlu0(16*i- 4)*x4)
120  x2= dlu0(16*i-10)*(x2 - dlu0(16*i- 8)*x4 - dlu0(16*i- 9)*x3 )
121  x1= dlu0(16*i-15)*(x1 - dlu0(16*i-12)*x4 - dlu0(16*i-13)*x3 - dlu0(16*i-14)*x2)
122 
123  ww(4*i-3)= x1
124  ww(4*i-2)= x2
125  ww(4*i-1)= x3
126  ww(4*i )= x4
127  enddo
128 
129  !C
130  !C-- BACKWARD
131 
132  do i= n, 1, -1
133  isu= inumfi1u(i-1) + 1
134  ieu= inumfi1u(i)
135  sw1= 0.d0
136  sw2= 0.d0
137  sw3= 0.d0
138  sw4= 0.d0
139  do j= ieu, isu, -1
140  k= fi1u(j)
141  x1= ww(4*k-3)
142  x2= ww(4*k-2)
143  x3= ww(4*k-1)
144  x4= ww(4*k )
145  sw1= sw1 + aulu0(16*j-15)*x1+aulu0(16*j-14)*x2+aulu0(16*j-13)*x3+aulu0(16*j-12)*x4
146  sw2= sw2 + aulu0(16*j-11)*x1+aulu0(16*j-10)*x2+aulu0(16*j- 9)*x3+aulu0(16*j- 8)*x4
147  sw3= sw3 + aulu0(16*j- 7)*x1+aulu0(16*j- 6)*x2+aulu0(16*j- 5)*x3+aulu0(16*j- 4)*x4
148  sw4= sw4 + aulu0(16*j- 3)*x1+aulu0(16*j- 2)*x2+aulu0(16*j- 1)*x3+aulu0(16*j )*x4
149  enddo
150  x1= sw1
151  x2= sw2
152  x3= sw3
153  x4= sw4
154  x2= x2 - dlu0(16*i-11)*x1
155  x3= x3 - dlu0(16*i- 7)*x1 - dlu0(16*i-6)*x2
156  x4= x4 - dlu0(16*i- 3)*x1 - dlu0(16*i-2)*x2 - dlu0(16*i-1)*x3
157  x4= dlu0(16*i )* x4
158  x3= dlu0(16*i- 5)*( x3 - dlu0(16*i- 4)*x4 )
159  x2= dlu0(16*i-10)*( x2 - dlu0(16*i- 8)*x4 - dlu0(16*i- 9)*x3 )
160  x1= dlu0(16*i-15)*( x1 - dlu0(16*i-12)*x4 - dlu0(16*i-13)*x3 - dlu0(16*i-14)*x2)
161  ww(4*i-3)= ww(4*i-3) - x1
162  ww(4*i-2)= ww(4*i-2) - x2
163  ww(4*i-1)= ww(4*i-1) - x3
164  ww(4*i )= ww(4*i ) - x4
165  enddo
166  end subroutine hecmw_precond_bilu_44_apply
167 
168  subroutine hecmw_precond_bilu_44_clear()
169  implicit none
170  if (associated(dlu0)) deallocate(dlu0)
171  if (associated(allu0)) deallocate(allu0)
172  if (associated(aulu0)) deallocate(aulu0)
173  if (associated(inumfi1l)) deallocate(inumfi1l)
174  if (associated(inumfi1u)) deallocate(inumfi1u)
175  if (associated(fi1l)) deallocate(fi1l)
176  if (associated(fi1u)) deallocate(fi1u)
177  nullify(dlu0)
178  nullify(allu0)
179  nullify(aulu0)
180  nullify(inumfi1l)
181  nullify(inumfi1u)
182  nullify(fi1l)
183  nullify(fi1u)
184  initialized = .false.
185  end subroutine hecmw_precond_bilu_44_clear
186 
187  !C
188  !C***
189  !C*** FORM_ILU0_44
190  !C***
191  !C
192  !C form ILU(0) matrix
193  !C
194  subroutine form_ilu0_44 &
195  & (n, np, npl, npu, d, al, inl, ial, au, inu, iau, &
196  & sigma_diag)
197  implicit none
198  integer(kind=kint ), intent(in):: n, np, npu, npl
199  real (kind=kreal), intent(in):: sigma_diag
200 
201  real(kind=kreal), dimension(16*NPL), intent(in):: al
202  real(kind=kreal), dimension(16*NPU), intent(in):: au
203  real(kind=kreal), dimension(16*NP ), intent(in):: d
204 
205  integer(kind=kint ), dimension(0:NP) ,intent(in) :: inu, inl
206  integer(kind=kint ), dimension( NPL),intent(in) :: ial
207  integer(kind=kint ), dimension( NPU),intent(in) :: iau
208 
209  integer(kind=kint), dimension(:), allocatable :: iw1, iw2
210  real (kind=kreal), dimension(4,4) :: rhs_aij, dkinv, aik, akj
211  integer(kind=kint) :: i,jj,ij0,kk
212  integer(kind=kint) :: j,k
213  allocate (iw1(np) , iw2(np))
214  allocate(dlu0(9*np), allu0(9*npl), aulu0(9*npu))
215  allocate(inumfi1l(0:np), inumfi1u(0:np), fi1l(npl), fi1u(npu))
216 
217  do i=1,16*np
218  dlu0(i) = d(i)
219  end do
220  do i=1,16*npl
221  allu0(i) = al(i)
222  end do
223  do i=1,16*npu
224  aulu0(i) = au(i)
225  end do
226  do i=0,np
227  inumfi1l(i) = inl(i)
228  inumfi1u(i) = inu(i)
229  end do
230  do i=1,npl
231  fi1l(i) = ial(i)
232  end do
233  do i=1,npu
234  fi1u(i) = iau(i)
235  end do
236 
237  !C
238  !C +----------------------+
239  !C | ILU(0) factorization |
240  !C +----------------------+
241  !C===
242  do i=1,np
243  dlu0(16*i-15)=dlu0(16*i-15)*sigma_diag
244  dlu0(16*i-10)=dlu0(16*i-10)*sigma_diag
245  dlu0(16*i- 5)=dlu0(16*i- 5)*sigma_diag
246  dlu0(16*i )=dlu0(16*i )*sigma_diag
247  enddo
248 
249  i = 1
250  call ilu1a44 (dkinv, &
251  dlu0(16*i-15), dlu0(16*i-14), dlu0(16*i-13), dlu0(16*i-12), &
252  dlu0(16*i-11), dlu0(16*i-10), dlu0(16*i- 9), dlu0(16*i- 8), &
253  dlu0(16*i- 7), dlu0(16*i- 6), dlu0(16*i- 5), dlu0(16*i- 4), &
254  dlu0(16*i- 3), dlu0(16*i- 2), dlu0(16*i- 1), dlu0(16*i ) )
255  dlu0(16*i-15)= dkinv(1,1)
256  dlu0(16*i-14)= dkinv(1,2)
257  dlu0(16*i-13)= dkinv(1,3)
258  dlu0(16*i-12)= dkinv(1,4)
259  dlu0(16*i-11)= dkinv(2,1)
260  dlu0(16*i-10)= dkinv(2,2)
261  dlu0(16*i- 9)= dkinv(2,3)
262  dlu0(16*i- 8)= dkinv(2,4)
263  dlu0(16*i- 7)= dkinv(3,1)
264  dlu0(16*i- 6)= dkinv(3,2)
265  dlu0(16*i- 5)= dkinv(3,3)
266  dlu0(16*i- 4)= dkinv(3,4)
267  dlu0(16*i- 3)= dkinv(4,1)
268  dlu0(16*i- 2)= dkinv(4,2)
269  dlu0(16*i- 1)= dkinv(4,3)
270  dlu0(16*i )= dkinv(4,4)
271 
272  do i= 2, np
273  iw1= 0
274  iw2= 0
275 
276  do k= inumfi1l(i-1)+1, inumfi1l(i)
277  iw1(fi1l(k))= k
278  enddo
279 
280  do k= inumfi1u(i-1)+1, inumfi1u(i)
281  iw2(fi1u(k))= k
282  enddo
283 
284  do kk= inl(i-1)+1, inl(i)
285  k= ial(kk)
286 
287  dkinv(1,1)= dlu0(16*k-15)
288  dkinv(1,2)= dlu0(16*k-14)
289  dkinv(1,3)= dlu0(16*k-13)
290  dkinv(1,4)= dlu0(16*k-12)
291  dkinv(2,1)= dlu0(16*k-11)
292  dkinv(2,2)= dlu0(16*k-10)
293  dkinv(2,3)= dlu0(16*k- 9)
294  dkinv(2,4)= dlu0(16*k- 8)
295  dkinv(3,1)= dlu0(16*k- 7)
296  dkinv(3,2)= dlu0(16*k- 6)
297  dkinv(3,3)= dlu0(16*k- 5)
298  dkinv(3,4)= dlu0(16*k- 4)
299  dkinv(4,1)= dlu0(16*k- 3)
300  dkinv(4,2)= dlu0(16*k- 2)
301  dkinv(4,3)= dlu0(16*k- 1)
302  dkinv(4,4)= dlu0(16*k )
303 
304  aik(1,1)= allu0(16*kk-15)
305  aik(1,2)= allu0(16*kk-14)
306  aik(1,3)= allu0(16*kk-13)
307  aik(1,4)= allu0(16*kk-12)
308  aik(2,1)= allu0(16*kk-11)
309  aik(2,2)= allu0(16*kk-10)
310  aik(2,3)= allu0(16*kk- 9)
311  aik(2,4)= allu0(16*kk- 8)
312  aik(3,1)= allu0(16*kk- 7)
313  aik(3,2)= allu0(16*kk- 6)
314  aik(3,3)= allu0(16*kk- 5)
315  aik(3,4)= allu0(16*kk- 4)
316  aik(4,1)= allu0(16*kk- 3)
317  aik(4,2)= allu0(16*kk- 2)
318  aik(4,3)= allu0(16*kk- 1)
319  aik(4,4)= allu0(16*kk )
320 
321  do jj= inu(k-1)+1, inu(k)
322  j= iau(jj)
323  if (iw1(j).eq.0.and.iw2(j).eq.0) cycle
324 
325  akj(1,1)= aulu0(16*jj-15)
326  akj(1,2)= aulu0(16*jj-14)
327  akj(1,3)= aulu0(16*jj-13)
328  akj(1,4)= aulu0(16*jj-12)
329  akj(2,1)= aulu0(16*jj-11)
330  akj(2,2)= aulu0(16*jj-10)
331  akj(2,3)= aulu0(16*jj- 9)
332  akj(2,4)= aulu0(16*jj- 8)
333  akj(3,1)= aulu0(16*jj- 7)
334  akj(3,2)= aulu0(16*jj- 6)
335  akj(3,3)= aulu0(16*jj- 5)
336  akj(3,4)= aulu0(16*jj- 4)
337  akj(4,1)= aulu0(16*jj- 3)
338  akj(4,2)= aulu0(16*jj- 2)
339  akj(4,3)= aulu0(16*jj- 1)
340  akj(4,4)= aulu0(16*jj )
341 
342  call ilu1b44 (rhs_aij, dkinv, aik, akj)
343 
344  if (j.eq.i) then
345  dlu0(16*i-15)= dlu0(16*i-15) - rhs_aij(1,1)
346  dlu0(16*i-14)= dlu0(16*i-14) - rhs_aij(1,2)
347  dlu0(16*i-13)= dlu0(16*i-13) - rhs_aij(1,3)
348  dlu0(16*i-12)= dlu0(16*i-12) - rhs_aij(1,4)
349  dlu0(16*i-11)= dlu0(16*i-11) - rhs_aij(2,1)
350  dlu0(16*i-10)= dlu0(16*i-10) - rhs_aij(2,2)
351  dlu0(16*i- 9)= dlu0(16*i- 9) - rhs_aij(2,3)
352  dlu0(16*i- 8)= dlu0(16*i- 8) - rhs_aij(2,4)
353  dlu0(16*i- 7)= dlu0(16*i- 7) - rhs_aij(3,1)
354  dlu0(16*i- 6)= dlu0(16*i- 6) - rhs_aij(3,2)
355  dlu0(16*i- 5)= dlu0(16*i- 5) - rhs_aij(3,3)
356  dlu0(16*i- 4)= dlu0(16*i- 4) - rhs_aij(3,4)
357  dlu0(16*i- 3)= dlu0(16*i- 3) - rhs_aij(4,1)
358  dlu0(16*i- 2)= dlu0(16*i- 2) - rhs_aij(4,2)
359  dlu0(16*i- 1)= dlu0(16*i- 1) - rhs_aij(4,3)
360  dlu0(16*i )= dlu0(16*i ) - rhs_aij(4,4)
361  endif
362 
363  if (j.lt.i) then
364  ij0= iw1(j)
365  allu0(16*ij0-15)= allu0(16*ij0-15) - rhs_aij(1,1)
366  allu0(16*ij0-14)= allu0(16*ij0-14) - rhs_aij(1,2)
367  allu0(16*ij0-13)= allu0(16*ij0-13) - rhs_aij(1,3)
368  allu0(16*ij0-12)= allu0(16*ij0-12) - rhs_aij(1,4)
369  allu0(16*ij0-11)= allu0(16*ij0-11) - rhs_aij(2,1)
370  allu0(16*ij0-10)= allu0(16*ij0-10) - rhs_aij(2,2)
371  allu0(16*ij0- 9)= allu0(16*ij0- 9) - rhs_aij(2,3)
372  allu0(16*ij0- 8)= allu0(16*ij0- 8) - rhs_aij(2,4)
373  allu0(16*ij0- 7)= allu0(16*ij0- 7) - rhs_aij(3,1)
374  allu0(16*ij0- 6)= allu0(16*ij0- 6) - rhs_aij(3,2)
375  allu0(16*ij0- 5)= allu0(16*ij0- 5) - rhs_aij(3,3)
376  allu0(16*ij0- 4)= allu0(16*ij0- 4) - rhs_aij(3,4)
377  allu0(16*ij0- 3)= allu0(16*ij0- 3) - rhs_aij(4,1)
378  allu0(16*ij0- 2)= allu0(16*ij0- 2) - rhs_aij(4,2)
379  allu0(16*ij0- 1)= allu0(16*ij0- 1) - rhs_aij(4,3)
380  allu0(16*ij0 )= allu0(16*ij0 ) - rhs_aij(4,4)
381  endif
382 
383  if (j.gt.i) then
384  ij0= iw2(j)
385  aulu0(16*ij0-15)= aulu0(16*ij0-15) - rhs_aij(1,1)
386  aulu0(16*ij0-14)= aulu0(16*ij0-14) - rhs_aij(1,2)
387  aulu0(16*ij0-13)= aulu0(16*ij0-13) - rhs_aij(1,3)
388  aulu0(16*ij0-12)= aulu0(16*ij0-12) - rhs_aij(1,4)
389  aulu0(16*ij0-11)= aulu0(16*ij0-11) - rhs_aij(2,1)
390  aulu0(16*ij0-10)= aulu0(16*ij0-10) - rhs_aij(2,2)
391  aulu0(16*ij0- 9)= aulu0(16*ij0- 9) - rhs_aij(2,3)
392  aulu0(16*ij0- 8)= aulu0(16*ij0- 8) - rhs_aij(2,4)
393  aulu0(16*ij0- 7)= aulu0(16*ij0- 7) - rhs_aij(3,1)
394  aulu0(16*ij0- 6)= aulu0(16*ij0- 6) - rhs_aij(3,2)
395  aulu0(16*ij0- 5)= aulu0(16*ij0- 5) - rhs_aij(3,3)
396  aulu0(16*ij0- 4)= aulu0(16*ij0- 4) - rhs_aij(3,4)
397  aulu0(16*ij0- 3)= aulu0(16*ij0- 3) - rhs_aij(4,1)
398  aulu0(16*ij0- 2)= aulu0(16*ij0- 2) - rhs_aij(4,2)
399  aulu0(16*ij0- 1)= aulu0(16*ij0- 1) - rhs_aij(4,3)
400  aulu0(16*ij0 )= aulu0(16*ij0 ) - rhs_aij(4,4)
401  endif
402 
403  enddo
404  enddo
405 
406  call ilu1a44 (dkinv, &
407  dlu0(16*i-15), dlu0(16*i-14), dlu0(16*i-13), dlu0(16*i-12), &
408  dlu0(16*i-11), dlu0(16*i-10), dlu0(16*i- 9), dlu0(16*i- 8), &
409  dlu0(16*i- 7), dlu0(16*i- 6), dlu0(16*i- 5), dlu0(16*i- 4), &
410  dlu0(16*i- 3), dlu0(16*i- 2), dlu0(16*i- 1), dlu0(16*i ) )
411  dlu0(16*i-15)= dkinv(1,1)
412  dlu0(16*i-14)= dkinv(1,2)
413  dlu0(16*i-13)= dkinv(1,3)
414  dlu0(16*i-12)= dkinv(1,4)
415  dlu0(16*i-11)= dkinv(2,1)
416  dlu0(16*i-10)= dkinv(2,2)
417  dlu0(16*i- 9)= dkinv(2,3)
418  dlu0(16*i- 8)= dkinv(2,4)
419  dlu0(16*i- 7)= dkinv(3,1)
420  dlu0(16*i- 6)= dkinv(3,2)
421  dlu0(16*i- 5)= dkinv(3,3)
422  dlu0(16*i- 4)= dkinv(3,4)
423  dlu0(16*i- 3)= dkinv(4,1)
424  dlu0(16*i- 2)= dkinv(4,2)
425  dlu0(16*i- 1)= dkinv(4,3)
426  dlu0(16*i )= dkinv(4,4)
427  enddo
428 
429  deallocate (iw1, iw2)
430  end subroutine form_ilu0_44
431 
432  !C
433  !C***
434  !C*** ILU1a44
435  !C***
436  !C
437  !C computes LU factorization of 4*4 Diagonal Block
438  !C
439  subroutine ilu1a44 (ALU, D11,D12,D13,D14,D21,D22,D23,D24,D31,D32,D33,D34,D41,D42,D43,D44)
440  use hecmw_util
441  implicit none
442  real(kind=kreal) :: alu(4,4), pw(4)
443  real(kind=kreal) :: d11,d12,d13,d14,d21,d22,d23,d24,d31,d32,d33,d34,d41,d42,d43,d44
444  integer(kind=kint) :: i,j,k
445 
446  alu(1,1)= d11
447  alu(1,2)= d12
448  alu(1,3)= d13
449  alu(1,4)= d14
450  alu(2,1)= d21
451  alu(2,2)= d22
452  alu(2,3)= d23
453  alu(2,4)= d24
454  alu(3,1)= d31
455  alu(3,2)= d32
456  alu(3,3)= d33
457  alu(3,4)= d34
458  alu(4,1)= d41
459  alu(4,2)= d42
460  alu(4,3)= d43
461  alu(4,4)= d44
462 
463  do k= 1, 4
464  alu(k,k)= 1.d0/alu(k,k)
465  do i= k+1, 4
466  alu(i,k)= alu(i,k) * alu(k,k)
467  do j= k+1, 4
468  pw(j)= alu(i,j) - alu(i,k)*alu(k,j)
469  enddo
470  do j= k+1, 4
471  alu(i,j)= pw(j)
472  enddo
473  enddo
474  enddo
475 
476  return
477  end subroutine ilu1a44
478 
479  !C
480  !C***
481  !C*** ILU1b44
482  !C***
483  !C
484  !C computes L_ik * D_k_INV * U_kj at ILU factorization
485  !C for 4*4 Block Type Matrix
486  !C
487  subroutine ilu1b44 (RHS_Aij, DkINV, Aik, Akj)
488  use hecmw_util
489  implicit none
490  real(kind=kreal) :: rhs_aij(4,4), dkinv(4,4), aik(4,4), akj(4,4)
491  real(kind=kreal) :: x1,x2,x3,x4
492 
493  !C
494  !C-- 1st Col.
495  x1= akj(1,1)
496  x2= akj(2,1)
497  x3= akj(3,1)
498  x4= akj(4,1)
499 
500  x2= x2 - dkinv(2,1)*x1
501  x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
502  x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
503 
504  x4= dkinv(4,4)* x4
505  x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
506  x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
507  x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
508 
509  rhs_aij(1,1)= aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
510  rhs_aij(2,1)= aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
511  rhs_aij(3,1)= aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
512  rhs_aij(4,1)= aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
513 
514  !C
515  !C-- 2nd Col.
516  x1= akj(1,2)
517  x2= akj(2,2)
518  x3= akj(3,2)
519  x4= akj(4,2)
520 
521  x2= x2 - dkinv(2,1)*x1
522  x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
523  x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
524 
525  x4= dkinv(4,4)* x4
526  x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
527  x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
528  x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
529 
530  rhs_aij(1,2)= aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
531  rhs_aij(2,2)= aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
532  rhs_aij(3,2)= aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
533  rhs_aij(4,2)= aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
534 
535  !C
536  !C-- 3rd Col.
537  x1= akj(1,3)
538  x2= akj(2,3)
539  x3= akj(3,3)
540  x4= akj(4,3)
541 
542  x2= x2 - dkinv(2,1)*x1
543  x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
544  x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
545 
546  x4= dkinv(4,4)* x4
547  x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
548  x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
549  x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
550 
551  rhs_aij(1,3)= aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
552  rhs_aij(2,3)= aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
553  rhs_aij(3,3)= aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
554  rhs_aij(4,3)= aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
555 
556  !C
557  !C-- 4th Col.
558  x1= akj(1,4)
559  x2= akj(2,4)
560  x3= akj(3,4)
561  x4= akj(4,4)
562 
563  x2= x2 - dkinv(2,1)*x1
564  x3= x3 - dkinv(3,1)*x1 - dkinv(3,2)*x2
565  x4= x4 - dkinv(4,1)*x1 - dkinv(4,2)*x2 - dkinv(4,3)*x3
566 
567  x4= dkinv(4,4)* x4
568  x3= dkinv(3,3)*( x3 - dkinv(3,4)*x4 )
569  x2= dkinv(2,2)*( x2 - dkinv(2,4)*x4 - dkinv(2,3)*x3 )
570  x1= dkinv(1,1)*( x1 - dkinv(1,4)*x4 - dkinv(1,3)*x3 - dkinv(1,2)*x2)
571 
572  rhs_aij(1,4)= aik(1,1)*x1 + aik(1,2)*x2 + aik(1,3)*x3 + aik(1,4)*x4
573  rhs_aij(2,4)= aik(2,1)*x1 + aik(2,2)*x2 + aik(2,3)*x3 + aik(2,4)*x4
574  rhs_aij(3,4)= aik(3,1)*x1 + aik(3,2)*x2 + aik(3,3)*x3 + aik(3,4)*x4
575  rhs_aij(4,4)= aik(4,1)*x1 + aik(4,2)*x2 + aik(4,3)*x3 + aik(4,4)*x4
576 
577  return
578  end subroutine ilu1b44
579 
580 end module hecmw_precond_bilu_44
581 
hecmw_matrix_misc
Definition: hecmw_matrix_misc.F90:6
hecmw_matrix_misc::hecmw_mat_get_sigma_diag
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
Definition: hecmw_matrix_misc.F90:811
hecmw_matrix_misc::hecmw_mat_get_precond
integer(kind=kint) function, public hecmw_mat_get_precond(hecMAT)
Definition: hecmw_matrix_misc.F90:381
hecmw_precond_bilu_44
Definition: hecmw_precond_BILU_44.f90:11
hecmw_precond_bilu_44::hecmw_precond_bilu_44_apply
subroutine, public hecmw_precond_bilu_44_apply(WW)
Definition: hecmw_precond_BILU_44.f90:85
hecmw_precond_bilu_44::hecmw_precond_bilu_44_clear
subroutine, public hecmw_precond_bilu_44_clear()
Definition: hecmw_precond_BILU_44.f90:169
hecmw_precond_bilu_44::hecmw_precond_bilu_44_setup
subroutine, public hecmw_precond_bilu_44_setup(hecMAT)
Definition: hecmw_precond_BILU_44.f90:35
hecmw_precond_bilu_44::ilu1b44
subroutine ilu1b44(RHS_Aij, DkINV, Aik, Akj)
Definition: hecmw_precond_BILU_44.f90:488
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:465