Creating a sparse matrix in CHOLMOD or SuiteSparseQR

Question

In SparseSuiteQR, all of the examples I can find use stdin or a file read to create a sparse matrix. Could someone provide a simple example of how to create one directly in C++?

Even better, in the CHOLMOD documentation, there is mention of a sparse2 function available in matlab, which behaves the same as the sparse. Can this be used in C++?

Answer 1

The data structures used by SuiteSparseQR (e.g. cholmod_sparse) are defined in the CHOLMOD library. You can find more information about it on the CHOLMOD documentation, which is much larger than the one from SuiteSparseQR.

Answer 2

I am assuming that you try to solve a linear system, see the CSparse package from Tim Davies, or boost matrix libraries which also have numeric bindings which interface umfpack and some lapack functions AFAIK...

Answer 3

CHOLMOD is a pretty awesome project - thanks Tim Davis :)

There is surprisingly a lot of code on GitHub that makes use of CHOLMOD, but you have to be logged into GitHub and know what you're looking for!

So, after crawling through CHOLMOD documentation and source code and then searching through GitHub for source code that uses CHOLMOD you would find out what to do.

But for most developers who want/need a quick example, here it is below.

*Note that your mileage might vary depending on how you compiled SuiteSparse. (You might need to use the cholmod_ variant (without the l), i.e. not cholmod_l_; and use int for indexing, not long int).

// example.cpp
#include "SuiteSparseQR.hpp"
#include "SuiteSparse_config.h"

int main (int argc, char **argv)
{
  cholmod_common Common, *cc;
  cholmod_sparse *A;
  cholmod_dense *X, *B;

  // start CHOLMOD
  cc = &Common;
  cholmod_l_start (cc);

  /* A =
      [
        1.1,  0.0, -0.5,  0.7
        0.0, -2.0,  0.0,  0.0
        0.0,  0.0,  0.9,  0.0
        0.0,  0.0,  0.0,  0.6
      ]
  */
  int m = 4;   // num rows in A
  int n = 4;   // num cols in A
  int nnz = 6; // num non-zero elements in A
  int unsymmetric = 0; // A is non-symmetric: see cholmod.h > search for `stype` for more details

  // In coordinate form (COO) a.k.a. triplet form (zero-based indexing)
  int i[nnz] = {0, 1, 0, 2, 0, 3}; // row indices
  int j[nnz] = {0, 1, 2, 2, 3, 3}; // col indices
  double x[nnz] = {1.1, -2.0, -0.5, 0.9, 0.7, 0.6}; // values

  // Set up the cholmod matrix in COO/triplet form
  cholmod_triplet *T = cholmod_l_allocate_triplet(m, n, nnz, unsymmetric, CHOLMOD_REAL, cc);
  T->nnz = nnz;

  for (int ind = 0; ind < nnz; ind++)
  {
    ((long int *) T->i)[ind] = i[ind];     // Notes:
    ((long int *) T->j)[ind] = j[ind];     // (1) casting necessary because these are void* (see cholmod.h)
    ((double *) T->x)[ind] = x[ind];       // (2) direct assignment will cause memory corruption
  }                                        // (3) long int for index pointers corresponds to usage of cholmod_l_* functions

  // convert COO/triplet to CSC (compressed sparse column) format
  A = (cholmod_sparse *) cholmod_l_triplet_to_sparse(T, nnz, cc);
  // note: if you already know CSC format you can skip the triplet allocation and instead use cholmod_allocate_sparse
  //       and assign the member variables: see cholmod.h > cholmod_sparse_struct definition

  // B = ones (size (A,1),1)
  B = cholmod_l_ones (A->nrow, 1, A->xtype, cc);

  // X = A\B
  X = SuiteSparseQR <double> (A, B, cc);

  // Print contents of X
  printf("X = [\n");
  for (int ind = 0; ind < n; ind++)
  {
    printf("%f\n", ((double *) X->x)[ind]);
  }
  printf("]\n");
  fflush(stdout);

  // free everything and finish CHOLMOD
  cholmod_l_free_triplet (&T, cc);
  cholmod_l_free_sparse (&A, cc);
  cholmod_l_free_dense (&X, cc);
  cholmod_l_free_dense (&B, cc);
  cholmod_l_finish (cc);
  return 0;
}

Supposing you have compiled SuiteSparse successfully and you have saved example.cpp in the base directory, then the following should work (on Linux):

gcc example.cpp -I./include -L./lib -lcholmod -lspqr -lsuitesparseconfig -o example
#Add SuiteSpare libraries to your `ld` search path if necessary
LD_LIBRARY_PATH=$(pwd)/lib
export LD_LIBRARY_PATH
./example

Output:

X = [
0.353535
-0.500000
1.111111
1.666667
]