hecmw_solver_las_66.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_66
  use hecmw_util
  implicit none
 
  private
 
  public :: hecmw_matvec_66
  public :: hecmw_matresid_66
  public :: hecmw_rel_resid_l2_66
 
contains
 
  !C
  !C***
  !C*** hecmw_matvec_66
  !C***
  !C
  subroutine hecmw_matvec_66 (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, js, je, in
    real(kind=kreal) :: yv1, yv2, yv3, x1, x2, x3
    real(kind=kreal) :: yv4, yv5, yv6, x4, x5, x6
 
    integer(kind=kint) :: n, np
    integer(kind=kint), pointer :: indexl(:), iteml(:), indexu(:), itemu(:), indexa(:), itema(:)
    real(kind=kreal), pointer :: al(:), au(:), d(:), a(:)
 
    ! added for tuning >>>
    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 tuning
 
    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
 
      ! added for tuning >>>
#ifndef _OPENACC
      if (isfirst .eqv. .true.) 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, 6, &
          sectorcachesize0, sectorcachesize1)
 
        isfirst = .false.
      endif
#endif
      ! <<< added for tuning
 
      start_time= hecmw_wtime()
      call hecmw_update_r (hecmesh, x, np, 6)
      end_time= hecmw_wtime()
      if (present(commtime)) commtime = commtime + end_time - start_time
 
      start_time = hecmw_wtime()
 
#ifdef _OPENACC
      !$acc kernels
      !$acc loop independent
      do i = 1, n
        x1= x(6*i-5)
        x2= x(6*i-4)
        x3= x(6*i-3)
        x4= x(6*i-2)
        x5= x(6*i-1)
        x6= x(6*i  )
        yv1= 0
        yv2= 0
        yv3= 0
        yv4= 0
        yv5= 0
        yv6= 0
        js= indexa(i-1) + 1
        je= indexa(i)
        do j= js, je
          in  = itema(j)
          x1= x(6*in-5)
          x2= x(6*in-4)
          x3= x(6*in-3)
          x4= x(6*in-2)
          x5= x(6*in-1)
          x6= x(6*in  )
          yv1= yv1 + a(36*j-35)*x1 + a(36*j-34)*x2 + a(36*j-33)*x3 + a(36*j-32)*x4 + a(36*j-31)*x5 + a(36*j-30)*x6
          yv2= yv2 + a(36*j-29)*x1 + a(36*j-28)*x2 + a(36*j-27)*x3 + a(36*j-26)*x4 + a(36*j-25)*x5 + a(36*j-24)*x6
          yv3= yv3 + a(36*j-23)*x1 + a(36*j-22)*x2 + a(36*j-21)*x3 + a(36*j-20)*x4 + a(36*j-19)*x5 + a(36*j-18)*x6
          yv4= yv4 + a(36*j-17)*x1 + a(36*j-16)*x2 + a(36*j-15)*x3 + a(36*j-14)*x4 + a(36*j-13)*x5 + a(36*j-12)*x6
          yv5= yv5 + a(36*j-11)*x1 + a(36*j-10)*x2 + a(36*j-9 )*x3 + a(36*j-8 )*x4 + a(36*j-7 )*x5 + a(36*j-6 )*x6
          yv6= yv6 + a(36*j-5 )*x1 + a(36*j-4 )*x2 + a(36*j-3 )*x3 + a(36*j-2 )*x4 + a(36*j-1 )*x5 + a(36*j   )*x6
        enddo
        y(6*i-5)= yv1
        y(6*i-4)= yv2
        y(6*i-3)= yv3
        y(6*i-2)= yv4
        y(6*i-1)= yv5
        y(6*i  )= yv6
      enddo
      !$acc end kernels
