hecmw_solver_las_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_solver_las_nn
  use hecmw_util
  implicit none
 
  private
 
  public :: hecmw_matvec_nn
  public :: hecmw_matvec_nn_set_async
  public :: hecmw_matvec_nn_unset_async
  public :: hecmw_matresid_nn
  public :: hecmw_rel_resid_l2_nn
  public :: hecmw_tvec_nn
  public :: hecmw_ttvec_nn
  public :: hecmw_ttmattvec_nn
  public :: hecmw_mat_diag_sr_nn
  public :: hecmw_mat_add_nn
  public :: hecmw_mat_multiple_nn
 
  ! ! for communication hiding in matvec
  ! integer(kind=kint), save, allocatable :: index_o(:), item_o(:)
  ! real(kind=kreal), save, allocatable :: A_o(:)
  logical, save :: async_matvec_flg = .false.
 
contains
 
  !C
  !C***
  !C*** hecmw_matvec_nn
  !C***
  !C
  subroutine hecmw_matvec_nn (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
    use hecmw_util
    use hecmw_matrix_misc
    implicit none
    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) :: time_ax
    real(kind=kreal), intent(inout), optional :: commtime
 
    real(kind=kreal) :: tcomm
    real(kind=kreal), allocatable :: wk(:)
 
    tcomm = 0.d0
 
    if (hecmw_mat_get_flag_mpcmatvec(hecmat) /= 0) then
      allocate(wk(hecmat%NP * hecmat%NDOF))
      call hecmw_ttmattvec_nn(hecmesh, hecmat, x, y, wk, time_ax, tcomm)
      deallocate(wk)
    else
      call hecmw_matvec_nn_inner(hecmesh, hecmat, x, y, time_ax, tcomm)
    endif
 
    if (present(commtime)) commtime = commtime + tcomm
  end subroutine hecmw_matvec_nn
 
  !C
  !C***
  !C*** hecmw_matvec_nn_set_async
  !C***
  !C
  subroutine hecmw_matvec_nn_set_async (hecMAT)
    use hecmw_util
    implicit none
    type (hecmwst_matrix), intent(in) :: hecmat
    ! integer(kind=kint) :: i, j, jS, jE, idx, in
 
    ! allocate(index_o(0:hecMAT%N))
    ! index_o(0) = 0
    ! do i = 1, hecMAT%N
    !   jS= hecMAT%indexU(i-1) + 1
    !   jE= hecMAT%indexU(i  )
    !   idx = index_o(i-1)
    !   do j= jS, jE
    !     in  = hecMAT%itemU(j)
    !     if (in <= hecMAT%N) cycle
    !     idx = idx + 1
    !   enddo
    !   index_o(i) = idx
    ! enddo
    ! allocate(item_o(idx))
    ! allocate(A_o(idx*9))
    ! do i = 1, hecMAT%N
    !   jS= hecMAT%indexU(i-1) + 1
    !   jE= hecMAT%indexU(i  )
    !   idx = index_o(i-1)
    !   do j= jS, jE
    !     in  = hecMAT%itemU(j)
    !     if (in <= hecMAT%N) cycle
    !     idx = idx + 1
    !     item_o(idx) = hecMAT%itemU(j) - hecMAT%N
    !     A_o(9*idx-8:9*idx) = hecMAT%AU(9*j-8:9*j)
    !   enddo
    ! enddo
    ! async_matvec_flg = .true.
  end subroutine hecmw_matvec_nn_set_async
 
  !C
  !C***
  !C*** hecmw_matvec_nn_unset_async
  !C***
  !C
  subroutine hecmw_matvec_nn_unset_async
    implicit none
    ! if (allocated(index_o)) deallocate(index_o)
    ! if (allocated(item_o)) deallocate(item_o)
    ! if (allocated(A_o)) deallocate(A_o)
    ! async_matvec_flg = .false.
  end subroutine hecmw_matvec_nn_unset_async
 
  !C
  !C***
  !C*** hecmw_matvec_nn_inner ( private subroutine )
  !C***
  !C
  subroutine hecmw_matvec_nn_inner (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
    use hecmw_matrix_misc
    use hecmw_jad_type
    use hecmw_tuning_fx
    !$ use omp_lib
 
    implicit none
    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) :: time_ax
    real(kind=kreal), intent(inout), optional :: commtime
 
    real(kind=kreal) :: start_time, end_time, tcomm
    integer(kind=kint) :: i, j, k, l, js, je, in
    real(kind=kreal) :: yv(hecmat%NDOF), xv(hecmat%NDOF)
 
    integer(kind=kint) :: n, np, ndof, ndof2
    integer(kind=kint), pointer :: indexl(:), iteml(:), indexu(:), itemu(:), indexa(:), itema(:)
    real(kind=kreal), pointer :: al(:), au(:), d(:), a(:)
 
    ! added for turning >>>
    integer, parameter :: numofblockperthread = 100
    logical, save :: isfirst = .true.
    integer, save :: numofthread = 1
    integer, save, allocatable :: startpos(:), endpos(:)
    integer(kind=kint), save :: sectorcachesize0, sectorcachesize1
    integer(kind=kint) :: threadnum, blocknum, numofblock
    integer(kind=kint) :: numofelement, elementcount, blockindex
    real(kind=kreal) :: numofelementperblock
    ! <<< added for turning
 
    if (hecmw_jad_is_initialized().ne.0) then
      tcomm = 0.d0
      start_time = hecmw_wtime()
      call hecmw_jad_matvec(hecmesh, hecmat, x, y, tcomm)
      end_time = hecmw_wtime()
      time_ax = time_ax + end_time - start_time - tcomm
      if (present(commtime)) commtime = commtime + tcomm
    else
 
      n = hecmat%N
      np = hecmat%NP
      indexl => hecmat%indexL
      indexu => hecmat%indexU
      indexa => hecmat%indexA
      iteml => hecmat%itemL
      itemu => hecmat%itemU
      itema => hecmat%itemA
      al => hecmat%AL
      au => hecmat%AU
      d => hecmat%D
      a => hecmat%A
      ndof =  hecmat%NDOF
      ndof2 = ndof*ndof
 
      ! added for turning >>>
