小灰在整理电脑文档的时候，有个就工程里用fortran写了点调用cuda的部分，就用一个用例子做了记录分享。cuda编程是Fortran中通过一个 C/CUDA 胶水层调用 CUDA Runtime、cuBLAS 与 cuSOLVER，在 GPU 上求解广义对称正定特征值问题 A x = λ B x（LAPACK DSYGVD 的 GPU 版本）。示例代码位于 test_dsygvd.F90 与 gpudriver.cu。

目录

  • 方法概览
  • 代码结构与数据流
  • 关键接口说明（含 Fortran bind(C) 示例）
  • 构建与运行
  • 常见问题与建议
  • 扩展方向

方法概览

  • 在 C/CUDA 侧以 extern "C" 暴露函数，并使用 F90(name)=name_ 适配 Fortran 默认的尾随下划线名约定。
  • Fortran 侧调用这些接口完成：设备初始化、内存分配与主机–设备拷贝、调用 cuSOLVER 的 cusolverDnDsygvd* 求解、回传结果并释放资源。
  • 字符选项（jobz/uplo）由 Fortran 以 ICHAR('v'/'u') 传递到 C 侧，C 侧以 ASCII 值比较后映射到 cuSOLVER/cuBLAS 的枚举。

代码结构与数据流

Fortran 驱动（test_dsygvd.F90）

  • 初始化参数与矩阵（uplo='u', jobz='v', itype=1, n=3）；
  • 初始化 CUDA 与库句柄：cudalib_init, set_gpu, init_test_gpus；
  • 设备内存分配：cuda_malloc(dev_a/dev_b/dev_w)；
  • 主机→设备拷贝：cuda_memcpy(a→dev_a), cuda_memcpy(b→dev_b)；
  • 求解：dsygvd_gpu(itype,uplo,jobz,n,dev_a,dev_b,dev_w)；
  • 设备→主机拷回结果：cuda_memcpy(dev_a→a), cuda_memcpy(dev_b→b), cuda_memcpy(dev_w→w)；
  • 释放设备内存：cuda_free(dev_a/dev_b/dev_w)。

CUDA 胶水层（gpudriver.cu）

  • 句柄管理：cudalib_init 创建 cublasHandle_t 和 cusolverDnHandle_t，cudalib_free 销毁；
  • 内存管理：cuda_malloc/cuda_free，数据拷贝：cuda_memcpy（flag=1: H2D, flag=2: D2H）；
  • 求解封装：dsygvd_gpu 完成工作区查询 cusolverDnDsygvd_bufferSize、分配临时缓冲、调用 cusolverDnDsygvd、同步并释放临时内存；
  • 调试辅助：cuda_print 通过设备核打印数组前若干元素。

关键接口说明（含 Fortran bind(C) 示例）

当前示例使用隐式接口（依赖符号名与默认传参约定）。为更稳妥的跨语言边界，这里给出等价的显式 bind(C) 接口示例，便于在实际项目中采用。

说明：示例基于双精度（double），指针以 type(c_ptr) 表示；字符选项以 c_char 或 c_int 传递，两种写法都提供。

C/CUDA 侧（简要节选）

// gpudriver.cu（节选）
extern "C" void cudalib_init(int *dev_id);
extern "C" void cudalib_free(void);
extern "C" void set_gpu(int *dev_id);
extern "C" void init_test_gpus(int *myid);

extern "C" void cuda_malloc(double **ptr, int *length, int *elemSize);
extern "C" void cuda_memcpy(double *host_ptr, double **device_ptr, int *length,
                            int *elemSize, int *flag);
extern "C" void cuda_free(long *ptr);

extern "C" void dsygvd_gpu(int *iitype, int *jobz_t, int *uplo_t, int *len,
                           double **dev_A, double **dev_B, double **dev_eigs);

Fortran 显式接口（推荐做法）

! fortran_bindings.f90（示例片段）
use iso_c_binding
implicit none

