hecmw_solver_las_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_solver_las_33
  use hecmw_util
  use hecmw_matrix_misc
  implicit none
 
  private
 
  public :: hecmw_matvec_33_setup
  public :: hecmw_matvec_33_teardown
  public :: hecmw_matvec_33
  public :: hecmw_matresid_33
  public :: hecmw_rel_resid_l2_33
  public :: hecmw_tvec_33
  public :: hecmw_ttvec_33
 
  ! Implementation in use between hecmw_matvec_33_setup and hecmw_matvec_33_teardown.
  ! Callers outside that interval (eigenvalue analysis, contact, heat conduction) find
  ! BSR here, which is the only implementation that needs no data of its own.
  integer(kind=kint), save :: matvec_impl = hecmw_matvec_impl_bsr
 
  ! a format may be selected on a build that carries no matvec for it; saying so once
  ! per run keeps the fallback from looking like the requested format took effect
  logical, save :: matvec_impl_missing_reported = .false.
#ifdef HECMW_ARCH_FX64
 
  integer(kind=kint), save :: N, NP, NZIN
  integer(kind=kint), allocatable, save :: index_i(:), item_i(:)
  real(kind=kreal), allocatable, save   :: a_i(:)
#endif
 
  ! 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
  ! <<< added for tuning
 
contains
 
  !C
  !C***
  !C*** hecmw_matvec_33_setup
  !C***
  !C
  subroutine hecmw_matvec_33_setup (hecMESH, hecMAT)
    use hecmw_util
    use hecmw_jad_type
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in), target :: hecmat
 
    ! JAD takes over the whole matvec, so the tuned implementations must stay out of it
    if (hecmw_jad_is_initialized().ne.0) return
 
    matvec_impl = hecmw_mat_get_matvec_impl(hecmat)
    if (matvec_impl == hecmw_matvec_impl_csr) call hecmw_matvec_setup_33(hecmesh, hecmat)
  end subroutine hecmw_matvec_33_setup
 
  subroutine hecmw_matvec_setup_33 (hecMESH, hecMAT)
    use hecmw_jad_type
    use hecmw_matrix_misc
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    type (hecmwST_matrix), intent(in), target :: hecMAT
#ifdef HECMW_ARCH_FX64
 
    integer(kind=kint) :: i, j, k
 
    if (hecmw_jad_is_initialized().ne.0) return
 
    n = hecmat%N
    np = hecmat%NP
    if( allocated(index_i) ) deallocate(index_i)
    if( allocated(item_i) ) deallocate(item_i)
    if( allocated(a_i) ) deallocate(a_i)
    !if( allocated(index_o) ) deallocate(index_o)
    !if( allocated(item_o) ) deallocate(item_o)
    !if( allocated(A_o) ) deallocate(A_o)
 
    !count nonzero elements
    nzin = hecmat%indexL(n) + n + hecmat%indexU(n)
    !NZOUT = 0
    !! upper
    !do i=1,N
    !  do k=hecMAT%indexU(i-1)+1,hecMAT%indexU(i)
    !    j = hecMAT%itemU(k)
    !    if( j > N ) NZOUT = NZOUT + 1
    !  end do
    !end do
    !NZIN = NZIN - NZOUT
 
    allocate(index_i(0:n),item_i(3*nzin),a_i(9*nzin))
    !allocate(index_o(0:N),item_o(3*NZOUT),A_o(9*NZOUT))
 
    ! set value
    nzin = 0
    !NZOUT = 0
    index_i(0) = 0
    !index_o(0) = 0
    do i=1,n
      ! lower
      do k=hecmat%indexL(i-1)+1,hecmat%indexL(i)
        j = hecmat%itemL(k)
        item_i(nzin+1) = 3*j-2
        item_i(nzin+2) = 3*j-1
        item_i(nzin+3) = 3*j
        nzin = nzin + 3
        a_i(3*nzin-8:3*nzin) = hecmat%AL(9*k-8:9*k)
        call hecmw_mat_block_transpose_33(a_i(3*nzin-8:3*nzin))
      end do
      ! diag
      item_i(nzin+1) = 3*i-2
      item_i(nzin+2) = 3*i-1
      item_i(nzin+3) = 3*i
      nzin = nzin + 3
      a_i(3*nzin-8:3*nzin) = hecmat%D(9*i-8:9*i)
      call hecmw_mat_block_transpose_33(a_i(3*nzin-8:3*nzin))
      ! upper
      do k=hecmat%indexU(i-1)+1,hecmat%indexU(i)
        j = hecmat%itemU(k)
        !if( j > N ) then
        !  item_o(NZOUT+1) = 3*j-2
        !  item_o(NZOUT+2) = 3*j-1
        !  item_o(NZOUT+3) = 3*j
        !  NZOUT = NZOUT + 3
        !  A_o(3*NZOUT-8:3*NZOUT) = hecMAT%AU(9*k-8:9*k)
        !  call hecmw_mat_block_transpose_33(A_o(3*NZOUT-8:3*NZOUT))
        !else
          item_i(nzin+1) = 3*j-2
          item_i(nzin+2) = 3*j-1
          item_i(nzin+3) = 3*j
          nzin = nzin + 3
          a_i(3*nzin-8:3*nzin) = hecmat%AU(9*k-8:9*k)
          call hecmw_mat_block_transpose_33(a_i(3*nzin-8:3*nzin))
        !end if
      end do
 
      index_i(i) = nzin
      !index_o(i) = NZOUT
    end do
 
