hecmw_precond_RIF_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_rif_nn
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
 
  private
 
  public:: hecmw_precond_rif_nn_setup
  public:: hecmw_precond_rif_nn_apply
  public:: hecmw_precond_rif_nn_clear
 
  integer(4),parameter :: krealp = 8
 
  integer(kind=kint) :: NPFIU, NPFIL
  integer(kind=kint) :: N
  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 :: sainvl(:) => null()
  real(kind=krealp), pointer :: sainvd(:) => null()
  real(kind=krealp), pointer :: rifu(:)   => null()
  real(kind=krealp), pointer :: rifl(:)   => null()
  real(kind=krealp), pointer :: rifd(:)   => null()
 
contains
 
  !C***
  !C*** hecmw_precond_nn_sainv_setup
  !C***
  subroutine hecmw_precond_rif_nn_setup(hecMAT)
    implicit none
    type(hecmwst_matrix) :: hecmat
 
    integer(kind=kint ) :: precond, ndof, ndof2
 
    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.21) call form_ilu1_rif_nn(hecmat)
 
    allocate (sainvd(ndof2*hecmat%NP))
    allocate (sainvl(ndof2*npfiu))
    allocate (rifd(ndof2*hecmat%NP))
    allocate (rifu(ndof2*npfiu))
    sainvd  = 0.0d0
    sainvl  = 0.0d0
    rifd  = 0.0d0
    rifu  = 0.0d0
 
    filter= hecmat%Rarray(5)
 
    write(*,"(a,F15.8)")"### RIF FILTER   :",filter
 
    call hecmw_rif_nn(hecmat)
 
    allocate (rifl(ndof2*npfiu))
    rifl  = 0.0d0
 
    call hecmw_rif_make_u_nn(hecmat)
 
  end subroutine hecmw_precond_rif_nn_setup
 
  subroutine hecmw_precond_rif_nn_apply(ZP, NDOF)
    implicit none
    real(kind=kreal), intent(inout)  :: zp(:)
    integer(kind=kint) :: in, i, j, isl, iel,ndof,ndof2,idof,jdof
    real(kind=kreal) :: sw(ndof),x(ndof)
    ndof2=ndof*ndof
    !C-- FORWARD
    do i= 1, n
      do idof = 1, ndof
        sw(idof) = zp(ndof*(i-1)+idof)
      end do
      isl= inumfi1l(i-1)+1
      iel= inumfi1l(i)
      do j= isl, iel
        in= fi1l(j)
        do idof = 1, ndof
          x(idof) = zp(ndof*(in-1)+idof)
        end do
        do idof = 1, ndof
          do jdof = 1, ndof
            sw(idof) = sw(idof) - rifl(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
          end do
        end do
      enddo
      x = sw
      do idof = 2,ndof
        do jdof = 1, idof-1
          x(idof) = x(idof) - rifd(ndof2*(i-1)+ndof*(idof-1)+jdof )*x(jdof)
        end do
      end do
      zp(ndof*(i-1)+1:ndof*(i-1)+ndof) = x(1:ndof)
    enddo
 
    do i=1, n
      do idof=1,ndof
        zp(ndof*(i-1)+idof)= zp(ndof*(i-1)+idof)*rifd(ndof2*(i-1)+(idof-1)*ndof+idof)
      end do
    enddo
    !C-- BACKWARD
    do i= n, 1, -1
      do idof = 1, ndof
        sw(idof) = zp(ndof*(i-1)+idof)
      end do
      isl= inumfi1u(i-1) + 1
      iel= inumfi1u(i)
      do j= iel, isl, -1
        in= fi1u(j)
        do idof = 1, ndof
          x(idof) = zp(ndof*(in-1)+idof)
        end do
        do idof = 1, ndof
          do jdof = 1, ndof
            sw(idof) = sw(idof) - rifu(ndof2*(j-1)+ndof*(idof-1)+jdof)*x(jdof)
          end do
        end do
      enddo
      x=sw
      do idof = ndof, 1, -1
        do jdof = ndof, idof+1, -1
          x(idof) = x(idof) - rifd(ndof*ndof*(i-1)+ndof*(jdof-1)+idof)*x(jdof)
        end do
      end do
      do idof = 1, ndof
        zp(ndof*(i-1)+idof) = x(idof)
      end do
    enddo
  end subroutine hecmw_precond_rif_nn_apply
 
 
  !C***
  !C*** hecmw_rif_nn
  !C***
  subroutine hecmw_rif_nn(hecMAT)
    implicit none
    type (hecmwst_matrix)     :: hecmat
    integer(kind=kint) :: i, j, k, js, je, in, itr, np, ndof, ndof2, 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)
 
    np = hecmat%NP
    ndof = hecmat%NDOF
    ndof2= ndof*ndof
    allocate(vv(ndof*np))
    allocate(zz(ndof*np))
 
    do itr=1,np
      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,np
          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,np
          do idof = 1, ndof
            sainvd(ndof2*(i-1)+ndof*(idof-1)+idof) = sainvd(ndof2*(i-1)+ndof*(idof-1)+idof)*dd
          end do
        enddo
 
        do idof = iitr, ndof
          rifd(ndof2*(itr-1)+ndof*(idof-1)+iitr) = sainvd(ndof2*(itr-1)+ndof*(idof-1)+idof) !RIF UPDATE
        end do
 
        js= inumfi1u(itr-1) + 1
        je= inumfi1u(itr  )
        do j= js, je
          do idof = 1, ndof
            rifu(ndof2*(j-1)+ndof*(iitr-1)+idof) = sainvd(ndof2*(fi1u(j)-1)+ndof*(idof-1)+idof)
          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,np
          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)
 
  end subroutine hecmw_rif_nn
 
  subroutine hecmw_rif_make_u_nn(hecMAT)
    implicit none
    type (hecmwst_matrix)     :: hecmat
    integer(kind=kint) i,j,k,n,m,o,idof,jdof,ndof,ndof2
    integer(kind=kint) is,ie,js,je
    ndof=hecmat%NDOF
    ndof2=ndof*ndof
    n = 1
    do i= 1, hecmat%NP
      is=inumfi1l(i-1) + 1
      ie=inumfi1l(i  )
      flag1:do k= is, ie
        m = fi1l(k)
        js=inumfi1u(m-1) + 1
        je=inumfi1u(m  )
        do j= js,je
          o = fi1u(j)
          if (o == i)then
            do idof = 1, ndof
              do jdof = 1, ndof
                rifl(ndof2*(n-1)+ndof*(idof-1)+jdof)=rifu(ndof2*(j-1)+ndof*(jdof-1)+idof)
              end do
            end do
            n = n + 1
            cycle flag1
          endif
        enddo
      enddo flag1
    enddo
  end subroutine hecmw_rif_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_rif_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 rif_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 rif_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_rif_nn
 
  !C
  !C***
  !C*** fill_in_S33_SORT
  !C***
  !C
  subroutine rif_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 rif_sort_nn
 
  subroutine hecmw_precond_rif_nn_clear()
    implicit none
 
    if (associated(sainvd)) deallocate(sainvd)
    if (associated(sainvl)) deallocate(sainvl)
    if (associated(rifd)) deallocate(rifd)
    if (associated(rifu)) deallocate(rifu)
    if (associated(rifl)) deallocate(rifl)
    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_rif_nn_clear
end module hecmw_precond_rif_nn