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