#endif
  end subroutine hecmw_matvec_setup_33
 
  !C
  !C***
  !C*** hecmw_matvec_33_teardown
  !C***
  !C
  subroutine hecmw_matvec_33_teardown
    implicit none
 
#ifdef HECMW_ARCH_FX64
    if( allocated(index_i) ) deallocate(index_i)
    if( allocated(item_i) ) deallocate(item_i)
    if( allocated(a_i) ) deallocate(a_i)
#endif
    matvec_impl = hecmw_matvec_impl_bsr
  end subroutine hecmw_matvec_33_teardown
 
  !C
  !C***
  !C*** hecmw_matvec_33
  !C***
  !C
  subroutine hecmw_matvec_33 (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
    use hecmw_util
    use hecmw_jad_type
    use m_hecmw_comm_f
 
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in), target :: hecmat
    real(kind=kreal), intent(inout) :: 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
 
    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
      select case (matvec_impl)
        case (hecmw_matvec_impl_bsr)
          call hecmw_matvec_33_generic(hecmesh, hecmat, x, y, time_ax, commtime)
#ifdef HECMW_ARCH_FX64
        case (hecmw_matvec_impl_csr)
          call hecmw_matvec_33_csr(hecmesh, hecmat, x, y, time_ax, commtime)
#endif
        case default
          if (.not. matvec_impl_missing_reported) then
            matvec_impl_missing_reported = .true.
            if (hecmw_comm_get_rank() == 0) write(*,'(a)') &
              '#### MATRIXFORMAT: this build has no matvec for the selected format '// &
              '-- running the BSR implementation'
          endif
          call hecmw_matvec_33_generic(hecmesh, hecmat, x, y, time_ax, commtime)
      end select
    endif
  end subroutine hecmw_matvec_33
 
  !C
  !C***
  !C*** hecmw_matvec_33_generic ( private subroutine )
  !C***
  !C
  subroutine hecmw_matvec_33_generic (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
    !$ use omp_lib
 
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    type (hecmwST_matrix), intent(in), target :: hecMAT
    real(kind=kreal), intent(inout) :: 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
    integer(kind=kint) :: i, j, jS, jE, in
    real(kind=kreal) :: yv1, yv2, yv3, x1, x2, x3
 
    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(kind=kint) :: threadNum, blockNum, numOfBlock
    integer(kind=kint) :: blockIndex
    ! <<< added for tuning
 
    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 (.not. isfirst) then
      numofblock = numofthread * numofblockperthread
      if (endpos(numofblock-1) .ne. n-1) then
        deallocate(startpos, endpos)
        isfirst = .true.
      endif
    endif
    if (isfirst) then
      call setup_tuning_parameters(n, np, indexl, indexu)
      isfirst = .false.
    endif
#endif
    ! <<< added for tuning
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, np, 3)
    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(3*i-2)
      x2= x(3*i-1)
      x3= x(3*i  )
      yv1= 0
      yv2= 0
      yv3= 0
      js= indexa(i-1) + 1
      je= indexa(i)
      do j = js, je
        in = itema(j)
        x1= x(3*in-2)
        x2= x(3*in-1)
        x3= x(3*in  )
        yv1= yv1 + a(9*j-8)*x1 + a(9*j-7)*x2 + a(9*j-6)*x3
        yv2= yv2 + a(9*j-5)*x1 + a(9*j-4)*x2 + a(9*j-3)*x3
        yv3= yv3 + a(9*j-2)*x1 + a(9*j-1)*x2 + a(9*j  )*x3
      enddo
      y(3*i-2)= yv1
      y(3*i-1)= yv2
      y(3*i  )= yv3
    enddo
    !$acc end kernels
