suanPan
SparseMatLis.hpp
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
29#ifndef SPARSEMATLIS_HPP
30#define SPARSEMATLIS_HPP
31
32#include <lis/lislib.h>
33#include "SparseMat.hpp"
34#include "csr_form.hpp"
35
36template<sp_d T> class SparseMatLis final : public SparseMat<T> {
37 LIS_SOLVER solver = nullptr;
38
39protected:
41 using SparseMat<T>::setting;
42
43 int direct_solve(Mat<T>&, const Mat<T>&) override;
44
45public:
46 SparseMatLis(const uword in_row, const uword in_col, const uword in_elem = 0)
47 : SparseMat<T>(in_row, in_col, in_elem) {
48 lis_solver_create(&solver);
49 lis_solver_set_option("-i fgmres", solver);
50 }
51
53 : SparseMat<T>(other) {
54 lis_solver_create(&solver);
55 lis_solver_set_option("-i fgmres", solver);
56 }
57
58 SparseMatLis(SparseMatLis&&) noexcept = delete;
59 SparseMatLis& operator=(const SparseMatLis&) = delete;
60 SparseMatLis& operator=(SparseMatLis&&) noexcept = delete;
61
62 ~SparseMatLis() override { lis_solver_destroy(solver); }
63
64 unique_ptr<MetaMat<T>> make_copy() override { return std::make_unique<SparseMatLis>(*this); }
65};
66
67template<sp_d T> int SparseMatLis<T>::direct_solve(Mat<T>& X, const Mat<T>& B) {
68 X.set_size(B.n_rows, B.n_cols);
69
70 csr_form<double, LIS_INT> csr_mat(this->triplet_mat, SparseBase::ZERO, true);
71
72 const sp_i auto n = csr_mat.n_rows;
73 const sp_i auto nnz = csr_mat.n_elem;
74
75 LIS_MATRIX A;
76 LIS_VECTOR b, x;
77
78 lis_matrix_create(0, &A);
79 lis_matrix_set_size(A, n, 0);
80 lis_matrix_set_csr(nnz, csr_mat.row_mem(), csr_mat.col_mem(), csr_mat.val_mem(), A);
81 lis_matrix_assemble(A);
82
83 lis_vector_create(0, &b);
84 lis_vector_create(0, &x);
85 lis_vector_set_size(b, n, 0);
86 lis_vector_set_size(x, n, 0);
87
88 lis_solver_set_option(setting.lis_options.c_str(), solver);
89
90 for(uword I = 0; I < B.n_cols; ++I) {
91 // ReSharper disable CppCStyleCast
92 lis_vector_set(b, (double*)B.colptr(I));
93 lis_vector_set(x, (double*)X.colptr(I));
94 // ReSharper restore CppCStyleCast
95
96 lis_solve(A, b, x, solver);
97 }
98
99 A->ptr = nullptr;
100 A->index = nullptr;
101 A->value = nullptr;
102 b->value = nullptr;
103 x->value = nullptr;
104
105 lis_matrix_destroy(A);
106 lis_vector_destroy(b);
107 lis_vector_destroy(x);
108
109 return SUANPAN_SUCCESS;
110}
111
112#endif
113
SolverSetting< T > setting
Definition: MetaMat.hpp:76
A SparseMat class that holds matrices.
Definition: SparseMat.hpp:34
A SparseMatLis class that holds matrices.
Definition: SparseMatLis.hpp:36
SparseMatLis(const SparseMatLis &other)
Definition: SparseMatLis.hpp:52
SparseMatLis(SparseMatLis &&) noexcept=delete
unique_ptr< MetaMat< T > > make_copy() override
Definition: SparseMatLis.hpp:64
SparseMatLis(const uword in_row, const uword in_col, const uword in_elem=0)
Definition: SparseMatLis.hpp:46
Definition: csr_form.hpp:25
const index_t * row_mem() const
Definition: csr_form.hpp:61
const index_t n_rows
Definition: csr_form.hpp:50
const data_t * val_mem() const
Definition: csr_form.hpp:65
const index_t * col_mem() const
Definition: csr_form.hpp:63
const index_t n_elem
Definition: csr_form.hpp:52
Definition: suanPan.h:331
int direct_solve(Mat< T > &, const Mat< T > &) override
Definition: SparseMatLis.hpp:67
constexpr auto SUANPAN_SUCCESS
Definition: suanPan.h:172