LoginSignup
0
0
@snowxcrash

Fortranで有限要素法解析 04 : cuSOLVER を使う 〜 其の2 高速化 〜

Last updated at Posted at 2024-06-02

前回，LAPACK の関数を cuSOLVER に置き換えて GPU オフロードを実現しました。

Fortranで有限要素法解析 03 : cuSOLVER を使う 〜 其の1 ナイーブな実装 〜

しかし，そこでは方程式を解いて新しい時間ステップのデータを得るたびに GPU から CPU にデータを戻し，次の解を計算する際に再び CPU から GPU にデータコピーを行うという，無駄に冗長なコードになっていました。今回，これを解消して高速化することを目指します。

LSVQR

前回作成した LSVQR（疎行列用の関数）を使用するプログラムで解説していきます。
CPU→GPU (Host-to-Device) と GPU→CPU (Device-to-Host) のデータコピー処理を，時間ループのイテレーションの外側に移動します。前者はイテレーションの前，後者はイテレーションの後で使用します。

OpenACC

fem1d_cusp_csr.f90
subroutine fem1d_to_device(n, nnz)
    integer, intent(in) :: n, nnz
    real(real64) :: csrValA(nnz)
    integer :: csrRowPtrA(n+1), csrColIndA(nnz)
    real(real64) :: bvec(n), xvec(n)

    !$acc enter data copyin(xvec) &
    !$acc            create(csrValA, csrRowPtrA, csrColIndA, bvec)

end subroutine

subroutine fem1d_to_host(n, nnz)
    integer, intent(in) :: n, nnz
    real(real64) :: csrValA(nnz)
    integer :: csrRowPtrA(n+1), csrColIndA(nnz)
    real(real64) :: bvec(n), xvec(n)

    !$acc exit data copyout(xvec) &
    !$acc           delete(csrValA, csrRowPtrA, csrColIndA, bvec)

end subroutine

OpenMP版

fem1d_cusp_csr.f90
subroutine fem1d_to_device(n, nnz)
    integer, intent(in) :: n, nnz
    real(real64) :: csrValA(nnz)
    integer :: csrRowPtrA(n+1), csrColIndA(nnz)
    real(real64) :: bvec(n), xvec(n)

    !$omp target enter data map (to:xvec) &
    !$omp                   map (alloc:csrValA, csrRowPtrA, csrColIndA, bvec)

end subroutine

subroutine fem1d_to_host(n, nnz)
    integer, intent(in) :: n, nnz
    real(real64) :: csrValA(nnz)
    integer :: csrRowPtrA(n+1), csrColIndA(nnz)
    real(real64) :: bvec(n), xvec(n)

    !$omp target exit data map (from:xvec) &
    !$omp                  map (release:csrValA, csrRowPtrA, csrColIndA, bvec)

end subroutine

左辺係数行列と右辺ベクトルの計算処理を GPU に移します。

OpenACC版

fem1d_cusp_csr.f90
subroutine assemble_matrix(n, nnz, h, csrValA, csrRowPtrA, csrColIndA)
    integer, intent(in) :: n, nnz
    real(real64), intent(in) :: h
    real(real64), intent(inout) :: csrValA(nnz)
    integer, intent(inout) :: csrRowPtrA(n+1), csrColIndA(nnz)
    real(real64), parameter :: one6th = 1.0_real64 / 6.0_real64
    real(real64), parameter :: two3rd = 2.0_real64 / 3.0_real64

    !$acc parallel present(csrValA, csrRowPtrA, csrColIndA)
    i = 1
    ind = 1
    csrValA(ind)   =  1.0_real64
    csrValA(ind+1) = -1.0_real64
    csrColIndA(ind)   = i
    csrColIndA(ind+1) = i+1
    csrRowPtrA(i) = ind
    !$acc end parallel

    !$acc parallel present(csrValA, csrRowPtrA, csrColIndA)
    !$acc loop
    do i = 2, n-1
      ind = 3 * (i - 1)
      csrValA(ind)   = one6th * h
      csrValA(ind+1) = two3rd * h
      csrValA(ind+2) = one6th * h
      csrColIndA(ind)   = i-1
      csrColIndA(ind+1) = i
      csrColIndA(ind+2) = i+1
      csrRowPtrA(i) = ind
    end do
    !$acc end parallel

    !$acc parallel present(csrValA, csrRowPtrA, csrColIndA)
    i = n
    ind = 3 * (n - 1)
    csrValA(ind)   = -1.0_real64
    csrValA(ind+1) =  1.0_real64
    csrColIndA(ind)   = i-1
    csrColIndA(ind+1) = i
    csrRowPtrA(i)   = ind
    csrRowPtrA(i+1) = ind + 2
    !$acc end parallel