#else
    !call fapp_start("loopInMatvec33", 1, 0)
    !call start_collection("loopInMatvec33")
 
    !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
    !OCL CACHE_SUBSECTOR_ASSIGN(X)
 
    !$OMP PARALLEL DEFAULT(NONE) &
      !$OMP&PRIVATE(i,X1,X2,X3,YV1,YV2,YV3,jS,jE,j,in,threadNum,blockNum,blockIndex) &
      !$OMP&SHARED(D,AL,AU,indexL,itemL,indexU,itemU,X,Y,startPos,endPos,numOfThread,N)
    threadnum = 0
    !$ threadNum = omp_get_thread_num()
    do blocknum = 0 , numofblockperthread - 1
      blockindex = blocknum * numofthread  + threadnum
      do i = startpos(blockindex), endpos(blockindex)
        x1= x(3*i-2)
        x2= x(3*i-1)
        x3= x(3*i  )
        yv1= d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
        yv2= d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
        yv3= d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i  )*x3
 
        js= indexl(i-1) + 1
        je= indexl(i  )
        do j= js, je
          in  = iteml(j)
          x1= x(3*in-2)
          x2= x(3*in-1)
          x3= x(3*in  )
          yv1= yv1 + al(9*j-8)*x1 + al(9*j-7)*x2 + al(9*j-6)*x3
          yv2= yv2 + al(9*j-5)*x1 + al(9*j-4)*x2 + al(9*j-3)*x3
          yv3= yv3 + al(9*j-2)*x1 + al(9*j-1)*x2 + al(9*j  )*x3
        enddo
        js= indexu(i-1) + 1
        je= indexu(i  )
        do j= js, je
          in  = itemu(j)
          x1= x(3*in-2)
          x2= x(3*in-1)
          x3= x(3*in  )
          yv1= yv1 + au(9*j-8)*x1 + au(9*j-7)*x2 + au(9*j-6)*x3
          yv2= yv2 + au(9*j-5)*x1 + au(9*j-4)*x2 + au(9*j-3)*x3
          yv3= yv3 + au(9*j-2)*x1 + au(9*j-1)*x2 + au(9*j  )*x3
        enddo
        y(3*i-2)= yv1
        y(3*i-1)= yv2
        y(3*i  )= yv3
      enddo
    enddo
    !$OMP END PARALLEL
 
    !OCL END_CACHE_SUBSECTOR
    !OCL END_CACHE_SECTOR_SIZE
 
    !call stop_collection("loopInMatvec33")
    !call fapp_stop("loopInMatvec33", 1, 0)
