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