dune-composites (2.5.1)

cg_fork.hh
1#ifndef CGSOLVERFORK
2#define CGSOLVERFORK
3
4#include <dune/istl/eigenvalue/arpackpp.hh>
5#include <dune/istl/scalarproducts.hh>
6#include <dune/common/timer.hh>
7
8namespace Dune {
9 namespace Geneo{
11 template<class X>
12 class CGSolverFork : public InverseOperator<X,X> {
13 public:
15 typedef X domain_type;
17 typedef X range_type;
19 typedef typename X::field_type field_type;
21 typedef typename FieldTraits<field_type>::real_type real_type;
22 //typedef typename FieldTraits<field_type>::real_type real_type;
23
24 // copy base class constructors
25 using InverseOperator<X,X>::InverseOperator;
26
32 template<class L, class P>
33 CGSolverFork (L& op, P& prec, real_type reduction, int maxit, int verbose) :
34 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
35 {
36 /*static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
37 "L and P must have the same category!");
38 static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential),
39 "L must be sequential!");*/
40 }
46 template<class L, class S, class P>
47 CGSolverFork (L& op, S& sp, P& prec, real_type reduction, int maxit, int verbose) :
48 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose)
49 {
50 /*static_assert(static_cast<int>(L::category) == static_cast<int>(P::category),
51 "L and P must have the same category!");
52 static_assert(static_cast<int>(L::category) == static_cast<int>(S::category),
53 "L and S must have the same category!");*/
54 }
55
56 CGSolverFork (LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
57 real_type reduction, int maxit, int verbose, bool condition_estimate) : CGSolverFork<X>(op, prec, reduction, maxit, verbose)
58 {
59 _condition_estimate = condition_estimate;
60 }
61
62 CGSolverFork (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
63 real_type reduction, int maxit, int verbose, bool condition_estimate) : CGSolverFork<X>(op, sp, prec, reduction, maxit, verbose)
64 {
65 _condition_estimate = condition_estimate;
66 }
67
68
80 virtual void apply (X& x, X& b, InverseOperatorResult& res)
81 {
82 using std::isfinite;
83
84 res.clear(); // clear solver statistics
85 _prec.pre(x,b); // prepare preconditioner
86 Timer watch; // start a timer
87 _op.applyscaleadd(-1,x,b); // overwrite b with defect
88
89 X p(x); // the search direction
90 X q(x); // a temporary vector
91
92 real_type def0 = _sp.norm(b); // compute norm
93
94 if (!isfinite(def0)) // check for inf or NaN
95 {
96 if (_verbose>0)
97 std::cout << "=== CGSolver: abort due to infinite or NaN initial defect"
98 << std::endl;
99 DUNE_THROW(SolverAbort, "CGSolver: initial defect=" << def0
100 << " is infinite or NaN");
101 }
102
103 if (def0<1E-30) // convergence check
104 {
105 res.converged = true;
106 res.iterations = 0; // fill statistics
107 res.reduction = 0;
108 res.conv_rate = 0;
109 res.elapsed=0;
110 if (_verbose>0) // final print
111 std::cout << "=== rate=" << res.conv_rate
112 << ", T=" << res.elapsed << ", TIT=" << res.elapsed
113 << ", IT=0" << std::endl;
114 return;
115 }
116
117 if (_verbose>0) // printing
118 {
119 std::cout << "=== CGSolver" << std::endl;
120 if (_verbose>1) {
121 this->printHeader(std::cout);
122 this->printOutput(std::cout,real_type(0),def0);
123 }
124 }
125
126 // some local variables
127 real_type def=def0; // loop variables
128 field_type rho,rholast,lambda,alpha,beta;
129
130 // determine initial search direction
131 p = 0; // clear correction
132 _prec.apply(p,b); // apply preconditioner
133 rholast = _sp.dot(p,b); // orthogonalization
134
135 // Remember alpha and beta values for condition estimate
136 std::vector<real_type> alphas(0);
137 std::vector<real_type> betas(0);
138 //std::vector<std::shared_ptr<X>> ps(0);
139
140 // the loop
141 int i=1;
142 for ( ; i<=_maxit; i++ )
143 {
144
145 /*if (configuration.get<bool>("OrthogonalizeCG",false)) {
146
147 for (auto previous_p : ps) {
148 X correction = x;
149 _op.apply (*previous_p, correction);
150
151 double scp = _sp.dot(p, correction) / _sp.dot(*previous_p, correction);
152 correction = *previous_p;
153 correction *= scp;
154 p -= correction;
155 }
156 ps.push_back(std::make_shared<X>(p));
157
158 }*/
159
160 // minimize in given search direction p
161 _op.apply(p,q); // q=Ap
162 alpha = _sp.dot(p,q); // scalar product
163 lambda = rholast/alpha; // minimization
164 alphas.push_back(lambda);
165 x.axpy(lambda,p); // update solution
166 b.axpy(-lambda,q); // update defect
167
168 // convergence test
169 real_type defnew=_sp.norm(b); // comp defect norm
170
171 if (_verbose>1) // print
172 this->printOutput(std::cout,real_type(i),defnew,def);
173
174 def = defnew; // update norm
175 if (!isfinite(def)) // check for inf or NaN
176 {
177 if (_verbose>0)
178 std::cout << "=== CGSolver: abort due to infinite or NaN defect"
179 << std::endl;
180 DUNE_THROW(SolverAbort,
181 "CGSolver: defect=" << def << " is infinite or NaN");
182 }
183
184 if (def<def0*_reduction || def<1E-30) // convergence check
185 {
186 res.converged = true;
187 break;
188 }
189
190 // determine new search direction
191 q = 0; // clear correction
192 _prec.apply(q,b); // apply preconditioner
193 rho = _sp.dot(q,b); // orthogonalization
194 beta = rho/rholast; // scaling factor
195 betas.push_back(beta);
196 p *= beta; // scale old search direction
197 p += q; // orthogonalization with correction
198 rholast = rho; // remember rho for recurrence
199 }
200
201 //correct i which is wrong if convergence was not achieved.
202 i=std::min(_maxit,i);
203
204 res.elapsed = watch.elapsed();
205
206 if (_condition_estimate) {
207#if HAVE_ARPACKPP
208
209 // Build T matrix
210 MAT T(i, i, MAT::row_wise);
211
212 for (auto row = T.createbegin(); row != T.createend(); ++row) {
213 if (row.index() > 0)
214 row.insert(row.index()-1);
215 row.insert(row.index());
216 if (row.index() < T.N() - 1)
217 row.insert(row.index()+1);
218 }
219 for (int row = 0; row < i; ++row) {
220 if (row > 0) {
221 T[row][row-1] = std::sqrt(betas[row-1]) / alphas[row-1];
222 }
223
224 T[row][row] = 1.0 / alphas[row];
225 if (row > 0) {
226 T[row][row] += betas[row-1] / alphas[row-1];
227 }
228
229 if (row < i - 1) {
230 T[row][row+1] = std::sqrt(betas[row]) / alphas[row];
231 }
232 }
233
234 // Compute largest and smallest eigenvalue of T matrix and return as estimate
235 Dune::ArPackPlusPlus_Algorithms<MAT, VEC> arpack(T);
236
237 double eps = 0.0;
238 VEC eigv;
239 double min_eigv, max_eigv;
240 arpack.computeSymMinMagnitude (eps, eigv, min_eigv);
241 arpack.computeSymMaxMagnitude (eps, eigv, max_eigv);
242
243 if (_verbose > 0) {
244 std::cout << "Min eigv: " << min_eigv << std::endl;
245 std::cout << "Max eigv: " << max_eigv << std::endl;
246 std::cout << "Condition: " << max_eigv / min_eigv << std::endl;
247 }
248
249#else
250 std::cerr << "WARNING: Condition estimate was requested. This requires ARPACK, but ARPACK was not found!" << std::endl;
251#endif
252 }
253
254 if (_verbose==1) // printing for non verbose
255 this->printOutput(std::cout,real_type(i),def);
256
257 _prec.post(x); // postprocess preconditioner
258 res.iterations = i; // fill statistics
259 res.reduction = static_cast<double>(def/def0);
260 res.conv_rate = static_cast<double>(pow(res.reduction,1.0/i));
261
262 if (_verbose>0) // final print
263 {
264 std::cout << "=== rate=" << res.conv_rate
265 << ", T=" << res.elapsed
266 << ", TIT=" << res.elapsed/i
267 << ", IT=" << i << std::endl;
268 }
269 }
270
282 virtual void apply (X& x, X& b, double reduction,
283 InverseOperatorResult& res)
284 {
285 real_type saved_reduction = _reduction;
286 _reduction = reduction;
287 (*this).apply(x,b,res);
288 _reduction = saved_reduction;
289 }
290
291 private:
292 typedef Dune::BCRSMatrix<Dune::FieldMatrix<real_type,1,1> > MAT;
293 typedef Dune::BlockVector<Dune::FieldVector<real_type,1> > VEC;
294 bool _condition_estimate;
295 private:
296 SeqScalarProduct<X> ssp;
297 LinearOperator<X,X>& _op;
298 Preconditioner<X,X>& _prec;
299 ScalarProduct<X>& _sp;
300 real_type _reduction;
301 int _maxit;
302 int _verbose;
303 /*SeqScalarProduct<X> ssp;
304 LinearOperator<X,X>& _op;
305 Preconditioner<X,X>& _prec;
306 ScalarProduct<X>& _sp;
307 real_type _reduction;
308 int _maxit;
309 int _verbose;*/
310 };
311
312 }
313}
314#endif
Dune::Geneo::CGSolverFork
conjugate gradient method
Definition: cg_fork.hh:12
Dune::Geneo::CGSolverFork::domain_type
X domain_type
The domain type of the operator to be inverted.
Definition: cg_fork.hh:15
Dune::Geneo::CGSolverFork::apply
virtual void apply(X &x, X &b, InverseOperatorResult &res)
Apply inverse operator.
Definition: cg_fork.hh:80
Dune::Geneo::CGSolverFork::CGSolverFork
CGSolverFork(L &op, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: cg_fork.hh:33
Dune::Geneo::CGSolverFork::CGSolverFork
CGSolverFork(L &op, S &sp, P &prec, real_type reduction, int maxit, int verbose)
Set up conjugate gradient solver.
Definition: cg_fork.hh:47
Dune::Geneo::CGSolverFork::range_type
X range_type
The range type of the operator to be inverted.
Definition: cg_fork.hh:17
Dune::Geneo::CGSolverFork::real_type
FieldTraits< field_type >::real_type real_type
The real type of the field type (is the same if using real numbers, but differs for std::complex)
Definition: cg_fork.hh:21
Dune::Geneo::CGSolverFork::apply
virtual void apply(X &x, X &b, double reduction, InverseOperatorResult &res)
Apply inverse operator with given reduction factor.
Definition: cg_fork.hh:282
Dune::Geneo::CGSolverFork::field_type
X::field_type field_type
The field type of the operator to be inverted.
Definition: cg_fork.hh:19
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 5, 22:35, 2025)