cuSolver 未返回正确的解
cuSolver doesn't return the correct solution
我正在尝试在cuSOLVER中使用QR线性系统求解器,这个
#include <cusparse_v2.h>
#include <stdio.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include "device_launch_parameters.h"
#include <iostream>
#include <cassert>
#include <cublas_v2.h>
#include <cusolverDn.h>
#include <cusolverSp.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
using namespace thrust;
#define CUSPARSE_CHECK(x) {cusparseStatus_t _c=x; if (_c != CUSPARSE_STATUS_SUCCESS) {printf("cusparse fail: %d, line: %dn", (int)_c, __LINE__); exit(-1);}}
void init_handlers_and_matr_descriptor(cusolverSpHandle_t& cusolverH,cusparseHandle_t& cusparseH,cusparseMatDescr_t& descrA) {
assert(CUSOLVER_STATUS_SUCCESS == cusolverSpCreate(&cusolverH));
assert(CUSPARSE_STATUS_SUCCESS == cusparseCreate(&cusparseH));
assert(CUSPARSE_STATUS_SUCCESS == cusparseCreateMatDescr(&descrA)); //CUSPARSE_INDEX_ZERO,CUSPARSE_MATRIX_TYPE_GENERAL
assert(CUSPARSE_STATUS_SUCCESS == cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ZERO));
assert(CUSPARSE_STATUS_SUCCESS == cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL));
}
int sparse_solver_test() {
//Init csr format A and b for solving Ax = b
/*
A =
[ 1.0 2.0 0.0 0.0 0.0 0.0 ]
[ 3.0 4.0 5.0 0.0 0.0 0.0 ]
[ 0.0 6.0 7.0 8.0 0.0 0.0 ]
[ 0.0 0.0 9.0 10.0 11.0 0.0 ]
[ 0.0 0.0 0.0 12.0 13.0 14.0 ]
[ 0.0 0.0 0.0 0.0 15.0 16.0 ]
b = [0.0,2.0,4.0,6.0,8.0,10.0]^T
*/
int nnz_A = 16, m = 6;
cusolverSpHandle_t cusolverH = NULL;
cusparseHandle_t cusparseH = NULL; //cuBLAS or cuSPARSE?
cusolverStatus_t cusolver_status = CUSOLVER_STATUS_SUCCESS;
cusparseMatDescr_t descrA;
cudaError_t cudaStat;
init_handlers_and_matr_descriptor(cusolverH, cusparseH, descrA);
host_vector<double> h_csrValA(nnz_A), h_dataB(m), h_dataX(m);
host_vector<int> h_csrRowPtrA(m + 1), h_csrColIndA(nnz_A);
device_vector<double> d_csrValA(nnz_A), d_dataB(m), d_dataX(m);
device_vector<int> d_csrRowPtrA(m + 1), d_csrColIndA(nnz_A);
//Init matrix
for (auto i = 0; i < nnz_A; ++i) h_csrValA[i] = static_cast<double>(i + 1);
h_csrRowPtrA[0] = 0;
h_csrRowPtrA[1] = 2;
for (auto i = 2; i < m; ++i) h_csrRowPtrA[i] = h_csrRowPtrA[i - 1] + 3;
h_csrRowPtrA[m] = nnz_A;
h_csrColIndA[0] = 0;
int v[] = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5 };
for (auto i = 0; i < nnz_A; ++i) h_csrColIndA[i] = v[i];
for (auto i = 0; i < m; ++i) h_dataB[i] = static_cast<double>(2 * i);
//device memory and descriptor A init
d_csrValA = h_csrValA;
d_csrRowPtrA = h_csrRowPtrA;
d_csrColIndA = h_csrColIndA;
d_dataB = h_dataB;
//step4, solve the linear system?
int singularity;
cusolver_status = cusolverSpDcsrlsvqr(
cusolverH, m, nnz_A, descrA,
d_csrValA.data().get(), d_csrRowPtrA.data().get(), d_csrColIndA.data().get(), d_dataB.data().get(),
0.0, 0, d_dataX.data().get(), &singularity);
std::cout << "singularity = " << singularity << std::endl;
assert(CUSOLVER_STATUS_SUCCESS == cusolver_status);
h_dataX = d_dataX;
std::cout << "x = (";
for (auto i = 0; i < m; ++i) {
std::cout << h_dataX[i];
if (i < m - 1) std::cout << ", ";
}
std::cout << std::endl;
if (cusparseH)
cusparseDestroy(cusparseH);
if (cusolverH)
cusolverSpDestroy(cusolverH);
if (descrA)
cusparseDestroyMatDescr(descrA);
cudaDeviceReset();
return 0;
}
int main(int argc, char** argv) {
sparse_solver_test();
return 0;
}
不确定我的函数设置是否错误,有人可以帮忙吗?
更新 我使用推力库简化了代码,虽然错误仍然相同,但至少我摆脱了所有的 malloc 等......
更新 按照建议更正了csrIndColA
(相应地更改了代码(数组。现在求解器工作了(即我不再得到我之前得到的错误(,尽管我得到的结果仍然是 0。
更新 随着我所做的所有更改,我也忘记了初始化h_dataB
,以及解决问题的csrIndColA
索引,完整的代码在上面供将来参考。
示例中
的csrColIndA
数组太短,因此 cuSOLVER 会尝试读取它的末尾。
根据 cuSOLVER 文档和通用约定,列索引数组与非零矩阵条目数组具有相同的长度,并存储每个非零元素的列索引(而不仅仅是每列中的第一个非零元素,如您的示例所示,这会将格式限制为所有非零元素垂直连续的稀疏模式(。
所以你的示例输出应该有
csrColIndA = {0, 1, 0, 1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5}
相关文章:
- 来自 std::list 的迭代器 .end() 按预期返回"0xcdcdcdcdcdcdcdcd"但 .begin()
- 什么时候在C++中返回常量引用是个好主意
- 你能重载对象变量名本身返回的内容吗
- 为什么 Serial.println(<char[]>);返回随机字符?
- C++映射:具有自定义类的运算符[]不起作用(总是返回0)
- 如何获取std::result_of函数的返回类型
- QueryWorkingSet总是返回false
- (C++)分析树以计算返回错误值的简单算术表达式
- 访问者访问变体并返回不同类型时出错
- 如何返回一个类的两个对象相加的结果
- OpenInventor从9.8升级到10.4.2后,GLSL纹理返回零
- lower_bound()返回最后一个元素
- "throw expression code" 1e7 >返回 d 是什么?投掷标准::overflow_error( "too big" ) : d;意味 着?
- 奇怪的(对我来说)返回声明 - 在谷歌上找不到任何关于它的信息
- 如何取消对nullptr的屏蔽,返回正确的对象
- 奇怪的结构&GCC&clang(void*返回类型)
- 架构决策:返回std::future还是提供回调
- 从python中调用C++函数并获取返回值
- 矩阵向量乘法(cublasDgemv)返回零
- cuSolver 未返回正确的解