#endif
 
    end_time = hecmw_wtime()
    time_ax = time_ax + end_time - start_time
 
  end subroutine hecmw_matvec_33_generic
 
  !C
  !C***
  !C*** hecmw_matvec_33_csr ( private subroutine )
  !C***
  !C
  subroutine hecmw_matvec_33_csr (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
 
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    type (hecmwST_matrix), intent(in), target :: hecMAT
    real(kind=kreal), intent(inout) :: x(:)
    real(kind=kreal), intent(out) :: y(:)
    real(kind=kreal), intent(inout) :: time_ax
    real(kind=kreal), intent(inout), optional :: commtime
#ifdef HECMW_ARCH_FX64
 
    real(kind=kreal) :: start_time, end_time
    integer(kind=kint) :: N, NP
    integer(kind=kint) :: numOfBlock
 
    n = hecmat%N
    np = hecmat%NP
 
    ! added for tuning >>>
    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
      call setup_tuning_parameters(n, np, hecmat%indexL, hecmat%indexU)
      isfirst = .false.
    endif
    numofblock = numofthread * numofblockperthread
    ! <<< added for tuning
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, np, 3)
    end_time= hecmw_wtime()
    if (present(commtime)) commtime = commtime + end_time - start_time
 
    start_time = hecmw_wtime()
 
    !OCL CACHE_SECTOR_SIZE(sectorCacheSize0,sectorCacheSize1)
    call hecmw_matvec_33_core(n,np,nzin,x,y,a_i,item_i,index_i, &
      numofblock,numofthread,startpos,endpos)
 
    !OCL END_CACHE_SECTOR_SIZE
 
    end_time = hecmw_wtime()
    time_ax = time_ax + end_time - start_time
 
#endif
  end subroutine hecmw_matvec_33_csr
 
  subroutine hecmw_matvec_33_core(N,NP,NNZ,X,Y,aval,item,idx,numOfBlock,numOfThread,startPos,endPos)
    !$ use omp_lib
    integer(kind=kint), intent(in)  :: N
    integer(kind=kint), intent(in)  :: NP
    integer(kind=kint), intent(in)  :: NNZ
    real(kind=kreal), intent(in)    :: x(3*np)
    real(kind=kreal), intent(inout) :: y(3*n)
    real(kind=kreal), intent(in)    :: aval(3*nnz)
    integer(kind=kint), intent(in)  :: item(NNZ)
    integer(kind=kint), intent(in)  :: idx(0:N)
    integer(kind=kint), intent(in)  :: numOfBlock
    integer(kind=kint), intent(in)  :: numOfThread
    integer(kind=kint), intent(in)  :: startPos(0:numOfBlock)
    integer(kind=kint), intent(in)  :: endPos(0:numOfBlock)
#ifdef HECMW_ARCH_FX64
 
    integer(kind=kint) :: i, j, jS, jE, in
    real(kind=kreal) :: yv1, yv2, yv3, x1
    integer(kind=kint) :: threadNum, blockNum, blockIndex
 
    !call fapp_start("loopInMatvec33", 1, 0)
    !call start_collection("loopInMatvec33")
 
    !OCL CACHE_SUBSECTOR_ASSIGN(X)
 
    if( numofthread > 1 ) then
 
    !$OMP PARALLEL DEFAULT(NONE) &
      !$OMP&PRIVATE(i,X1,YV1,YV2,YV3,j,threadNum,blockNum,blockIndex) &
      !$OMP&SHARED(aval,item,idx,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)
        yv1= 0.d0
        yv2= 0.d0
        yv3= 0.d0
 
        do j= idx(i-1)+1, idx(i)
          x1= x(item(j))
          yv1= yv1 + aval(3*j-2)*x1
          yv2= yv2 + aval(3*j-1)*x1
          yv3= yv3 + aval(3*j  )*x1
        enddo
 
        y(3*i-2) = yv1
        y(3*i-1) = yv2
        y(3*i  ) = yv3
      enddo
    enddo
    !$OMP END PARALLEL
 
    else
 
      do i = 1, n
        yv1= 0.d0
        yv2= 0.d0
        yv3= 0.d0
 
        do j= idx(i-1)+1, idx(i)
          x1= x(item(j))
          yv1= yv1 + aval(3*j-2)*x1
          yv2= yv2 + aval(3*j-1)*x1
          yv3= yv3 + aval(3*j  )*x1
        enddo
 
        y(3*i-2)= yv1
        y(3*i-1)= yv2
        y(3*i  )= yv3
      enddo
 
    end if
 
 
    !OCL END_CACHE_SUBSECTOR
 
    !call stop_collection("loopInMatvec33")
    !call fapp_stop("loopInMatvec33", 1, 0)
 
