rheolef  6.6
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_type& opt,
24  const solver_option_type& 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_type& opt,
46  const solver_option_type& 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
65  if (_opt.iterative == solver_option_type::decide) {
66  _opt.iterative = (_a.pattern_dimension() > 2);
67  }
68 
69 #if defined(_RHEOLEF_HAVE_QD) || defined(_RHEOLEF_HAVE_FLOAT128)
70  _opt.iterative = true;
71 #endif // _RHEOLEF_HAVE_QD
72 
73  if (_opt.iterative) {
74 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
75  _sa = solver_basic<T,M>(_a, _sub_opt),
76  _smp = solver_basic<T,M>(_mp,_sub_opt);
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  }
87  _need_constraint = (fabs(z) <= std::numeric_limits<T>::epsilon());
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  }
104  _sA = solver_basic<T,M>(A,_opt);
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  size_type max_iter = _opt.max_iter;
117  T tol = _opt.tol;
118  if (_c.dis_nnz() == 0) {
119  pcg_abtb (_a, _b, u, p, f, g, _smp, _sa, max_iter, tol, &derr);
120  } else {
121  pcg_abtbc (_a, _b, _c, u, p, f, g, _smp, _sa, max_iter, tol, &derr);
122  }
123  check_macro (tol <= _opt.tol,
124  "solver: precision "<<_opt.tol<<" not reached: get "<<tol
125  << " after " << max_iter << " iterations");
126 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
127  return;
128  }
129  vec<T,M> L;
130  if (_need_constraint) {
131  L = { f, g, T(0) };
132  } else {
133  L = { f, g };
134  }
135  vec<T,M> U = _sA.solve (L);
136  u = U [range(0,u.size())];
137  p = U [range(u.size(),u.size()+p.size())];
138 
139  if (_need_constraint) {
140  T lambda = (U.size() == u.size()+p.size()+1) ? U [u.size()+p.size()] : 0;
141 #ifdef _RHEOLEF_HAVE_MPI
142  mpi::broadcast (U.comm(), lambda, U.comm().size() - 1);
143 #endif // _RHEOLEF_HAVE_MPI
144  }
145 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
146 }
148 
149 #ifdef _RHEOLEF_HAVE_MPI
151 #endif // _RHEOLEF_HAVE_MPI
152 
153 } // namespace rheolef
solver_option_type - options for direct or interative solvers
int pcg_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, Size &max_iter, Real &tol, odiststream *p_derr=0, std::string label="pcg_abtb")
Definition: mixed_solver.h:133
boost::proto::result_of::make_expr< boost::proto::tag::function, vec_domain, const vec_detail::fabs_, const vec< T, M > & >::type fabs(const vec< T, M > &x)
solver_abtb – direct or iterative solver iterface for mixed linear systems
Definition: solver_abtb.h:58
csr< T, M >::size_type size_type
Definition: solver_abtb.h:62
irheostream, orheostream - large data streams
Definition: compiler.h:7
static const long int decide
#define check_macro(ok_condition, message)
Definition: compiler.h:104
size_t d
int pcg_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, Size &max_iter, Real &tol, odiststream *p_derr=0, std::string label="pcg_abtbc")
pcg_abtb, pcg_abtbc, pminres_abtb, pminres_abtbc – solvers for mixed linear problems ...
Definition: mixed_solver.h:119
Definition: rotating-hill.h:1
void g(boost::multi_array< double, 3 >::reference subarray, size_t jmax, size_t kmax)
Float epsilon
T dot(const vec< T, M > &x, const int &y)
vec - vector in distributed environment
Definition: vec.h:44
odiststream derr(cerr)
Definition: diststream.h:340
solver - direct or interative solver interface
Definition: solver.h:8
csr< T, sequential > trans(const csr< T, sequential > &a)
Definition: csr.h:383
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
ublas::vector_range< ublas::vector< double > > f()
csr - compressed sparse row matrix
Definition: asr.h:8