18 integer(4),
parameter :: krealp = 8
20 integer(kind=kint) :: NPFIU, NPFIL
21 integer(kind=kint) :: N
22 integer(kind=kint) :: NDOF, NDOF2
23 integer(kind=kint),
pointer :: inumFI1L(:) => null()
24 integer(kind=kint),
pointer :: inumFI1U(:) => null()
25 integer(kind=kint),
pointer :: FI1L(:) => null()
26 integer(kind=kint),
pointer :: FI1U(:) => null()
28 integer(kind=kint),
pointer :: indexL(:) => null()
29 integer(kind=kint),
pointer :: indexU(:) => null()
30 integer(kind=kint),
pointer :: itemL(:) => null()
31 integer(kind=kint),
pointer :: itemU(:) => null()
32 real(kind=
kreal),
pointer :: d(:) => null()
33 real(kind=
kreal),
pointer :: al(:) => null()
34 real(kind=
kreal),
pointer :: au(:) => null()
36 real(kind=krealp),
pointer :: sainvu(:) => null()
37 real(kind=krealp),
pointer :: sainvl(:) => null()
38 real(kind=krealp),
pointer :: sainvd(:) => null()
39 real(kind=
kreal),
pointer :: t(:) => null()
50 integer(kind=kint ) :: precond
51 real(kind=krealp) :: filter
61 indexl => hecmat%indexL
62 indexu => hecmat%indexU
66 if (precond.eq.20)
call form_ilu1_sainv_nn(hecmat)
68 allocate (sainvd(ndof2*hecmat%NP))
69 allocate (sainvl(ndof2*npfiu))
70 allocate (t(ndof*hecmat%NP))
75 filter= hecmat%Rarray(5)
77 write(*,
"(a,F15.8)")
"### SAINV FILTER :",filter
79 call hecmw_sainv_nn(hecmat)
81 allocate (sainvu(ndof2*npfiu))
84 call hecmw_sainv_make_u_nn(hecmat)
90 real(kind=
kreal),
intent(inout) :: zp(:)
91 real(kind=
kreal),
intent(in) :: r(:)
92 integer(kind=kint) :: in, i, j, isl, iel, isu, ieu,idof,jdof
93 real(kind=
kreal) :: sw(ndof),x(ndof)
105 x(idof) = r(ndof*(in-1)+idof)
109 sw(idof) = sw(idof) + sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
114 x(idof) = r(ndof*(i-1)+idof)
115 t(ndof*(i-1)+idof)=x(idof)+sw(idof)
119 t(ndof*(i-1)+idof)=t(ndof*(i-1)+idof)+sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)*x(jdof)
121 t(ndof*(i-1)+idof)=t(ndof*(i-1)+idof)*sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
131 isu= inumfi1u(i-1) + 1
136 x(idof) = t(ndof*(in-1)+idof)
140 sw(idof) = sw(idof) + sainvu(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
146 x(idof) = t(ndof*(i-1)+idof)
149 zp(ndof*(i-1)+idof) = x(idof) + sw(idof)
150 do jdof = ndof, idof+1, -1
151 zp(ndof*(i-1)+idof) = zp(ndof*(i-1)+idof)+sainvd(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(jdof)
162 subroutine hecmw_sainv_nn(hecMAT)
166 integer(kind=kint) :: i, j, k, js, je, in, itr, idof, jdof, iitr
167 real(kind=krealp) :: dd, dtmp(hecmat%NDOF), x(hecmat%NDOF)
169 real(kind=krealp) :: filter
170 real(kind=krealp),
allocatable :: zz(:), vv(:)
172 filter= hecmat%Rarray(5)
174 allocate (vv(ndof*hecmat%NP))
175 allocate (zz(ndof*hecmat%NP))
182 zz(ndof*(itr-1)+idof)= sainvd(ndof2*(itr-1)+ndof*(idof-1)+iitr)
184 zz(ndof*(itr-1)+iitr)= 1.0d0
186 js= inumfi1l(itr-1) + 1
191 zz(ndof*(in-1)+idof)=sainvl(ndof2*(j-1)+ndof*(iitr-1)+idof)
197 x(idof)=zz(ndof*(i-1)+idof)
201 vv(ndof*(i-1)+idof) = vv(ndof*(i-1)+idof) + d(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(jdof)
211 vv(ndof*(in-1)+idof) = vv(ndof*(in-1)+idof) + al(ndof2*(j-1)+ndof*(jdof-1)+idof)*x(jdof)
221 vv(ndof*(in-1)+idof) = vv(ndof*(in-1)+idof) + au(ndof2*(j-1)+ndof*(jdof-1)+idof)*x(jdof)
229 dtmp(idof)= sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof)
234 sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = vv(ndof*(i-1)+idof)
236 sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) + &
237 & sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)*vv(ndof*(i-1)+jdof)
240 js= inumfi1l(i-1) + 1
245 x(idof)=vv(ndof*(in-1)+idof)
249 sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) + &
250 & sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
257 dd = 1.0d0/sainvd(ndof2*(itr-1)+ndof*(iitr-1)+iitr)
259 sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof) = dtmp(idof)
262 do idof = iitr+1, ndof
263 sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof)*dd
268 sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)*dd
274 dd = sainvd(ndof2*(itr-1)+ndof*(k-1)+k)
275 if(abs(dd) > filter)
then
277 sainvd(ndof2*(itr-1)+ndof*(jdof-1)+k)= sainvd(ndof2*(itr-1)+ndof*(jdof-1)+k) - dd*zz(ndof*(itr-1)+jdof)
279 js= inumfi1l(itr-1) + 1
284 sainvl(ndof2*(j-1)+ndof*(k-1)+idof) = sainvl(ndof2*(j-1)+ndof*(k-1)+idof)-dd*zz(ndof*(in-1)+idof)
291 js= inumfi1l(i-1) + 1
294 dd = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
295 if(abs(dd) > filter)
then
300 sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)=sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)-dd*zz(ndof*(in-1)+jdof)
313 sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)=1/sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
314 do jdof = idof+1, ndof
315 sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)=sainvd(ndof2*(i-1)+ndof*(idof-1)+jdof)
319 end subroutine hecmw_sainv_nn
321 subroutine hecmw_sainv_make_u_nn(hecMAT)
324 integer(kind=kint) i,j,k,n,m,o,idof,jdof
325 integer(kind=kint) is,ie,js,je
340 sainvu(ndof2*(n-1)+ndof*(jdof-1)+idof)=sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)
349 end subroutine hecmw_sainv_make_u_nn
354 subroutine form_ilu0_sainv_nn(hecMAT)
358 allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
359 allocate (fi1l(hecmat%NPL), fi1u(hecmat%NPU))
366 inumfi1l = hecmat%indexL
367 inumfi1u = hecmat%indexU
374 end subroutine form_ilu0_sainv_nn
379 subroutine form_ilu1_sainv_nn(hecMAT)
383 integer(kind=kint),
allocatable :: iwsl(:), iwsu(:), iw1(:), iw2(:)
384 integer(kind=kint) :: nplf1,npuf1
385 integer(kind=kint) :: i,jj,kk,l,isk,iek,isj,iej
386 integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
387 integer(kind=kint) :: j,k,isl,isu
396 allocate (iw1(hecmat%NP) , iw2(hecmat%NP))
397 allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
408 do l= indexl(i-1)+1, indexl(i)
411 do l= indexu(i-1)+1, indexu(i)
419 isj= indexu(kk-1) + 1
423 if (iw1(jj).eq.0 .and. jj.lt.i)
then
424 inumfi1l(i)= inumfi1l(i)+1
427 if (iw1(jj).eq.0 .and. jj.gt.i)
then
428 inumfi1u(i)= inumfi1u(i)+1
433 nplf1= nplf1 + inumfi1l(i)
434 npuf1= npuf1 + inumfi1u(i)
439 allocate (iwsl(0:hecmat%NP), iwsu(0:hecmat%NP))
440 allocate (fi1l(hecmat%NPL+nplf1), fi1u(hecmat%NPU+npuf1))
442 npfiu = hecmat%NPU+npuf1
443 npfil = hecmat%NPL+nplf1
451 iwsl(i)= indexl(i)-indexl(i-1) + inumfi1l(i) + iwsl(i-1)
452 iwsu(i)= indexu(i)-indexu(i-1) + inumfi1u(i) + iwsu(i-1)
458 inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
459 inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
463 do l= indexl(i-1)+1, indexl(i)
466 do l= indexu(i-1)+1, indexu(i)
474 isj= indexu(kk-1) + 1
478 if (iw1(jj).eq.0 .and. jj.lt.i)
then
480 fi1l(icoul+iwsl(i-1)+indexl(i)-indexl(i-1))= jj
483 if (iw1(jj).eq.0 .and. jj.gt.i)
then
485 fi1u(icouu+iwsu(i-1)+indexu(i)-indexu(i-1))= jj
495 icoul1= indexl(i) - indexl(i-1)
496 icoul2= inumfi1l(i) - inumfi1l(i-1)
497 icoul3= icoul1 + icoul2
498 icouu1= indexu(i) - indexu(i-1)
499 icouu2= inumfi1u(i) - inumfi1u(i-1)
500 icouu3= icouu1 + icouu2
504 do k= indexl(i-1)+1, indexl(i)
509 do k= inumfi1l(i-1)+1, inumfi1l(i)
511 iw1(icou0)= fi1l(icou0+iwsl(i-1))
517 call sainv_sort_nn (iw1, iw2, icoul3, hecmat%NP)
525 do k= indexu(i-1)+1, indexu(i)
530 do k= inumfi1u(i-1)+1, inumfi1u(i)
532 iw1(icou0)= fi1u(icou0+iwsu(i-1))
538 call sainv_sort_nn (iw1, iw2, icouu3, hecmat%NP)
554 deallocate (iw1, iw2)
555 deallocate (iwsl, iwsu)
557 end subroutine form_ilu1_sainv_nn
564 subroutine sainv_sort_nn(STEM, INUM, N, NP)
567 integer(kind=kint) :: n, np
568 integer(kind=kint),
dimension(NP) :: stem
569 integer(kind=kint),
dimension(NP) :: inum
570 integer(kind=kint),
dimension(:),
allocatable :: istack
571 integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
573 allocate (istack(-np:+np))
592 if (stem(i).le.ss)
goto 2
603 if (jstack.eq.0)
then
622 if (stem(l+1).gt.stem(ir))
then
631 if (stem(l).gt.stem(ir))
then
640 if (stem(l+1).gt.stem(l))
then
657 if (stem(i).lt.ss)
goto 3
661 if (stem(j).gt.ss)
goto 4
684 if (jstack.gt.nstack)
then
685 write (*,*)
'NSTACK overflow'
689 if (ir-i+1.ge.j-1)
then
703 end subroutine sainv_sort_nn
708 if (
associated(sainvd))
deallocate(sainvd)
709 if (
associated(sainvl))
deallocate(sainvl)
710 if (
associated(sainvu))
deallocate(sainvu)
711 if (
associated(inumfi1l))
deallocate(inumfi1l)
712 if (
associated(inumfi1u))
deallocate(inumfi1u)
713 if (
associated(fi1l))
deallocate(fi1l)
714 if (
associated(fi1u))
deallocate(fi1u)
integer(kind=kint) function, public hecmw_mat_get_precond(hecMAT)
subroutine, public hecmw_precond_nn_sainv_clear()
subroutine, public hecmw_precond_nn_sainv_setup(hecMAT)
subroutine, public hecmw_precond_nn_sainv_apply(R, ZP)
integer(kind=4), parameter kreal