FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_precond_SAINV_nn.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5 
6 module hecmw_precond_sainv_nn
7  use hecmw_util
8  use m_hecmw_comm_f
9  use hecmw_matrix_misc
10 
11  private
12 
13  public:: hecmw_precond_sainv_nn_setup
14  public:: hecmw_precond_sainv_nn_apply
15  public:: hecmw_precond_sainv_nn_clear
16 
17  integer(4),parameter :: krealp = 8
18 
19  integer(kind=kint) :: NPFIU, NPFIL
20  integer(kind=kint) :: N
21  integer(kind=kint) :: NDOF, NDOF2
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 :: sainvu(:) => null()
36  real(kind=krealp), pointer :: sainvl(:) => null()
37  real(kind=krealp), pointer :: sainvd(:) => null()
38  real(kind=kreal), pointer :: t(:) => null()
39 
40 contains
41 
42  !C***
43  !C*** hecmw_precond_nn_sainv_setup
44  !C***
45  subroutine hecmw_precond_sainv_nn_setup(hecMAT)
46  implicit none
47  type(hecmwst_matrix) :: hecmat
48 
49  integer(kind=kint ) :: precond
50  real(kind=krealp) :: filter
51 
52  n = hecmat%N
53  ndof = hecmat%NDOF
54  ndof2 = ndof*ndof
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_nn(hecmat)
66 
67  allocate (sainvd(ndof2*hecmat%NP))
68  allocate (sainvl(ndof2*npfiu))
69  allocate (t(ndof*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_nn(hecmat)
79 
80  allocate (sainvu(ndof2*npfiu))
81  sainvu = 0.0d0
82 
83  call hecmw_sainv_make_u_nn(hecmat)
84 
85  end subroutine hecmw_precond_sainv_nn_setup
86 
87  subroutine hecmw_precond_sainv_nn_apply(R, ZP)
88  implicit none
89  real(kind=kreal), intent(inout) :: zp(:)
90  real(kind=kreal), intent(in) :: r(:)
91  integer(kind=kint) :: in, i, j, isl, iel, isu, ieu,idof,jdof
92  real(kind=kreal) :: sw(ndof),x(ndof)
93 
94  !$OMP PARALLEL DEFAULT(NONE) &
95  !$OMP&PRIVATE(i,X,SW,j,in,isL,ieL,isU,ieU,idof,jdof) &
96  !$OMP&SHARED(N,SAINVD,SAINVL,SAINVU,inumFI1U,FI1U,inumFI1L,FI1L,R,T,ZP,NDOF,NDOF2)
97  !$OMP DO
98  !C-- FORWARD
99  do i= 1, n
100  do idof = 1, ndof
101  sw(idof) = 0.0d0
102  end do
103  isl= inumfi1l(i-1)+1
104  iel= inumfi1l(i)
105  do j= isl, iel
106  in= fi1l(j)
107  do idof = 1, ndof
108  x(idof) = r(ndof*(in-1)+idof)
109  end do
110  do idof = 1, ndof
111  do jdof = 1, ndof
112  sw(idof) = sw(idof) + sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
113  end do
114  end do
115  enddo
116  do idof = 1, ndof
117  x(idof) = r(ndof*(i-1)+idof)
118  t(ndof*(i-1)+idof)=x(idof)+sw(idof)
119  end do
120  do idof = 1, ndof
121  do jdof = 1, idof-1
122  t(ndof*(i-1)+idof)=t(ndof*(i-1)+idof)+sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)*x(jdof)
123  end do
124  t(ndof*(i-1)+idof)=t(ndof*(i-1)+idof)*sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
125  end do
126  enddo
127  !$OMP END DO
128  !$OMP DO
129  !C-- BACKWARD
130  do i= 1, n
131  do idof = 1, ndof
132  sw(idof) = 0.0d0
133  end do
134 
135  isu= inumfi1u(i-1) + 1
136  ieu= inumfi1u(i)
137  do j= isu, ieu
138  in= fi1u(j)
139  do idof = 1, ndof
140  x(idof) = t(ndof*(in-1)+idof)
141  end do
142  do idof = 1, ndof
143  do jdof = 1, ndof
144  sw(idof) = sw(idof) + sainvu(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
145  end do
146  end do
147  enddo
148 
149  do idof = 1, ndof
150  x(idof) = t(ndof*(i-1)+idof)
151  end do
152  do idof = 1, ndof
153  zp(ndof*(i-1)+idof) = x(idof) + sw(idof)
154  do jdof = ndof, idof+1, -1
155  zp(ndof*(i-1)+idof) = zp(ndof*(i-1)+idof)+sainvd(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(jdof)
156  end do
157  end do
158  enddo
159  !$OMP END DO
160  !$OMP END PARALLEL
161 
162  end subroutine hecmw_precond_sainv_nn_apply
163 
164 
165  !C***
166  !C*** hecmw_rif_nn
167  !C***
168  subroutine hecmw_sainv_nn(hecMAT)
169  implicit none
170  type (hecmwst_matrix) :: hecmat
171 
172  integer(kind=kint) :: i, j, k, js, je, in, itr, idof, jdof, iitr
173  real(kind=krealp) :: dd, dtmp(hecmat%NDOF), x(hecmat%NDOF)
174 
175  real(kind=krealp) :: filter
176  real(kind=krealp), allocatable :: zz(:), vv(:)
177 
178  filter= hecmat%Rarray(5)
179 
180  allocate (vv(ndof*hecmat%NP))
181  allocate (zz(ndof*hecmat%NP))
182  do itr=1,n
183  do iitr=1,ndof
184  zz(:) = 0.0d0
185  vv(:) = 0.0d0
186  !{v}=[A]{zi}
187  do idof = 1,ndof
188  zz(ndof*(itr-1)+idof)= sainvd(ndof2*(itr-1)+ndof*(idof-1)+iitr)
189  end do
190  zz(ndof*(itr-1)+iitr)= 1.0d0
191 
192  js= inumfi1l(itr-1) + 1
193  je= inumfi1l(itr )
194  do j= js, je
195  in = fi1l(j)
196  do idof = 1, ndof
197  zz(ndof*(in-1)+idof)=sainvl(ndof2*(j-1)+ndof*(iitr-1)+idof)
198  end do
199  enddo
200 
201  do i= 1, itr
202  do idof = 1,ndof
203  x(idof)=zz(ndof*(i-1)+idof)
204  end do
205  do idof = 1, ndof
206  do jdof = 1, ndof
207  vv(ndof*(i-1)+idof) = vv(ndof*(i-1)+idof) + d(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(jdof)
208  end do
209  end do
210 
211  js= indexl(i-1) + 1
212  je= indexl(i )
213  do j=js,je
214  in = iteml(j)
215  do idof = 1, ndof
216  do jdof = 1, ndof
217  vv(ndof*(in-1)+idof) = vv(ndof*(in-1)+idof) + al(ndof2*(j-1)+ndof*(jdof-1)+idof)*x(jdof)
218  end do
219  end do
220  enddo
221  js= indexu(i-1) + 1
222  je= indexu(i )
223  do j= js, je
224  in = itemu(j)
225  do idof = 1, ndof
226  do jdof = 1, ndof
227  vv(ndof*(in-1)+idof) = vv(ndof*(in-1)+idof) + au(ndof2*(j-1)+ndof*(jdof-1)+idof)*x(jdof)
228  end do
229  end do
230  enddo
231  enddo
232 
233  !{d}={v^t}{z_j}
234  do idof = 1, ndof
235  dtmp(idof)= sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof)
236  end do
237 
238  do i= itr,n
239  do idof = 1,ndof
240  sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = vv(ndof*(i-1)+idof)
241  do jdof = 1, idof-1
242  sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) + &
243  & sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)*vv(ndof*(i-1)+jdof)
244  end do
245  end do
246  js= inumfi1l(i-1) + 1
247  je= inumfi1l(i )
248  do j= js, je
249  in = fi1l(j)
250  do idof = 1,ndof
251  x(idof)=vv(ndof*(in-1)+idof)
252  end do
253  do idof = 1, ndof
254  do jdof = 1, ndof
255  sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) + &
256  & sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
257  end do
258  end do
259  enddo
260  enddo
261 
262  !Update D
263  dd = 1.0d0/sainvd(ndof2*(itr-1)+ndof*(iitr-1)+iitr)
264  do idof=1,iitr-1
265  sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof) = dtmp(idof)
266  end do
267  ! SAINVD(NDOF2*(itr-1)+NDOF*(iitr-1)+iitr)=dd
268  do idof = iitr+1, ndof
269  sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof)*dd
270  end do
271 
272  do i =itr+1,n
273  do idof = 1, ndof
274  sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)*dd
275  end do
276  enddo
277 
278  !Update Z
279  do k=iitr+1,ndof
280  dd = sainvd(ndof2*(itr-1)+ndof*(k-1)+k)
281  if(abs(dd) > filter)then
282  do jdof = 1, iitr
283  sainvd(ndof2*(itr-1)+ndof*(jdof-1)+k)= sainvd(ndof2*(itr-1)+ndof*(jdof-1)+k) - dd*zz(ndof*(itr-1)+jdof)
284  end do
285  js= inumfi1l(itr-1) + 1
286  je= inumfi1l(itr )
287  do j= js, je
288  in = fi1l(j)
289  do idof = 1, ndof
290  sainvl(ndof2*(j-1)+ndof*(k-1)+idof) = sainvl(ndof2*(j-1)+ndof*(k-1)+idof)-dd*zz(ndof*(in-1)+idof)
291  end do
292  enddo
293  endif
294  end do
295 
296  do i= itr +1,n
297  js= inumfi1l(i-1) + 1
298  je= inumfi1l(i )
299  do idof = 1, ndof
300  dd = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
301  if(abs(dd) > filter)then
302  do j= js, je
303  in = fi1l(j)
304  if (in > itr) exit
305  do jdof=1,ndof
306  sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)=sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)-dd*zz(ndof*(in-1)+jdof)
307  end do
308  enddo
309  endif
310  end do
311  enddo
312  end do
313  enddo
314  deallocate(vv)
315  deallocate(zz)
316 
317  do i =1,n
318  do idof = 1, ndof
319  sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)=1/sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
320  do jdof = idof+1, ndof
321  sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)=sainvd(ndof2*(i-1)+ndof*(idof-1)+jdof)
322  end do
323  end do
324  enddo
325  end subroutine hecmw_sainv_nn
326 
327  subroutine hecmw_sainv_make_u_nn(hecMAT)
328  implicit none
329  type (hecmwst_matrix) :: hecmat
330  integer(kind=kint) i,j,k,n,m,o,idof,jdof
331  integer(kind=kint) is,ie,js,je
332 
333  n = 1
334  do i= 1, hecmat%NP
335  is=inumfi1u(i-1) + 1
336  ie=inumfi1u(i )
337  flag1:do k= is, ie
338  m = fi1u(k)
339  js=inumfi1l(m-1) + 1
340  je=inumfi1l(m )
341  do j= js,je
342  o = fi1l(j)
343  if (o == i)then
344  do idof = 1, ndof
345  do jdof = 1, ndof
346  sainvu(ndof2*(n-1)+ndof*(jdof-1)+idof)=sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)
347  end do
348  end do
349  n = n + 1
350  cycle flag1
351  endif
352  enddo
353  enddo flag1
354  enddo
355  end subroutine hecmw_sainv_make_u_nn
356 
357  !C***
358  !C*** FORM_ILU1_nn
359  !C*** form ILU(1) matrix
360  subroutine form_ilu0_sainv_nn(hecMAT)
361  implicit none
362  type(hecmwst_matrix) :: hecmat
363 
364  allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
365  allocate (fi1l(hecmat%NPL), fi1u(hecmat%NPU))
366 
367  inumfi1l = 0
368  inumfi1u = 0
369  fi1l = 0
370  fi1u = 0
371 
372  inumfi1l = hecmat%indexL
373  inumfi1u = hecmat%indexU
374  fi1l = hecmat%itemL
375  fi1u = hecmat%itemU
376 
377  npfiu = hecmat%NPU
378  npfil = hecmat%NPL
379 
380  end subroutine form_ilu0_sainv_nn
381 
382  !C***
383  !C*** FORM_ILU1_nn
384  !C*** form ILU(1) matrix
385  subroutine form_ilu1_sainv_nn(hecMAT)
386  implicit none
387  type(hecmwst_matrix) :: hecmat
388 
389  integer(kind=kint),allocatable :: iwsl(:), iwsu(:), iw1(:), iw2(:)
390  integer(kind=kint) :: nplf1,npuf1
391  integer(kind=kint) :: i,jj,kk,l,isk,iek,isj,iej
392  integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
393  integer(kind=kint) :: j,k,isl,isu
394  !C
395  !C +--------------+
396  !C | find fill-in |
397  !C +--------------+
398  !C===
399 
400  !C
401  !C-- count fill-in
402  allocate (iw1(hecmat%NP) , iw2(hecmat%NP))
403  allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
404 
405  inumfi1l= 0
406  inumfi1u= 0
407 
408  nplf1= 0
409  npuf1= 0
410  do i= 2, hecmat%NP
411  icou= 0
412  iw1= 0
413  iw1(i)= 1
414  do l= indexl(i-1)+1, indexl(i)
415  iw1(iteml(l))= 1
416  enddo
417  do l= indexu(i-1)+1, indexu(i)
418  iw1(itemu(l))= 1
419  enddo
420 
421  isk= indexl(i-1) + 1
422  iek= indexl(i)
423  do k= isk, iek
424  kk= iteml(k)
425  isj= indexu(kk-1) + 1
426  iej= indexu(kk )
427  do j= isj, iej
428  jj= itemu(j)
429  if (iw1(jj).eq.0 .and. jj.lt.i) then
430  inumfi1l(i)= inumfi1l(i)+1
431  iw1(jj)= 1
432  endif
433  if (iw1(jj).eq.0 .and. jj.gt.i) then
434  inumfi1u(i)= inumfi1u(i)+1
435  iw1(jj)= 1
436  endif
437  enddo
438  enddo
439  nplf1= nplf1 + inumfi1l(i)
440  npuf1= npuf1 + inumfi1u(i)
441  enddo
442 
443  !C
444  !C-- specify fill-in
445  allocate (iwsl(0:hecmat%NP), iwsu(0:hecmat%NP))
446  allocate (fi1l(hecmat%NPL+nplf1), fi1u(hecmat%NPU+npuf1))
447 
448  npfiu = hecmat%NPU+npuf1
449  npfil = hecmat%NPL+nplf1
450 
451  fi1l= 0
452  fi1u= 0
453 
454  iwsl= 0
455  iwsu= 0
456  do i= 1, hecmat%NP
457  iwsl(i)= indexl(i)-indexl(i-1) + inumfi1l(i) + iwsl(i-1)
458  iwsu(i)= indexu(i)-indexu(i-1) + inumfi1u(i) + iwsu(i-1)
459  enddo
460 
461  do i= 2, hecmat%NP
462  icoul= 0
463  icouu= 0
464  inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
465  inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
466  icou= 0
467  iw1= 0
468  iw1(i)= 1
469  do l= indexl(i-1)+1, indexl(i)
470  iw1(iteml(l))= 1
471  enddo
472  do l= indexu(i-1)+1, indexu(i)
473  iw1(itemu(l))= 1
474  enddo
475 
476  isk= indexl(i-1) + 1
477  iek= indexl(i)
478  do k= isk, iek
479  kk= iteml(k)
480  isj= indexu(kk-1) + 1
481  iej= indexu(kk )
482  do j= isj, iej
483  jj= itemu(j)
484  if (iw1(jj).eq.0 .and. jj.lt.i) then
485  icoul = icoul + 1
486  fi1l(icoul+iwsl(i-1)+indexl(i)-indexl(i-1))= jj
487  iw1(jj) = 1
488  endif
489  if (iw1(jj).eq.0 .and. jj.gt.i) then
490  icouu = icouu + 1
491  fi1u(icouu+iwsu(i-1)+indexu(i)-indexu(i-1))= jj
492  iw1(jj) = 1
493  endif
494  enddo
495  enddo
496  enddo
497 
498  isl = 0
499  isu = 0
500  do i= 1, hecmat%NP
501  icoul1= indexl(i) - indexl(i-1)
502  icoul2= inumfi1l(i) - inumfi1l(i-1)
503  icoul3= icoul1 + icoul2
504  icouu1= indexu(i) - indexu(i-1)
505  icouu2= inumfi1u(i) - inumfi1u(i-1)
506  icouu3= icouu1 + icouu2
507  !C
508  !C-- LOWER part
509  icou0= 0
510  do k= indexl(i-1)+1, indexl(i)
511  icou0 = icou0 + 1
512  iw1(icou0)= iteml(k)
513  enddo
514 
515  do k= inumfi1l(i-1)+1, inumfi1l(i)
516  icou0 = icou0 + 1
517  iw1(icou0)= fi1l(icou0+iwsl(i-1))
518  enddo
519 
520  do k= 1, icoul3
521  iw2(k)= k
522  enddo
523  call sainv_sort_nn (iw1, iw2, icoul3, hecmat%NP)
524 
525  do k= 1, icoul3
526  fi1l(k+isl)= iw1(k)
527  enddo
528  !C
529  !C-- UPPER part
530  icou0= 0
531  do k= indexu(i-1)+1, indexu(i)
532  icou0 = icou0 + 1
533  iw1(icou0)= itemu(k)
534  enddo
535 
536  do k= inumfi1u(i-1)+1, inumfi1u(i)
537  icou0 = icou0 + 1
538  iw1(icou0)= fi1u(icou0+iwsu(i-1))
539  enddo
540 
541  do k= 1, icouu3
542  iw2(k)= k
543  enddo
544  call sainv_sort_nn (iw1, iw2, icouu3, hecmat%NP)
545 
546  do k= 1, icouu3
547  fi1u(k+isu)= iw1(k)
548  enddo
549 
550  isl= isl + icoul3
551  isu= isu + icouu3
552  enddo
553 
554  !C===
555  do i= 1, hecmat%NP
556  inumfi1l(i)= iwsl(i)
557  inumfi1u(i)= iwsu(i)
558  enddo
559 
560  deallocate (iw1, iw2)
561  deallocate (iwsl, iwsu)
562  !C===
563  end subroutine form_ilu1_sainv_nn
564 
565  !C
566  !C***
567  !C*** fill_in_S33_SORT
568  !C***
569  !C
570  subroutine sainv_sort_nn(STEM, INUM, N, NP)
571  use hecmw_util
572  implicit none
573  integer(kind=kint) :: n, np
574  integer(kind=kint), dimension(NP) :: stem
575  integer(kind=kint), dimension(NP) :: inum
576  integer(kind=kint), dimension(:), allocatable :: istack
577  integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
578 
579  allocate (istack(-np:+np))
580 
581  m = 100
582  nstack= np
583 
584  jstack= 0
585  l = 1
586  ir = n
587 
588  ip= 0
589  1 continue
590  ip= ip + 1
591 
592  if (ir-l.lt.m) then
593  do j= l+1, ir
594  ss= stem(j)
595  ii= inum(j)
596 
597  do i= j-1,1,-1
598  if (stem(i).le.ss) goto 2
599  stem(i+1)= stem(i)
600  inum(i+1)= inum(i)
601  end do
602  i= 0
603 
604  2 continue
605  stem(i+1)= ss
606  inum(i+1)= ii
607  end do
608 
609  if (jstack.eq.0) then
610  deallocate (istack)
611  return
612  endif
613 
614  ir = istack(jstack)
615  l = istack(jstack-1)
616  jstack= jstack - 2
617  else
618 
619  k= (l+ir) / 2
620  temp = stem(k)
621  stem(k) = stem(l+1)
622  stem(l+1)= temp
623 
624  it = inum(k)
625  inum(k) = inum(l+1)
626  inum(l+1)= it
627 
628  if (stem(l+1).gt.stem(ir)) then
629  temp = stem(l+1)
630  stem(l+1)= stem(ir)
631  stem(ir )= temp
632  it = inum(l+1)
633  inum(l+1)= inum(ir)
634  inum(ir )= it
635  endif
636 
637  if (stem(l).gt.stem(ir)) then
638  temp = stem(l)
639  stem(l )= stem(ir)
640  stem(ir)= temp
641  it = inum(l)
642  inum(l )= inum(ir)
643  inum(ir)= it
644  endif
645 
646  if (stem(l+1).gt.stem(l)) then
647  temp = stem(l+1)
648  stem(l+1)= stem(l)
649  stem(l )= temp
650  it = inum(l+1)
651  inum(l+1)= inum(l)
652  inum(l )= it
653  endif
654 
655  i= l + 1
656  j= ir
657 
658  ss= stem(l)
659  ii= inum(l)
660 
661  3 continue
662  i= i + 1
663  if (stem(i).lt.ss) goto 3
664 
665  4 continue
666  j= j - 1
667  if (stem(j).gt.ss) goto 4
668 
669  if (j.lt.i) goto 5
670 
671  temp = stem(i)
672  stem(i)= stem(j)
673  stem(j)= temp
674 
675  it = inum(i)
676  inum(i)= inum(j)
677  inum(j)= it
678 
679  goto 3
680 
681  5 continue
682 
683  stem(l)= stem(j)
684  stem(j)= ss
685  inum(l)= inum(j)
686  inum(j)= ii
687 
688  jstack= jstack + 2
689 
690  if (jstack.gt.nstack) then
691  write (*,*) 'NSTACK overflow'
692  stop
693  endif
694 
695  if (ir-i+1.ge.j-1) then
696  istack(jstack )= ir
697  istack(jstack-1)= i
698  ir= j-1
699  else
700  istack(jstack )= j-1
701  istack(jstack-1)= l
702  l= i
703  endif
704 
705  endif
706 
707  goto 1
708 
709  end subroutine sainv_sort_nn
710 
711  subroutine hecmw_precond_sainv_nn_clear()
712  implicit none
713 
714  if (associated(sainvd)) deallocate(sainvd)
715  if (associated(sainvl)) deallocate(sainvl)
716  if (associated(sainvu)) deallocate(sainvu)
717  if (associated(inumfi1l)) deallocate(inumfi1l)
718  if (associated(inumfi1u)) deallocate(inumfi1u)
719  if (associated(fi1l)) deallocate(fi1l)
720  if (associated(fi1u)) deallocate(fi1u)
721  nullify(inumfi1l)
722  nullify(inumfi1u)
723  nullify(fi1l)
724  nullify(fi1u)
725  nullify(d)
726  nullify(al)
727  nullify(au)
728  nullify(indexl)
729  nullify(indexu)
730  nullify(iteml)
731  nullify(itemu)
732 
733  end subroutine hecmw_precond_sainv_nn_clear
734 end module hecmw_precond_sainv_nn
hecmw_precond_sainv_nn::hecmw_precond_sainv_nn_setup
subroutine, public hecmw_precond_sainv_nn_setup(hecMAT)
Definition: hecmw_precond_SAINV_nn.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_sainv_nn::hecmw_precond_sainv_nn_apply
subroutine, public hecmw_precond_sainv_nn_apply(R, ZP)
Definition: hecmw_precond_SAINV_nn.f90:88
hecmw_precond_sainv_nn
Definition: hecmw_precond_SAINV_nn.f90:6
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_matrix_misc
Definition: hecmw_matrix_misc.f90:6
hecmw_precond_sainv_nn::hecmw_precond_sainv_nn_clear
subroutine, public hecmw_precond_sainv_nn_clear()
Definition: hecmw_precond_SAINV_nn.f90:712
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444