#ifndef _OPENACC
      if (.not. isfirst) then
        numofblock = numofthread * numofblockperthread
        if (endpos(numofblock-1) .ne. n-1) then
          deallocate(startpos, endpos)
          isfirst = .true.
        endif
      endif
      if (isfirst) then
        !$ numOfThread = omp_get_max_threads()
        numofblock = numofthread * numofblockperthread
        allocate (startpos(0 : numofblock - 1), endpos(0 : numofblock - 1))
        numofelement = n + indexl(n) + indexu(n)
        numofelementperblock = dble(numofelement) / numofblock
        blocknum = 0
        elementcount = 0
        startpos(blocknum) = 1
        do i= 1, n
          elementcount = elementcount + 1
          elementcount = elementcount + (indexl(i) - indexl(i-1))
          elementcount = elementcount + (indexu(i) - indexu(i-1))
          if (elementcount > (blocknum + 1) * numofelementperblock) then
            endpos(blocknum) = i
            ! write(9000+hecMESH%my_rank,*) mod(blockNum, numOfThread), &
              !      startPos(blockNum), endPos(blockNum)
            blocknum = blocknum + 1
            startpos(blocknum) = i + 1
            if (blocknum == (numofblock - 1)) exit
          endif
        enddo
        endpos(blocknum) = n
        ! write(9000+hecMESH%my_rank,*) mod(blockNum, numOfThread), &
          !      startPos(blockNum), endPos(blockNum)
        ! for irregular data
        do i= blocknum+1, numofblock-1
          startpos(i) = n
          endpos(i) = n-1
          ! write(9000+hecMESH%my_rank,*) mod(i, numOfThread), &
            !      startPos(i), endPos(i)
        end do
 
        call hecmw_tuning_fx_calc_sector_cache(np, ndof, &
          sectorcachesize0, sectorcachesize1)
 
        isfirst = .false.
      endif
