Skip to content
Snippets Groups Projects
solver_cholmod.hpp 12 KiB
Newer Older
  • Learn to ignore specific revisions
  • #pragma once
    
    
    Jakub Homola's avatar
    Jakub Homola committed
    #include "solver.h"
    
    Jakub Homola's avatar
    Jakub Homola committed
    #include <suitesparse/cholmod.h>
    
    Jakub Homola's avatar
    Jakub Homola committed
    const char * solver_str = "cholmod";
    
    
    
    
    
    
    template<typename T, typename I>
    
    Jakub Homola's avatar
    Jakub Homola committed
    struct direct_sparse_solver_data
    
        cholmod_common * cm_common = nullptr;
    
        cholmod_factor * cm_factor_super = nullptr;
        cholmod_factor * cm_factor_simpl = nullptr;
    
        cholmod_sparse * cm_Kreg_view = nullptr;
    
        MatrixDense<I,I> map_simpl_super;
    
    Jakub Homola's avatar
    Jakub Homola committed
    template<typename I>
    static void my_cholmod_start(cholmod_common * common)
    {
        if constexpr(std::is_same_v<I,int32_t>) cholmod_start(common);
        if constexpr(std::is_same_v<I,int64_t>) cholmod_l_start(common);
    }
    template<typename I>
    static cholmod_factor * my_cholmod_analyze(cholmod_sparse * A, cholmod_common * common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_analyze(A, common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_analyze(A, common);
    }
    template<typename I>
    static void my_cholmod_factorize(cholmod_sparse * A, cholmod_factor * F, cholmod_common * common)
    {
        if constexpr(std::is_same_v<I,int32_t>) cholmod_factorize(A, F, common);
        if constexpr(std::is_same_v<I,int64_t>) cholmod_l_factorize(A, F, common);
    }
    template<typename I>
    static cholmod_sparse * my_cholmod_factor_to_sparse(cholmod_factor * F, cholmod_common * common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_factor_to_sparse(F, common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_factor_to_sparse(F, common);
    }
    template<typename I>
    static cholmod_dense * my_cholmod_solve(int sys, cholmod_factor * L, cholmod_dense * B, cholmod_common * common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_solve(sys, L, B, common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_solve(sys, L, B, common);
    }
    template<typename I>
    static void my_cholmod_free_sparse(cholmod_sparse ** A, cholmod_common * common)
    {
        if constexpr(std::is_same_v<I,int32_t>) cholmod_free_sparse(A, common);
        if constexpr(std::is_same_v<I,int64_t>) cholmod_l_free_sparse(A, common);
    }
    template<typename I>
    static void my_cholmod_free_factor(cholmod_factor ** F, cholmod_common * common)
    {
        if constexpr(std::is_same_v<I,int32_t>) cholmod_free_factor(F, common);
        if constexpr(std::is_same_v<I,int64_t>) cholmod_l_free_factor(F, common);
    }
    template<typename I>
    static int my_cholmod_sdmult(cholmod_sparse * A, int transpose, double alpha[2], double beta[2], cholmod_dense * X, cholmod_dense * Y, cholmod_common * Common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_sdmult(A, transpose, alpha, beta, X, Y, Common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_sdmult(A, transpose, alpha, beta, X, Y, Common);
    }
    template<typename I>
    static cholmod_dense * my_cholmod_sparse_to_dense(cholmod_sparse * A, cholmod_common * Common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_sparse_to_dense(A, Common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_sparse_to_dense(A, Common);
    }
    template<typename I>
    static cholmod_dense * my_cholmod_allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *Common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_allocate_dense(nrow, ncol, d, xtype, Common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_allocate_dense(nrow, ncol, d, xtype, Common);
    }
    template<typename I>
    static int my_cholmod_free_dense(cholmod_dense **X, cholmod_common *Common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_free_dense(X, Common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_free_dense(X, Common);
    }
    template<typename I>
    static cholmod_factor * my_cholmod_copy_factor(cholmod_factor * F, cholmod_common * Common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_copy_factor(F, Common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_copy_factor(F, Common);
    }
    template<typename I>
    static int my_cholmod_change_factor(int to_xtype, int to_ll, int to_super, int to_packed, int to_monotonic, cholmod_factor *L, cholmod_common *Common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_change_factor(to_xtype, to_ll, to_super, to_packed, to_monotonic, L, Common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_change_factor(to_xtype, to_ll, to_super, to_packed, to_monotonic, L, Common);
    }
    template<typename I>
    static int my_cholmod_resymbol(cholmod_sparse *A, I *fset, size_t fsize, int pack, cholmod_factor *L, cholmod_common *Common)
    {
        if constexpr(std::is_same_v<I,int32_t>) return cholmod_resymbol(A, fset, fsize, pack, L, Common);
        if constexpr(std::is_same_v<I,int64_t>) return cholmod_l_resymbol(A, fset, fsize, pack, L, Common);
    }
    template<typename T>
    static int my_cholmod_dtype()
    {
        if constexpr(std::is_same_v<T,double>) return CHOLMOD_DOUBLE;
        if constexpr(std::is_same_v<T,float>)  return CHOLMOD_SINGLE;
    }
    template<typename I>
    static int my_cholmod_itype()
    {
        if constexpr(std::is_same_v<I,int32_t>) return CHOLMOD_INT;
        if constexpr(std::is_same_v<I,int64_t>) return CHOLMOD_LONG;
    }
    
    
    
    
    
    
    Jakub Homola's avatar
    Jakub Homola committed
    template<typename T, typename I>
    char direct_sparse_solver<T,I>::get_factors_symmetry()
    {
        return 'U';
    }
    
    
    
    
    template<typename T, typename I>
    
    Jakub Homola's avatar
    Jakub Homola committed
    direct_sparse_solver<T,I>::direct_sparse_solver()
    
    template<typename T, typename I>
    direct_sparse_solver<T,I>::~direct_sparse_solver()
    {
    
    Jakub Homola's avatar
    Jakub Homola committed
        destroy();
    
    template<typename T, typename I>
    
    Jakub Homola's avatar
    Jakub Homola committed
    void direct_sparse_solver<T,I>::init(const char * magicstring)
    
    Jakub Homola's avatar
    Jakub Homola committed
        {
    
        if(stage != 0) throw std::runtime_error("init: invalid order of operations in solver");
    
    Jakub Homola's avatar
    Jakub Homola committed
        data = std::make_unique<direct_sparse_solver_data<T,I>>();
    
    Jakub Homola's avatar
    Jakub Homola committed
    
        check_magicstringsolver_length(magicstring);
    
    Jakub Homola's avatar
    Jakub Homola committed
        data->zerodrop = mss_get_zerodrop(magicstring, true);
    
    Jakub Homola's avatar
    Jakub Homola committed
        data->cm_common = new cholmod_common();
        my_cholmod_start<I>(data->cm_common);
        data->cm_common->final_ll = 1;
        data->cm_common->nthreads_max = 1;
        data->cm_common->nmethods = 1;
        data->cm_common->method[0].ordering = CHOLMOD_METIS;
        data->cm_common->itype = my_cholmod_itype<I>();
        data->cm_common->supernodal = CHOLMOD_SUPERNODAL;
    
    
        stage = 1;
    
    Jakub Homola's avatar
    Jakub Homola committed
        }
    
    }
    
    
    
    template<typename T, typename I>
    
    Jakub Homola's avatar
    Jakub Homola committed
    template<typename A>
    
    void direct_sparse_solver<T,I>::set_matrix(const MatrixCSR<T,I,A> & Kreg)
    {
    
        if(stage < 1) throw std::runtime_error("set_matrix: invalid order of operations in solver");
    
    Jakub Homola's avatar
    Jakub Homola committed
        if(data->cm_Kreg_view == nullptr)
    
    Jakub Homola's avatar
    Jakub Homola committed
            data->cm_Kreg_view = new cholmod_sparse();
            data->cm_Kreg_view->nrow = Kreg.ncols;
            data->cm_Kreg_view->ncol = Kreg.nrows;
            data->cm_Kreg_view->nzmax = Kreg.nvals;
            data->cm_Kreg_view->nz = nullptr;
            data->cm_Kreg_view->z = nullptr;
            data->cm_Kreg_view->stype = -1;
            data->cm_Kreg_view->itype = my_cholmod_itype<I>();
            data->cm_Kreg_view->xtype = CHOLMOD_REAL;
            data->cm_Kreg_view->dtype = my_cholmod_dtype<T>();
            data->cm_Kreg_view->sorted = 1;
            data->cm_Kreg_view->packed = 1;
    
    Jakub Homola's avatar
    Jakub Homola committed
        data->cm_Kreg_view->p = Kreg.rowptrs;
        data->cm_Kreg_view->i = Kreg.colidxs;
        data->cm_Kreg_view->x = Kreg.vals;
    
        if(stage == 1) stage = 2;
        if(stage == 4) stage = 3;
    
    template<typename T, typename I>
    MatrixCSR<T,I> direct_sparse_solver<T,I>::get_matrix_view()
    {
        if(stage < 2) throw std::runtime_error("get_matrix_view: invalid order of operations in solver");
    
        MatrixCSR<T,I> view(data->cm_Kreg_view->nrow, data->cm_Kreg_view->ncol, data->cm_Kreg_view->nzmax, false);
        view.rowptrs = reinterpret_cast<I*>(data->cm_Kreg_view->p);
        view.colidxs = reinterpret_cast<I*>(data->cm_Kreg_view->i);
        view.vals = reinterpret_cast<T*>(data->cm_Kreg_view->x);
        return view;
    }
    
    
    
    
    template<typename T, typename I>
    void direct_sparse_solver<T,I>::factorize_symbolic()
    {
    
        if(stage != 2) throw std::runtime_error("factorize_symbolic: invalid order of operations in solver");
    
    Jakub Homola's avatar
    Jakub Homola committed
        data->cm_factor_super = my_cholmod_analyze<I>(data->cm_Kreg_view, data->cm_common);
    
        if(data->cm_factor_super->xsize > get_max_val_no_precision_loss_in_fp<T>()) throw std::runtime_error("factorize_symbolic: factor nnz too large for my super->simpl map");
    
    Jakub Homola's avatar
    Jakub Homola committed
        data->cm_factor_simpl = my_cholmod_copy_factor<I>(data->cm_factor_super, data->cm_common);
    
        my_cholmod_change_factor<I>(CHOLMOD_REAL, 1, 1, 1, 1, data->cm_factor_simpl, data->cm_common);
    
    Jakub Homola's avatar
    Jakub Homola committed
        for(size_t i = 0; i < data->cm_factor_simpl->xsize; i++) reinterpret_cast<T*>(data->cm_factor_simpl->x)[i] = static_cast<T>(i);
        my_cholmod_change_factor<I>(CHOLMOD_REAL, 1, 0, 1, 1, data->cm_factor_simpl, data->cm_common);
        if(data->zerodrop == 'D') my_cholmod_resymbol<I>(data->cm_Kreg_view, nullptr, 0, 1, data->cm_factor_simpl, data->cm_common);
        data->map_simpl_super.resize(data->cm_factor_simpl->nzmax, 1, -1, true);
        for(I i = 0; i < data->map_simpl_super.nrows; i++) data->map_simpl_super.vals[i] = static_cast<I>(reinterpret_cast<T*>(data->cm_factor_simpl->x)[i]);
    
        
        stage = 3;
    
    }
    
    
    
    template<typename T, typename I>
    void direct_sparse_solver<T,I>::factorize_numeric()
    {
    
        if(stage < 3) throw std::runtime_error("factorize_numeric: invalid order of operations in solver");
    
        my_cholmod_factorize<I>(data->cm_Kreg_view, data->cm_factor_super, data->cm_common);
    
        stage = 4;
    
    }
    
    
    
    template<typename T, typename I>
    
    Jakub Homola's avatar
    Jakub Homola committed
    template<typename A>
    void direct_sparse_solver<T,I>::get_factor_L(MatrixCSR<T,I,A> & /*L*/, bool /*copy_pattern*/, bool /*copy_values*/) const
    {
    
        throw std::runtime_error("L factor is not provided");
    
    Jakub Homola's avatar
    Jakub Homola committed
    }
    
    
    
    template<typename T, typename I>
    template<typename A>
    void direct_sparse_solver<T,I>::get_factor_U(MatrixCSR<T,I,A> & U, bool copy_pattern, bool copy_values) const
    
        if(stage < 3) throw std::runtime_error("get_factor_U: invalid order of operations in solver");
        if(copy_values && stage < 4) throw std::runtime_error("get_factor_U: invalid order of operations in solver");
        if((size_t)U.nrows != data->cm_factor_simpl->n || (size_t)U.ncols != data->cm_factor_simpl->n) throw std::runtime_error("get_factor_U: output matrix has wrong dimensions");
        if((size_t)U.nvals != data->cm_factor_simpl->nzmax) throw std::runtime_error("get_factor_U: output matrix has wrong nnz");
    
    Jakub Homola's avatar
    Jakub Homola committed
        if(copy_pattern) std::copy_n(static_cast<I*>(data->cm_factor_simpl->p), data->cm_factor_simpl->n+1, U.rowptrs);
        if(copy_pattern) std::copy_n(static_cast<I*>(data->cm_factor_simpl->i), data->cm_factor_simpl->nzmax, U.colidxs);
        if(copy_values) for(I i = 0; i < data->map_simpl_super.nrows; i++) U.vals[i] = reinterpret_cast<T*>(data->cm_factor_super->x)[data->map_simpl_super.vals[i]];
    
    }
    
    
    
    template<typename T, typename I>
    
    Jakub Homola's avatar
    Jakub Homola committed
    template<typename A>
    
    void direct_sparse_solver<T,I>::get_permutation(Permutation<I,A> & perm) const
    {
    
        if(stage < 3) throw std::runtime_error("get_permutation: invalid order of operations in solver");
    
    Jakub Homola's avatar
    Jakub Homola committed
        if((size_t)perm.size != data->cm_factor_simpl->n) MY_ABORT("get_permutation: output permutation has wrong size");
    
    Jakub Homola's avatar
    Jakub Homola committed
        std::copy_n(static_cast<I*>(data->cm_factor_simpl->Perm), data->cm_factor_simpl->n, perm.map_dst_to_src);
    
    Jakub Homola's avatar
    Jakub Homola committed
        if(data->cm_factor_simpl->IPerm != nullptr) std::copy_n(static_cast<I*>(data->cm_factor_simpl->IPerm), data->cm_factor_simpl->n, perm.map_dst_to_src);
    
    Jakub Homola's avatar
    Jakub Homola committed
        else for(I i = 0; i < perm.size; i++) perm.map_src_to_dst[perm.map_dst_to_src[i]] = i;
    
    }
    
    
    
    template<typename T, typename I>
    I direct_sparse_solver<T,I>::get_nnz_factor() const
    {
    
        if(stage < 3) throw std::runtime_error("get_nnz_factor: invalid order of operations in solver");
    
    
        // https://github.com/DrTimothyAldenDavis/SuiteSparse/issues/523
    
    
    Jakub Homola's avatar
    Jakub Homola committed
        return data->cm_factor_simpl->nzmax;
    }
    
    
    
    template<typename T, typename I>
    void direct_sparse_solver<T,I>::destroy()
    {
        if(stage == 0) return;
    
        delete data->cm_Kreg_view;
        data->cm_Kreg_view = nullptr;
    
        if(data->cm_factor_super != nullptr) my_cholmod_free_factor<I>(&data->cm_factor_super, data->cm_common);
        data->cm_factor_super = nullptr;
    
        if(data->cm_factor_super != nullptr) my_cholmod_free_factor<I>(&data->cm_factor_simpl, data->cm_common);
        data->cm_factor_simpl = nullptr;
    
        if(data->cm_common != nullptr) cholmod_finish(data->cm_common);
        delete data->cm_common;
        data->cm_common = nullptr;
    
        data->map_simpl_super.deallocate();
    
        stage = 0;