hecmw_precond_SAINV_nn.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module hecmw_precond_sainv_nn
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
 
  private
 
  public:: hecmw_precond_sainv_nn_setup
  public:: hecmw_precond_sainv_nn_apply
  public:: hecmw_precond_sainv_nn_clear
 
  integer(4),parameter :: krealp = 8
 
  integer(kind=kint) :: NPFIU, NPFIL
  integer(kind=kint) :: N
  integer(kind=kint) :: NDOF, NDOF2
  integer(kind=kint), pointer :: inumFI1L(:) => null()
  integer(kind=kint), pointer :: inumFI1U(:) => null()
  integer(kind=kint), pointer :: FI1L(:) => null()
  integer(kind=kint), pointer :: FI1U(:) => null()
 
  integer(kind=kint), pointer :: indexL(:) => null()
  integer(kind=kint), pointer :: indexU(:) => null()
  integer(kind=kint), pointer :: itemL(:)  => null()
  integer(kind=kint), pointer :: itemU(:)  => null()
  real(kind=kreal), pointer :: d(:)  => null()
  real(kind=kreal), pointer :: al(:) => null()
  real(kind=kreal), pointer :: au(:) => null()
 
  real(kind=krealp), pointer :: sainvu(:) => null()
  real(kind=krealp), pointer :: sainvl(:) => null()
  real(kind=krealp), pointer :: sainvd(:) => null()
  real(kind=kreal),  pointer :: t(:) => null()
 
contains
 
  !C***
  !C*** hecmw_precond_nn_sainv_setup
  !C***
  subroutine hecmw_precond_sainv_nn_setup(hecMAT)
    implicit none
    type(hecmwst_matrix) :: hecmat
 
    integer(kind=kint ) :: precond
    real(kind=krealp) :: filter
 
    n = hecmat%N
    ndof = hecmat%NDOF
    ndof2 = ndof*ndof
    precond = hecmw_mat_get_precond(hecmat)
 
    d => hecmat%D
    au=> hecmat%AU
    al=> hecmat%AL
    indexl => hecmat%indexL
    indexu => hecmat%indexU
    iteml => hecmat%itemL
    itemu => hecmat%itemU
 
    if (precond.eq.20) call form_ilu1_sainv_nn(hecmat)
 
    allocate (sainvd(ndof2*hecmat%NP))
    allocate (sainvl(ndof2*npfiu))
    allocate (t(ndof*hecmat%NP))
    sainvd  = 0.0d0
    sainvl  = 0.0d0
    t = 0.0d0
 
    filter= hecmat%Rarray(5)
 
    write(*,"(a,F15.8)")"### SAINV FILTER   :",filter
 
    call hecmw_sainv_nn(hecmat)
 
    allocate (sainvu(ndof2*npfiu))
    sainvu  = 0.0d0
 
    call hecmw_sainv_make_u_nn(hecmat)
 
  end subroutine hecmw_precond_sainv_nn_setup
 
  subroutine hecmw_precond_sainv_nn_apply(R, ZP)
    implicit none
    real(kind=kreal), intent(inout)  :: zp(:)
    real(kind=kreal), intent(in)  :: r(:)
    integer(kind=kint) :: in, i, j, isl, iel, isu, ieu,idof,jdof
    real(kind=kreal) :: sw(ndof),x(ndof)
 
#ifndef _OPENACC
    !$OMP PARALLEL DEFAULT(NONE) &
      !$OMP&PRIVATE(i,X,SW,j,in,isL,ieL,isU,ieU,idof,jdof) &
      !$OMP&SHARED(N,SAINVD,SAINVL,SAINVU,inumFI1U,FI1U,inumFI1L,FI1L,R,T,ZP,NDOF,NDOF2)
#endif
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
#else
    !$OMP DO
