I am using the open-source CHLOMOD (as here http://faculty.cse.tamu.edu/davis/suitesparse.html) in order to solve a linear system Ax=b (performing A/b=x) in my domain decomposition code but I am unsure how to use it in my C serial code on a regular CPU rather than a GPU. Could someone give me an example of how to import this into some example of code that solves for a system of equations? I am fairly new to C and am unsure how to do this. Thank you so much!
-
1$\begingroup$ Did you look at "Simple example program" on page 10 of the CHOLMOD User Guide? It is only 33 lines long so I doubt it someone can come up with a much simpler example. What is it you didn't understand about this example? Using CHOLMOD on a CPU is the more typical and simpler usage; using it on a GPU would be more challenging. $\endgroup$ – Bill Greene Jul 29 '16 at 18:03
-
$\begingroup$ Thank you taking the time to respond. I was trying to locate a guide of sorts, I read all the read me files but I could not find an example. I am sorry if this seems like a silly question, I can't seem to find it, where did you see this user guide? I do not see it within the tar file nor on the webpage. $\endgroup$ – user20973 Jul 29 '16 at 18:09
-
$\begingroup$ The doc is in SuiteSparse/CHOLMOD/Doc/UserGuide.pdf after SuiteSparse is built. I also found copies of it on the web using a simple search. $\endgroup$ – Bill Greene Jul 29 '16 at 18:20
$\begingroup$
$\endgroup$
Here is how I call CHOLMOD from a C program. It uses as input my own C++ structure for a sparse matrix (in standard Compressed Row Storage format). There are several gotchas: understanding the conventions for lower and upper triangular (stype parameter) was what took me some time (CHOLMOD silently fails if this is not correct).
#include <suitesparse/cholmod.h>
#include <suitesparse/cholmod_check.h>
namespace {
using namespace OGF;
/**
* \brief Solves a linear system with CHOLMOD.
* \param[in] M a reference to a matrix. The matrix is de-initialized
* on exit.
* \param[out] x_out a pointer to the solution vector
* \param[in] b_in a pointer to the right hand side
* \retval true if the linear system could be solved
* \retval false otherwise
*/
bool solve_linear_system_with_cholmod(
SparseMatrix& M, double* x_out, const double* b_in
) {
ogf_assert(M.m() == M.n());
index_t n = M.n();
NLSparseMatrix& MM = *M.implementation();
// Compute required nnz (works if MM in symmetric and
// MM is in full storage mode).
index_t nnz = 0;
for(index_t i=0; i<n; ++i) {
NLRowColumn& Ri = MM.row[i];
for(index_t jj=0; jj<Ri.size; ++jj) {
if(Ri.coeff[jj].index <= i) {
++nnz;
}
}
}
// Step 1: initialize CHOLMOD library
//----------------------------------------------------
static cholmod_common c ;
static bool cholmod_initialized = false;
if(!cholmod_initialized) {
cholmod_start(&c) ;
cholmod_initialized = true;
}
// Step 2: translate sparse matrix into cholmod format
//----------------------------------------------------
cholmod_sparse* A = cholmod_allocate_sparse(
n, n, nnz, // Dimensions and number of non-zeros
false, // Sorted = false
true, // Packed = true
1, // stype (-1 = lower triangular, 1 = upper triangular)
CHOLMOD_REAL, // Entries are real numbers
&c
);
int* colptr = (int*)A->p;
int* rowind = (int*)A->i;
double* val = (double*)A->x;
// Convert Geogram Matrix into CHOLMOD Matrix
index_t count = 0 ;
for(index_t i=0; i<n; ++i) {
colptr[i] = int(count);
NLRowColumn& Ri = MM.row[i];
for(index_t jj=0; jj<Ri.size; ++jj) {
const NLCoeff& C = Ri.coeff[jj];
index_t j = C.index;
if(j <= i) {
val[count] = C.value;
rowind[count] = int(j);
++count;
}
}
}
geo_assert(count == nnz);
colptr[n] = int(nnz);
/*
geo_assert(cholmod_check_sparse(A,&c) != 0);
if(n < 10) {
cholmod_write_sparse(stdout,A,NULL,NULL,&c);
}
cholmod_print_sparse(A,"A",&c);
*/
// Step 2: construct right-hand side
cholmod_dense* b = cholmod_allocate_dense(n, 1, n, CHOLMOD_REAL, &c) ;
Memory::copy(b->x, b_in, n * sizeof(double)) ;
// Step 3: factorize and solve
cholmod_factor* L = cholmod_analyze(A, &c) ;
geo_debug_assert(cholmod_check_factor(L,&c) != 0);
if(!cholmod_factorize(A, L, &c)) {
std::cerr << "COULD NOT FACTORIZE !!!" << std::endl;
}
cholmod_dense* x = cholmod_solve(CHOLMOD_A, L, b, &c) ;
Memory::copy(x_out, x->x, n * sizeof(double)) ;
// Step 4: cleanup
cholmod_free_factor(&L, &c) ;
cholmod_free_sparse(&A, &c) ;
cholmod_free_dense(&x, &c) ;
cholmod_free_dense(&b, &c) ;
// To be called at exit.
// cholmod_finish(&c) ;
// Commented-out sanity check
/*
vector<double> check(n, 0.0);
M.mult(x_out,check.data());
double err = 0.0;
for(index_t i=0; i<n; ++i) {
err += geo_sqr(b_in[i] - check[i]);
}
err = ::sqrt(err);
std::cerr << " CHOLMOD || Ax - b || = " << err << std::endl;
*/
return true;
}
}