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