#endif
    !C-- FORWARD
    do i= 1, n
      do idof = 1, ndof
        sw(idof) = 0.0d0
      end do
      isl= inumfi1l(i-1)+1
      iel= inumfi1l(i)
      do j= isl, iel
        in= fi1l(j)
        do idof = 1, ndof
          x(idof) = r(ndof*(in-1)+idof)
        end do
        do idof = 1, ndof
          do jdof = 1, ndof
            sw(idof) = sw(idof) + sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
          end do
        end do
      enddo
      do idof = 1, ndof
        x(idof) = r(ndof*(i-1)+idof)
        t(ndof*(i-1)+idof)=x(idof)+sw(idof)
      end do
      do idof = 1, ndof
        do jdof = 1, idof-1
          t(ndof*(i-1)+idof)=t(ndof*(i-1)+idof)+sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)*x(jdof)
        end do
        t(ndof*(i-1)+idof)=t(ndof*(i-1)+idof)*sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
      end do
    enddo
#ifdef _OPENACC
    !$acc end kernels
#else
    !$OMP END DO
#endif
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
#else
    !$OMP DO
#endif
    !C-- BACKWARD
    do i= 1, n
      do idof = 1, ndof
        sw(idof) = 0.0d0
      end do
 
      isu= inumfi1u(i-1) + 1
      ieu= inumfi1u(i)
      do j= isu, ieu
        in= fi1u(j)
        do idof = 1, ndof
          x(idof) = t(ndof*(in-1)+idof)
        end do
        do idof = 1, ndof
          do jdof = 1, ndof
            sw(idof) = sw(idof) + sainvu(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
          end do
        end do
      enddo
 
      do idof = 1, ndof
        x(idof) = t(ndof*(i-1)+idof)
      end do
      do idof = 1, ndof
        zp(ndof*(i-1)+idof) = x(idof) + sw(idof)
        do jdof = ndof, idof+1, -1
          zp(ndof*(i-1)+idof) = zp(ndof*(i-1)+idof)+sainvd(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(jdof)
        end do
      end do
    enddo
#ifdef _OPENACC
    !$acc end kernels
#else
    !$OMP END DO
#endif
 
#ifndef _OPENACC
    !$OMP END PARALLEL
#endif
 
  end subroutine hecmw_precond_sainv_nn_apply
 
 
  !C***
  !C*** hecmw_rif_nn
  !C***
  subroutine hecmw_sainv_nn(hecMAT)
    implicit none
    type (hecmwst_matrix)     :: hecmat
 
    integer(kind=kint) :: i, j, k, js, je, in, itr, idof, jdof, iitr
    real(kind=krealp) :: dd, dtmp(hecmat%NDOF), x(hecmat%NDOF)
 
    real(kind=krealp) :: filter
    real(kind=krealp), allocatable :: zz(:), vv(:)
 
    filter= hecmat%Rarray(5)
 
    allocate (vv(ndof*hecmat%NP))
    allocate (zz(ndof*hecmat%NP))
    do itr=1,n
      do iitr=1,ndof
        zz(:) = 0.0d0
        vv(:) = 0.0d0
        !{v}=[A]{zi}
        do idof = 1,ndof
          zz(ndof*(itr-1)+idof)= sainvd(ndof2*(itr-1)+ndof*(idof-1)+iitr)
        end do
        zz(ndof*(itr-1)+iitr)= 1.0d0
 
        js= inumfi1l(itr-1) + 1
        je= inumfi1l(itr  )
        do j= js, je
          in  = fi1l(j)
          do idof = 1, ndof
            zz(ndof*(in-1)+idof)=sainvl(ndof2*(j-1)+ndof*(iitr-1)+idof)
          end do
        enddo
 
        do i= 1, itr
          do idof = 1,ndof
            x(idof)=zz(ndof*(i-1)+idof)
          end do
          do idof = 1, ndof
            do jdof = 1, ndof
              vv(ndof*(i-1)+idof) = vv(ndof*(i-1)+idof) + d(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(jdof)
            end do
          end do
 
          js= indexl(i-1) + 1
          je= indexl(i  )
          do j=js,je
            in = iteml(j)
            do idof = 1, ndof
              do jdof = 1, ndof
                vv(ndof*(in-1)+idof) = vv(ndof*(in-1)+idof) + al(ndof2*(j-1)+ndof*(jdof-1)+idof)*x(jdof)
              end do
            end do
          enddo
          js= indexu(i-1) + 1
          je= indexu(i  )
          do j= js, je
            in = itemu(j)
            do idof = 1, ndof
              do jdof = 1, ndof
                vv(ndof*(in-1)+idof) = vv(ndof*(in-1)+idof) + au(ndof2*(j-1)+ndof*(jdof-1)+idof)*x(jdof)
              end do
            end do
          enddo
        enddo
 
        !{d}={v^t}{z_j}
        do idof = 1, ndof
          dtmp(idof)= sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof)
        end do
 
        do i= itr,n
          do idof = 1,ndof
            sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = vv(ndof*(i-1)+idof)
            do jdof = 1, idof-1
              sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) + &
                & sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)*vv(ndof*(i-1)+jdof)
            end do
          end do
          js= inumfi1l(i-1) + 1
          je= inumfi1l(i  )
          do j= js, je
            in  = fi1l(j)
            do idof = 1,ndof
              x(idof)=vv(ndof*(in-1)+idof)
            end do
            do idof = 1, ndof
              do jdof = 1, ndof
                sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) + &
                  & sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
              end do
            end do
          enddo
        enddo
 
        !Update D
        dd = 1.0d0/sainvd(ndof2*(itr-1)+ndof*(iitr-1)+iitr)
        do idof=1,iitr-1
          sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof) = dtmp(idof)
        end do
        !        SAINVD(NDOF2*(itr-1)+NDOF*(iitr-1)+iitr)=dd
        do idof = iitr+1, ndof
          sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof)*dd
        end do
 
        do i =itr+1,n
          do idof = 1, ndof
            sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)*dd
          end do
        enddo
 
        !Update Z
        do k=iitr+1,ndof
          dd = sainvd(ndof2*(itr-1)+ndof*(k-1)+k)
          if(abs(dd) > filter)then
            do jdof = 1, iitr
              sainvd(ndof2*(itr-1)+ndof*(jdof-1)+k)= sainvd(ndof2*(itr-1)+ndof*(jdof-1)+k) - dd*zz(ndof*(itr-1)+jdof)
            end do
            js= inumfi1l(itr-1) + 1
            je= inumfi1l(itr  )
            do j= js, je
              in  = fi1l(j)
              do idof = 1, ndof
                sainvl(ndof2*(j-1)+ndof*(k-1)+idof) = sainvl(ndof2*(j-1)+ndof*(k-1)+idof)-dd*zz(ndof*(in-1)+idof)
              end do
            enddo
          endif
        end do
 
        do i= itr +1,n
          js= inumfi1l(i-1) + 1
          je= inumfi1l(i  )
          do idof = 1, ndof
            dd = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
            if(abs(dd) > filter)then
              do j= js, je
                in  = fi1l(j)
                if (in > itr) exit
                do jdof=1,ndof
                  sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)=sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)-dd*zz(ndof*(in-1)+jdof)
                end do
              enddo
            endif
          end do
        enddo
      end do
    enddo
    deallocate(vv)
    deallocate(zz)
 
    do i =1,n
      do idof = 1, ndof
        sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)=1/sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)
        do jdof = idof+1, ndof
          sainvd(ndof2*(i-1)+ndof*(jdof-1)+idof)=sainvd(ndof2*(i-1)+ndof*(idof-1)+jdof)
        end do
      end do
    enddo
  end subroutine hecmw_sainv_nn
 
  subroutine hecmw_sainv_make_u_nn(hecMAT)
    implicit none
    type (hecmwst_matrix)     :: hecmat
    integer(kind=kint) i,j,k,n,m,o,idof,jdof
    integer(kind=kint) is,ie,js,je
 
    n = 1
    do i= 1, hecmat%NP
      is=inumfi1u(i-1) + 1
      ie=inumfi1u(i  )
      flag1:do k= is, ie
        m = fi1u(k)
        js=inumfi1l(m-1) + 1
        je=inumfi1l(m  )
        do j= js,je
          o = fi1l(j)
          if (o == i)then
            do idof = 1, ndof
              do jdof = 1, ndof
                sainvu(ndof2*(n-1)+ndof*(jdof-1)+idof)=sainvl(ndof2*(j-1)+ndof*(idof-1)+jdof)
              end do
            end do
            n = n + 1
            cycle flag1
          endif
        enddo
      enddo flag1
    enddo
  end subroutine hecmw_sainv_make_u_nn
 
  !C***
  !C*** FORM_ILU1_nn
  !C*** form ILU(1) matrix
 
  !C***
  !C*** FORM_ILU1_nn
  !C*** form ILU(1) matrix
  subroutine form_ilu1_sainv_nn(hecMAT)
    implicit none
    type(hecmwst_matrix) :: hecmat
 
    integer(kind=kint),allocatable :: iwsl(:), iwsu(:), iw1(:), iw2(:)
    integer(kind=kint) :: nplf1,npuf1
    integer(kind=kint) :: i,jj,kk,l,isk,iek,isj,iej
    integer(kind=kint) :: icou,icou0,icouu,icouu1,icouu2,icouu3,icoul,icoul1,icoul2,icoul3
    integer(kind=kint) :: j,k,isl,isu
    !C
    !C +--------------+
    !C | find fill-in |
    !C +--------------+
    !C===
 
    !C
    !C-- count fill-in
    allocate (iw1(hecmat%NP) , iw2(hecmat%NP))
    allocate (inumfi1l(0:hecmat%NP), inumfi1u(0:hecmat%NP))
 
    inumfi1l= 0
    inumfi1u= 0
 
    nplf1= 0
    npuf1= 0
    do i= 2, hecmat%NP
      icou= 0
      iw1= 0
      iw1(i)= 1
      do l= indexl(i-1)+1, indexl(i)
        iw1(iteml(l))= 1
      enddo
      do l= indexu(i-1)+1, indexu(i)
        iw1(itemu(l))= 1
      enddo
 
      isk= indexl(i-1) + 1
      iek= indexl(i)
      do k= isk, iek
        kk= iteml(k)
        isj= indexu(kk-1) + 1
        iej= indexu(kk  )
        do j= isj, iej
          jj= itemu(j)
          if (iw1(jj).eq.0 .and. jj.lt.i) then
            inumfi1l(i)= inumfi1l(i)+1
            iw1(jj)= 1
          endif
          if (iw1(jj).eq.0 .and. jj.gt.i) then
            inumfi1u(i)= inumfi1u(i)+1
            iw1(jj)= 1
          endif
        enddo
      enddo
      nplf1= nplf1 + inumfi1l(i)
      npuf1= npuf1 + inumfi1u(i)
    enddo
 
    !C
    !C-- specify fill-in
    allocate (iwsl(0:hecmat%NP), iwsu(0:hecmat%NP))
    allocate (fi1l(hecmat%NPL+nplf1), fi1u(hecmat%NPU+npuf1))
 
    npfiu = hecmat%NPU+npuf1
    npfil = hecmat%NPL+nplf1
 
    fi1l= 0
    fi1u= 0
 
    iwsl= 0
    iwsu= 0
    do i= 1, hecmat%NP
      iwsl(i)= indexl(i)-indexl(i-1) + inumfi1l(i) + iwsl(i-1)
      iwsu(i)= indexu(i)-indexu(i-1) + inumfi1u(i) + iwsu(i-1)
    enddo
 
    do i= 2, hecmat%NP
      icoul= 0
      icouu= 0
      inumfi1l(i)= inumfi1l(i-1) + inumfi1l(i)
      inumfi1u(i)= inumfi1u(i-1) + inumfi1u(i)
      icou= 0
      iw1= 0
      iw1(i)= 1
      do l= indexl(i-1)+1, indexl(i)
        iw1(iteml(l))= 1
      enddo
      do l= indexu(i-1)+1, indexu(i)
        iw1(itemu(l))= 1
      enddo
 
      isk= indexl(i-1) + 1
      iek= indexl(i)
      do k= isk, iek
        kk= iteml(k)
        isj= indexu(kk-1) + 1
        iej= indexu(kk  )
        do j= isj, iej
          jj= itemu(j)
          if (iw1(jj).eq.0 .and. jj.lt.i) then
            icoul           = icoul + 1
            fi1l(icoul+iwsl(i-1)+indexl(i)-indexl(i-1))= jj
            iw1(jj)          = 1
          endif
          if (iw1(jj).eq.0 .and. jj.gt.i) then
            icouu           = icouu + 1
            fi1u(icouu+iwsu(i-1)+indexu(i)-indexu(i-1))= jj
            iw1(jj)          = 1
          endif
        enddo
      enddo
    enddo
 
    isl  = 0
    isu  = 0
    do i= 1, hecmat%NP
      icoul1= indexl(i) - indexl(i-1)
      icoul2= inumfi1l(i) - inumfi1l(i-1)
      icoul3= icoul1 + icoul2
      icouu1= indexu(i) - indexu(i-1)
      icouu2= inumfi1u(i) - inumfi1u(i-1)
      icouu3= icouu1 + icouu2
      !C
      !C-- LOWER part
      icou0= 0
      do k= indexl(i-1)+1, indexl(i)
        icou0 = icou0 + 1
        iw1(icou0)= iteml(k)
      enddo
 
      do k= inumfi1l(i-1)+1, inumfi1l(i)
        icou0 = icou0 + 1
        iw1(icou0)= fi1l(icou0+iwsl(i-1))
      enddo
 
      do k= 1, icoul3
        iw2(k)= k
      enddo
      call sainv_sort_nn (iw1, iw2, icoul3, hecmat%NP)
 
      do k= 1, icoul3
        fi1l(k+isl)= iw1(k)
      enddo
      !C
      !C-- UPPER part
      icou0= 0
      do k= indexu(i-1)+1, indexu(i)
        icou0 = icou0 + 1
        iw1(icou0)= itemu(k)
      enddo
 
      do k= inumfi1u(i-1)+1, inumfi1u(i)
        icou0 = icou0 + 1
        iw1(icou0)= fi1u(icou0+iwsu(i-1))
      enddo
 
      do k= 1, icouu3
        iw2(k)= k
      enddo
      call sainv_sort_nn (iw1, iw2, icouu3, hecmat%NP)
 
      do k= 1, icouu3
        fi1u(k+isu)= iw1(k)
      enddo
 
      isl= isl + icoul3
      isu= isu + icouu3
    enddo
 
    !C===
    do i= 1, hecmat%NP
      inumfi1l(i)= iwsl(i)
      inumfi1u(i)= iwsu(i)
    enddo
 
    deallocate (iw1, iw2)
    deallocate (iwsl, iwsu)
    !C===
  end subroutine form_ilu1_sainv_nn
 
  !C
  !C***
  !C*** fill_in_S33_SORT
  !C***
  !C
  subroutine sainv_sort_nn(STEM, INUM, N, NP)
    use hecmw_util
    implicit none
    integer(kind=kint) :: n, np
    integer(kind=kint), dimension(NP) :: stem
    integer(kind=kint), dimension(NP) :: inum
    integer(kind=kint), dimension(:), allocatable :: istack
    integer(kind=kint) :: m,nstack,jstack,l,ir,ip,i,j,k,ss,ii,temp,it
 
    allocate (istack(-np:+np))
 
    m     = 100
    nstack= np
 
    jstack= 0
    l     = 1
    ir    = n
 
    ip= 0
    1   continue
    ip= ip + 1
 
    if (ir-l.lt.m) then
      do j= l+1, ir
        ss= stem(j)
        ii= inum(j)
 
        do i= j-1,1,-1
          if (stem(i).le.ss) goto 2
          stem(i+1)= stem(i)
          inum(i+1)= inum(i)
        end do
        i= 0
 
        2       continue
        stem(i+1)= ss
        inum(i+1)= ii
      end do
 
      if (jstack.eq.0) then
        deallocate (istack)
        return
      endif
 
      ir = istack(jstack)
      l = istack(jstack-1)
      jstack= jstack - 2
    else
 
      k= (l+ir) / 2
      temp = stem(k)
      stem(k)  = stem(l+1)
      stem(l+1)= temp
 
      it = inum(k)
      inum(k)  = inum(l+1)
      inum(l+1)= it
 
      if (stem(l+1).gt.stem(ir)) then
        temp = stem(l+1)
        stem(l+1)= stem(ir)
        stem(ir )= temp
        it = inum(l+1)
        inum(l+1)= inum(ir)
        inum(ir )= it
      endif
 
      if (stem(l).gt.stem(ir)) then
        temp = stem(l)
        stem(l )= stem(ir)
        stem(ir)= temp
        it = inum(l)
        inum(l )= inum(ir)
        inum(ir)= it
      endif
 
      if (stem(l+1).gt.stem(l)) then
        temp = stem(l+1)
        stem(l+1)= stem(l)
        stem(l  )= temp
        it = inum(l+1)
        inum(l+1)= inum(l)
        inum(l  )= it
      endif
 
      i= l + 1
      j= ir
 
      ss= stem(l)
      ii= inum(l)
 
      3     continue
      i= i + 1
      if (stem(i).lt.ss) goto 3
 
      4     continue
      j= j - 1
      if (stem(j).gt.ss) goto 4
 
      if (j.lt.i)        goto 5
 
      temp   = stem(i)
      stem(i)= stem(j)
      stem(j)= temp
 
      it     = inum(i)
      inum(i)= inum(j)
      inum(j)= it
 
      goto 3
 
      5     continue
 
      stem(l)= stem(j)
      stem(j)= ss
      inum(l)= inum(j)
      inum(j)= ii
 
      jstack= jstack + 2
 
      if (jstack.gt.nstack) then
        write (*,*) 'NSTACK overflow'
        stop
      endif
 
      if (ir-i+1.ge.j-1) then
        istack(jstack  )= ir
        istack(jstack-1)= i
        ir= j-1
      else
        istack(jstack  )= j-1
        istack(jstack-1)= l
        l= i
      endif
 
    endif
 
    goto 1
 
  end subroutine sainv_sort_nn
 
  subroutine hecmw_precond_sainv_nn_clear()
    implicit none
 
    if (associated(sainvd)) deallocate(sainvd)
    if (associated(sainvl)) deallocate(sainvl)
    if (associated(sainvu)) deallocate(sainvu)
    if (associated(inumfi1l)) deallocate(inumfi1l)
    if (associated(inumfi1u)) deallocate(inumfi1u)
    if (associated(fi1l)) deallocate(fi1l)
    if (associated(fi1u)) deallocate(fi1u)
    nullify(inumfi1l)
    nullify(inumfi1u)
    nullify(fi1l)
    nullify(fi1u)
    nullify(d)
    nullify(al)
    nullify(au)
    nullify(indexl)
    nullify(indexu)
    nullify(iteml)
    nullify(itemu)
 
  end subroutine hecmw_precond_sainv_nn_clear
end module hecmw_precond_sainv_nn