如何在 Fortran 中将 OpenACC 与 cublasDgetrfBatched 接口？

Posted 2023-03-16

【中文标题】如何在 Fortran 中将 OpenACC 与 cublasDgetrfBatched 接口？

【英文标题】：How to interface OpenACC with cublasDgetrfBatched in Fortran?

【发布时间】：2014-05-24 13:33:59

【问题描述】：

我一直在研究 Fortran 代码，它使用 cuBLAS 批处理 LU 和 cuSPARSE 批处理三对角求解器作为带有 ADI 预条件器的 BiCG 迭代求解器的一部分。我使用的是计算能力为 3.5 和 CUDA 5.5 的 Kepler K20X。我在没有 PGI 的 CUDA Fortran 的情况下这样做，所以我正在编写自己的接口：

FUNCTION cublasDgetrfBatched(handle, n, dA, ldda, dP, dInfo, nbatch) BIND(C, NAME="cublasDgetrfBatched")
    USE, INTRINSIC :: ISO_C_BINDING
    INTEGER(KIND(CUBLAS_STATUS_SUCCESS)) :: cublasDgetrfBatched
    TYPE(C_PTR), VALUE :: handle
    INTEGER(C_INT), VALUE :: n
    TYPE(C_PTR), VALUE :: dA
    INTEGER(C_INT), VALUE :: ldda
    TYPE(C_PTR), VALUE :: dP
    TYPE(C_PTR), VALUE :: dInfo
    INTEGER(C_INT), VALUE :: nbatch
END FUNCTION cublasDgetrfBatched

我使用 cudaHostAlloc 在主机上分配固定内存，为矩阵和包含指向矩阵的设备指针的设备数组分配设备内存，将每个矩阵异步复制到设备，执行操作，然后异步复制分解的矩阵并返回到主机以使用单个右侧执行反向替换：

REAL(8), POINTER, DIMENSION(:,:,:) :: A
INTEGER, DIMENSION(:,:), POINTER :: ipiv
TYPE(C_PTR) :: cPtr_A, cPtr_ipiv
TYPE(C_PTR), ALLOCATABLE, DIMENSION(:), TARGET :: dPtr_A
TYPE(C_PTR) :: dPtr_ipiv, dPtr_A_d, dPtr_info
INTEGER(C_SIZE_T) :: sizeof_A, sizeof_ipiv

...

stat = cudaHostAlloc(cPtr_A, sizeof_A, cudaHostAllocDefault)
CALL C_F_POINTER(cPtr_A, A, (/m,m,nbatch/))
stat = cudaHostAlloc(cPtr_ipiv, sizeof_ipiv, cudaHostAllocDefault)
CALL C_F_POINTER(cPtr_ipiv, ipiv, (/m,nbatch/))

ALLOCATE(dPtr_A(nbatch))
DO ibatch=1,nbatch
  stat = cudaMalloc(dPtr_A(ibatch), m*m*sizeof_double)
END DO
stat = cudaMalloc(dPtr_A_d, nbatch*sizeof_cptr)
stat = cublasSetVector(nbatch, sizeof_cptr, C_LOC(dPtr_A(1)), 1, dPtr_A_d, 1)
stat = cudaMalloc(dPtr_ipiv, m*nbatch*sizeof_cint)
stat = cudaMalloc(dPtr_info, nbatch*sizeof_cint)

...

!$OMP PARALLEL DEFAULT(shared) PRIVATE( stat, ibatch )
!$OMP DO
DO ibatch = 1,nbatch
  stat = cublasSetMatrixAsync(m, m, sizeof_double, C_LOC(A(1,1,ibatch)), m, dPtr_A(ibatch), m, mystream)
END DO
!$OMP END DO
!$OMP END PARALLEL

...

stat = cublasDgetrfBatched(cublas_handle, m, dPtr_A_d, m, dPtr_ipiv, dPtr_info, nbatch)

...

stat = cublasGetMatrixAsync(m, nbatch, sizeof_cint, dPtr_ipiv, m, C_LOC(ipiv(1,1)), m, mystream)

