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