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

前々回，LAPACK を使用して計算するプログラムを作成しました。これを cuSOLVER を使用するコードに書き換えて，GPU オフロード処理を実現してみようと思います。

使用する関数は，係数行列が三重対角行列のときに使用可能な GTSV2 と，疎行列として取り扱いながら解く LSVQR とします。

GTSV2

まずは，GTSV2（三重対角行列用の関数）を使用するプログラムからコーディングしていきます。OpenACC と OpenMP の両方で書いてみます。オフロードするデータは，左辺行列を構成する dl, d, du と右辺ベクトルの bvec です。

変更のメインは，LAPACK の dgetrf と dgetrs の両関数を，cuSPARSE の cusparseDgtsv2 に置き換えることです。そのために，cusparseCreate と cusparseDgtsv2_bufferSizeExt で事前に関数の使用準備をしています。

OpenACC版

fem1d_cusp_tridiag_naive.f90

subroutine solve_cusparse(n, dl, d, du, bvec) 
    use openacc
    use cusparse 
    integer, intent(in) :: n
    real(real64), intent(inout) :: dl(n), d(n), du(n)
    real(real64), intent(inout) :: bvec(n)
    integer :: nrhs, lda, ldb, istat, lwork, d_info
    integer(c_size_t) :: buffsize
    integer(c_intptr_t), allocatable :: buff(:)
    type(cusparseHandle), save :: handle
    logical, save :: is_initialized = .false.
    
    nrhs = 1
    ldb = n

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

    !$acc data copyin(dl, d, du) &
    !$acc      copy(bvec) &
    !$acc      create(buff)

    !$acc host_data use_device(dl, d, du, bvec)
    istat = cusparseDgtsv2_bufferSizeExt(handle, n, nrhs, dl, d, du, bvec, ldb, buffsize)
    !$acc end host_data

    if (.not. allocated(buff)) allocate(buff(buffsize))

    !$acc enter data create(buff)

    !$acc host_data use_device(dl, d, du, bvec)
    istat = cusparseDgtsv2(handle, n, nrhs, dl, d, du, bvec, ldb, buff)
    !$acc end host_data

    !$acc exit data delete(buff)

    !$acc end data

end subroutine

OpenMP版

fem1d_cusp_tridiag_naive.f90

subroutine solve_cusparse(n, dl, d, du, bvec)
    use cusparse 
    integer, intent(in) :: n
    real(real64), intent(inout) :: dl(n), d(n), du(n)
    real(real64), intent(inout) :: bvec(n)
    integer :: nrhs, lda, ldb, istat, lwork, d_info
    integer(c_size_t) :: buffsize
    integer(c_intptr_t), allocatable :: buff(:)
    type(cusparseHandle), save :: handle
    logical, save :: is_initialized = .false.
    
    nrhs = 1
    ldb = n

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

    !$omp target data map (to:dl, d, du) &
    !$omp             map (tofrom:bvec) &
    !$omp             map (alloc:buff)

    !$omp target data use_device_ptr(dl, d, du, bvec)
    istat = cusparseDgtsv2_bufferSizeExt(handle, n, nrhs, dl, d, du, bvec, ldb, buffsize)
    !$omp end target data

    if (.not. allocated(buff)) allocate(buff(buffsize))

    !$omp target enter data map (alloc:buff)

    !$omp target data use_device_ptr(dl, d, du, bvec)
    istat = cusparseDgtsv2(handle, n, nrhs, dl, d, du, bvec, ldb, buff)
    !$omp end target data

    !$omp target exit data map (release:buff)

    !$omp end target data

end subroutine

LSVQR

LSVQR（疎行列用の関数）を使用するプログラムも同様にします。
こちらも変更のメインは，LAPACK の dgetrf / dgetrs を，cuSOLVER の cusolverSpDcsrlsvqr に置き換えることです。左辺係数行列は所定の疎行列の形式で csrValA, csrRowPtrA, csrColIndA に与えられているものとします。

OpenACC版

fem1d_cusp_csr_naive.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(out) :: 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 copyin(csrValA, csrRowPtrA, csrColIndA, bvec) &
    !$acc      copyout(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_naive.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(out) :: 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 (to:csrValA, csrRowPtrA, csrColIndA, bvec) &
    !$omp             map (from: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.3
LSVQR（疎行列 cuSOLVER）	GPU w/ OpenMP	5.3

コンパイラは NVIDIA HPC SDK 23.5 を使用しています。オプションは次のとおりとしました。また，CPU/GPU はそれぞれ Intel Xeon W-1270，NVIDIA T400 を使用しています。

	コンパイルオプション	リンクオプション
OpenACC	-O2 -acc=gpu -gpu=cc75 -Minfo=accel -Mpreprocess	-cudalib= cusolver,cusparse
OpenMP	-O2 -mp=gpu -gpu=cc75 -Minfo=accel,mp -Mpreprocess	-cudalib= cusolver,cusparse

前回 CMake の話をしましたが，バージョン 3.24 ではまだ，nvfortran の OpenACC/OpenMP offloading に対して対応が途上のようです。コンパイルでエラーが出るか、そうでなかったとしても実行時にエラーが出て計算ができません。仕方がないので，Makefile を使ってコンパイルしました。

まとめ

LAPACK を使用して解く有限要素法コードを GPU にオフロードして解きました。この段階では LAPACK の関数を cuSOLVER の関数に置き換えただけなので，関数を呼び出すたびに CPU-GPU 間のデータ転送をしており無駄の多いコードとなっています。次回，その改善に取り組みます。