为何cusolverDnDsytri在Fortran中无法计算对称矩阵逆？

在Fortran项目中使用CUDA求解非正定对称矩阵的逆，因矩阵非正定无法使用Cholesky分解，故采用cusolverDnDsytrf执行Bunch-Kaufman分解，再调用cusolverDnDsytri求逆。编译使用nvhpc 24.9版本的nvfortran，所有错误返回值均被捕获且无报错，分解步骤结果符合预期（已用scipy验证），但cusolverDnDsytri未按预期将矩阵更新为逆矩阵，复制上三角到下三角的步骤正确，但求逆后矩阵无变化。

最小示例代码：

program symmetric_inverse
      use cusolverDn
      use cublas
      integer n,lwork_factor,lwork_solve,error,col,row
      parameter(n=3)
      double precision A(n,n), A_inv(n,n), I(n,n), prnt(n,n)
      double precision, device :: dA_inv(n,n), dI(n,n), dA(n,n)
      double precision, allocatable, device :: dwork_factor(:)
      double precision, allocatable, device :: dwork_solve(:)
      integer, device :: dpivots(n)
      integer, device :: devinfo
      integer local_devinfo
      type(cusolverDnHandle) :: h

      A = reshape((/ -1, 2, 0, 2, -5, 0, 0, 0, -1 /), shape(A))

C     send A to device so A can be inverted
      dA_inv = A

C     get the cusolver handle
      error = cusolverDnCreate(h)
      if (error .ne. CUBLAS_STATUS_SUCCESS) then
            print*, 'Could not create handle', error
            stop
      endif

C     get the buffer size needed for LDL^T calculations
      error = cusolverDnDsytrf_buffersize(h,n,dA_inv,n,lwork_factor)
      if (error .ne. CUBLAS_STATUS_SUCCESS) then
            print*, 'Could not get get LDL^T buffer size', error
            stop
      endif

C     allocate the workspace
      allocate(dwork_factor(lwork_factor))

      prnt = dA_inv
      print*, "Before LDL"
      print*, prnt

C     Compute the LDL^T factorization
      error = cusolverDnDsytrf(h,CUBLAS_FILL_MODE_UPPER,
     &                         n,dA_inv,n,dpivots,dwork_factor,
     &                         lwork_factor,devinfo)
      local_devinfo = devinfo
      if ((error.ne.CUBLAS_STATUS_SUCCESS).or.(local_devinfo.ne.0)) then
            print*, 'Could not find the LDL^T factorization', error
            print*, 'SYTRF Info:', local_devinfo
            stop
      endif

      prnt = dA_inv
      print*, "After LDL"
      print*, prnt

C     get the buffer size needed for inverse calculations 
      error = cusolverDnDsytri_buffersize(h,CUBLAS_FILL_MODE_UPPER,n,
     &                                    dA_inv,n,dpivots,lwork_solve)
      if(error.ne.CUBLAS_STATUS_SUCCESS) then
            print*, 'Could not get inverse buffer size', error
            stop
      endif

      prnt = dA_inv
      print*, "Before INV"
      print*, prnt
 
C     Find the inverse of the matrix 
      allocate(dwork_solve(lwork_solve))
      error = cusolverDnDsytri(h,CUBLAS_FILL_MODE_UPPER,
     &                         n,dA_inv,n,dpivots,dwork_solve,
     &                         lwork_solve,devinfo)
      local_devinfo = devinfo
      if((error.ne.CUBLAS_STATUS_SUCCESS).or.(local_devinfo.ne.0)) then
            print*, 'Could not compute the inverse from LDL', error
            print*, 'SYTRI Info:', local_devinfo
            stop
      endif

      prnt = dA_inv
      print*, "After INV"
      print*, prnt

C     get the inverse from device
      A_inv = dA_inv

      print*, "Before Copy"
      print*, A_inv

C     fill in the lower triangle of A
C     since upper is used for results in previous blas calls
      do row=1, n
            do col=row+1, n
                  A_inv(col,row) = A_inv(row,col)
            enddo
      enddo

      print*, "After Copy"
      print*, A_inv

      deallocate(dwork_factor)
      deallocate(dwork_solve)
      error = cusolverDnDestroy(h)
      end 

解决建议 #

1. 显式同步设备操作 #

虽然未使用异步流，但nvfortran的隐式设备数据拷贝可能存在延迟，导致在设备求逆操作完成前就读取了旧数据。在调用cusolverDnDsytri之后、拷贝数据到主机之前，添加显式同步：

error = cusolverDnDsytri(h,CUBLAS_FILL_MODE_UPPER,
     &                         n,dA_inv,n,dpivots,dwork_solve,
     &                         lwork_solve,devinfo)
      local_devinfo = devinfo
      if((error.ne.CUBLAS_STATUS_SUCCESS).or.(local_devinfo.ne.0)) then
            print*, 'Could not compute the inverse from LDL', error
            print*, 'SYTRI Info:', local_devinfo
            stop
      endif

C     显式同步设备，确保求逆操作完成
      call cudaDeviceSynchronize()

      prnt = dA_inv
      print*, "After INV"
      print*, prnt

2. 验证dpivots数组的正确性 #

Bunch-Kaufman分解的dpivots是求逆的关键输入，确认分解后的dpivots没有被意外修改。可以在分解后将其拷贝到主机打印，与scipy的Bunch-Kaufman分解结果对比：

C     分解后拷贝dpivots到主机
      integer pivots_host(n)
      pivots_host = dpivots
      print*, "Pivots after SYTRF:"
      print*, pivots_host

3. 切换填充模式测试 #

尝试将分解和求逆的填充模式从CUBLAS_FILL_MODE_UPPER改为CUBLAS_FILL_MODE_LOWER，确认是否是填充模式的处理逻辑问题：

C     Compute the LDL^T factorization
      error = cusolverDnDsytrf(h,CUBLAS_FILL_MODE_LOWER,
     &                         n,dA_inv,n,dpivots,dwork_factor,
     &                         lwork_factor,devinfo)
      ...
C     get the buffer size needed for inverse calculations 
      error = cusolverDnDsytri_buffersize(h,CUBLAS_FILL_MODE_LOWER,n,
     &                                    dA_inv,n,dpivots,lwork_solve)
      ...
C     Find the inverse of the matrix 
      error = cusolverDnDsytri(h,CUBLAS_FILL_MODE_LOWER,
     &                         n,dA_inv,n,dpivots,dwork_solve,
     &                         lwork_solve,devinfo)

4. 验证矩阵乘积是否为单位矩阵 #

通过cublas计算原矩阵与结果矩阵的乘积，确认是否真的未完成求逆，还是打印/拷贝环节的问题：

C     计算原矩阵与逆矩阵的乘积，验证是否为单位矩阵
      double precision, device :: dA_copy(n,n), dProd(n,n)
      dA_copy = A  ! 保存原矩阵的设备副本
      call cublasDgemm(h, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, 1.0d0,
     &                 dA_copy, n, dA_inv, n, 0.0d0, dProd, n)
      call cudaDeviceSynchronize()
      double precision prod_host(n,n)
      prod_host = dProd
      print*, "Product of original and inverse matrix:"
      print*, prod_host

5. 排查编译器/库版本问题 #

尝试更新nvhpc到最新稳定版本，或者回退到更早的版本（如24.7），确认是否是当前版本的cusolver库存在bug。

内容的提问来源于stack exchange，提问作者Ethan