#endif
  end subroutine
 
  subroutine setup_tuning_parameters(N, NP, indexL, indexU)
    use hecmw_tuning_fx
    !$ use omp_lib
    implicit none
    integer(kind=kint), intent(in) :: N, NP
    integer(kind=kint), pointer, intent(in) :: indexL(:), indexU(:)
    integer(kind=kint) :: i, blockNum, numOfBlock
    integer(kind=kint) :: numOfElement, elementCount
    real(kind=kreal) :: numofelementperblock
 
    !$ 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+hecmw_comm_get_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+hecmw_comm_get_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+hecmw_comm_get_rank(),*) mod(i, numOfThread), &
        !      startPos(i), endPos(i)
    end do
 
    call hecmw_tuning_fx_calc_sector_cache(np, 3, &
      sectorcachesize0, sectorcachesize1)
  end subroutine setup_tuning_parameters
 
  !C
  !C***
  !C*** hecmw_matresid_33
  !C***
  !C
  subroutine hecmw_matresid_33 (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) :: b(:)
    real(kind=kreal), intent(inout) :: x(:)
    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_33 (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 * 3
      r(i) = b(i) - r(i)
    enddo
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
  end subroutine hecmw_matresid_33
 
  !C
  !C***
  !C*** hecmw_rel_resid_L2_33
  !C***
  !C
  function hecmw_rel_resid_l2_33 (hecMESH, hecMAT, time_Ax, COMMtime)
    use hecmw_util
    use hecmw_solver_misc
    implicit none
    real(kind=kreal) :: hecmw_rel_resid_l2_33
    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_33(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_33 = sqrt(rnorm2 / bnorm2)
 
    if (present(commtime)) commtime = commtime + tcomm
 
    deallocate(r)
  end function hecmw_rel_resid_l2_33
 
  !C
  !C***
  !C*** hecmw_Tvec_33
  !C***
  !C
  subroutine hecmw_tvec_33 (hecMESH, X, Y, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    real(kind=kreal), intent(inout) :: 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, 3)
    end_time= hecmw_wtime()
    commtime = commtime + end_time - start_time
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent private(i)
#else
    !$omp parallel default(none),private(i,k,kk,j,jj),shared(hecMESH,X,Y)
    !$omp do
#endif
    do i= 1, hecmesh%nn_internal * hecmesh%n_dof
      y(i)= x(i)
    enddo
#ifndef _OPENACC
    !$omp end do
#endif
 
#ifdef _OPENACC
    !$acc loop independent private(i,k,kk,j,jj)
#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) > 3) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1) + 1
      kk = 3 * (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 = 3 * (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_33
 
  !C
  !C***
  !C*** hecmw_Ttvec_33
  !C***
  !C
  subroutine hecmw_ttvec_33 (hecMESH, X, Y, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    real(kind=kreal), intent(inout) :: x(:)
    real(kind=kreal), intent(inout) :: 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, 3)
    end_time= hecmw_wtime()
    commtime = commtime + end_time - start_time
 
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent private(i)
#else
    !$omp parallel default(none),private(i,k,kk,j,jj),shared(hecMESH,X,Y)
    !$omp do
#endif
    do i= 1, hecmesh%nn_internal * hecmesh%n_dof
      y(i)= x(i)
    enddo
#ifndef _OPENACC
    !$omp end do
#endif
 
#ifdef _OPENACC
    !$acc loop independent private(i,k,kk,j,jj)
#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) > 3) cycle outer
      enddo
      k = hecmesh%mpc%mpc_index(i-1) + 1
      kk = 3 * (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 = 3 * (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)
        !$acc end atomic
      enddo
    enddo outer
#ifdef _OPENACC
    !$acc end kernels
#else
    !$omp end do
    !$omp end parallel
#endif
 
  end subroutine hecmw_ttvec_33
 
end module hecmw_solver_las_33