rheolef  7.0
solver_abtb.cc
Go to the documentation of this file.
1 #include "rheolef/solver_abtb.h"
2 #include "rheolef/skit.h"
3 
4 namespace rheolef {
5 
6 template <class T, class M>
8  : _opt(),
9  _sub_opt(),
10  _a(),
11  _b(),
12  _c(),
13  _mp(),
14  _sA(),
15  _need_constraint(false)
16 {
17 }
18 template <class T, class M>
20  const csr<T,M>& a,
21  const csr<T,M>& b,
22  const csr<T,M>& mp,
23  const solver_option& opt,
24  const solver_option& sub_opt)
25  : _opt(opt),
26  _sub_opt(sub_opt),
27  _a(a),
28  _b(b),
29  _c(),
30  _mp(mp),
31  _sA(),
32  _sa(),
33  _smp(),
34  _need_constraint(false)
35 {
36  _c.resize (_mp.row_ownership(), _mp.col_ownership());
37  init();
38 }
39 template <class T, class M>
41  const csr<T,M>& a,
42  const csr<T,M>& b,
43  const csr<T,M>& c,
44  const csr<T,M>& mp,
45  const solver_option& opt,
46  const solver_option& sub_opt)
47  : _opt(opt),
48  _sub_opt(sub_opt),
49  _a(a),
50  _b(b),
51  _c(c),
52  _mp(mp),
53  _sA(),
54  _sa(),
55  _smp(),
56  _need_constraint(false)
57 {
58  init();
59 }
60 template <class T, class M>
61 void
63 {
64 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
66  _opt.iterative = (_a.pattern_dimension() > 2);
67  }
68 
69 #if defined(_RHEOLEF_HAVE_FLOAT128)
70  _opt.iterative = true;
71 #endif // _RHEOLEF_HAVE_FLOAT128
72 
73  if (_opt.iterative) {
74 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
77 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
78  } else {
79  csr<T,M> A;
80  vec<T,M> one (_b.row_ownership(), 1);
81  vec<T,M> r =_b.trans_mult (one);
82  T z = dot(r,r);
83  if (_c.dis_nnz() != 0) {
84  vec<T,M> r2 =_c*one;
85  z += dot(r2,r2);
86  }
88  if (_need_constraint) {
89  vec<T,M> zu (_a.col_ownership(), 0);
90  vec<T,M> d = _mp*one;
91  A = { { _a, trans(_b), zu },
92  { _b, - _c, d },
93  { trans(zu), trans(d), T(0)} };
94  } else {
95  A = { { _a, trans(_b)},
96  { _b, - _c } };
97  }
98  A.set_pattern_dimension (_a.pattern_dimension());
99  if (_c.dis_nnz() == 0) {
100  A.set_symmetry (_a.is_symmetric());
101  } else {
102  A.set_symmetry (_a.is_symmetric() && _c.is_symmetric());
103  }
105  }
106 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
107 }
108 template <class T, class M>
109 void
111 {
112 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
113  if (_opt.iterative) {
114 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
115  check_macro (_a.is_symmetric(), "solver_abtb: iterative for unsymmetric system: not yet (HINT: use direct solver)");
116 #ifdef TO_CLEAN
117  size_type max_iter = _opt.max_iter;
118  T tol = _opt.tol;
119 #endif // TO_CLEAN
120  if (_c.dis_nnz() == 0) {
121  cg_abtb (_a, _b, u, p, f, g, _smp, _sa, option());
122  } else {
123  cg_abtbc (_a, _b, _c, u, p, f, g, _smp, _sa, option());
124  }
125  check_macro (option().residue <= option().tol,
126  "solver: precision "<<option().tol<<" not reached: get "<<option().residue
127  << " after " << option().max_iter << " iterations");
128 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
129  return;
130  }
131  vec<T,M> L;
132  if (_need_constraint) {
133  L = { f, g, T(0) };
134  } else {
135  L = { f, g };
136  }
137  vec<T,M> U = _sA.solve (L);
138  u = U [range(0,u.size())];
139  p = U [range(u.size(),u.size()+p.size())];
140 
141  if (_need_constraint) {
142  T lambda = (U.size() == u.size()+p.size()+1) ? U [u.size()+p.size()] : 0;
143 #ifdef _RHEOLEF_HAVE_MPI
144  mpi::broadcast (U.comm(), lambda, U.comm().size() - 1);
145 #endif // _RHEOLEF_HAVE_MPI
146  }
147 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
148 }
150 
151 #ifdef _RHEOLEF_HAVE_MPI
153 #endif // _RHEOLEF_HAVE_MPI
154 
155 } // namespace rheolef
void solve(const vec< T, M > &f, const vec< T, M > &g, vec< T, M > &u, vec< T, M > &p) const
Definition: solver_abtb.cc:110
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition: vec_expr_v2.h:335
solver_abtb – direct or iterative solver interface for mixed linear systems
Definition: solver_abtb.h:58
solver_basic< T, M > _sa
Definition: solver_abtb.h:86
csr< T, M >::size_type size_type
Definition: solver_abtb.h:62
irheostream, orheostream - large data streams
Definition: compiler.h:7
field residue(Float p, const field &uh)
#define check_macro(ok_condition, message)
Definition: compiler.h:104
size_t d
solver - direct or interative solver interface
Definition: solver.h:8
Definition: rotating-hill.h:1
solver_basic< T, M > _sA
Definition: solver_abtb.h:85
void g(boost::multi_array< double, 3 >::reference subarray, size_t jmax, size_t kmax)
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())
cg_abtb, cg_abtbc, minres_abtb, minres_abtbc – solvers for mixed linear problems ...
Definition: mixed_solver.h:119
Float epsilon
vec - vector in distributed environment
Definition: vec.h:48
static const long int decide
csr< T, sequential > trans(const csr< T, sequential > &a)
Definition: csr.h:381
solver_basic< T, M > _smp
Definition: solver_abtb.h:87
ublas::vector_range< ublas::vector< double > > f()
csr - compressed sparse row matrix
Definition: asr.h:8
const solver_option & option() const
Definition: solver_abtb.h:75
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())
Definition: mixed_solver.h:133