0
$\begingroup$

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.

$\endgroup$

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Thank you, it worked! But, unfortunately, the result is wrong. Do you know what I am doing wrong? $\endgroup$ Oct 31, 2022 at 19:46
  • $\begingroup$ Do you know other sparse solvers that are more new-fashioned? $\endgroup$ Oct 31, 2022 at 19:48
  • $\begingroup$ I really appreciated your edit, Charlie. Thank you so much! If I am not asking too much, what do you mean by saying "Usually the math is far more expensive than cheap, fast copies like this". I am working with the Finite Element Method, and big sparse matrices. The faster the better. $\endgroup$ Nov 3, 2022 at 11:46
  • $\begingroup$ "[...] but once you get better at c++ you can optimize all this stuff". Could you indicate some c++ resources so I can get better to optimize this (when necessary, should I say)? $\endgroup$ Nov 3, 2022 at 11:47
  • 1
    $\begingroup$ As you have experienced, trying to avoid a copy can lead to an incorrect result. Perhaps its avoidable, or perhaps you will spend hours/days trying to "fix" something that can't be fixed due to the way someone else designed their linear algebra library. I think its best to focus on getting correct results first. Then, if its absolutely paramount that you reduce the run time of your simulation by 0.001%, then you can replace copy_n with something more clever. $\endgroup$
    – Charlie S
    Nov 3, 2022 at 14:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.