#endif
      ! <<< added for turning
 
      start_time= hecmw_wtime()
      ! if (async_matvec_flg) then
      !   call hecmw_update_3_R_async (hecMESH, X, NP, ireq)
      ! else
      call hecmw_update_r (hecmesh, x, np, ndof)
      ! endif
      end_time= hecmw_wtime()
      if (present(commtime)) commtime = commtime + end_time - start_time
 
      start_time = hecmw_wtime()
 
      !call fapp_start("loopInMatvec33", 1, 0)
      !call start_collection("loopInMatvec33")
 
      !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
      !OCL CACHE_SUBSECTOR_ASSIGN(X)
 
#ifdef _OPENACC
      !$acc kernels
      !$acc loop independent
      do i = 1, n
        do k=1,ndof
          xv(k) = x(ndof*(i-1)+k)
        end do
        yv(:)=0.0d0
        js= indexa(i) + 1
        je= indexa(i+1)
        do j= js, je
          in  = itema(j)
          do k=1,ndof
            xv(k) = x(ndof*(in-1)+k)
          end do
          do k=1,ndof
            do l=1,ndof
              yv(k)=yv(k)+a(ndof2*(j-1)+(k-1)*ndof+l)*xv(l)
            end do
          end do
        enddo
        do k=1,ndof
          y(ndof*(i-1)+k) = yv(k)
        end do
      enddo
      !$acc end kernels
#else
      !$OMP PARALLEL DEFAULT(NONE) &
        !$OMP&PRIVATE(i,XV,YV,jS,jE,j,k,l,in,threadNum,blockNum,blockIndex) &
        !$OMP&SHARED(D,AL,AU,indexL,itemL,indexU,itemU,X,Y,startPos,endPos,numOfThread,N,NDOF,NDOF2,async_matvec_flg)
      threadnum = 0
      !$ threadNum = omp_get_thread_num()
      do blocknum = 0 , numofblockperthread - 1
        blockindex = blocknum * numofthread  + threadnum
        do i = startpos(blockindex), endpos(blockindex)
          do k=1,ndof
            xv(k) = x(ndof*(i-1)+k)
          end do
          yv(:)=0.0d0
          do k=1,ndof
            do l=1,ndof
              yv(k)=yv(k)+d(ndof2*(i-1)+(k-1)*ndof+l)*xv(l)
            end do
          end do
          js= indexl(i-1) + 1
          je= indexl(i  )
          do j= js, je
            in  = iteml(j)
            do k=1,ndof
              xv(k) = x(ndof*(in-1)+k)
            end do
            do k=1,ndof
              do l=1,ndof
                yv(k)=yv(k)+al(ndof2*(j-1)+(k-1)*ndof+l)*xv(l)
              end do
            end do
          enddo
          js= indexu(i-1) + 1
          je= indexu(i  )
          do j= js, je
            in  = itemu(j)
            ! if (async_matvec_flg .and. in > N) cycle
            do k=1,ndof
              xv(k) = x(ndof*(in-1)+k)
            end do
            do k=1,ndof
              do l=1,ndof
                yv(k)=yv(k)+au(ndof2*(j-1)+(k-1)*ndof+l)*xv(l)
              end do
            end do
          enddo
          do k=1,ndof
            y(ndof*(i-1)+k) = yv(k)
          end do
        enddo
      enddo
      !$OMP END PARALLEL
#endif
 
      !OCL END_CACHE_SUBSECTOR
      !OCL END_CACHE_SECTOR_SIZE
 
      !call stop_collection("loopInMatvec33")
      !call fapp_stop("loopInMatvec33", 1, 0)
 
      end_time = hecmw_wtime()
      time_ax = time_ax + end_time - start_time
    endif
 
  end subroutine hecmw_matvec_nn_inner
 
 
 
  !C
  !C***
  !C*** hecmw_matresid_nn
  !C***
  !C
  subroutine hecmw_matresid_nn (hecMESH, hecMAT, X, B, R, time_Ax, COMMtime)
    use hecmw_util
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in)     :: hecmat
    real(kind=kreal), intent(in) :: x(:), b(:)
    real(kind=kreal), intent(out) :: r(:)
    real(kind=kreal), intent(inout) :: time_ax
    real(kind=kreal), intent(inout), optional :: commtime
 
    integer(kind=kint) :: i
    real(kind=kreal) :: tcomm
 
    tcomm = 0.d0
    call hecmw_matvec_nn (hecmesh, hecmat, x, r, time_ax, tcomm)
    if (present(commtime)) commtime = commtime + tcomm
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
#else
    !$omp parallel default(none),private(i),shared(hecMAT,R,B)
    !$omp do
