Assign and print the results of CHOLMOD package

Question

I am trying to solve a simple working example, a linear system $Ax=b$, where $A$ is sparse SPD and $b$ is dense, using CHOLMOD.

#include <cmath>
#include <iostream>
#include <vector>
#include "suitesparse/cholmod.h"

int main()
{
    cholmod_sparse *a;
    cholmod_dense *x, *b;
    cholmod_factor *L;
    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 = sqrt(A.size());

    // Dense vector B

    b = cholmod_allocate_dense(n, 1, n, CHOLMOD_REAL, &c);
    b->x = &*B.begin();

    cholmod_print_dense(b, "b", &c);

    // Sparse matrix A

    cholmod_triplet *T;

    std::vector<int> Ti, Tj;
    std::vector<double> Tx;

    size_t index;
    double value;

    for (size_t i = 0; i < n; ++i)
    {
        for (size_t j = 0; j < n; ++j)
        {
            index = n * j + i;

            value = A.at(index);

            if (abs(value) > 0)
            {
                Ti.push_back(i);
                Tj.push_back(j);
                Tx.push_back(value);
            }
        }
    }

    T = cholmod_allocate_triplet(n, n, Tx.size(), 0, CHOLMOD_REAL, &c);
    T->i = &*Ti.begin();
    T->j = &*Tj.begin();
    T->x = &*Tx.begin();
    T->nnz = Tx.size();

    cholmod_print_triplet(T, "T", &c);

    a = cholmod_triplet_to_sparse(T, Tx.size(), &c);

    cholmod_print_sparse(a, "A", &c);

    L = cholmod_analyze(a, &c);
    cholmod_factorize(a, L, &c);
    x = cholmod_solve(CHOLMOD_A, L, b, &c);

    cholmod_free_factor(&L, &c);
    cholmod_free_sparse(&a, &c);
    cholmod_free_dense(&x, &c);
    cholmod_free_dense(&r, &c);
    cholmod_free_dense(&b, &c);

    cholmod_finish(&c);

    return 0;
}

Does someone know how to print the values of x = cholmod_solve(CHOLMOD_A, L, b, &c) and how to assign his values to a std::vector? As I do not know how to print and assign the values of the solution, I even know whether my code is working properly.

For sake of clarity, cholmod_dense is a Struct:

typedef struct cholmod_dense_struct
{
    size_t nrow ;   /* the matrix is nrow-by-ncol */
    size_t ncol ;
    size_t nzmax ;  /* maximum number of entries in the matrix */
    size_t d ;      /* leading dimension (d >= nrow must hold) */
    void *x ;       /* size nzmax or 2*nzmax, if present */
    void *z ;       /* size nzmax, if present */
    int xtype ;     /* pattern, real, complex, or zomplex */
    int dtype ;     /* x and z double or float */

} cholmod_dense ;

I would like to see the values of void *x of this struct.

Answer 1

Cholmod is a bit old-fashioned in that the data arrays are void*, which could be anything. However, you know that the type is double (or you could deduce it from xtype and dtype, but lets just assume you know its double).

Since you named the variable x, and cholmod just happens to also name its data member x, the data array is accessed by invoking static_cast<double*>(x->x), which is a bit awkward!

#include <algorithm> // in case you haven't already for std::copy_n
// the number of values you need to copy
int size = x->nrow*x->ncol; 
// create a vector with pre-allocated capacity
std::vector<double> x_vector(size);
// tell the compiler that x->x is an array of doubles, and copy the values into the vector's array
std::copy_n(static_cast<double*>(x->x),size,x_vector.data()); array

Or, just use a simple loop

int size = x->nrow*x->ncol; 
std::vector<double> x_vector(size);
for(int i = 0; i < size; i++)
  x_vector[i] = static_cast<double*>(x->x)[i];

Edit: Few things to consider. Doing things like b->x = &*B.begin() is risky: while you know that you are simply telling b->x to point to B's data array, cholmod doesnt know you're doing that. So when you call free, all hell breaks loose. For now, just copy everything. Its not efficient, but once you get better at c++ you can optimize all this stuff -- if its even necessary. Usually the math is far more expensive than cheap, fast copies like this.

Secondly, cholmod prefers if you only give it the upper or lower half of a symmetric matrix. Many sparse libraries behave this way. For cholmod, use stype == 1, and only entries where i>=j. You can experiment with other settings, but this seems to work.

int main()
{
    cholmod_sparse* a;
    cholmod_dense* x, * b;
    cholmod_factor* L;
    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 = sqrt(A.size()); // this seems needlessly complicated!

    // Dense vector B

    b = cholmod_allocate_dense(n, 1, n, CHOLMOD_REAL, &c);
    std::copy_n(B.begin(), n, static_cast<double*>(b->x));

    cholmod_print_dense(b, "b", &c);

    // Sparse matrix A

    cholmod_triplet* T;

    std::vector<int> Ti, Tj;
    std::vector<double> Tx;

    size_t index;
    double value;
    size_t nnz = 0;
  
    for (size_t i = 0; i < n; ++i)
    {
        for (size_t j = 0; j < n; ++j)
        {
            index = n * j + i;

            value = A.at(index);
            //only add "lower" triangle
            if (i >= j && abs(value) > 0)
            {
                Ti.push_back(i);
                Tj.push_back(j);
                Tx.push_back(value);
                nnz++;
            }
        }
    }
    //stype 1 refers to "lower" triangle
    T = cholmod_allocate_triplet(n, n, nnz , 1, CHOLMOD_REAL, &c);
    std::copy_n(Ti.begin(), nnz, static_cast<int*>(T->i));
    std::copy_n(Tj.begin(), nnz, static_cast<int*>(T->j));
    std::copy_n(Tx.begin(), nnz, static_cast<double*>(T->x));
    T->nnz = nnz;

    cholmod_print_triplet(T, "T", &c);

    a = cholmod_triplet_to_sparse(T, Tx.size(), &c);

    cholmod_print_sparse(a, "A", &c);
    cholmod_check_sparse(a, &c);

    L = cholmod_analyze(a, &c);
    cholmod_factorize(a, L, &c);
    x = cholmod_solve(CHOLMOD_A, L, b, &c);
    std::vector<double> sol(n);
    std::copy_n(static_cast<double*>(x->x), n, sol.data());
    for (int i = 0; i < n; i++)
        std::cout << sol[i] << "\n";

    cholmod_free_factor(&L, &c);
    cholmod_free_sparse(&a, &c);
    cholmod_free_dense(&x, &c);
    cholmod_free_dense(&b, &c);

    cholmod_finish(&c);

    return 0;
}

While there are easier sparse libraries to work with, cholmod is one of the best. Its not the prettiest or newest, but its algorithms are top notch. If you can stomach it, I'd stick with it.

For something easier, though considerably slower, try Eigen. I try to contribute to Eigen when I can, though I doubt I can ever integrate all the thought that went into cholmod.