interface
  subroutine cudalib_init(dev_id) bind(C,name="cudalib_init")
    import :: c_int
    integer(c_int), value :: dev_id
  end subroutine

  subroutine cudalib_free() bind(C,name="cudalib_free")
  end subroutine

  subroutine set_gpu(dev_id) bind(C,name="set_gpu")
    import :: c_int
    integer(c_int), value :: dev_id
  end subroutine

  subroutine init_test_gpus(myid) bind(C,name="init_test_gpus")
    import :: c_int
    integer(c_int), value :: myid
  end subroutine

  subroutine cuda_malloc(ptr, length, elemSize) bind(C,name="cuda_malloc")
    import :: c_int, c_double, c_ptr
    type(c_ptr) :: ptr
    integer(c_int), value :: length, elemSize
  end subroutine

  subroutine cuda_memcpy(host_ptr, device_ptr, length, elemSize, flag) &
    bind(C,name="cuda_memcpy")
    import :: c_int, c_double, c_ptr
    real(c_double) :: host_ptr(*)
    type(c_ptr)    :: device_ptr
    integer(c_int), value :: length, elemSize, flag  ! 1:H2D, 2:D2H
  end subroutine

  subroutine cuda_free(ptr) bind(C,name="cuda_free")
    import :: c_long
    integer(c_long), value :: ptr   ! 与 C 侧 long* 对齐
  end subroutine

  subroutine dsygvd_gpu(itype, jobz_t, uplo_t, n, dev_A, dev_B, dev_W) &
    bind(C,name="dsygvd_gpu")
    import :: c_int, c_ptr
    integer(c_int), value :: itype, jobz_t, uplo_t, n
    type(c_ptr) :: dev_A, dev_B, dev_W
  end subroutine
end interface

传递字符选项的两种方式：

  1. 整数 ASCII：与现有代码一致（jobz_t = ichar('V') 等）；
  2. 字符绑定：定义 character(c_char) :: jobz,uplo，在 C 侧以 char 接收。

使用示例（Fortran 主程序片段）

program test
  use iso_c_binding
  implicit none
  integer(c_int) :: n, itype
  integer(c_int) :: CPU2GPU, GPU2CPU
  integer(c_int) :: jobz_t, uplo_t
  type(c_ptr)    :: dev_A, dev_B, dev_W
  real(c_double) :: A(3,3), B(3,3), W(3)

  CPU2GPU = 1; GPU2CPU = 2
  n = 3; itype = 1
  jobz_t = ichar('V'); uplo_t = ichar('U')

  call cudalib_init(0_c_int)
  call set_gpu(0_c_int)
  call init_test_gpus(0_c_int)

  call cuda_malloc(dev_A, n*n, int(c_sizeof(0.0d0),c_int))
  call cuda_malloc(dev_B, n*n, int(c_sizeof(0.0d0),c_int))
  call cuda_malloc(dev_W, n,   int(c_sizeof(0.0d0),c_int))

  ! 初始化 A/B ...
  call cuda_memcpy(A, dev_A, n*n, int(c_sizeof(0.0d0),c_int), CPU2GPU)
  call cuda_memcpy(B, dev_B, n*n, int(c_sizeof(0.0d0),c_int), CPU2GPU)

  call dsygvd_gpu(itype, jobz_t, uplo_t, n, dev_A, dev_B, dev_W)

  call cuda_memcpy(W, dev_W, n, int(c_sizeof(0.0d0),c_int), GPU2CPU)

  ! ... 使用 W
end program

构建与运行

  1. 编译 CUDA 胶水层
nvcc -O2 -Xcompiler -fPIC -c gpudriver.cu -o gpudriver.o
  1. 链接 Fortran 与 CUDA 库
gfortran -O2 test_dsygvd.F90 gpudriver.o \
  -L/usr/local/cuda/lib64 -lcudart -lcublas -lcusolver \
  -o test_dsygvd

如需自定义 include 或库路径，请根据安装位置调整 -I 与 -L。

  1. 运行前设置库搜索路径（按需）
export LD_LIBRARY_PATH=/usr/local/cuda/lib64:${LD_LIBRARY_PATH}
  1. 运行
./test_dsygvd