an efficient C++ finite element environment
|
|
Go to the documentation of this file. 1 #include "rheolef/solver_abtb.h"
22 #include "rheolef/linalg.h"
26 template <
class T,
class M>
36 _need_constraint(false)
39 template <
class T,
class M>
53 _need_constraint(false)
55 _c.resize (
_mp.row_ownership(),
_mp.col_ownership());
58 template <
class T,
class M>
73 _need_constraint(false)
77 template <
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) {
95 if (_need_constraint) {
98 A = { { _a,
trans(_b), zu },
102 A = { { _a,
trans(_b)},
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());
115 template <
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);
151 #endif // _RHEOLEF_HAVE_MPI
154 template <
class T,
class M>
158 return (_a.dis_nnz() != 0);
165 #ifdef _RHEOLEF_HAVE_MPI
167 #endif // _RHEOLEF_HAVE_MPI
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
field residue(Float p, const field &uh)
see the range page for the full documentation
see the vec page for the full documentation
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())
see the csr page for the full documentation
This file is part of Rheolef.
void solve(const vec< T, M > &f, const vec< T, M > &g, vec< T, M > &u, vec< T, M > &p) const
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
static const long int decide
see the solver_option page for the full documentation
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())