end subroutine

subroutine assemble_rhs(n, h, kappa, dt, u, bvec)
    integer, intent(in) :: n
    real(real64), intent(in) :: h, kappa, dt
    real(real64), intent(in) :: u(n)
    real(real64), intent(inout) :: bvec(n)
    real(real64), parameter :: one6th = 1.0_real64 / 6.0_real64
    real(real64), parameter :: two3rd = 2.0_real64 / 3.0_real64

    !$acc parallel present(bvec, u)
    !$acc loop
    do i = 2, n-1
      bvec(i) = (one6th * h +              kappa * dt / h) * u(i-1) &
              + (two3rd * h - 2.0_real64 * kappa * dt / h) * u(i)   &
              + (one6th * h +              kappa * dt / h) * u(i+1)
    end do
    !$acc end parallel

    !$acc parallel present(bvec)
    bvec(1) = 0.0_real64
    bvec(n) = 0.0_real64
    !$acc end parallel

end subroutine

OpenMP版

fem1d_cusp_csr.f90
subroutine assemble_matrix(n, nnz, h, csrValA, csrRowPtrA, csrColIndA)
    integer, intent(in) :: n, nnz
    real(real64), intent(in) :: h
    real(real64), intent(inout) :: csrValA(nnz)
    integer, intent(inout) :: csrRowPtrA(n+1), csrColIndA(nnz)
    real(real64), parameter :: one6th = 1.0_real64 / 6.0_real64
    real(real64), parameter :: two3rd = 2.0_real64 / 3.0_real64

    !$omp target teams map (alloc:csrValA, csrRowPtrA, csrColIndA)
    i = 1
    ind = 1
    csrValA(ind)   =  1.0_real64
    csrValA(ind+1) = -1.0_real64
    csrColIndA(ind)   = i
    csrColIndA(ind+1) = i+1 
    csrRowPtrA(i) = ind
    !$omp end target teams

    !$omp target teams map (alloc:csrValA, csrRowPtrA, csrColIndA)
    !$omp distribute parallel do
    do i = 2, n-1
      ind = 3 * (i - 1)
      csrValA(ind)   = one6th * h
      csrValA(ind+1) = two3rd * h
      csrValA(ind+2) = one6th * h
      csrColIndA(ind)   = i-1
      csrColIndA(ind+1) = i
      csrColIndA(ind+2) = i+1
      csrRowPtrA(i) = ind
    end do
    !$omp end distribute parallel do
    !$omp end target teams

    !$omp target teams map (alloc:csrValA, csrRowPtrA, csrColIndA)
    i = n
    ind = 3 * (n - 1)
    csrValA(ind)   = -1.0_real64
    csrValA(ind+1) =  1.0_real64
    csrColIndA(ind)   = i-1
    csrColIndA(ind+1) = i
    csrRowPtrA(i)   = ind
    csrRowPtrA(i+1) = ind + 2
    !$omp end target teams

end subroutine

subroutine assemble_rhs(n, h, kappa, dt, u, bvec)
    integer, intent(in) :: n
    real(real64), intent(in) :: h, kappa, dt
    real(real64), intent(in) :: u(n)
    real(real64), intent(inout) :: bvec(n)
    real(real64), parameter :: one6th = 1.0_real64 / 6.0_real64
    real(real64), parameter :: two3rd = 2.0_real64 / 3.0_real64

    !$omp target teams map (alloc:bvec, u)
    !$omp distribute parallel do
    do i = 2, n-1
      bvec(i) = (one6th * h +              kappa * dt / h) * u(i-1) &
              + (two3rd * h - 2.0_real64 * kappa * dt / h) * u(i)   &
              + (one6th * h +              kappa * dt / h) * u(i+1)
    end do
    !$omp end distribute parallel do
    !$omp end target teams

    !$omp target teams map (alloc:bvec)
    bvec(1) = 0.0_real64
    bvec(n) = 0.0_real64
    !$omp end target teams

end subroutine

cuSOLVER の cusolverSpDcsrlsvqr 関数を使用して方程式を解きます。上に示したように，係数行列や右辺ベクトル等は GPU メモリ上に移してあるのでコピーの必要がありません。

OpenACC版

