37template<
typename T,
typename U>
concept ArmaContainer = std::is_floating_point_v<U> && (std::is_convertible_v<T, Mat<U>> || std::is_convertible_v<T, SpMat<U>>) ;
44 shared_ptr<MetaMat<T>> mat_a, mat_b;
69 , bracket(std::forward<
op_add<
T>>(B)) {}
90 template<std::invocable<fmat&> F>
int mixed_trs(mat& X, mat&& B, F trs) {
93 X = arma::zeros(size(B));
95 auto multiplier =
norm(B);
101 auto residual = conv_to<fmat>::from(B / multiplier);
103 if(0 != (INFO = trs(residual)))
break;
105 const mat incre = multiplier * conv_to<mat>::from(residual);
109 suanpan_debug(
"Mixed precision algorithm multiplier: {:.5E}.\n", multiplier =
arma::norm(B -= this->
operator*(incre)));
122 MetaMat(
const uword in_rows,
const uword in_cols,
const uword in_elem)
147 this->
at(K,
K) =
T(1);
152 [[nodiscard]]
virtual T max()
const = 0;
153 [[nodiscard]]
virtual Col<T>
diag()
const = 0;
164 virtual T&
unsafe_at(
const uword I,
const uword J) {
return this->
at(I, J); }
170 virtual T&
at(uword, uword) = 0;
172 [[nodiscard]]
virtual const T*
memptr()
const = 0;
183 const auto& bracket =
M.bracket;
184 if(
nullptr != bracket.mat_a) this->
scale_accu(M.scalar, bracket.mat_a);
185 if(
nullptr != bracket.mat_b) this->
scale_accu(M.scalar, bracket.mat_b);
189 const auto& bracket =
M.bracket;
190 if(
nullptr != bracket.mat_a) this->
scale_accu(-M.scalar, bracket.mat_a);
191 if(
nullptr != bracket.mat_b) this->
scale_accu(-M.scalar, bracket.mat_b);
204 template<ArmaContainer<T> C> Mat<T>
solve(C&& B) {
207 if(
SUANPAN_SUCCESS != this->
solve(X, std::forward<C>(B)))
throw std::runtime_error(
"fail to solve the system");
215 if(!
to_mat(*this).save(name, raw_ascii))
227 this->csc_condense();
229 X.zeros(arma::size(B));
231 unique_ptr<Preconditioner<T>> preconditioner;
232 if(
PreconditionerType::JACOBI == this->setting.preconditioner_type) preconditioner = std::make_unique<Jacobi<T>>(this->diag());
233#ifndef SUANPAN_SUPERLUMT
235 if(this->triplet_mat.is_empty()) preconditioner = std::make_unique<
ILU<T>>(to_triplet_form<T, int>(
this));
236 else preconditioner = std::make_unique<ILU<T>>(this->triplet_mat);
239 else if(
PreconditionerType::NONE == this->setting.preconditioner_type) preconditioner = std::make_unique<UnityPreconditioner<T>>();
243 this->setting.preconditioner = preconditioner.get();
245 std::atomic_int code = 0;
249 Col<T> sub_x(X.colptr(I), X.n_rows, false, true);
250 const Col<T> sub_b(B.colptr(I), B.n_rows);
251 auto col_setting = setting;
252 code += GMRES(this, sub_x, sub_b, col_setting);
256 Col<T> sub_x(X.colptr(I), X.n_rows, false, true);
257 const Col<T> sub_b(B.colptr(I), B.n_rows);
258 auto col_setting = setting;
259 code += BiCGSTAB(this, sub_x, sub_b, col_setting);
261 else throw invalid_argument(
"no proper iterative solver assigned but somehow iterative solving is called");
268 for(uword J = 0; J < in_mat.
n_cols; ++J)
for(uword I = 0; I < in_mat.
n_rows; ++I) out_mat(I, J) = in_mat(I, J);
275 Mat<data_t> out_mat(in_mat.
n_rows, in_mat.
n_cols, fill::zeros);
276 for(index_t I = 0; I < in_mat.
n_elem; ++I) out_mat(in_mat.
row(I), in_mat.
col(I)) += in_mat.
val(I);
281 Mat<data_t> out_mat(in_mat.
n_rows, in_mat.
n_cols, fill::zeros);
284 for(index_t I = 0; I < in_mat.
n_elem; ++I) {
285 if(I >= in_mat.
row_mem()[c_idx]) ++c_idx;
293 Mat<data_t> out_mat(in_mat.
n_rows, in_mat.
n_cols, fill::zeros);
296 for(index_t I = 0; I < in_mat.
n_elem; ++I) {
297 if(I >= in_mat.
col_mem()[c_idx]) ++c_idx;
307 const sp_i auto n_rows = index_t(in_mat->
n_rows);
308 const sp_i auto n_cols = index_t(in_mat->
n_cols);
309 const sp_i auto n_elem = index_t(in_mat->
n_elem);
312 for(index_t J = 0; J < n_cols; ++J)
for(index_t I = 0; I < n_rows; ++I) out_mat.
at(I, J) = in_mat->operator()(I, J);
A ILU class.
Definition: ILU.hpp:40
Definition: MetaMat.hpp:41
op_add(const shared_ptr< MetaMat< T > > &A, const shared_ptr< MetaMat< T > > &B)
Definition: MetaMat.hpp:51
op_add(const shared_ptr< MetaMat< T > > &A)
Definition: MetaMat.hpp:47
Definition: MetaMat.hpp:56
op_scale(const T A, const shared_ptr< MetaMat< T > > &B)
Definition: MetaMat.hpp:63
op_scale(const T A, op_add< T > &&B)
Definition: MetaMat.hpp:67
Definition: MetaMat.hpp:37
Definition: suanPan.h:331
void for_each(const IT start, const IT end, F &&FN)
Definition: utility.h:28
double norm(const vec &)
Definition: tensor.cpp:370
unsigned iterative_refinement
Definition: SolverSetting.hpp:44
data_t tolerance
Definition: SolverSetting.hpp:43
IterativeSolver iterative_solver
Definition: SolverSetting.hpp:46
#define suanpan_debug(...)
Definition: suanPan.h:307
constexpr auto SUANPAN_SUCCESS
Definition: suanPan.h:172
constexpr auto SUANPAN_FAIL
Definition: suanPan.h:173
#define suanpan_error(...)
Definition: suanPan.h:309