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