fem1d_cusp_csr.f90
subroutine solve_cusolver_sp(n, nnz, csrValA, csrRowPtrA, csrColIndA, bvec, xvec)
    use openacc
    use m_cusolversp_interface
    integer, intent(in) :: n
    integer, intent(in) :: nnz
    real(real64), intent(in) :: csrValA(nnz)
    integer, intent(in) :: csrRowPtrA(n+1)
    integer, intent(in) :: csrColIndA(nnz)
    real(real64), intent(in) :: bvec(n)
    real(real64), intent(inout) :: xvec(n)
    integer :: istat
    integer(c_int) :: reorder, singularity
    real(c_double) :: tol
    type(cusparseMatDescr) :: descrA
    type(cusolverSpHandle), save :: handle
    logical, save :: is_initialized = .false.

    tol = 1.0e-12_real64
    reorder = 0

    if (.not. is_initialized) then
      istat = cusolverSpCreate(handle)
      is_initialized = .true.
    end if

    istat = cusparseCreateMatDescr(descrA)
    istat = cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL)
    istat = cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ONE)
    
    !$acc data present(csrValA, csrRowPtrA, csrColIndA, bvec, xvec)

    !$acc host_data use_device(csrValA, csrRowPtrA, csrColIndA, bvec, xvec)
    istat = cusolverSpDcsrlsvqr(handle, n, nnz, &
                                descrA, csrValA, csrRowPtrA, csrColIndA, bvec, &
                                tol, reorder, xvec, singularity)
    !$acc end host_data

    !$acc end data

end subroutine

OpenMP版

fem1d_cusp_csr.f90
subroutine solve_cusolver_sp(n, nnz, csrValA, csrRowPtrA, csrColIndA, bvec, xvec)
    use m_cusolversp_interface
    integer, intent(in) :: n
    integer, intent(in) :: nnz
    real(real64), intent(in) :: csrValA(nnz)
    integer, intent(in) :: csrRowPtrA(n+1)
    integer, intent(in) :: csrColIndA(nnz)
    real(real64), intent(in) :: bvec(n)
    real(real64), intent(inout) :: xvec(n)
    integer :: istat
    integer(c_int) :: reorder, singularity
    real(c_double) :: tol
    type(cusparseMatDescr) :: descrA
    type(cusolverSpHandle), save :: handle
    logical, save :: is_initialized = .false.

    tol = 1.0e-12_real64
    reorder = 0

    if (.not. is_initialized) then
      istat = cusolverSpCreate(handle)
      is_initialized = .true.
    end if

    istat = cusparseCreateMatDescr(descrA)
    istat = cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL)
    istat = cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ONE)

    !$omp target data map (alloc:csrValA, csrRowPtrA, csrColIndA, bvec, xvec)

    !$omp target data use_device_ptr(csrValA, csrRowPtrA, csrColIndA, bvec, xvec)
    istat = cusolverSpDcsrlsvqr(handle, n, nnz, &
                                descrA, csrValA, csrRowPtrA, csrColIndA, bvec, &
                                tol, reorder, xvec, singularity)
    !$omp end target data

    !$omp end target data

end subroutine

ベンチマーク

メッシュ数 4096，計算ステップ数 100 の条件で，改めて LAPACK 使用コードと cuSOLVER 使用コードの速度比較を行いました。結果は前回とほとんど変わりません。データの規模が小さいためと思われます。

ソルバ デバイス 処理時間 (秒)
GBTRF/GBTRS（帯行列 LAPACK） CPU 2.9
GTSV2（三重対角行列 cuSOLVER） GPU w/ OpenACC 0.13
GTSV2（三重対角行列 cuSOLVER） GPU w/ OpenMP 0.13
LSVQR（疎行列 cuSOLVER） GPU w/ OpenACC 5.2
LSVQR（疎行列 cuSOLVER） GPU w/ OpenMP 5.2

メッシュ数（nx）を多くして比較してみます。

ソルバ デバイス nx=4,096 nx=8,192 nx=16,384
GBTRF/GBTRS（帯行列 LAPACK） CPU 2.9 11.6 45.6
GTSV2（三重対角行列 cuSOLVER） GPU 0.13 0.17 0.24
LSVQR（疎行列 cuSOLVER） GPU 5.2 10.2 20.3

大規模計算では GPU を使用した cuSOLVER の方が速くなります。特に，三重対角行列ソルバは非常に高速です。疎行列ソルバはそれほどでもありませんが，LAPACK の帯行列ソルバがメッシュ数に対して 2 次のオーダーで依存するのに対して，cuSOLVER の疎行列ソルバが 1 次のオーダーでしかないことには注目すべきでしょう。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0