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