21#include "rheolef/solver_abtb.h"
22#include "rheolef/linalg.h"
26template <
class T,
class M>
36 _need_constraint(false)
39template <
class T,
class M>
53 _need_constraint(false)
55 _c.resize (
_mp.row_ownership(),
_mp.col_ownership());
58template <
class T,
class M>
73 _need_constraint(false)
77template <
class T,
class M>
82 _opt.iterative = (_a.pattern_dimension() > 2);
84 if (! _opt.iterative) {
87 vec<T,M> one (_b.row_ownership(), 1);
90 if (_c.dis_nnz() != 0) {
94 _need_constraint = (fabs(z) <= std::numeric_limits<T>::epsilon());
95 if (_need_constraint) {
98 A = { { _a,
trans(_b), zu },
105 A.set_pattern_dimension (_a.pattern_dimension());
106 A.set_definite_positive (
false);
107 if (_c.dis_nnz() == 0) {
108 A.set_symmetry (_a.is_symmetric());
110 A.set_symmetry (_a.is_symmetric() && _c.is_symmetric());
115template <
class T,
class M>
119 if (_opt.iterative) {
124 check_macro (_a.is_symmetric(),
"solver_abtb: iterative for unsymmetric system: not yet (HINT: use direct solver)");
125 if (_c.dis_nnz() == 0) {
131 "solver: precision "<<option().tol<<
" not reached: get "<<option().
residue
132 <<
" after " << option().max_iter <<
" iterations");
137 if (_need_constraint) {
142 vec<T,M> U = _sA.solve (L);
143 u = U [range(0,
u.size())];
146 if (_need_constraint) {
148 T lambda = (U.size() ==
u.size()+
p.size()+1) ? U [
u.size()+
p.size()] : 0;
149#ifdef _RHEOLEF_HAVE_MPI
150 mpi::broadcast (U.comm(),
lambda, U.comm().size() - 1);
154template <
class T,
class M>
158 return (_a.dis_nnz() != 0);
165#ifdef _RHEOLEF_HAVE_MPI
see the csr page for the full documentation
void solve(const vec< T, M > &f, const vec< T, M > &g, vec< T, M > &u, vec< T, M > &p) const
see the solver_option page for the full documentation
static const long int decide
see the vec page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
int cg_abtbc(const Matrix &A, const Matrix &B, const Matrix &C, Vector &u, Vector &p, const VectorExpr1 &Mf, const VectorExpr2 &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const solver_option &sopt=solver_option())
int cg_abtb(const Matrix &A, const Matrix &B, Vector &u, Vector &p, const VectorExpr1 &Mf, const VectorExpr2 &Mg, const Preconditioner &S1, const Solver &inner_solver_A, const solver_option &sopt=solver_option())
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
field residue(Float p, const field &uh)
see the range page for the full documentation