Newer
Older
cholmod_factor * cm_factor_super = nullptr;
cholmod_factor * cm_factor_simpl = nullptr;
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
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;
}
template<typename T, typename I>
char direct_sparse_solver<T,I>::get_factors_symmetry()
{
return 'U';
}
direct_sparse_solver<T,I>::direct_sparse_solver()
template<typename T, typename I>
direct_sparse_solver<T,I>::~direct_sparse_solver()
{
void direct_sparse_solver<T,I>::init(const char * magicstring)
if(stage != 0) throw std::runtime_error("init: invalid order of operations in solver");
data = std::make_unique<direct_sparse_solver_data<T,I>>();
check_magicstringsolver_length(magicstring);
data->zerodrop = mss_get_zerodrop(magicstring, true);
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;
}
template<typename T, typename I>
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");
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;
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");
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");
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);
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]);
}
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);
}
template<typename T, typename I>
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");
}
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");
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>
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");
if((size_t)perm.size != data->cm_factor_simpl->n) MY_ABORT("get_permutation: output permutation has wrong size");
std::copy_n(static_cast<I*>(data->cm_factor_simpl->Perm), data->cm_factor_simpl->n, perm.map_dst_to_src);
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);
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
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
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;