#endif
    do i = 1, hecmat%N * hecmat%NDOF
      r(i) = b(i) - r(i)
    enddo
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
  end subroutine hecmw_matresid_nn
 
  !C
  !C***
  !C*** hecmw_rel_resid_L2_nn
  !C***
  !C
  function hecmw_rel_resid_l2_nn (hecMESH, hecMAT, time_Ax, COMMtime)
    use hecmw_util
    use hecmw_solver_misc
    implicit none
    real(kind=kreal) :: hecmw_rel_resid_l2_nn
    type ( hecmwst_local_mesh ), intent(in) :: hecmesh
    type ( hecmwst_matrix     ), intent(in) :: hecmat
    real(kind=kreal), intent(inout) :: time_ax
    real(kind=kreal), intent(inout), optional :: commtime
 
    real(kind=kreal), allocatable :: r(:)
    real(kind=kreal) :: bnorm2, rnorm2
    real(kind=kreal) :: tcomm
 
    allocate(r(hecmat%NDOF*hecmat%NP))
 
    tcomm = 0.d0
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, &
      hecmat%B, hecmat%B, bnorm2, tcomm)
    if (bnorm2 == 0.d0) then
      bnorm2 = 1.d0
    endif
    call hecmw_matresid_nn(hecmesh, hecmat, hecmat%X, hecmat%B, r, time_ax, tcomm)
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, r, r, rnorm2, tcomm)
    hecmw_rel_resid_l2_nn = sqrt(rnorm2 / bnorm2)
 
    if (present(commtime)) commtime = commtime + tcomm
 
    deallocate(r)
  end function hecmw_rel_resid_l2_nn
 
  !C
  !C***
  !C*** hecmw_Tvec_nn
  !C***
  !C
  subroutine hecmw_tvec_nn (hecMESH, ndof, X, Y, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: ndof
    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
    integer(kind=kint) :: i, j, jj, k, kk
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmesh%n_node, ndof)
    end_time= hecmw_wtime()
    commtime = commtime + end_time - start_time
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
#else
    !$omp parallel default(none),private(i,k,kk,j,jj),shared(hecMESH,X,Y),firstprivate(ndof)
    !$omp do
#endif
    do i= 1, hecmesh%nn_internal * ndof
      y(i)= x(i)
    enddo
#ifndef _OPENACC
    !$omp end do
#endif
 
#ifdef _OPENACC
    !$acc loop independent
#else
    !$omp do
#endif
    outer: do i= 1, hecmesh%mpc%n_mpc
      do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
        if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1) + 1
      kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
      y(kk) = 0.d0
      do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
        jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
        y(kk) = y(kk) - hecmesh%mpc%mpc_val(j) * x(jj)
      enddo
    enddo outer
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
 
  end subroutine hecmw_tvec_nn
 
  !C
  !C***
  !C*** hecmw_Ttvec_nn
  !C***
  !C
  subroutine hecmw_ttvec_nn (hecMESH, ndof, X, Y, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: ndof
    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
    integer(kind=kint) :: i, j, jj, k, kk
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmesh%n_node,ndof)
    end_time= hecmw_wtime()
    commtime = commtime + end_time - start_time
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
#else
    !$omp parallel default(none),private(i,k,kk,j,jj),shared(hecMESH,X,Y),firstprivate(ndof)
    !$omp do
#endif
    do i= 1, hecmesh%nn_internal * ndof
      y(i)= x(i)
    enddo
#ifndef _OPENACC
    !$omp end do
#endif
 
#ifdef _OPENACC
    !$acc loop independent
#else
    !$omp do
#endif
    outer: do i= 1, hecmesh%mpc%n_mpc
      do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
        if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1) + 1
      kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
      y(kk) = 0.d0
      do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
        jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
#ifdef _OPENACC
        !$acc atomic update
#else
        !$omp atomic
