FrontISTR  5.9.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 
648  !C***
649  !C*** FORM_ILU1_33
650  !C*** form ILU(1) matrix
651  subroutine form_ilu1_rif_33(hecMAT)
652  implicit none
653  type(hecmwst_matrix) :: hecmat
654 
655  integer(kind=kint),allocatable :: iwsl(:), iwsu(:), iw1(:), iw2(:)
656  integer(kind=kint) :: nplf1,npuf1
657  integer(kind=kint) :: i,jj,kk,l,isk,iek,isj,iej
658  integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
659  integer(kind=kint) :: j,k,isl,isu
660  !C
661  !C +--------------+
662  !C | find fill-in |
663  !C +--------------+
664  !C===
665 
666  !C
667  !C-- count fill-in
668  allocate (iw1(hecmat%NP) , iw2(hecmat%NP))
669  allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
670 
671  inumfi1l= 0
672  inumfi1u= 0
673 
674  nplf1= 0
675  npuf1= 0
676  do i= 2, hecmat%NP
677  icou= 0
678  iw1= 0
679  iw1(i)= 1
680  do l= indexl(i-1)+1, indexl(i)
681  iw1(iteml(l))= 1
682  enddo
683  do l= indexu(i-1)+1, indexu(i)
684  iw1(itemu(l))= 1
685  enddo
686 
687  isk= indexl(i-1) + 1
688  iek= indexl(i)
689  do k= isk, iek
690  kk= iteml(k)
691  isj= indexu(kk-1) + 1
692  iej= indexu(kk )
693  do j= isj, iej
694  jj= itemu(j)
695  if (iw1(jj).eq.0 .and. jj.lt.i) then
696  inumfi1l(i)= inumfi1l(i)+1
697  iw1(jj)= 1
698  endif
699  if (iw1(jj).eq.0 .and. jj.gt.i) then
700  inumfi1u(i)= inumfi1u(i)+1
701  iw1(jj)= 1
702  endif
703  enddo
704  enddo
705  nplf1= nplf1 + inumfi1l(i)
706  npuf1= npuf1 + inumfi1u(i)
707  enddo
708 
709  !C
710  !C-- specify fill-in
711  allocate (iwsl(0:hecmat%NP), iwsu(0:hecmat%NP))
712  allocate (fi1l(hecmat%NPL+nplf1), fi1u(hecmat%NPU+npuf1))
713 
714  npfiu = hecmat%NPU+npuf1
715  npfil = hecmat%NPL+nplf1
716 
717  fi1l= 0
718  fi1u= 0
719 
720  iwsl= 0
721  iwsu= 0
722  do i= 1, hecmat%NP
723  iwsl(i)= indexl(i)-indexl(i-1) + inumfi1l(i) + iwsl(i-1)
724  iwsu(i)= indexu(i)-indexu(i-1) + inumfi1u(i) + iwsu(i-1)
725  enddo
726 
727  do i= 2, hecmat%NP
728  icoul= 0
729  icouu= 0
730  inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
731  inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
732  icou= 0
733  iw1= 0
734  iw1(i)= 1
735  do l= indexl(i-1)+1, indexl(i)
736  iw1(iteml(l))= 1
737  enddo
738  do l= indexu(i-1)+1, indexu(i)
739  iw1(itemu(l))= 1
740  enddo
741 
742  isk= indexl(i-1) + 1
743  iek= indexl(i)
744  do k= isk, iek
745  kk= iteml(k)
746  isj= indexu(kk-1) + 1
747  iej= indexu(kk )
748  do j= isj, iej
749  jj= itemu(j)
750  if (iw1(jj).eq.0 .and. jj.lt.i) then
751  icoul = icoul + 1
752  fi1l(icoul+iwsl(i-1)+indexl(i)-indexl(i-1))= jj
753  iw1(jj) = 1
754  endif
755  if (iw1(jj).eq.0 .and. jj.gt.i) then
756  icouu = icouu + 1
757  fi1u(icouu+iwsu(i-1)+indexu(i)-indexu(i-1))= jj
758  iw1(jj) = 1
759  endif
760  enddo
761  enddo
762  enddo
763 
764  isl = 0
765  isu = 0
766  do i= 1, hecmat%NP
767  icoul1= indexl(i) - indexl(i-1)
768  icoul2= inumfi1l(i) - inumfi1l(i-1)
769  icoul3= icoul1 + icoul2
770  icouu1= indexu(i) - indexu(i-1)
771  icouu2= inumfi1u(i) - inumfi1u(i-1)
772  icouu3= icouu1 + icouu2
773  !C
774  !C-- LOWER part
775  icou0= 0
776  do k= indexl(i-1)+1, indexl(i)
777  icou0 = icou0 + 1
778  iw1(icou0)= iteml(k)
779  enddo
780 
781  do k= inumfi1l(i-1)+1, inumfi1l(i)
782  icou0 = icou0 + 1
783  iw1(icou0)= fi1l(icou0+iwsl(i-1))
784  enddo
785 
786  do k= 1, icoul3
787  iw2(k)= k
788  enddo
789  call rif_sort_33 (iw1, iw2, icoul3, hecmat%NP)
790 
791  do k= 1, icoul3
792  fi1l(k+isl)= iw1(k)
793  enddo
794  !C
795  !C-- UPPER part
796  icou0= 0
797  do k= indexu(i-1)+1, indexu(i)
798  icou0 = icou0 + 1
799  iw1(icou0)= itemu(k)
800  enddo
801 
802  do k= inumfi1u(i-1)+1, inumfi1u(i)
803  icou0 = icou0 + 1
804  iw1(icou0)= fi1u(icou0+iwsu(i-1))
805  enddo
806 
807  do k= 1, icouu3
808  iw2(k)= k
809  enddo
810  call rif_sort_33 (iw1, iw2, icouu3, hecmat%NP)
811 
812  do k= 1, icouu3
813  fi1u(k+isu)= iw1(k)
814  enddo
815 
816  isl= isl + icoul3
817  isu= isu + icouu3
818  enddo
819 
820  !C===
821  do i= 1, hecmat%NP
822  inumfi1l(i)= iwsl(i)
823  inumfi1u(i)= iwsu(i)
824  enddo
825 
826  deallocate (iw1, iw2)
827  deallocate (iwsl, iwsu)
828  !C===
829  end subroutine form_ilu1_rif_33
830 
831  !C
832  !C***
833  !C*** fill_in_S33_SORT
834  !C***
835  !C
836  subroutine rif_sort_33(STEM, INUM, N, NP)
837  use hecmw_util
838  implicit none
839  integer(kind=kint) :: n, np
840  integer(kind=kint), dimension(NP) :: stem
841  integer(kind=kint), dimension(NP) :: inum
842  integer(kind=kint), dimension(:), allocatable :: istack
843  integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
844 
845  allocate (istack(-np:+np))
846 
847  m = 100
848  nstack= np
849 
850  jstack= 0
851  l = 1
852  ir = n
853 
854  ip= 0
855  1 continue
856  ip= ip + 1
857 
858  if (ir-l.lt.m) then
859  do j= l+1, ir
860  ss= stem(j)
861  ii= inum(j)
862 
863  do i= j-1,1,-1
864  if (stem(i).le.ss) goto 2
865  stem(i+1)= stem(i)
866  inum(i+1)= inum(i)
867  end do
868  i= 0
869 
870  2 continue
871  stem(i+1)= ss
872  inum(i+1)= ii
873  end do
874 
875  if (jstack.eq.0) then
876  deallocate (istack)
877  return
878  endif
879 
880  ir = istack(jstack)
881  l = istack(jstack-1)
882  jstack= jstack - 2
883  else
884 
885  k= (l+ir) / 2
886  temp = stem(k)
887  stem(k) = stem(l+1)
888  stem(l+1)= temp
889 
890  it = inum(k)
891  inum(k) = inum(l+1)
892  inum(l+1)= it
893 
894  if (stem(l+1).gt.stem(ir)) then
895  temp = stem(l+1)
896  stem(l+1)= stem(ir)
897  stem(ir )= temp
898  it = inum(l+1)
899  inum(l+1)= inum(ir)
900  inum(ir )= it
901  endif
902 
903  if (stem(l).gt.stem(ir)) then
904  temp = stem(l)
905  stem(l )= stem(ir)
906  stem(ir)= temp
907  it = inum(l)
908  inum(l )= inum(ir)
909  inum(ir)= it
910  endif
911 
912  if (stem(l+1).gt.stem(l)) then
913  temp = stem(l+1)
914  stem(l+1)= stem(l)
915  stem(l )= temp
916  it = inum(l+1)
917  inum(l+1)= inum(l)
918  inum(l )= it
919  endif
920 
921  i= l + 1
922  j= ir
923 
924  ss= stem(l)
925  ii= inum(l)
926 
927  3 continue
928  i= i + 1
929  if (stem(i).lt.ss) goto 3
930 
931  4 continue
932  j= j - 1
933  if (stem(j).gt.ss) goto 4
934 
935  if (j.lt.i) goto 5
936 
937  temp = stem(i)
938  stem(i)= stem(j)
939  stem(j)= temp
940 
941  it = inum(i)
942  inum(i)= inum(j)
943  inum(j)= it
944 
945  goto 3
946 
947  5 continue
948 
949  stem(l)= stem(j)
950  stem(j)= ss
951  inum(l)= inum(j)
952  inum(j)= ii
953 
954  jstack= jstack + 2
955 
956  if (jstack.gt.nstack) then
957  write (*,*) 'NSTACK overflow'
958  stop
959  endif
960 
961  if (ir-i+1.ge.j-1) then
962  istack(jstack )= ir
963  istack(jstack-1)= i
964  ir= j-1
965  else
966  istack(jstack )= j-1
967  istack(jstack-1)= l
968  l= i
969  endif
970 
971  endif
972 
973  goto 1
974 
975  end subroutine rif_sort_33
976 
977  subroutine hecmw_precond_rif_33_clear()
978  implicit none
979 
980  if (associated(sainvd)) deallocate(sainvd)
981  if (associated(sainvl)) deallocate(sainvl)
982  if (associated(rifd)) deallocate(rifd)
983  if (associated(rifu)) deallocate(rifu)
984  if (associated(rifl)) deallocate(rifl)
985  if (associated(inumfi1l)) deallocate(inumfi1l)
986  if (associated(inumfi1u)) deallocate(inumfi1u)
987  if (associated(fi1l)) deallocate(fi1l)
988  if (associated(fi1u)) deallocate(fi1u)
989  nullify(inumfi1l)
990  nullify(inumfi1u)
991  nullify(fi1l)
992  nullify(fi1u)
993  nullify(d)
994  nullify(al)
995  nullify(au)
996  nullify(indexl)
997  nullify(indexu)
998  nullify(iteml)
999  nullify(itemu)
1000 
1001  end subroutine hecmw_precond_rif_33_clear
1002 end module hecmw_precond_rif_33
hecmw_matrix_misc
Definition: hecmw_matrix_misc.F90:6
hecmw_matrix_misc::hecmw_mat_get_precond
integer(kind=kint) function, public hecmw_mat_get_precond(hecMAT)
Definition: hecmw_matrix_misc.F90:379
hecmw_precond_rif_33
Definition: hecmw_precond_RIF_33.f90:6
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_precond_rif_33::hecmw_precond_rif_33_setup
subroutine, public hecmw_precond_rif_33_setup(hecMAT)
Definition: hecmw_precond_RIF_33.f90:46
hecmw_precond_rif_33::hecmw_precond_rif_33_clear
subroutine, public hecmw_precond_rif_33_clear()
Definition: hecmw_precond_RIF_33.f90:978
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_hecmw_comm_f
Definition: hecmw_comm_f.F90:6
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:467