在 C 中使用 CHLOMOD 进行线性求解

计算科学 矩阵 稀疏矩阵 C
2021-11-29 12:18:52

我正在使用开源 CHLOMOD(如这里http://faculty.cse.tamu.edu/davis/suitesparse.html)来解决我的域中的线性系统 Ax=b(执行 A/b=x)分解代码,但我不确定如何在常规 CPU 而不是 GPU 上的 C 串行代码中使用它。有人可以给我一个示例,说明如何将其导入到求解方程组的代码示例中吗?我对C相当陌生,不确定如何做到这一点。太感谢了!

1个回答

这是我从 C 程序中调用 CHOLMOD 的方法。它使用我自己的 C++ 结构作为稀疏矩阵的输入(采用标准压缩行存储格式)。有几个问题:理解下三角和上三角(stype 参数)的约定让我花了一些时间(如果这不正确,CHOLMOD 会默默地失败)。

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

}