使用cusolverDnDpotrfBatched计算Cholesky分解下三角无预期零值如何解决

现象说明 #

你观察到的下三角区域外没有出现预期零值的现象完全正常，属于cusolverDnDpotrfBatched接口的原生设计行为，并非代码逻辑错误。

底层原因 #

cuSolver的Cholesky分解系列接口完全兼容BLAS/LAPACK的通用约定：

• 当你通过uplo参数指定CUBLAS_FILL_MODE_LOWER（下三角模式）时，接口只会计算并覆写输入矩阵的下三角元素（包含主对角线），上三角区域的原始输入值会被完整保留，接口不会主动将其置零
• 反之如果指定CUBLAS_FILL_MODE_UPPER（上三角模式），则只会处理上三角区域，下三角区域的原始值保持不变
  该设计的核心目的是减少不必要的内存写入操作，降低计算开销，用户如果需要完整的三角矩阵，自行对不需要的区域做置零处理即可。

遗漏的操作/配置 #

你的代码主要存在三处可补充的点：

1. 缺少结果后处理：从设备拷贝分解结果回主机后，需要手动将上三角区域的元素置为0，才能得到标准的下三角Cholesky矩阵
2. 缺少错误校验逻辑：没有读取infoArray的返回值判断分解是否成功，infoArray[i]返回0代表第i个矩阵分解成功，返回大于0的值代表第i个矩阵的第infoArray[i]个主元为0，矩阵非正定无法完成分解
3. 缺少资源释放逻辑：代码中申请的cusolver句柄、设备侧内存都没有主动释放，会造成显存泄漏和资源占用

修正后的参考代码 #

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <cuda_runtime.h>
#include <cusolverDn.h>
#include <iostream>

void printMatrix(int m, int n, const double*A, int lda, const char* name)
{
    for(int row = 0 ; row < m ; row++){
        for(int col = 0 ; col < n ; col++){
            double Areg = A[row + col*lda];
            printf("%s(%d,%d) = %f\n", name, row+1, col+1, Areg);
        }
    }
}

int main(){
    cusolverDnHandle_t handle = NULL;
    cusolverDnCreate(&handle);

    const cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
    const int batchSize = 1;
    const int m = 3;
    const int lda = m;
    double A0[lda*m] = { 1.0, 2.0, 3.0, 2.0, 5.0, 5.0, 3.0, 5.0, 12.0 };
    int infoArray[batchSize];

    double L0[lda*m];
    double *Aarray[batchSize];
    double **d_Aarray = NULL;
    int *d_infoArray = NULL;

    for(int j = 0 ; j < batchSize ; j++){
        cudaMalloc((void**)&Aarray[j], sizeof(double) * lda * m);
    }
    cudaMalloc((void**)&d_infoArray, sizeof(int)*batchSize);
    cudaMalloc((void**)&d_Aarray, sizeof(double*) * batchSize);

    cudaMemcpy(Aarray[0], A0, sizeof(double) * lda * m, cudaMemcpyHostToDevice);
    cudaMemcpy(d_Aarray, Aarray, sizeof(double*)*batchSize, cudaMemcpyHostToDevice);

    cusolverDnDpotrfBatched(handle, uplo, m, d_Aarray, lda, d_infoArray, batchSize);
    cudaDeviceSynchronize();

    cudaMemcpy(infoArray, d_infoArray, sizeof(int)*batchSize, cudaMemcpyDeviceToHost);
    // 校验分解结果
    if(infoArray[0] != 0){
        printf("Cholesky分解失败，错误码：%d\n", infoArray[0]);
        return 1;
    }
    cudaMemcpy(L0, Aarray[0], sizeof(double) * lda * m, cudaMemcpyDeviceToHost);

    // 手动将上三角区域置零
    for(int row = 0; row < m; row++){
        for(int col = row + 1; col < m; col++){
            L0[row + col * lda] = 0.0;
        }
    }

    printMatrix(m, m, L0, lda, "L0");

    // 资源释放
    for(int j = 0; j < batchSize; j++){
        cudaFree(Aarray[j]);
    }
    cudaFree(d_Aarray);
    cudaFree(d_infoArray);
    cusolverDnDestroy(handle);

    return 0;
}

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