如何在C++中将CHOLMOD求解结果打印并赋值给std::vector

核心问题解决方法 #

CHOLMOD的cholmod_solve返回cholmod_dense*类型的结果，其内部数据通过x字段存储（类型为void*），需要先强制转换为对应数值类型指针（这里是double*，因为使用了CHOLMOD_REAL），再进行打印或转换。

1. 打印求解结果x

将x->x转换为double*后，按列优先顺序遍历元素（对于n行1列的向量，元素是连续存储的）：

// 假设x是cholmod_solve返回的cholmod_dense*指针
double* x_data = static_cast<double*>(x->x);
for (size_t i = 0; i < n; ++i) {
    // 列索引为0，行索引i，对应位置为 0 * x->d + i（因x->d = n，直接取x_data[i]即可）
    std::cout << "x[" << i << "] = " << x_data[i] << std::endl;
}

2. 将结果转换为std::vector<double>

利用vector的迭代器构造或循环赋值，直接复制连续存储的元素：

// 方法1：迭代器构造（最简洁）
std::vector<double> x_vec(static_cast<double*>(x->x), static_cast<double*>(x->x) + n);

// 方法2：循环赋值
std::vector<double> x_vec(n);
for (size_t i = 0; i < n; ++i) {
    x_vec[i] = static_cast<double*>(x->x)[i];
}

修正后的完整代码 #

原代码存在两处关键问题：直接赋值b->x会导致野指针风险、未填充稀疏矩阵的triplet结构，以下是修正后的完整代码：

#include <iostream>
#include <vector>
#include <cmath>
#include <cstring>
#include "cholmod.h"

int main()
{
    // 初始化CHOLMOD环境
    cholmod_common c;
    cholmod_start(&c);

    std::vector<double> A = {4.0, 0.0, 1.0,
                             0.0, 3.0, 1.0,
                             1.0, 1.0, 2.0};
    std::vector<double> B = {1.0, 1.0, 1.0};
    size_t n = static_cast<size_t>(sqrt(A.size()));    

    // 正确初始化稠密向量b：拷贝数据而非直接赋值指针
    cholmod_dense* b = cholmod_allocate_dense(n, 1, n, CHOLMOD_REAL, &c);
    std::memcpy(b->x, B.data(), n * sizeof(double));

    // 将稠密矩阵A转换为CHOLMOD triplet格式（存储非零元素）
    std::vector<int> Ti, Tj;
    std::vector<double> Tx;
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            double val = A[i * n + j];
            if (val != 0.0) {
                Ti.push_back(static_cast<int>(i));
                Tj.push_back(static_cast<int>(j));
                Tx.push_back(val);
            }
        }
    }
    cholmod_triplet* T = cholmod_allocate_triplet(n, n, Tx.size(), CHOLMOD_REAL, CHOLMOD_SYMMETRIC, &c);
    std::memcpy(T->i, Ti.data(), Ti.size() * sizeof(int));
    std::memcpy(T->j, Tj.data(), Tj.size() * sizeof(int));
    std::memcpy(T->x, Tx.data(), Tx.size() * sizeof(double));
    T->nnz = static_cast<int>(Tx.size());

    // 转换为稀疏矩阵并进行Cholesky分解
    cholmod_sparse* a = cholmod_triplet_to_sparse(T, static_cast<int>(Tx.size()), &c);
    cholmod_factor* L = cholmod_analyze(a, &c);
    cholmod_factorize(a, L, &c);

    // 求解Ax=b
    cholmod_dense* x = cholmod_solve(CHOLMOD_A, L, b, &c);

    // 打印求解结果
    std::cout << "求解结果x：" << std::endl;
    double* x_data = static_cast<double*>(x->x);
    for (size_t i = 0; i < n; ++i) {
        std::cout << "x[" << i << "] = " << x_data[i] << std::endl;
    }

    // 转换为std::vector<double>
    std::vector<double> x_vec(static_cast<double*>(x->x), static_cast<double*>(x->x) + n);
    std::cout << "\n转换后的std::vector：" << std::endl;
    for (size_t i = 0; i < x_vec.size(); ++i) {
        std::cout << "x_vec[" << i << "] = " << x_vec[i] << std::endl;
    }

    // 释放CHOLMOD资源
    cholmod_free_dense(&x, &c);
    cholmod_free_factor(&L, &c);
    cholmod_free_sparse(&a, &c);
    cholmod_free_triplet(&T, &c);
    cholmod_free_dense(&b, &c);
    cholmod_finish(&c);

    return 0;
}

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