!$OMP PARALLEL DEFAULT(shared) PRIVATE( ibatch, stat )
!$OMP DO
DO ibatch = 1,nbatch
  stat = cublasGetMatrixAsync(m, m, sizeof_double, dPtr_A(ibatch), m, C_LOC(A(1,1,ibatch)), m, mystream)
END DO
!$OMP END DO
!$OMP END PARALLEL

...

!$OMP PARALLEL DEFAULT(shared) PRIVATE( ibatch, x, stat )
!$OMP DO
DO ibatch = 1,nbatch
  x = rhs(:,ibatch)
  CALL dgetrs( 'N', m, 1, A(1,1,ibatch), m, ipiv(1,ibatch), x(1), m, info )
  rhs(:,ibatch) = x
END DO
!$OMP END DO
!$OMP END PARALLEL

...

我宁愿不必在最后一步执行此操作，但 cublasDtrsmBatched 例程将矩阵大小限制为 32，而我的大小为 80（批处理 Dtrsv 会更好，但这不存在）。启动多个单独的 cublasDtrsv 内核的成本使得在设备上执行 back-sub 变得难以为继。

在对 cublasDgetrfBatched 和 cusparseDgtsvStridedBatch 的调用之间，我还需要执行其他操作。其中大部分当前在主机上执行，OpenMP 用于在批处理级别并行化循环。一些操作，例如每个被分解的矩阵的矩阵向量乘法，正在使用 OpenACC 的设备上进行计算：

!$ACC DATA COPYIN(A) COPYIN(x) COPYOUT(Ax)

...

!$ACC KERNELS
  DO ibatch = 1,nbatch
    DO i = 1,m
      Ax(i,ibatch) = zero
    END DO
    DO j = 1,m
      DO i = 1,m
        Ax(i,ibatch) = Ax(i,ibatch) + A(i,j,ibatch)*x(j,ibatch)
      END DO
    END DO
  END DO
!$ACC END KERNELS

...

!$ACC END DATA

我想通过 OpenACC 将更多的计算放在 GPU 上，但要做到这一点，我需要能够将两者连接起来。类似于以下内容：

!$ACC DATA COPYIN(A) CREATE(info,A_d) COPYOUT(ipiv)

!$ACC HOST_DATA USE_DEVICE(A)
DO ibatch = 1,nbatch
  A_d(ibatch) = acc_deviceptr(A(1,1,ibatch))
END DO
!$ACC END HOST_DATA

...

!$ACC HOST_DATA USE_DEVICE(ipiv,info)
stat = cublasDgetrfBatched(cublas_handle, m, A_d, m, ipiv, info, nbatch)
!$ACC END HOST_DATA

...

!$ACC END DATA

我知道在大多数情况下，带有 host_device 子句的 host_data 构造是合适的，但由于我实际上需要向 cuBLAS 传递一个包含指向设备上矩阵的指针的设备数组，所以我不确定如何继续。

谁能提供任何见解？

谢谢

【参考方案1】：

！！把所有东西都放在设备上 !$ACC 数据复制（A）创建（信息，A_d）复制（ipiv）

！！填充设备 A_d 数组 !$ACC 并行循环 DO ibatch = 1,nbatch A_d(ibatch) = A(1,1,ibatch) 结束做 !$ACC 并行结束

...

！！将A_d的设备地址发送给设备 !$ACC HOST_DATA USE_DEVICE(A_d,ipiv,info) stat = cublasDgetrfBatched(cublas_handle, m, A_d, m, ipiv, info, nbatch) !$ACC END HOST_DATA

...

!$ACC 结束数据

---

                                       or

---

！！将除 A_d 之外的所有内容都放在设备上 !$ACC 数据复制（A）创建（信息）复制（ipiv）

！！填充主机 A_d 数组 DO ibatch = 1,nbatch A_d(ibatch) = acc_deviceptr( A(1,1,ibatch) ) 结束做

！！将 A_d 复制到设备 !$acc 数据复制（A_d） ...

！！将A_d等的设备地址发送给设备 !$ACC HOST_DATA USE_DEVICE(A_d,ipiv,info) stat = cublasDgetrfBatched(cublas_handle, m, A_d, m, ipiv, info, nbatch) !$ACC END HOST_DATA

... !$acc 结束数据

!$ACC 结束数据

【讨论】：

介意解释一下吗？