#else
      !call fapp_start("loopInMatvec66", 1, 0)
      !call start_collection("loopInMatvec66")
 
      !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
      !OCL CACHE_SUBSECTOR_ASSIGN(X)
 
      !$OMP PARALLEL DEFAULT(NONE) &
        !$OMP&PRIVATE(i,X1,X2,X3,X4,X5,X6,YV1,YV2,YV3,YV4,YV5,YV6,jS,jE,j,in,threadNum,blockNum,blockIndex) &
        !$OMP&SHARED(D,AL,AU,indexL,itemL,indexU,itemU,X,Y,startPos,endPos,numOfThread)
      threadnum = 0
      !$ threadNum = omp_get_thread_num()
      do blocknum = 0 , numofblockperthread - 1
        blockindex = blocknum * numofthread  + threadnum
        do i = startpos(blockindex), endpos(blockindex)
          x1= x(6*i-5)
          x2= x(6*i-4)
          x3= x(6*i-3)
          x4= x(6*i-2)
          x5= x(6*i-1)
          x6= x(6*i  )
          yv1= d(36*i-35)*x1 + d(36*i-34)*x2 + d(36*i-33)*x3 + d(36*i-32)*x4 + d(36*i-31)*x5 + d(36*i-30)*x6
          yv2= d(36*i-29)*x1 + d(36*i-28)*x2 + d(36*i-27)*x3 + d(36*i-26)*x4 + d(36*i-25)*x5 + d(36*i-24)*x6
          yv3= d(36*i-23)*x1 + d(36*i-22)*x2 + d(36*i-21)*x3 + d(36*i-20)*x4 + d(36*i-19)*x5 + d(36*i-18)*x6
          yv4= d(36*i-17)*x1 + d(36*i-16)*x2 + d(36*i-15)*x3 + d(36*i-14)*x4 + d(36*i-13)*x5 + d(36*i-12)*x6
          yv5= d(36*i-11)*x1 + d(36*i-10)*x2 + d(36*i-9 )*x3 + d(36*i-8 )*x4 + d(36*i-7 )*x5 + d(36*i-6 )*x6
          yv6= d(36*i-5 )*x1 + d(36*i-4 )*x2 + d(36*i-3 )*x3 + d(36*i-2 )*x4 + d(36*i-1 )*x5 + d(36*i   )*x6
 
          js= indexl(i-1) + 1
          je= indexl(i  )
          do j= js, je
            in  = iteml(j)
            x1= x(6*in-5)
            x2= x(6*in-4)
            x3= x(6*in-3)
            x4= x(6*in-2)
            x5= x(6*in-1)
            x6= x(6*in  )
            yv1= yv1 + al(36*j-35)*x1 + al(36*j-34)*x2 + al(36*j-33)*x3 + al(36*j-32)*x4 + al(36*j-31)*x5 + al(36*j-30)*x6
            yv2= yv2 + al(36*j-29)*x1 + al(36*j-28)*x2 + al(36*j-27)*x3 + al(36*j-26)*x4 + al(36*j-25)*x5 + al(36*j-24)*x6
            yv3= yv3 + al(36*j-23)*x1 + al(36*j-22)*x2 + al(36*j-21)*x3 + al(36*j-20)*x4 + al(36*j-19)*x5 + al(36*j-18)*x6
            yv4= yv4 + al(36*j-17)*x1 + al(36*j-16)*x2 + al(36*j-15)*x3 + al(36*j-14)*x4 + al(36*j-13)*x5 + al(36*j-12)*x6
            yv5= yv5 + al(36*j-11)*x1 + al(36*j-10)*x2 + al(36*j-9 )*x3 + al(36*j-8 )*x4 + al(36*j-7 )*x5 + al(36*j-6 )*x6
            yv6= yv6 + al(36*j-5 )*x1 + al(36*j-4 )*x2 + al(36*j-3 )*x3 + al(36*j-2 )*x4 + al(36*j-1 )*x5 + al(36*j   )*x6
          enddo
          js= indexu(i-1) + 1
          je= indexu(i  )
          do j= js, je
            in  = itemu(j)
            x1= x(6*in-5)
            x2= x(6*in-4)
            x3= x(6*in-3)
            x4= x(6*in-2)
            x5= x(6*in-1)
            x6= x(6*in  )
            yv1= yv1 + au(36*j-35)*x1 + au(36*j-34)*x2 + au(36*j-33)*x3 + au(36*j-32)*x4 + au(36*j-31)*x5 + au(36*j-30)*x6
            yv2= yv2 + au(36*j-29)*x1 + au(36*j-28)*x2 + au(36*j-27)*x3 + au(36*j-26)*x4 + au(36*j-25)*x5 + au(36*j-24)*x6
            yv3= yv3 + au(36*j-23)*x1 + au(36*j-22)*x2 + au(36*j-21)*x3 + au(36*j-20)*x4 + au(36*j-19)*x5 + au(36*j-18)*x6
            yv4= yv4 + au(36*j-17)*x1 + au(36*j-16)*x2 + au(36*j-15)*x3 + au(36*j-14)*x4 + au(36*j-13)*x5 + au(36*j-12)*x6
            yv5= yv5 + au(36*j-11)*x1 + au(36*j-10)*x2 + au(36*j-9 )*x3 + au(36*j-8 )*x4 + au(36*j-7 )*x5 + au(36*j-6 )*x6
            yv6= yv6 + au(36*j-5 )*x1 + au(36*j-4 )*x2 + au(36*j-3 )*x3 + au(36*j-2 )*x4 + au(36*j-1 )*x5 + au(36*j   )*x6
          enddo
          y(6*i-5)= yv1
          y(6*i-4)= yv2
          y(6*i-3)= yv3
          y(6*i-2)= yv4
          y(6*i-1)= yv5
          y(6*i  )= yv6
        enddo
      enddo
      !$OMP END PARALLEL
 
      !OCL END_CACHE_SUBSECTOR
      !OCL END_CACHE_SECTOR_SIZE
 
      !call stop_collection("loopInMatvec66")
      !call fapp_stop("loopInMatvec66", 1, 0)
#endif
 
      end_time = hecmw_wtime()
      time_ax = time_ax + end_time - start_time
 
    endif
  end subroutine hecmw_matvec_66
 
  !C
  !C***
  !C*** hecmw_matresid_66
  !C***
  !C
  subroutine hecmw_matresid_66 (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_66 (hecmesh, hecmat, x, r, time_ax, tcomm)
    if (present(commtime)) commtime = commtime + tcomm
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
#endif
    do i = 1, hecmat%N * 6
      r(i) = b(i) - r(i)
    enddo
#ifdef _OPENACC
    !$acc end kernels
#endif
 
  end subroutine hecmw_matresid_66
 
  !C
  !C***
  !C*** hecmw_rel_resid_L2_66
  !C***
  !C
  function hecmw_rel_resid_l2_66 (hecMESH, hecMAT, time_Ax, COMMtime)
    use hecmw_util
    use hecmw_solver_misc
    implicit none
    real(kind=kreal) :: hecmw_rel_resid_l2_66
    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_66(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_66 = sqrt(rnorm2 / bnorm2)
 
    if (present(commtime)) commtime = commtime + tcomm
 
    deallocate(r)
  end function hecmw_rel_resid_l2_66
 
end module hecmw_solver_las_66