hecmw_jadm_33.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_33
  use hecmw_util
  use m_hecmw_comm_f
  implicit none
 
  private
 
  public :: hecmw_jad_init_33
  public :: hecmw_jad_finalize_33
  public :: hecmw_jad_matvec_33
 
  !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(:)
  integer(kind=kint) :: INITIALIZED = 0
 
contains
 
  subroutine hecmw_jad_init_33(hecMAT)
    type(hecmwst_matrix) :: hecmat
    allocate(wp1(hecmat%NP), wp2(hecmat%NP), wp3(hecmat%NP))
    allocate(ajad((hecmat%NPL+hecmat%NPU)*9))
    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_33
 
  subroutine hecmw_jad_finalize_33()
    deallocate(ajad)
    deallocate(jajad)
    deallocate(jadord)
    deallocate(iajad)
    deallocate(wp1,wp2,wp3)
    initialized = 0
  end subroutine hecmw_jad_finalize_33
 
  subroutine hecmw_jad_matvec_33(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
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmat%NP, 3)
    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(3*i-2)
      x2= x(3*i-1)
      x3= x(3*i  )
      y(3*i -2)= d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
      y(3*i -1)= d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
      y(3*i   )= d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i  )*x3
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
 
    call matjad(hecmat%N, mjad, iajad, jajad, ajad, jadord, x, y, wp1, wp2, wp3)
  end subroutine hecmw_jad_matvec_33
 
  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*9-8:jc*9) = hecmat%AL(9*j-8:9*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*9-8:jc*9) = hecmat%AU(9*j-8:9*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)
    use hecmw_util
    integer(kind=kint) :: N, MJAD
    integer(kind=kint) :: IAJAD(*), JAJAD(*), JADORD(*)
    real(kind=kreal)   :: ajad(*), x(*), y(*), w1(*), w2(*), w3(*)
 
    integer(kind=kint) :: I, K, NZ, IXX
    real(kind=kreal)   :: x1, x2, x3
 
    !$OMP PARALLEL PRIVATE(I,K,X1,X2,X3,IXX)
    !$OMP DO
    do i=1,n
      w1(i)=0.d0
      w2(i)=0.d0
      w3(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)*3-2)
        x2=x(jajad(k)*3-1)
        x3=x(jajad(k)*3  )
        ixx = k-iajad(nz)+1
        w1(ixx)=w1(ixx) + ajad(k*9-8)*x1 + ajad(k*9-7)*x2 + ajad(k*9-6)*x3
        w2(ixx)=w2(ixx) + ajad(k*9-5)*x1 + ajad(k*9-4)*x2 + ajad(k*9-3)*x3
        w3(ixx)=w3(ixx) + ajad(k*9-2)*x1 + ajad(k*9-1)*x2 + ajad(k*9-0)*x3
      enddo
      !$OMP END DO
    enddo
 
    !$OMP DO
    do i=1,n
      y(3*i-2)=y(3*i-2)+w1(jadord(i))
      y(3*i-1)=y(3*i-1)+w2(jadord(i))
      y(3*i  )=y(3*i  )+w3(jadord(i))
    enddo
    !$OMP END DO
    !$OMP END PARALLEL
  end subroutine matjad
 
end module hecmw_jad_type_33