Linear solve using CHLOMOD in C

Question

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!

Answer 1

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

}