使用cusolverDnDpotrfBatched得到不正确的结果



我想使用cusolverDnDpotrfBatched找到3x3矩阵的Cholesky分解,但我没有得到应该存在于较低三角形矩阵中的零。这是我要计算cholesky分解的矩阵[1 2 3;2 5 5;[5]。应该是这样吗?我错过了什么?我知道这个帖子基于CUDA的Cholesky分解我的代码如下:

#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) = %fn", 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 nrhs = 1;
const int m = 3;
const int lda = m;
//const int ldb = 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]; /* host copy of error info */
double L0[lda*m]; /* cholesky factor of A0 */
double *Aarray[batchSize];
//double *Barray[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);
//assert(cudaSuccess == cudaStat1);
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);
cudaMemcpy(L0, Aarray[0], sizeof(double) * lda * m, cudaMemcpyDeviceToHost);
for(int i =0; i<9;i++)std::cout<<L0[i]<<std::endl;
//printMatrix(m, m, L0, lda, "L0");
//printf("=====n");
}

我没有得到应该存在于下三角矩阵中的零

也许您希望阅读文档:

如果输入参数uplo为CUBLAS_FILL_MODE_LOWER,则只处理A的下三角部分,并将其替换为下三角choolesky因子l

备注:A的另一部分用作工作区。例如,uplo为CUBLAS_FILL_MODE_UPPER,则A的上三角形包含choolesky因子U, potrfBatched后A的下三角形被破坏。

所以不期望矩阵的另一部分是0

相关内容

  • 没有找到相关文章

最新更新