0
$\begingroup$

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!

$\endgroup$
  • 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
2
$\begingroup$

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;
    }

}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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