#endif
        y(jj) = y(jj) - hecmesh%mpc%mpc_val(j) * x(kk)
      enddo
    enddo outer
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
 
 
  end subroutine hecmw_ttvec_nn
 
  !C
  !C***
  !C*** hecmw_TtmatTvec_nn
  !C***
  !C
  subroutine hecmw_ttmattvec_nn (hecMESH, hecMAT, X, Y, W, time_Ax, COMMtime)
    use hecmw_util
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in)     :: hecmat
    real(kind=kreal), intent(in) :: x(:)
    real(kind=kreal), intent(out) :: y(:), w(:)
    real(kind=kreal), intent(inout) :: time_ax
    real(kind=kreal), intent(inout) :: commtime
 
    call hecmw_tvec_nn(hecmesh, hecmat%NDOF, x, y, commtime)
    call hecmw_matvec_nn_inner(hecmesh, hecmat, y, w, time_ax, commtime)
    call hecmw_ttvec_nn(hecmesh, hecmat%NDOF, w, y, commtime)
 
  end subroutine hecmw_ttmattvec_nn
 
  !C
  !C***
  !C*** hecmw_mat_diag_sr_nn
  !C***
  !C
  subroutine hecmw_mat_diag_sr_nn(hecMESH, hecMAT, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(inout), target :: hecmat
    real(kind=kreal), intent(inout), optional :: commtime
    real(kind=kreal), allocatable :: w(:,:)
    real(kind=kreal), pointer :: d(:)
    integer(kind=kint) :: ip, ndof, i, j
    real(kind=kreal) :: start_time, end_time
    ndof = hecmat%NDOF
    allocate(w(ndof*hecmat%NP,ndof))
    d => hecmat%D
    do ip= 1, hecmat%N
      do i=1,ndof
        do j=1,ndof
          w(ndof*(ip-1)+i,j) = d(ndof*ndof*(ip-1)+(i-1)*ndof+j)
        end do
      end do
    enddo
    start_time= hecmw_wtime()
    do i=1,ndof
      call hecmw_update_r (hecmesh, w(:,i), hecmat%NP, ndof)
    end do
    end_time= hecmw_wtime()
    if (present(commtime)) commtime = commtime + end_time - start_time
    do ip= hecmat%N+1, hecmat%NP
      do i=1,ndof
        do j=1,ndof
          d(ndof*ndof*(ip-1)+(i-1)*ndof+j) = w(ndof*(ip-1)+i,j)
        end do
      end do
    enddo
    deallocate(w)
  end subroutine hecmw_mat_diag_sr_nn
 
  subroutine hecmw_mat_add_nn(hecMAT1, hecMAT2, hecMAT3)
    use hecmw_util
    implicit none
    type (hecmwst_matrix)     :: hecmat1, hecmat2, hecmat3
    integer(kind=kint) :: i
 
    do i = 1, hecmat1%NP*hecmat1%NDOF*hecmat1%NDOF
      hecmat3%D(i) = hecmat1%D(i) + hecmat2%D(i)
    enddo
 
    do i = 1, hecmat1%NPU*hecmat1%NDOF*hecmat1%NDOF
      hecmat3%AU(i) = hecmat1%AU(i) + hecmat2%AU(i)
    enddo
 
    do i = 1, hecmat1%NPL*hecmat1%NDOF*hecmat1%NDOF
      hecmat3%AL(i) = hecmat1%AL(i) + hecmat2%AL(i)
    enddo
  end subroutine hecmw_mat_add_nn
 
  subroutine hecmw_mat_multiple_nn(hecMAT, alpha)
    use hecmw_util
    implicit none
    type (hecmwst_matrix)     :: hecmat
    real(kind=kreal), intent(in) :: alpha
    integer(kind=kint) :: i
 
    do i = 1, hecmat%NP*hecmat%NDOF*hecmat%NDOF
      hecmat%D(i) = alpha*hecmat%D(i)
    enddo
 
    do i = 1, hecmat%NPU*hecmat%NDOF*hecmat%NDOF
      hecmat%AU(i) = alpha*hecmat%AU(i)
    enddo
 
    do i = 1, hecmat%NPL*hecmat%NDOF*hecmat%NDOF
      hecmat%AL(i) = alpha*hecmat%AL(i)
    enddo
  end subroutine hecmw_mat_multiple_nn
end module hecmw_solver_las_nn