hecmw_jadm_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_jad_type_nn
  use hecmw_util
  use m_hecmw_comm_f
  implicit none
 
  private
 
  public :: hecmw_jad_init_nn
  public :: hecmw_jad_finalize_nn
  public :: hecmw_jad_matvec_nn
 
  !C---------------------- AU&AL
  real(kind=kreal), allocatable      :: ajad(:)
  integer(kind=kint), allocatable    :: JAJAD(:)
  integer(kind=kint), allocatable    :: JADORD(:)
  integer(kind=kint), allocatable    :: IAJAD(:)
  integer(kind=kint) :: MJAD
  real(kind=kreal), allocatable  :: wp(:,:)
  integer(kind=kint) :: INITIALIZED = 0
 
contains
 
  subroutine hecmw_jad_init_nn(hecMAT)
    type(hecmwst_matrix) :: hecmat
    allocate(wp(hecmat%NDOF,hecmat%NP))
    allocate(ajad((hecmat%NPL+hecmat%NPU)*hecmat%NDOF*hecmat%NDOF))
    allocate(jajad(hecmat%NPL+hecmat%NPU))
    allocate(jadord(hecmat%NP))
    allocate(iajad(hecmat%NP+1))
    call repack(hecmat%N, hecmat, mjad, ajad, jajad, iajad, jadord)
    initialized = 1
  end subroutine hecmw_jad_init_nn
 
  subroutine hecmw_jad_finalize_nn()
    deallocate(ajad)
    deallocate(jajad)
    deallocate(jadord)
    deallocate(iajad)
    deallocate(wp)
    initialized = 0
  end subroutine hecmw_jad_finalize_nn
 
  subroutine hecmw_jad_matvec_nn(hecMESH, hecMAT, X, Y, COMMtime)
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(hecmwst_matrix), intent(in), target :: hecmat
    real(kind=kreal), intent(in) :: x(:)
    real(kind=kreal), intent(out) :: y(:)
    real(kind=kreal), intent(inout) :: commtime
    real(kind=kreal) :: start_time, end_time
    real(kind=kreal), pointer :: d(:)
    integer(kind=kint) :: i,idof,jdof,ndof,ndof2
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmat%NP,hecmat%NDOF)
    end_time= hecmw_wtime()
    commtime = commtime + end_time - start_time
 
    d => hecmat%D
    ndof = hecmat%NDOF
    ndof2 = ndof*ndof
    y=0.0d0
 
    !$OMP PARALLEL PRIVATE(i, idof, jdof)
    !$OMP DO
    do i= 1, hecmat%N
      do idof=1,hecmat%NDOF
        do jdof=1,hecmat%NDOF
          y(ndof*(i-1)+idof) = y(ndof*(i-1)+idof) + d(ndof2*(i-1)+ndof*(idof-1)+jdof)*x(ndof*(i-1)+jdof)
        end do
      end do
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
    call matjad(hecmat%N,hecmat%NDOF, mjad, iajad, jajad, ajad, jadord, x, y, wp)
  end subroutine hecmw_jad_matvec_nn
 
  subroutine repack(N, hecMAT, MJAD,  AJAD, JAJAD, IAJAD, JADORD)
    use hecmw_util
    !C---------------------------------
    type (hecmwst_matrix)     :: hecmat
    !C----------------------
    integer(kind = kint) :: n, mjad
    real(kind = kreal), dimension(*) :: ajad
    integer(kind = kint), dimension(*) :: jajad
    integer(kind = kint), dimension(*) :: iajad
    integer(kind = kint), dimension(*) :: jadord
 
    integer(kind = kint) :: ijad, maxnz, minnz,ndof,ndof2
    integer(kind = kint) :: i, j, js, je, in, jc
    integer(kind = kint), allocatable :: len(:), lenz(:), jadreord(:)
    ndof = hecmat%NDOF;ndof2=ndof*ndof
    allocate(len(n))
    allocate(jadreord(n))
    do i=1,n
      len(i)= hecmat%indexL(i) - hecmat%indexL(i-1) &
        &        + hecmat%indexU(i) - hecmat%indexU(i-1)
    end do
    maxnz=maxval(len(1:n))
    minnz=minval(len(1:n))
    mjad =maxnz
    allocate(lenz(0:mjad))
    lenz = 0
    do i=1,n
      lenz(len(i))=lenz(len(i))+1
    enddo
    do i=maxnz-1,minnz,-1
      lenz(i)=lenz(i)+lenz(i+1)
    enddo
    do i=1,n
      jadord(i)=lenz(len(i))
      lenz(len(i))=lenz(len(i))-1
    enddo
    do i=1,n
      jadreord(jadord(i))=i
    enddo
    do i=1,n
      lenz(len(jadreord(i)))=i
    enddo
    do i=maxnz-1,1,-1
      lenz(i)=max(lenz(i+1),lenz(i))
    enddo
    iajad(1)=1
    do i=1,maxnz
      iajad(i+1)=iajad(i)+lenz(i)
    enddo
    len=0
    do i= 1, n
      ijad=jadord(i)
      js= hecmat%indexL(i-1) + 1
      je= hecmat%indexL(i  )
      do j=js,je
        in  = hecmat%itemL(j)
        len(ijad)=len(ijad)+1
        jc=iajad(len(ijad))+ijad-1
        ajad(ndof2*(jc-1)+1:ndof2*(jc)) = hecmat%AL(ndof2*(j-1)+1:ndof2*(j))
        jajad(jc) = in
      end do
    end do
    do i= 1, n
      ijad=jadord(i)
      js= hecmat%indexU(i-1) + 1
      je= hecmat%indexU(i  )
      do j=js,je
        in  = hecmat%itemU(j)
        len(ijad)=len(ijad)+1
        jc=iajad(len(ijad))+ijad-1
        ajad(ndof2*(jc-1)+1:ndof2*(jc)) = hecmat%AU(ndof2*(j-1)+1:ndof2*(j))
        jajad(jc) = in
      end do
    end do
    deallocate(len)
    deallocate(jadreord)
    deallocate(lenz)
  end subroutine repack
 
  subroutine matjad(N, NDOF,  MJAD, IAJAD, JAJAD, AJAD, JADORD, X, Y, W)
    use hecmw_util
    integer(kind=kint) :: N,NDOF, MJAD,NDOF2
    integer(kind=kint) :: IAJAD(*), JAJAD(*), JADORD(*)
    real(kind=kreal)   :: ajad(*), x(*), y(*),  w(ndof,n)
    integer(kind=kint) :: I, K, NZ, IXX, idof,jdof
    ndof2=ndof*ndof
 
    w=0.0d0
 
    do nz=1,mjad
      !$OMP PARALLEL PRIVATE(K,IXX,idof,jdof)
      !$OMP DO
      do k=iajad(nz),iajad(nz+1)-1
        ixx = k-iajad(nz)+1
        do idof = 1, ndof
          do jdof = 1, ndof
            w(idof,ixx)=w(idof,ixx)+ajad(ndof2*(k-1)+ndof*(idof-1)+jdof)*x(ndof*(jajad(k)-1)+jdof)
          end do
        end do
      enddo
      !$OMP END DO
      !$OMP END PARALLEL
    enddo
 
    !$OMP PARALLEL PRIVATE(I,idof)
    !$OMP DO
    do i=1,n
      do idof = 1, ndof
        y(ndof*(i-1)+idof)=y(ndof*(i-1)+idof)+w(idof,jadord(i))
      end do
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
  end subroutine matjad
 
end module hecmw_jad_type_nn