rheolef  6.5
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  _a(),
10  _b(),
11  _c(),
12  _mp(),
13  _sA(),
14  _need_constraint(false)
15 {
16 }
17 template <class T, class M>
19  const csr<T,M>& a,
20  const csr<T,M>& b,
21  const csr<T,M>& mp,
22  const solver_option_type& opt)
23  : _opt(opt),
24  _a(a),
25  _b(b),
26  _c(),
27  _mp(mp),
28  _sA(),
29  _sa(),
30  _smp(),
31  _need_constraint(false)
32 {
33  _c.resize (_mp.row_ownership(), _mp.col_ownership());
34  init();
35 }
36 template <class T, class M>
38  const csr<T,M>& a,
39  const csr<T,M>& b,
40  const csr<T,M>& c,
41  const csr<T,M>& mp,
42  const solver_option_type& opt)
43  : _opt(opt),
44  _a(a),
45  _b(b),
46  _c(c),
47  _mp(mp),
48  _sA(),
49  _sa(),
50  _smp(),
51  _need_constraint(false)
52 {
53  init();
54 }
55 template <class T, class M>
56 void
58 {
59 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
60  if (_opt.iterative == solver_option_type::decide) {
61  _opt.iterative = (_a.pattern_dimension() > 2);
62  }
63  if (_opt.iterative) {
64 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
65  _sa = solver_basic<T,M>(_a, _opt),
66  _smp = solver_basic<T,M>(_mp,_opt);
67 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
68  } else {
69  csr<T,M> A;
70  vec<T,M> one (_b.row_ownership(), 1);
71  vec<T,M> r =_b.trans_mult (one);
72  T z = dot(r,r);
73  if (_c.dis_nnz() != 0) {
74  vec<T,M> r2 =_c*one;
75  z += dot(r2,r2);
76  }
77  _need_constraint = (fabs(z) + 1 == 1);
78  if (_need_constraint) {
79  vec<T,M> zu (_a.col_ownership(), 0);
80  vec<T,M> d = _mp*one;
81  A = { { _a, trans(_b), zu },
82  { _b, - _c, d },
83  { trans(zu), trans(d), 0 } };
84  } else {
85  A = { { _a, trans(_b)},
86  { _b, - _c } };
87  }
88  A.set_pattern_dimension (_a.pattern_dimension());
89  if (_c.dis_nnz() == 0) {
90  A.set_symmetry (_a.is_symmetric());
91  } else {
92  A.set_symmetry (_a.is_symmetric() && _c.is_symmetric());
93  }
94  _sA = solver_basic<T,M>(A,_opt);
95  }
96 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
97 }
98 template <class T, class M>
99 void
101 {
102 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
103  if (_opt.iterative) {
104 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
105  check_macro (_a.is_symmetric(), "solver_abtb: iterative for unsymmetric system: not yet (HINT: use direct solver)");
106  size_type max_iter = _opt.max_iter;
107  T tol = _opt.tol;
108  if (_c.dis_nnz() == 0) {
109  pcg_abtb (_a, _b, u, p, f, g, _smp, _sa, max_iter, tol, &derr);
110  } else {
111  pcg_abtbc (_a, _b, _c, u, p, f, g, _smp, _sa, max_iter, tol, &derr);
112  }
113  check_macro (tol <= _opt.tol,
114  "solver: precision "<<_opt.tol<<" not reached: get "<<tol
115  << " after " << max_iter << " iterations");
116 #ifdef _RHEOLEF_HAVE_STD_INITIALIZER_LIST
117  return;
118  }
119  vec<T,M> L;
120  if (_need_constraint) {
121  L = { f, g, 0 };
122  } else {
123  L = { f, g };
124  }
125  vec<T,M> U = _sA.solve (L);
126  u = U [range(0,u.size())];
127  p = U [range(u.size(),u.size()+p.size())];
128 
129  if (_need_constraint) {
130  T lambda = (U.size() == u.size()+p.size()+1) ? U [u.size()+p.size()] : 0;
131 #ifdef _RHEOLEF_HAVE_MPI
132  mpi::broadcast (U.comm(), lambda, U.comm().size() - 1);
133 #endif // _RHEOLEF_HAVE_MPI
134  }
135 #endif // _RHEOLEF_HAVE_STD_INITIALIZER_LIST
136 }
138 
139 #ifdef _RHEOLEF_HAVE_MPI
141 #endif // _RHEOLEF_HAVE_MPI
142 
143 } // namespace rheolef
144