将 Eigen::SparseMatrix 转换为 cuSparse,反之亦然
Convert Eigen::SparseMatrix to cuSparse and vice versa
由于在线文档和示例很少,我无法弄清楚如何将 Eigen::SparseMatrix 转换为 cuSparse。对于密集矩阵,将 cublas 的 Eigen 转换为 CUDA 相当简单
Eigen::MatrixXd A = Eigen::MatrixXd::Identity(3,3);
double *d_A;
cudaMalloc(reinterpret_cast<void **>(&d_A), 3 * 3 * sizeof(double));
cudaMemcpy(d_A, A.data(), sizeof(double) * 3 * 3, cudaMemcpyHostToDevice);
// do cublas operations on d_A
如何对稀疏矩阵进行等价?
std::vector<Eigen::Triplet<double>> trip;
trip.emplace_back(0, 0, 1);
trip.emplace_back(1, 1, 1);
trip.emplace_back(2, 2, 1);
Eigen::SparseMatrix<double> A(3, 3);
A.setFromTriplets(trip.begin(), trip.end());
double *d_A;
// cudaMalloc?
// cudaMemcpy? some conversion?
// do cusparse operations
以防万一有人感兴趣,我弄明白了。棘手的部分是 Eigen 的稀疏矩阵是 CSC 格式,而 cuSparse 是 CSR 格式。幸运的是,只需将 CSC 转换为 CSR 即可完成转换。
void EigenSparseToCuSparseTranspose(
const Eigen::SparseMatrix<double> &mat, int *row, int *col, double *val)
{
const int num_non0 = mat.nonZeros();
const int num_outer = mat.cols() + 1;
cudaMemcpy(row,
mat.outerIndexPtr(),
sizeof(int) * num_outer,
cudaMemcpyHostToDevice);
cudaMemcpy(
col, mat.innerIndexPtr(), sizeof(int) * num_non0, cudaMemcpyHostToDevice);
cudaMemcpy(
val, mat.valuePtr(), sizeof(double) * num_non0, cudaMemcpyHostToDevice);
}
void CuSparseTransposeToEigenSparse(
const int *row,
const int *col,
const double *val,
const int num_non0,
const int mat_row,
const int mat_col,
Eigen::SparseMatrix<double> &mat)
{
std::vector<int> outer(mat_col + 1);
std::vector<int> inner(num_non0);
std::vector<double> value(num_non0);
cudaMemcpy(
outer.data(), row, sizeof(int) * (mat_col + 1), cudaMemcpyDeviceToHost);
cudaMemcpy(inner.data(), col, sizeof(int) * num_non0, cudaMemcpyDeviceToHost);
cudaMemcpy(
value.data(), val, sizeof(double) * num_non0, cudaMemcpyDeviceToHost);
Eigen::Map<Eigen::SparseMatrix<double>> mat_map(
mat_row, mat_col, num_non0, outer.data(), inner.data(), value.data());
mat = mat_map.eval();
}
由于在线文档和示例很少,我无法弄清楚如何将 Eigen::SparseMatrix 转换为 cuSparse。对于密集矩阵,将 cublas 的 Eigen 转换为 CUDA 相当简单
Eigen::MatrixXd A = Eigen::MatrixXd::Identity(3,3);
double *d_A;
cudaMalloc(reinterpret_cast<void **>(&d_A), 3 * 3 * sizeof(double));
cudaMemcpy(d_A, A.data(), sizeof(double) * 3 * 3, cudaMemcpyHostToDevice);
// do cublas operations on d_A
如何对稀疏矩阵进行等价?
std::vector<Eigen::Triplet<double>> trip;
trip.emplace_back(0, 0, 1);
trip.emplace_back(1, 1, 1);
trip.emplace_back(2, 2, 1);
Eigen::SparseMatrix<double> A(3, 3);
A.setFromTriplets(trip.begin(), trip.end());
double *d_A;
// cudaMalloc?
// cudaMemcpy? some conversion?
// do cusparse operations
以防万一有人感兴趣,我弄明白了。棘手的部分是 Eigen 的稀疏矩阵是 CSC 格式,而 cuSparse 是 CSR 格式。幸运的是,只需将 CSC 转换为 CSR 即可完成转换。
void EigenSparseToCuSparseTranspose(
const Eigen::SparseMatrix<double> &mat, int *row, int *col, double *val)
{
const int num_non0 = mat.nonZeros();
const int num_outer = mat.cols() + 1;
cudaMemcpy(row,
mat.outerIndexPtr(),
sizeof(int) * num_outer,
cudaMemcpyHostToDevice);
cudaMemcpy(
col, mat.innerIndexPtr(), sizeof(int) * num_non0, cudaMemcpyHostToDevice);
cudaMemcpy(
val, mat.valuePtr(), sizeof(double) * num_non0, cudaMemcpyHostToDevice);
}
void CuSparseTransposeToEigenSparse(
const int *row,
const int *col,
const double *val,
const int num_non0,
const int mat_row,
const int mat_col,
Eigen::SparseMatrix<double> &mat)
{
std::vector<int> outer(mat_col + 1);
std::vector<int> inner(num_non0);
std::vector<double> value(num_non0);
cudaMemcpy(
outer.data(), row, sizeof(int) * (mat_col + 1), cudaMemcpyDeviceToHost);
cudaMemcpy(inner.data(), col, sizeof(int) * num_non0, cudaMemcpyDeviceToHost);
cudaMemcpy(
value.data(), val, sizeof(double) * num_non0, cudaMemcpyDeviceToHost);
Eigen::Map<Eigen::SparseMatrix<double>> mat_map(
mat_row, mat_col, num_non0, outer.data(), inner.data(), value.data());
mat = mat_map.eval();
}