hecmw_jadm_44.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_44
  use hecmw_util
  use m_hecmw_comm_f
  implicit none
 
  private
 
  public :: hecmw_jad_init_44
  public :: hecmw_jad_finalize_44
  public :: hecmw_jad_matvec_44
 
  !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  :: wp1(:), wp2(:), wp3(:), wp4(:)
  integer(kind=kint) :: INITIALIZED = 0
 
contains
 
  subroutine hecmw_jad_init_44(hecMAT)
    type(hecmwst_matrix) :: hecmat
    allocate(wp1(hecmat%NP), wp2(hecmat%NP), wp3(hecmat%NP), wp4(hecmat%NP))
    allocate(ajad((hecmat%NPL+hecmat%NPU)*16))
    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_44
 
  subroutine hecmw_jad_finalize_44()
    deallocate(ajad)
    deallocate(jajad)
    deallocate(jadord)
    deallocate(iajad)
    deallocate(wp1,wp2,wp3,wp4)
    initialized = 0
  end subroutine hecmw_jad_finalize_44
 
  subroutine hecmw_jad_matvec_44(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
    real(kind=kreal) :: x1, x2, x3, x4
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmat%NP, 4)
    end_time= hecmw_wtime()
    commtime = commtime + end_time - start_time
 
    d => hecmat%D
 
    !$OMP PARALLEL PRIVATE(i)
    !$OMP DO
    do i= 1, hecmat%N
      x1= x(4*i-3)
      x2= x(4*i-2)
      x3= x(4*i-1)
      x4= x(4*i  )
      y(4*i-3)= d(16*i-15)*x1 + d(16*i-14)*x2 + d(16*i-13)*x3 + d(16*i-12)*x4
      y(4*i-2)= d(16*i-11)*x1 + d(16*i-10)*x2 + d(16*i- 9)*x3 + d(16*i- 8)*x4
      y(4*i-1)= d(16*i- 7)*x1 + d(16*i- 6)*x2 + d(16*i- 5)*x3 + d(16*i- 4)*x4
      y(4*i  )= d(16*i- 3)*x1 + d(16*i- 2)*x2 + d(16*i- 1)*x3 + d(16*i- 0)*x4
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
    call matjad(hecmat%N, mjad, iajad, jajad, ajad, jadord, x, y, wp1, wp2, wp3, wp4)
  end subroutine hecmw_jad_matvec_44
 
  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
    integer(kind = kint) :: i, j, js, je, in, jc
    integer(kind = kint), allocatable :: len(:), lenz(:), jadreord(:)
 
    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(jc*16-15:jc*16) = hecmat%AL(16*j-15:16*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(jc*16-15:jc*16) = hecmat%AU(16*j-15:16*j)
        jajad(jc) = in
      end do
    end do
    deallocate(len)
    deallocate(jadreord)
    deallocate(lenz)
  end subroutine repack
 
  subroutine matjad(N, MJAD, IAJAD, JAJAD, AJAD, JADORD, X, Y, W1, W2, W3, W4)
    use hecmw_util
    integer(kind=kint) :: N, MJAD
    integer(kind=kint) :: IAJAD(*), JAJAD(*), JADORD(*)
    real(kind=kreal)   :: ajad(*), x(*), y(*), w1(*), w2(*), w3(*), w4(*)
 
    integer(kind=kint) :: I, K, NZ, IXX
    real(kind=kreal)   :: x1, x2, x3, x4
 
    !$OMP PARALLEL PRIVATE(I,K,X1,X2,X3,X4,IXX)
    !$OMP DO
    do i=1,n
      w1(i)=0.d0
      w2(i)=0.d0
      w3(i)=0.d0
      w4(i)=0.d0
    enddo
    !$OMP END DO
 
    do nz=1,mjad
      !$OMP DO
      do k=iajad(nz),iajad(nz+1)-1
        x1=x(jajad(k)*4-3)
        x2=x(jajad(k)*4-2)
        x3=x(jajad(k)*4-1)
        x4=x(jajad(k)*4  )
        ixx = k-iajad(nz)+1
        w1(ixx)=w1(ixx) + ajad(k*16-15)*x1 + ajad(k*16-14)*x2 + ajad(k*16-13)*x3 + ajad(k*16-12)*x4
        w2(ixx)=w2(ixx) + ajad(k*16-11)*x1 + ajad(k*16-10)*x2 + ajad(k*16- 9)*x3 + ajad(k*16- 8)*x4
        w3(ixx)=w3(ixx) + ajad(k*16- 7)*x1 + ajad(k*16- 6)*x2 + ajad(k*16- 5)*x3 + ajad(k*16- 4)*x4
        w4(ixx)=w4(ixx) + ajad(k*16- 3)*x1 + ajad(k*16- 2)*x2 + ajad(k*16- 1)*x3 + ajad(k*16- 0)*x4
      enddo
      !$OMP END DO
    enddo
 
    !$OMP DO
    do i=1,n
      y(4*i-3)=y(4*i-3)+w1(jadord(i))
      y(4*i-2)=y(4*i-2)+w2(jadord(i))
      y(4*i-1)=y(4*i-1)+w3(jadord(i))
      y(4*i  )=y(4*i  )+w4(jadord(i))
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
  end subroutine matjad
 
end module hecmw_jad_type_44