FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_precond_SAINV_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 module hecmw_precond_sainv_33
7  use hecmw_util
8  use m_hecmw_comm_f
9  use hecmw_matrix_misc
10  !$ use omp_lib
11 
12  private
13 
14  public:: hecmw_precond_sainv_33_setup
15  public:: hecmw_precond_sainv_33_apply
16  public:: hecmw_precond_sainv_33_clear
17 
18  integer(4),parameter :: krealp = 8
19 
20  integer(kind=kint) :: NPFIU, NPFIL
21  integer(kind=kint) :: N
22  integer(kind=kint), pointer :: inumFI1L(:) => null()
23  integer(kind=kint), pointer :: inumFI1U(:) => null()
24  integer(kind=kint), pointer :: FI1L(:) => null()
25  integer(kind=kint), pointer :: FI1U(:) => null()
26 
27  integer(kind=kint), pointer :: indexL(:) => null()
28  integer(kind=kint), pointer :: indexU(:) => null()
29  integer(kind=kint), pointer :: itemL(:) => null()
30  integer(kind=kint), pointer :: itemU(:) => null()
31  real(kind=kreal), pointer :: d(:) => null()
32  real(kind=kreal), pointer :: al(:) => null()
33  real(kind=kreal), pointer :: au(:) => null()
34 
35  real(kind=krealp), pointer :: sainvu(:) => null()
36  real(kind=krealp), pointer :: sainvl(:) => null()
37  real(kind=krealp), pointer :: sainvd(:) => null()
38  real(kind=kreal), pointer :: t(:) => null()
39 
40 contains
41 
42  !C***
43  !C*** hecmw_precond_33_sainv_setup
44  !C***
45  subroutine hecmw_precond_sainv_33_setup(hecMAT)
46  implicit none
47  type(hecmwst_matrix) :: hecmat
48 
49  integer(kind=kint ) :: precond
50 
51  real(kind=krealp) :: filter
52 
53  n = hecmat%N
54  precond = hecmw_mat_get_precond(hecmat)
55 
56  d => hecmat%D
57  au=> hecmat%AU
58  al=> hecmat%AL
59  indexl => hecmat%indexL
60  indexu => hecmat%indexU
61  iteml => hecmat%itemL
62  itemu => hecmat%itemU
63 
64  if (precond.eq.20) call form_ilu1_sainv_33(hecmat)
65 
66  allocate (sainvd(9*hecmat%NP))
67  allocate (sainvl(9*npfiu))
68  allocate (t(3*hecmat%NP))
69  sainvd = 0.0d0
70  sainvl = 0.0d0
71  t = 0.0d0
72 
73  filter= hecmat%Rarray(5)
74 
75  write(*,"(a,F15.8)")"### SAINV FILTER :",filter
76 
77  call hecmw_sainv_33(hecmat)
78 
79  allocate (sainvu(9*npfiu))
80  sainvu = 0.0d0
81 
82  call hecmw_sainv_make_u_33(hecmat)
83 
84  end subroutine hecmw_precond_sainv_33_setup
85 
86  subroutine hecmw_sainv_lu_33()
87  implicit none
88  integer(kind=kint) :: i,j,js,je,in
89  real(kind=kreal) :: x1, x2, x3
90 
91  do i=1, n
92  sainvd(9*i-5) = sainvd(9*i-5)*sainvd(9*i-4)
93  sainvd(9*i-2) = sainvd(9*i-2)*sainvd(9*i )
94  sainvd(9*i-1) = sainvd(9*i-1)*sainvd(9*i )
95  enddo
96 
97  do i=1, n
98  js = inumfi1l(i-1)+1
99  je = inumfi1l(i)
100  do j= js,je
101  in= fi1l(j)
102  x1= sainvd(9*i-8)
103  x2= sainvd(9*i-4)
104  x3= sainvd(9*i )
105  sainvl(9*j-8) = sainvl(9*j-8)*x1
106  sainvl(9*j-7) = sainvl(9*j-7)*x1
107  sainvl(9*j-6) = sainvl(9*j-6)*x1
108  sainvl(9*j-5) = sainvl(9*j-5)*x2
109  sainvl(9*j-4) = sainvl(9*j-4)*x2
110  sainvl(9*j-3) = sainvl(9*j-3)*x2
111  sainvl(9*j-2) = sainvl(9*j-2)*x3
112  sainvl(9*j-1) = sainvl(9*j-1)*x3
113  sainvl(9*j ) = sainvl(9*j )*x3
114  enddo
115  enddo
116 
117  end subroutine hecmw_sainv_lu_33
118 
119  subroutine hecmw_precond_sainv_33_apply(R, ZP)
120  implicit none
121  real(kind=kreal), intent(inout) :: zp(:)
122  real(kind=kreal), intent(in) :: r(:)
123  integer(kind=kint) :: in, i, j, isl, iel, isu, ieu
124  real(kind=kreal) :: sw1, sw2, sw3, x1, x2, x3
125 
126  !$OMP PARALLEL DEFAULT(NONE) &
127  !$OMP&PRIVATE(i,X1,X2,X3,SW1,SW2,SW3,j,in,isL,ieL,isU,ieU) &
128  !$OMP&SHARED(N,SAINVD,SAINVL,SAINVU,inumFI1U,FI1U,inumFI1L,FI1L,R,T,ZP)
129  !$OMP DO
130  !C-- FORWARD
131  do i= 1, n
132  sw1= 0.0d0
133  sw2= 0.0d0
134  sw3= 0.0d0
135 
136  isl= inumfi1l(i-1)+1
137  iel= inumfi1l(i)
138  do j= isl, iel
139  in= fi1l(j)
140  x1= r(3*in-2)
141  x2= r(3*in-1)
142  x3= r(3*in )
143  sw1= sw1 + sainvl(9*j-8)*x1 + sainvl(9*j-7)*x2 + sainvl(9*j-6)*x3
144  sw2= sw2 + sainvl(9*j-5)*x1 + sainvl(9*j-4)*x2 + sainvl(9*j-3)*x3
145  sw3= sw3 + sainvl(9*j-2)*x1 + sainvl(9*j-1)*x2 + sainvl(9*j )*x3
146  enddo
147 
148  x1= r(3*i-2)
149  x2= r(3*i-1)
150  x3= r(3*i )
151 
152  t(3*i-2)= (x1 + sw1)*sainvd(9*i-8)
153  t(3*i-1)= (x2 + sainvd(9*i-7)*x1 + sw2)*sainvd(9*i-4)
154  t(3*i )= (x3 + sainvd(9*i-6)*x1 + sainvd(9*i-3)*x2 + sw3)*sainvd(9*i )
155  enddo
156  !$OMP END DO
157  !$OMP DO
158  !C-- BACKWARD
159  do i= 1, n
160  sw1= 0.0d0
161  sw2= 0.0d0
162  sw3= 0.0d0
163 
164  isu= inumfi1u(i-1) + 1
165  ieu= inumfi1u(i)
166  do j= isu, ieu
167  in= fi1u(j)
168  x1= t(3*in-2)
169  x2= t(3*in-1)
170  x3= t(3*in )
171  sw1= sw1 + sainvu(9*j-8)*x1 + sainvu(9*j-7)*x2 + sainvu(9*j-6)*x3
172  sw2= sw2 + sainvu(9*j-5)*x1 + sainvu(9*j-4)*x2 + sainvu(9*j-3)*x3
173  sw3= sw3 + sainvu(9*j-2)*x1 + sainvu(9*j-1)*x2 + sainvu(9*j )*x3
174  enddo
175 
176  x1= t(3*i-2)
177  x2= t(3*i-1)
178  x3= t(3*i )
179 
180  zp(3*i-2)= x1 + sw1 + sainvd(9*i-7)*x2 + sainvd(9*i-6)*x3
181  zp(3*i-1)= x2 + sw2 + sainvd(9*i-3)*x3
182  zp(3*i )= x3 + sw3
183  enddo
184  !$OMP END DO
185  !$OMP END PARALLEL
186 
187  end subroutine hecmw_precond_sainv_33_apply
188 
189 
190  !C***
191  !C*** hecmw_rif_33
192  !C***
193  subroutine hecmw_sainv_33(hecMAT)
194  implicit none
195  type (hecmwst_matrix) :: hecmat
196 
197  integer(kind=kint) :: i, j, js, je, in, itr
198  real(kind=krealp) :: x1, x2, x3, dd, dd1, dd2, dd3, dtemp(3)
199  real(kind=krealp) :: filter
200  real(kind=krealp), allocatable :: zz(:), vv(:)
201 
202  filter= hecmat%Rarray(5)
203 
204  allocate (vv(3*hecmat%NP))
205  allocate (zz(3*hecmat%NP))
206  do itr=1,n
207 
208  !------------------------------ iitr = 1 ----------------------------------------
209 
210  zz(:) = 0.0d0
211  vv(:) = 0.0d0
212 
213  !{v}=[A]{zi}
214 
215  zz(3*itr-2)= sainvd(9*itr-8)
216  zz(3*itr-1)= sainvd(9*itr-5)
217  zz(3*itr )= sainvd(9*itr-2)
218 
219  zz(3*itr-2)= 1.0d0! * SIGMA_DIAG
220 
221  js= inumfi1l(itr-1) + 1
222  je= inumfi1l(itr )
223  do j= js, je
224  in = fi1l(j)
225  zz(3*in-2)= sainvl(9*j-8)
226  zz(3*in-1)= sainvl(9*j-7)
227  zz(3*in )= sainvl(9*j-6)
228  enddo
229 
230  do i= 1, itr
231  x1= zz(3*i-2)
232  x2= zz(3*i-1)
233  x3= zz(3*i )
234  vv(3*i-2) = vv(3*i-2) + d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
235  vv(3*i-1) = vv(3*i-1) + d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
236  vv(3*i ) = vv(3*i ) + d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i )*x3
237 
238  js= indexl(i-1) + 1
239  je= indexl(i )
240  do j=js,je
241  in = iteml(j)
242  vv(3*in-2)= vv(3*in-2) + al(9*j-8)*x1 + al(9*j-5)*x2 + al(9*j-2)*x3
243  vv(3*in-1)= vv(3*in-1) + al(9*j-7)*x1 + al(9*j-4)*x2 + al(9*j-1)*x3
244  vv(3*in )= vv(3*in ) + al(9*j-6)*x1 + al(9*j-3)*x2 + al(9*j )*x3
245  enddo
246 
247  js= indexu(i-1) + 1
248  je= indexu(i )
249  do j= js, je
250  in = itemu(j)
251  vv(3*in-2)= vv(3*in-2) + au(9*j-8)*x1 + au(9*j-5)*x2 + au(9*j-2)*x3
252  vv(3*in-1)= vv(3*in-1) + au(9*j-7)*x1 + au(9*j-4)*x2 + au(9*j-1)*x3
253  vv(3*in )= vv(3*in ) + au(9*j-6)*x1 + au(9*j-3)*x2 + au(9*j )*x3
254  enddo
255  enddo
256 
257  !{d}={v^t}{z_j}
258 
259  !dtemp(1) = SAINVD(9*itr-8)
260  !dtemp(2) = SAINVD(9*itr-4)
261 
262  !$OMP PARALLEL DEFAULT(NONE) &
263  !$OMP&PRIVATE(i,j,jS,jE,in,X1,X2,X3) &
264  !$OMP&FIRSTPRIVATE(vv) &
265  !$OMP&SHARED(N,itr,SAINVD,SAINVL,inumFI1L,FI1L)
266  !$OMP DO
267  do i=itr,n
268  sainvd(9*i-8) = vv(3*i-2)
269  sainvd(9*i-4) = vv(3*i-2)*sainvd(9*i-7) + vv(3*i-1)
270  sainvd(9*i ) = vv(3*i-2)*sainvd(9*i-6) + vv(3*i-1)*sainvd(9*i-3) + vv(3*i)
271  js= inumfi1l(i-1) + 1
272  je= inumfi1l(i )
273  do j= js, je
274  in = fi1l(j)
275  x1= vv(3*in-2)
276  x2= vv(3*in-1)
277  x3= vv(3*in )
278  sainvd(9*i-8)= sainvd(9*i-8) + x1*sainvl(9*j-8) + x2*sainvl(9*j-7) + x3*sainvl(9*j-6)
279  sainvd(9*i-4)= sainvd(9*i-4) + x1*sainvl(9*j-5) + x2*sainvl(9*j-4) + x3*sainvl(9*j-3)
280  sainvd(9*i )= sainvd(9*i ) + x1*sainvl(9*j-2) + x2*sainvl(9*j-1) + x3*sainvl(9*j )
281  enddo
282  enddo
283  !$OMP END DO
284  !$OMP END PARALLEL
285 
286  !Update D
287  dd = 1.0d0/sainvd(9*itr-8)
288 
289  sainvd(9*itr-4) =sainvd(9*itr-4)*dd
290  sainvd(9*itr ) =sainvd(9*itr )*dd
291 
292  do i =itr+1,n
293  sainvd(9*i-8) = sainvd(9*i-8)*dd
294  sainvd(9*i-4) = sainvd(9*i-4)*dd
295  sainvd(9*i ) = sainvd(9*i )*dd
296  enddo
297 
298  !Update Z
299 
300  dd2=sainvd(9*itr-4)
301  if(dabs(dd2) > filter)then
302  sainvd(9*itr-7)= sainvd(9*itr-7) - dd2*zz(3*itr-2)
303  js= inumfi1l(itr-1) + 1
304  je= inumfi1l(itr )
305  do j= js, je
306  in = fi1l(j)
307  sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
308  sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
309  sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
310  enddo
311  endif
312 
313  dd3=sainvd(9*itr )
314  if(dabs(dd3) > filter)then
315  sainvd(9*itr-6)= sainvd(9*itr-6) - dd3*zz(3*itr-2)
316  js= inumfi1l(itr-1) + 1
317  je= inumfi1l(itr )
318  do j= js, je
319  in = fi1l(j)
320  sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
321  sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
322  sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
323  enddo
324  endif
325 
326  do i= itr +1,n
327  js= inumfi1l(i-1) + 1
328  je= inumfi1l(i )
329  dd1=sainvd(9*i-8)
330  if(dabs(dd1) > filter)then
331  do j= js, je
332  in = fi1l(j)
333  if (in > itr) exit
334  sainvl(9*j-8) = sainvl(9*j-8)-dd1*zz(3*in-2)
335  sainvl(9*j-7) = sainvl(9*j-7)-dd1*zz(3*in-1)
336  sainvl(9*j-6) = sainvl(9*j-6)-dd1*zz(3*in )
337  enddo
338  endif
339  dd2=sainvd(9*i-4)
340  if(dabs(dd2) > filter)then
341  do j= js, je
342  in = fi1l(j)
343  if (in > itr) exit
344  sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
345  sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
346  sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
347  enddo
348  endif
349  dd3=sainvd(9*i )
350  if(dabs(dd3) > filter)then
351  do j= js, je
352  in = fi1l(j)
353  if (in > itr) exit
354  sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
355  sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
356  sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
357  enddo
358  endif
359  enddo
360 
361  !------------------------------ iitr = 1 ----------------------------------------
362 
363  zz(:) = 0.0d0
364  vv(:) = 0.0d0
365 
366  !{v}=[A]{zi}
367 
368  zz(3*itr-2)= sainvd(9*itr-7)
369  zz(3*itr-1)= sainvd(9*itr-4)
370  zz(3*itr )= sainvd(9*itr-1)
371 
372  zz(3*itr-1)= 1.0d0
373 
374  js= inumfi1l(itr-1) + 1
375  je= inumfi1l(itr )
376  do j= js, je
377  in = fi1l(j)
378  zz(3*in-2)= sainvl(9*j-5)
379  zz(3*in-1)= sainvl(9*j-4)
380  zz(3*in )= sainvl(9*j-3)
381  enddo
382 
383  do i= 1, itr
384  x1= zz(3*i-2)
385  x2= zz(3*i-1)
386  x3= zz(3*i )
387  vv(3*i-2) = vv(3*i-2) + d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
388  vv(3*i-1) = vv(3*i-1) + d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
389  vv(3*i ) = vv(3*i ) + d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i )*x3
390 
391  js= indexl(i-1) + 1
392  je= indexl(i )
393  do j=js,je
394  in = iteml(j)
395  vv(3*in-2)= vv(3*in-2) + al(9*j-8)*x1 + al(9*j-5)*x2 + al(9*j-2)*x3
396  vv(3*in-1)= vv(3*in-1) + al(9*j-7)*x1 + al(9*j-4)*x2 + al(9*j-1)*x3
397  vv(3*in )= vv(3*in ) + al(9*j-6)*x1 + al(9*j-3)*x2 + al(9*j )*x3
398  enddo
399 
400  js= indexu(i-1) + 1
401  je= indexu(i )
402  do j= js, je
403  in = itemu(j)
404  vv(3*in-2)= vv(3*in-2) + au(9*j-8)*x1 + au(9*j-5)*x2 + au(9*j-2)*x3
405  vv(3*in-1)= vv(3*in-1) + au(9*j-7)*x1 + au(9*j-4)*x2 + au(9*j-1)*x3
406  vv(3*in )= vv(3*in ) + au(9*j-6)*x1 + au(9*j-3)*x2 + au(9*j )*x3
407  enddo
408  enddo
409 
410  !{d}={v^t}{z_j}
411  dtemp(1) = sainvd(9*itr-8)
412 
413  !$OMP PARALLEL DEFAULT(NONE) &
414  !$OMP&PRIVATE(i,j,jS,jE,in,X1,X2,X3) &
415  !$OMP&FIRSTPRIVATE(vv) &
416  !$OMP&SHARED(N,itr,SAINVD,SAINVL,inumFI1L,FI1L)
417  !$OMP DO
418  do i=itr,n
419  sainvd(9*i-8) = vv(3*i-2)
420  sainvd(9*i-4) = vv(3*i-2)*sainvd(9*i-7) + vv(3*i-1)
421  sainvd(9*i ) = vv(3*i-2)*sainvd(9*i-6) + vv(3*i-1)*sainvd(9*i-3) + vv(3*i)
422  js= inumfi1l(i-1) + 1
423  je= inumfi1l(i )
424  do j= js, je
425  in = fi1l(j)
426  x1= vv(3*in-2)
427  x2= vv(3*in-1)
428  x3= vv(3*in )
429  sainvd(9*i-8)= sainvd(9*i-8) + x1*sainvl(9*j-8) + x2*sainvl(9*j-7) + x3*sainvl(9*j-6)
430  sainvd(9*i-4)= sainvd(9*i-4) + x1*sainvl(9*j-5) + x2*sainvl(9*j-4) + x3*sainvl(9*j-3)
431  sainvd(9*i )= sainvd(9*i ) + x1*sainvl(9*j-2) + x2*sainvl(9*j-1) + x3*sainvl(9*j )
432  enddo
433  enddo
434  !$OMP END DO
435  !$OMP END PARALLEL
436 
437  !Update D
438  dd = 1.0d0/sainvd(9*itr-4)
439 
440  sainvd(9*itr-8) = dtemp(1)
441  sainvd(9*itr ) =sainvd(9*itr )*dd
442 
443  do i =itr+1,n
444  sainvd(9*i-8) = sainvd(9*i-8)*dd
445  sainvd(9*i-4) = sainvd(9*i-4)*dd
446  sainvd(9*i ) = sainvd(9*i )*dd
447  enddo
448 
449  !Update Z
450  dd3=sainvd(9*itr )
451  if(dabs(dd3) > filter)then
452  sainvd(9*itr-6)= sainvd(9*itr-6) - dd3*zz(3*itr-2)
453  sainvd(9*itr-3)= sainvd(9*itr-3) - dd3*zz(3*itr-1)
454 
455  js= inumfi1l(itr-1) + 1
456  je= inumfi1l(itr )
457  do j= js, je
458  in = fi1l(j)
459  sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
460  sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
461  sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
462  enddo
463  endif
464 
465  do i= itr +1,n
466  js= inumfi1l(i-1) + 1
467  je= inumfi1l(i )
468  dd1=sainvd(9*i-8)
469  if(dabs(dd1) > filter)then
470  do j= js, je
471  in = fi1l(j)
472  if (in > itr) exit
473  sainvl(9*j-8) = sainvl(9*j-8)-dd1*zz(3*in-2)
474  sainvl(9*j-7) = sainvl(9*j-7)-dd1*zz(3*in-1)
475  sainvl(9*j-6) = sainvl(9*j-6)-dd1*zz(3*in )
476  enddo
477  endif
478  dd2=sainvd(9*i-4)
479  if(dabs(dd2) > filter)then
480  do j= js, je
481  in = fi1l(j)
482  if (in > itr) exit
483  sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
484  sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
485  sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
486  enddo
487  endif
488  dd3=sainvd(9*i )
489  if(dabs(dd3) > filter)then
490  do j= js, je
491  in = fi1l(j)
492  if (in > itr) exit
493  sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
494  sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
495  sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
496  enddo
497  endif
498  enddo
499 
500 
501  !------------------------------ iitr = 1 ----------------------------------------
502 
503  zz(:) = 0.0d0
504  vv(:) = 0.0d0
505 
506  !{v}=[A]{zi}
507 
508  zz(3*itr-2)= sainvd(9*itr-6)
509  zz(3*itr-1)= sainvd(9*itr-3)
510  zz(3*itr )= sainvd(9*itr )
511 
512  zz(3*itr )= 1.0d0
513 
514  js= inumfi1l(itr-1) + 1
515  je= inumfi1l(itr )
516  do j= js, je
517  in = fi1l(j)
518  zz(3*in-2)= sainvl(9*j-2)
519  zz(3*in-1)= sainvl(9*j-1)
520  zz(3*in )= sainvl(9*j )
521  enddo
522 
523  do i= 1, itr
524  x1= zz(3*i-2)
525  x2= zz(3*i-1)
526  x3= zz(3*i )
527  vv(3*i-2) = vv(3*i-2) + d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
528  vv(3*i-1) = vv(3*i-1) + d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
529  vv(3*i ) = vv(3*i ) + d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i )*x3
530 
531  js= indexl(i-1) + 1
532  je= indexl(i )
533  do j=js,je
534  in = iteml(j)
535  vv(3*in-2)= vv(3*in-2) + al(9*j-8)*x1 + al(9*j-5)*x2 + al(9*j-2)*x3
536  vv(3*in-1)= vv(3*in-1) + al(9*j-7)*x1 + al(9*j-4)*x2 + al(9*j-1)*x3
537  vv(3*in )= vv(3*in ) + al(9*j-6)*x1 + al(9*j-3)*x2 + al(9*j )*x3
538  enddo
539 
540  js= indexu(i-1) + 1
541  je= indexu(i )
542  do j= js, je
543  in = itemu(j)
544  vv(3*in-2)= vv(3*in-2) + au(9*j-8)*x1 + au(9*j-5)*x2 + au(9*j-2)*x3
545  vv(3*in-1)= vv(3*in-1) + au(9*j-7)*x1 + au(9*j-4)*x2 + au(9*j-1)*x3
546  vv(3*in )= vv(3*in ) + au(9*j-6)*x1 + au(9*j-3)*x2 + au(9*j )*x3
547  enddo
548  enddo
549 
550  !{d}={v^t}{z_j}
551  dtemp(1) = sainvd(9*itr-8)
552  dtemp(2) = sainvd(9*itr-4)
553 
554  !$OMP PARALLEL DEFAULT(NONE) &
555  !$OMP&PRIVATE(i,j,jS,jE,in,X1,X2,X3) &
556  !$OMP&FIRSTPRIVATE(vv) &
557  !$OMP&SHARED(N,itr,SAINVD,SAINVL,inumFI1L,FI1L)
558  !$OMP DO
559  do i=itr,n
560  sainvd(9*i-8) = vv(3*i-2)
561  sainvd(9*i-4) = vv(3*i-2)*sainvd(9*i-7) + vv(3*i-1)
562  sainvd(9*i ) = vv(3*i-2)*sainvd(9*i-6) + vv(3*i-1)*sainvd(9*i-3) + vv(3*i)
563  js= inumfi1l(i-1) + 1
564  je= inumfi1l(i )
565  do j= js, je
566  in = fi1l(j)
567  x1= vv(3*in-2)
568  x2= vv(3*in-1)
569  x3= vv(3*in )
570  sainvd(9*i-8)= sainvd(9*i-8) + x1*sainvl(9*j-8) + x2*sainvl(9*j-7) + x3*sainvl(9*j-6)
571  sainvd(9*i-4)= sainvd(9*i-4) + x1*sainvl(9*j-5) + x2*sainvl(9*j-4) + x3*sainvl(9*j-3)
572  sainvd(9*i )= sainvd(9*i ) + x1*sainvl(9*j-2) + x2*sainvl(9*j-1) + x3*sainvl(9*j )
573  enddo
574  enddo
575  !$OMP END DO
576  !$OMP END PARALLEL
577 
578  !Update D
579  dd = 1.0d0/sainvd(9*itr )
580 
581  sainvd(9*itr-8) = dtemp(1)
582  sainvd(9*itr-4) = dtemp(2)
583 
584  do i =itr+1,n
585  sainvd(9*i-8) = sainvd(9*i-8)*dd
586  sainvd(9*i-4) = sainvd(9*i-4)*dd
587  sainvd(9*i ) = sainvd(9*i )*dd
588  enddo
589 
590  !Update Z
591  do i= itr +1,n
592  js= inumfi1l(i-1) + 1
593  je= inumfi1l(i )
594  dd1=sainvd(9*i-8)
595  if(dabs(dd1) > filter)then
596  do j= js, je
597  in = fi1l(j)
598  if (in > itr) exit
599  sainvl(9*j-8) = sainvl(9*j-8)-dd1*zz(3*in-2)
600  sainvl(9*j-7) = sainvl(9*j-7)-dd1*zz(3*in-1)
601  sainvl(9*j-6) = sainvl(9*j-6)-dd1*zz(3*in )
602  enddo
603  endif
604  dd2=sainvd(9*i-4)
605  if(dabs(dd2) > filter)then
606  do j= js, je
607  in = fi1l(j)
608  if (in > itr) exit
609  sainvl(9*j-5) = sainvl(9*j-5)-dd2*zz(3*in-2)
610  sainvl(9*j-4) = sainvl(9*j-4)-dd2*zz(3*in-1)
611  sainvl(9*j-3) = sainvl(9*j-3)-dd2*zz(3*in )
612  enddo
613  endif
614  dd3=sainvd(9*i )
615  if(dabs(dd3) > filter)then
616  do j= js, je
617  in = fi1l(j)
618  if (in > itr) exit
619  sainvl(9*j-2) = sainvl(9*j-2)-dd3*zz(3*in-2)
620  sainvl(9*j-1) = sainvl(9*j-1)-dd3*zz(3*in-1)
621  sainvl(9*j ) = sainvl(9*j )-dd3*zz(3*in )
622  enddo
623  endif
624  enddo
625  enddo
626  deallocate(vv)
627  deallocate(zz)
628 
629  do i =1,n
630  sainvd(9*i-8) = 1.0d0/sainvd(9*i-8)
631  sainvd(9*i-4) = 1.0d0/sainvd(9*i-4)
632  sainvd(9*i ) = 1.0d0/sainvd(9*i )
633  sainvd(9*i-5) = sainvd(9*i-7)
634  sainvd(9*i-2) = sainvd(9*i-6)
635  sainvd(9*i-1) = sainvd(9*i-3)
636  enddo
637 
638  end subroutine hecmw_sainv_33
639 
640  subroutine hecmw_sainv_make_u_33(hecMAT)
641  implicit none
642  type (hecmwst_matrix) :: hecmat
643  integer(kind=kint) i,j,k,n,m,o
644  integer(kind=kint) is,ie,js,je
645 
646  n = 1
647  do i= 1, hecmat%NP
648  is=inumfi1u(i-1) + 1
649  ie=inumfi1u(i )
650  flag1:do k= is, ie
651  m = fi1u(k)
652  js=inumfi1l(m-1) + 1
653  je=inumfi1l(m )
654  do j= js,je
655  o = fi1l(j)
656  if (o == i)then
657  sainvu(9*n-8)=sainvl(9*j-8)
658  sainvu(9*n-7)=sainvl(9*j-5)
659  sainvu(9*n-6)=sainvl(9*j-2)
660  sainvu(9*n-5)=sainvl(9*j-7)
661  sainvu(9*n-4)=sainvl(9*j-4)
662  sainvu(9*n-3)=sainvl(9*j-1)
663  sainvu(9*n-2)=sainvl(9*j-6)
664  sainvu(9*n-1)=sainvl(9*j-3)
665  sainvu(9*n )=sainvl(9*j )
666  n = n + 1
667  cycle flag1
668  endif
669  enddo
670  enddo flag1
671  enddo
672  end subroutine hecmw_sainv_make_u_33
673 
674  !C***
675  !C*** FORM_ILU1_33
676  !C*** form ILU(1) matrix
677  subroutine form_ilu0_sainv_33(hecMAT)
678  implicit none
679  type(hecmwst_matrix) :: hecmat
680 
681  allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
682  allocate (fi1l(hecmat%NPL), fi1u(hecmat%NPU))
683 
684  inumfi1l = 0
685  inumfi1u = 0
686  fi1l = 0
687  fi1u = 0
688 
689  inumfi1l = hecmat%indexL
690  inumfi1u = hecmat%indexU
691  fi1l = hecmat%itemL
692  fi1u = hecmat%itemU
693 
694  npfiu = hecmat%NPU
695  npfil = hecmat%NPL
696 
697  end subroutine form_ilu0_sainv_33
698 
699  !C***
700  !C*** FORM_ILU1_33
701  !C*** form ILU(1) matrix
702  subroutine form_ilu1_sainv_33(hecMAT)
703  implicit none
704  type(hecmwst_matrix) :: hecmat
705 
706  integer(kind=kint),allocatable :: iwsl(:), iwsu(:), iw1(:), iw2(:)
707  integer(kind=kint) :: nplf1,npuf1
708  integer(kind=kint) :: i,jj,kk,l,isk,iek,isj,iej
709  integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
710  integer(kind=kint) :: j,k,isl,isu
711  !C
712  !C +--------------+
713  !C | find fill-in |
714  !C +--------------+
715  !C===
716 
717  !C
718  !C-- count fill-in
719  allocate (iw1(hecmat%NP) , iw2(hecmat%NP))
720  allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
721 
722  inumfi1l= 0
723  inumfi1u= 0
724 
725  nplf1= 0
726  npuf1= 0
727  do i= 2, hecmat%NP
728  icou= 0
729  iw1= 0
730  iw1(i)= 1
731  do l= indexl(i-1)+1, indexl(i)
732  iw1(iteml(l))= 1
733  enddo
734  do l= indexu(i-1)+1, indexu(i)
735  iw1(itemu(l))= 1
736  enddo
737 
738  isk= indexl(i-1) + 1
739  iek= indexl(i)
740  do k= isk, iek
741  kk= iteml(k)
742  isj= indexu(kk-1) + 1
743  iej= indexu(kk )
744  do j= isj, iej
745  jj= itemu(j)
746  if (iw1(jj).eq.0 .and. jj.lt.i) then
747  inumfi1l(i)= inumfi1l(i)+1
748  iw1(jj)= 1
749  endif
750  if (iw1(jj).eq.0 .and. jj.gt.i) then
751  inumfi1u(i)= inumfi1u(i)+1
752  iw1(jj)= 1
753  endif
754  enddo
755  enddo
756  nplf1= nplf1 + inumfi1l(i)
757  npuf1= npuf1 + inumfi1u(i)
758  enddo
759 
760  !C
761  !C-- specify fill-in
762  allocate (iwsl(0:hecmat%NP), iwsu(0:hecmat%NP))
763  allocate (fi1l(hecmat%NPL+nplf1), fi1u(hecmat%NPU+npuf1))
764 
765  npfiu = hecmat%NPU+npuf1
766  npfil = hecmat%NPL+nplf1
767 
768  fi1l= 0
769  fi1u= 0
770 
771  iwsl= 0
772  iwsu= 0
773  do i= 1, hecmat%NP
774  iwsl(i)= indexl(i)-indexl(i-1) + inumfi1l(i) + iwsl(i-1)
775  iwsu(i)= indexu(i)-indexu(i-1) + inumfi1u(i) + iwsu(i-1)
776  enddo
777 
778  do i= 2, hecmat%NP
779  icoul= 0
780  icouu= 0
781  inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
782  inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
783  icou= 0
784  iw1= 0
785  iw1(i)= 1
786  do l= indexl(i-1)+1, indexl(i)
787  iw1(iteml(l))= 1
788  enddo
789  do l= indexu(i-1)+1, indexu(i)
790  iw1(itemu(l))= 1
791  enddo
792 
793  isk= indexl(i-1) + 1
794  iek= indexl(i)
795  do k= isk, iek
796  kk= iteml(k)
797  isj= indexu(kk-1) + 1
798  iej= indexu(kk )
799  do j= isj, iej
800  jj= itemu(j)
801  if (iw1(jj).eq.0 .and. jj.lt.i) then
802  icoul = icoul + 1
803  fi1l(icoul+iwsl(i-1)+indexl(i)-indexl(i-1))= jj
804  iw1(jj) = 1
805  endif
806  if (iw1(jj).eq.0 .and. jj.gt.i) then
807  icouu = icouu + 1
808  fi1u(icouu+iwsu(i-1)+indexu(i)-indexu(i-1))= jj
809  iw1(jj) = 1
810  endif
811  enddo
812  enddo
813  enddo
814 
815  isl = 0
816  isu = 0
817  do i= 1, hecmat%NP
818  icoul1= indexl(i) - indexl(i-1)
819  icoul2= inumfi1l(i) - inumfi1l(i-1)
820  icoul3= icoul1 + icoul2
821  icouu1= indexu(i) - indexu(i-1)
822  icouu2= inumfi1u(i) - inumfi1u(i-1)
823  icouu3= icouu1 + icouu2
824  !C
825  !C-- LOWER part
826  icou0= 0
827  do k= indexl(i-1)+1, indexl(i)
828  icou0 = icou0 + 1
829  iw1(icou0)= iteml(k)
830  enddo
831 
832  do k= inumfi1l(i-1)+1, inumfi1l(i)
833  icou0 = icou0 + 1
834  iw1(icou0)= fi1l(icou0+iwsl(i-1))
835  enddo
836 
837  do k= 1, icoul3
838  iw2(k)= k
839  enddo
840  call sainv_sort_33 (iw1, iw2, icoul3, hecmat%NP)
841 
842  do k= 1, icoul3
843  fi1l(k+isl)= iw1(k)
844  enddo
845  !C
846  !C-- UPPER part
847  icou0= 0
848  do k= indexu(i-1)+1, indexu(i)
849  icou0 = icou0 + 1
850  iw1(icou0)= itemu(k)
851  enddo
852 
853  do k= inumfi1u(i-1)+1, inumfi1u(i)
854  icou0 = icou0 + 1
855  iw1(icou0)= fi1u(icou0+iwsu(i-1))
856  enddo
857 
858  do k= 1, icouu3
859  iw2(k)= k
860  enddo
861  call sainv_sort_33 (iw1, iw2, icouu3, hecmat%NP)
862 
863  do k= 1, icouu3
864  fi1u(k+isu)= iw1(k)
865  enddo
866 
867  isl= isl + icoul3
868  isu= isu + icouu3
869  enddo
870 
871  !C===
872  do i= 1, hecmat%NP
873  inumfi1l(i)= iwsl(i)
874  inumfi1u(i)= iwsu(i)
875  enddo
876 
877  deallocate (iw1, iw2)
878  deallocate (iwsl, iwsu)
879  !C===
880  end subroutine form_ilu1_sainv_33
881 
882  !C
883  !C***
884  !C*** fill_in_S33_SORT
885  !C***
886  !C
887  subroutine sainv_sort_33(STEM, INUM, N, NP)
888  use hecmw_util
889  implicit none
890  integer(kind=kint) :: n, np
891  integer(kind=kint), dimension(NP) :: stem
892  integer(kind=kint), dimension(NP) :: inum
893  integer(kind=kint), dimension(:), allocatable :: istack
894  integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
895 
896  allocate (istack(-np:+np))
897 
898  m = 100
899  nstack= np
900 
901  jstack= 0
902  l = 1
903  ir = n
904 
905  ip= 0
906  1 continue
907  ip= ip + 1
908 
909  if (ir-l.lt.m) then
910  do j= l+1, ir
911  ss= stem(j)
912  ii= inum(j)
913 
914  do i= j-1,1,-1
915  if (stem(i).le.ss) goto 2
916  stem(i+1)= stem(i)
917  inum(i+1)= inum(i)
918  end do
919  i= 0
920 
921  2 continue
922  stem(i+1)= ss
923  inum(i+1)= ii
924  end do
925 
926  if (jstack.eq.0) then
927  deallocate (istack)
928  return
929  endif
930 
931  ir = istack(jstack)
932  l = istack(jstack-1)
933  jstack= jstack - 2
934  else
935 
936  k= (l+ir) / 2
937  temp = stem(k)
938  stem(k) = stem(l+1)
939  stem(l+1)= temp
940 
941  it = inum(k)
942  inum(k) = inum(l+1)
943  inum(l+1)= it
944 
945  if (stem(l+1).gt.stem(ir)) then
946  temp = stem(l+1)
947  stem(l+1)= stem(ir)
948  stem(ir )= temp
949  it = inum(l+1)
950  inum(l+1)= inum(ir)
951  inum(ir )= it
952  endif
953 
954  if (stem(l).gt.stem(ir)) then
955  temp = stem(l)
956  stem(l )= stem(ir)
957  stem(ir)= temp
958  it = inum(l)
959  inum(l )= inum(ir)
960  inum(ir)= it
961  endif
962 
963  if (stem(l+1).gt.stem(l)) then
964  temp = stem(l+1)
965  stem(l+1)= stem(l)
966  stem(l )= temp
967  it = inum(l+1)
968  inum(l+1)= inum(l)
969  inum(l )= it
970  endif
971 
972  i= l + 1
973  j= ir
974 
975  ss= stem(l)
976  ii= inum(l)
977 
978  3 continue
979  i= i + 1
980  if (stem(i).lt.ss) goto 3
981 
982  4 continue
983  j= j - 1
984  if (stem(j).gt.ss) goto 4
985 
986  if (j.lt.i) goto 5
987 
988  temp = stem(i)
989  stem(i)= stem(j)
990  stem(j)= temp
991 
992  it = inum(i)
993  inum(i)= inum(j)
994  inum(j)= it
995 
996  goto 3
997 
998  5 continue
999 
1000  stem(l)= stem(j)
1001  stem(j)= ss
1002  inum(l)= inum(j)
1003  inum(j)= ii
1004 
1005  jstack= jstack + 2
1006 
1007  if (jstack.gt.nstack) then
1008  write (*,*) 'NSTACK overflow'
1009  stop
1010  endif
1011 
1012  if (ir-i+1.ge.j-1) then
1013  istack(jstack )= ir
1014  istack(jstack-1)= i
1015  ir= j-1
1016  else
1017  istack(jstack )= j-1
1018  istack(jstack-1)= l
1019  l= i
1020  endif
1021 
1022  endif
1023 
1024  goto 1
1025 
1026  end subroutine sainv_sort_33
1027 
1028  subroutine hecmw_precond_sainv_33_clear()
1029  implicit none
1030 
1031  if (associated(sainvd)) deallocate(sainvd)
1032  if (associated(sainvl)) deallocate(sainvl)
1033  if (associated(sainvu)) deallocate(sainvu)
1034  if (associated(inumfi1l)) deallocate(inumfi1l)
1035  if (associated(inumfi1u)) deallocate(inumfi1u)
1036  if (associated(fi1l)) deallocate(fi1l)
1037  if (associated(fi1u)) deallocate(fi1u)
1038  nullify(inumfi1l)
1039  nullify(inumfi1u)
1040  nullify(fi1l)
1041  nullify(fi1u)
1042  nullify(d)
1043  nullify(al)
1044  nullify(au)
1045  nullify(indexl)
1046  nullify(indexu)
1047  nullify(iteml)
1048  nullify(itemu)
1049 
1050  end subroutine hecmw_precond_sainv_33_clear
1051 end module hecmw_precond_sainv_33
hecmw_precond_sainv_33::hecmw_precond_sainv_33_apply
subroutine, public hecmw_precond_sainv_33_apply(R, ZP)
Definition: hecmw_precond_SAINV_33.f90:120
hecmw_matrix_misc::hecmw_mat_get_precond
integer(kind=kint) function, public hecmw_mat_get_precond(hecMAT)
Definition: hecmw_matrix_misc.f90:321
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
m_hecmw_comm_f
Definition: hecmw_comm_f.F90:6
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
hecmw_precond_sainv_33
Definition: hecmw_precond_SAINV_33.f90:6
hecmw_matrix_misc
Definition: hecmw_matrix_misc.f90:6
hecmw_precond_sainv_33::hecmw_precond_sainv_33_clear
subroutine, public hecmw_precond_sainv_33_clear()
Definition: hecmw_precond_SAINV_33.f90:1029
hecmw_precond_sainv_33::hecmw_precond_sainv_33_setup
subroutine, public hecmw_precond_sainv_33_setup(hecMAT)
Definition: hecmw_precond_SAINV_33.f90:46
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444