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