Branch data Line data Source code
1 : : //
2 : : // General pre-conditioners for the pre-conditioned Conjugate Gradient Solver
3 : : // such as Jacobi, Polynomial Newton, Polynomial Chebyshev, Richardson, Gauss-Seidel
4 : : // provides a function operator() that returns the preconditioned field
5 : : //
6 : :
7 : : #ifndef IPPL_PRECONDITIONER_H
8 : : #define IPPL_PRECONDITIONER_H
9 : :
10 : : #include "Expression/IpplOperations.h" // get the function apply()
11 : :
12 : : // Expands to a lambda that acts as a wrapper for a differential operator
13 : : // fun: the function for which to create the wrapper, such as ippl::laplace
14 : : // type: the argument type, which should match the LHS type for the solver
15 : : #define IPPL_SOLVER_OPERATOR_WRAPPER(fun, type) \
16 : : [](type arg) { \
17 : : return fun(arg); \
18 : : }
19 : :
20 : : namespace ippl {
21 : : template <typename Field>
22 : : struct preconditioner {
23 : : constexpr static unsigned Dim = Field::dim;
24 : : using mesh_type = typename Field::Mesh_t;
25 : : using layout_type = typename Field::Layout_t;
26 : :
27 : 0 : preconditioner()
28 [ # # ]: 0 : : type_m("Identity") {}
29 : :
30 : 0 : preconditioner(std::string name)
31 : 0 : : type_m(name) {}
32 : :
33 : 0 : virtual ~preconditioner() = default;
34 : :
35 : : // Placeholder for the function operator, actually implemented in the derived classes
36 : 0 : virtual Field operator()(Field& u) {
37 : 0 : Field res = u.deepCopy();
38 : 0 : return res;
39 : : }
40 : :
41 : : // Placeholder for setting additional fields, actually implemented in the derived classes
42 : 0 : virtual void init_fields(Field& b) {
43 [ # # ]: 0 : Field res = b.deepCopy();
44 : 0 : return;
45 : 0 : }
46 : :
47 : : std::string get_type() { return type_m; };
48 : :
49 : : protected:
50 : : std::string type_m;
51 : : };
52 : :
53 : : /*!
54 : : * Jacobi preconditioner
55 : : * M = w*diag{A} // w is a damping factor
56 : : */
57 : : template <typename Field, typename InvDiagF>
58 : : struct jacobi_preconditioner : public preconditioner<Field> {
59 : : constexpr static unsigned Dim = Field::dim;
60 : : using mesh_type = typename Field::Mesh_t;
61 : : using layout_type = typename Field::Layout_t;
62 : :
63 : 0 : jacobi_preconditioner(InvDiagF&& inverse_diagonal, double w = 1.0)
64 : : : preconditioner<Field>("jacobi")
65 [ # # # # ]: 0 : , w_m(w) {
66 : 0 : inverse_diagonal_m = std::move(inverse_diagonal);
67 : 0 : }
68 : :
69 : 0 : Field operator()(Field& u) override {
70 : 0 : mesh_type& mesh = u.get_mesh();
71 : 0 : layout_type& layout = u.getLayout();
72 : 0 : Field res(mesh, layout);
73 : :
74 [ # # # # : 0 : res = inverse_diagonal_m(u);
# # ]
75 [ # # # # ]: 0 : res = w_m * res;
76 : 0 : return res;
77 : 0 : }
78 : :
79 : : protected:
80 : : InvDiagF inverse_diagonal_m;
81 : : double w_m; // Damping factor
82 : : };
83 : :
84 : : /*!
85 : : * Polynomial Newton Preconditioner
86 : : * Computes iteratively approximations for A^-1
87 : : */
88 : : template <typename Field, typename OperatorF>
89 : : struct polynomial_newton_preconditioner : public preconditioner<Field> {
90 : : constexpr static unsigned Dim = Field::dim;
91 : : using mesh_type = typename Field::Mesh_t;
92 : : using layout_type = typename Field::Layout_t;
93 : :
94 : 0 : polynomial_newton_preconditioner(OperatorF&& op, double alpha, double beta,
95 : : unsigned int max_level = 6, double zeta = 1e-3,
96 : : double* eta = nullptr)
97 : : : preconditioner<Field>("polynomial_newton")
98 : 0 : , alpha_m(alpha)
99 : 0 : , beta_m(beta)
100 : 0 : , level_m(max_level)
101 : 0 : , zeta_m(zeta)
102 [ # # # # ]: 0 : , eta_m(eta) {
103 : 0 : op_m = std::move(op);
104 : 0 : }
105 : :
106 : : // Memory management is needed because of runtime defined eta_m
107 : 0 : ~polynomial_newton_preconditioner() {
108 : 0 : if (eta_m != nullptr) {
109 [ # # ]: 0 : delete[] eta_m;
110 : 0 : eta_m = nullptr;
111 : : }
112 [ # # ]: 0 : }
113 : :
114 : : polynomial_newton_preconditioner(const polynomial_newton_preconditioner& other)
115 : : : preconditioner<Field>("polynomial_newton")
116 : : , level_m(other.level_m)
117 : : , alpha_m(other.alpha_m)
118 : : , beta_m(other.beta_m)
119 : : , zeta_m(other.zeta_m)
120 : : , eta_m(other.eta_m) {
121 : : op_m = std::move(other.op_m);
122 : : }
123 : :
124 : : polynomial_newton_preconditioner& operator=(const polynomial_newton_preconditioner& other) {
125 : : return *this = polynomial_newton_preconditioner(other);
126 : : }
127 : :
128 : 0 : Field recursive_preconditioner(Field& u, unsigned int level) {
129 : 0 : mesh_type& mesh = u.get_mesh();
130 [ # # ]: 0 : layout_type& layout = u.getLayout();
131 : : // Define etas if not defined yet
132 [ # # ]: 0 : if (eta_m == nullptr) {
133 : : // Precompute the etas for later use
134 [ # # # # ]: 0 : eta_m = new double[level_m + 1];
135 : 0 : eta_m[0] = 2.0 / ((alpha_m + beta_m) * (1.0 + zeta_m));
136 [ # # ]: 0 : if (level_m > 0) {
137 : 0 : eta_m[1] =
138 : : 2.0
139 : 0 : / (1.0 + 2 * alpha_m * eta_m[0] - alpha_m * eta_m[0] * alpha_m * eta_m[0]);
140 : : }
141 [ # # ]: 0 : for (unsigned int i = 2; i < level_m + 1; ++i) {
142 : 0 : eta_m[i] = 2.0 / (1.0 + 2 * eta_m[i - 1] - eta_m[i - 1] * eta_m[i - 1]);
143 : : }
144 : : }
145 : :
146 [ # # ]: 0 : Field res(mesh, layout);
147 : : // Base case
148 [ # # ]: 0 : if (level == 0) {
149 [ # # # # ]: 0 : res = eta_m[0] * u;
150 : 0 : return res;
151 : : }
152 : : // Recursive case
153 [ # # ]: 0 : Field PAPr(mesh, layout);
154 [ # # ]: 0 : Field Pr(mesh, layout);
155 : :
156 [ # # # # ]: 0 : Pr = recursive_preconditioner(u, level - 1);
157 [ # # # # : 0 : PAPr = op_m(Pr);
# # ]
158 [ # # # # ]: 0 : PAPr = recursive_preconditioner(PAPr, level - 1);
159 [ # # # # : 0 : res = eta_m[level] * (2.0 * Pr - PAPr);
# # # # ]
160 : 0 : return res;
161 : 0 : }
162 : :
163 : 0 : Field operator()(Field& u) override { return recursive_preconditioner(u, level_m); }
164 : :
165 : : protected:
166 : : OperatorF op_m; // Operator to be preconditioned
167 : : double alpha_m; // Smallest Eigenvalue
168 : : double beta_m; // Largest Eigenvalue
169 : : unsigned int level_m; // Number of recursive calls
170 : : double zeta_m; // smallest (alpha + beta) is multiplied by (1+zeta) to avoid clustering of
171 : : // Eigenvalues
172 : : double* eta_m = nullptr; // Size is determined at runtime
173 : : };
174 : :
175 : : /*!
176 : : * Polynomial Chebyshev Preconditioner
177 : : * Computes iteratively approximations for A^-1
178 : : */
179 : : template <typename Field, typename OperatorF>
180 : : struct polynomial_chebyshev_preconditioner : public preconditioner<Field> {
181 : : constexpr static unsigned Dim = Field::dim;
182 : : using mesh_type = typename Field::Mesh_t;
183 : : using layout_type = typename Field::Layout_t;
184 : :
185 : 0 : polynomial_chebyshev_preconditioner(OperatorF&& op, double alpha, double beta,
186 : : unsigned int degree = 63, double zeta = 1e-3)
187 : : : preconditioner<Field>("polynomial_chebyshev")
188 : 0 : , alpha_m(alpha)
189 : 0 : , beta_m(beta)
190 : 0 : , degree_m(degree)
191 : 0 : , zeta_m(zeta)
192 [ # # # # ]: 0 : , rho_m(nullptr) {
193 : 0 : op_m = std::move(op);
194 : 0 : }
195 : :
196 : : // Memory management is needed because of runtime defined rho_m
197 : 0 : ~polynomial_chebyshev_preconditioner() {
198 : 0 : if (rho_m != nullptr) {
199 [ # # ]: 0 : delete[] rho_m;
200 : 0 : rho_m = nullptr;
201 : : }
202 [ # # ]: 0 : }
203 : :
204 : : polynomial_chebyshev_preconditioner(const polynomial_chebyshev_preconditioner& other)
205 : : : preconditioner<Field>("polynomial_chebyshev")
206 : : , degree_m(other.degree_m)
207 : : , theta_m(other.theta_m)
208 : : , sigma_m(other.sigma_m)
209 : : , delta_m(other.delta_m)
210 : : , alpha_m(other.delta_m)
211 : : , beta_m(other.delta_m)
212 : : , zeta_m(other.zeta_m)
213 : : , rho_m(other.rho_m) {
214 : : op_m = std::move(other.op_m);
215 : : }
216 : :
217 : : polynomial_chebyshev_preconditioner& operator=(
218 : : const polynomial_chebyshev_preconditioner& other) {
219 : : return *this = polynomial_chebyshev_preconditioner(other);
220 : : }
221 : :
222 : 0 : Field operator()(Field& r) override {
223 : 0 : mesh_type& mesh = r.get_mesh();
224 [ # # ]: 0 : layout_type& layout = r.getLayout();
225 : :
226 [ # # ]: 0 : Field res(mesh, layout);
227 [ # # ]: 0 : Field x(mesh, layout);
228 [ # # ]: 0 : Field x_old(mesh, layout);
229 [ # # ]: 0 : Field A(mesh, layout);
230 [ # # ]: 0 : Field z(mesh, layout);
231 : :
232 : : // Precompute the coefficients if not done yet
233 [ # # ]: 0 : if (rho_m == nullptr) {
234 : : // Start precomputing the coefficients
235 : 0 : theta_m = (beta_m + alpha_m) / 2.0 * (1.0 + zeta_m);
236 : 0 : delta_m = (beta_m - alpha_m) / 2.0;
237 : 0 : sigma_m = theta_m / delta_m;
238 : :
239 [ # # # # ]: 0 : rho_m = new double[degree_m + 1];
240 : 0 : rho_m[0] = 1.0 / sigma_m;
241 [ # # ]: 0 : for (unsigned int i = 1; i < degree_m + 1; ++i) {
242 : 0 : rho_m[i] = 1.0 / (2.0 * sigma_m - rho_m[i - 1]);
243 : : }
244 : : } // End of precomputing the coefficients
245 : :
246 [ # # # # ]: 0 : res = r.deepCopy();
247 : :
248 [ # # # # ]: 0 : x_old = r / theta_m;
249 [ # # # # : 0 : A = op_m(r);
# # ]
250 [ # # # # : 0 : x = 2.0 * rho_m[1] / delta_m * (2.0 * r - A / theta_m);
# # # # #
# ]
251 : :
252 [ # # ]: 0 : if (degree_m == 0) {
253 [ # # ]: 0 : return x_old;
254 : : }
255 : :
256 [ # # ]: 0 : if (degree_m == 1) {
257 [ # # ]: 0 : return x;
258 : : }
259 [ # # ]: 0 : for (unsigned int i = 2; i < degree_m + 1; ++i) {
260 [ # # # # : 0 : A = op_m(x);
# # ]
261 [ # # # # : 0 : z = 2.0 / delta_m * (r - A);
# # ]
262 [ # # # # : 0 : res = rho_m[i] * (2 * sigma_m * x - rho_m[i - 1] * x_old + z);
# # # # #
# # # ]
263 [ # # # # ]: 0 : x_old = x.deepCopy();
264 [ # # # # ]: 0 : x = res.deepCopy();
265 : : }
266 [ # # ]: 0 : return res;
267 : 0 : }
268 : :
269 : : protected:
270 : : OperatorF op_m;
271 : : double alpha_m;
272 : : double beta_m;
273 : : double delta_m;
274 : : double theta_m;
275 : : double sigma_m;
276 : : unsigned degree_m;
277 : : double zeta_m;
278 : : double* rho_m = nullptr; // Size is determined at runtime
279 : : };
280 : :
281 : : /*!
282 : : * Richardson preconditioner
283 : : */
284 : : template <typename Field, typename UpperAndLowerF, typename InvDiagF>
285 : : struct richardson_preconditioner : public preconditioner<Field> {
286 : : constexpr static unsigned Dim = Field::dim;
287 : : using mesh_type = typename Field::Mesh_t;
288 : : using layout_type = typename Field::Layout_t;
289 : :
290 : 0 : richardson_preconditioner(UpperAndLowerF&& upper_and_lower, InvDiagF&& inverse_diagonal,
291 : : unsigned innerloops = 5)
292 : : : preconditioner<Field>("Richardson")
293 [ # # # # : 0 : , innerloops_m(innerloops) {
# # ]
294 : 0 : upper_and_lower_m = std::move(upper_and_lower);
295 : 0 : inverse_diagonal_m = std::move(inverse_diagonal);
296 : 0 : }
297 : :
298 : 0 : Field operator()(Field& r) override {
299 : 0 : mesh_type& mesh = r.get_mesh();
300 : 0 : layout_type& layout = r.getLayout();
301 : 0 : Field g(mesh, layout);
302 : :
303 [ # # ]: 0 : g = 0;
304 [ # # ]: 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
305 [ # # # # : 0 : ULg_m = upper_and_lower_m(g);
# # ]
306 [ # # # # ]: 0 : g = r - ULg_m;
307 [ # # # # : 0 : g = inverse_diagonal_m(g) * g;
# # # # ]
308 : : }
309 : 0 : return g;
310 : 0 : }
311 : :
312 : 0 : void init_fields(Field& b) override {
313 : 0 : layout_type& layout = b.getLayout();
314 : 0 : mesh_type& mesh = b.get_mesh();
315 : :
316 [ # # # # ]: 0 : ULg_m = Field(mesh, layout);
317 : 0 : }
318 : :
319 : : protected:
320 : : UpperAndLowerF upper_and_lower_m;
321 : : InvDiagF inverse_diagonal_m;
322 : : unsigned innerloops_m;
323 : : Field ULg_m;
324 : : };
325 : :
326 : : /*!
327 : : * 2-step Gauss-Seidel preconditioner
328 : : */
329 : : template <typename Field, typename LowerF, typename UpperF, typename InvDiagF>
330 : : struct gs_preconditioner : public preconditioner<Field> {
331 : : constexpr static unsigned Dim = Field::dim;
332 : : using mesh_type = typename Field::Mesh_t;
333 : : using layout_type = typename Field::Layout_t;
334 : :
335 : 0 : gs_preconditioner(LowerF&& lower, UpperF&& upper, InvDiagF&& inverse_diagonal,
336 : : unsigned innerloops, unsigned outerloops)
337 : : : preconditioner<Field>("Gauss-Seidel")
338 : 0 : , innerloops_m(innerloops)
339 [ # # # # : 0 : , outerloops_m(outerloops) {
# # # # ]
340 : 0 : lower_m = std::move(lower);
341 : 0 : upper_m = std::move(upper);
342 : 0 : inverse_diagonal_m = std::move(inverse_diagonal);
343 : 0 : }
344 : :
345 : 0 : Field operator()(Field& b) override {
346 : 0 : layout_type& layout = b.getLayout();
347 : 0 : mesh_type& mesh = b.get_mesh();
348 : :
349 : 0 : Field x(mesh, layout);
350 : :
351 [ # # ]: 0 : x = 0; // Initial guess
352 : :
353 [ # # ]: 0 : for (unsigned int k = 0; k < outerloops_m; ++k) {
354 [ # # # # : 0 : UL_m = upper_m(x);
# # ]
355 [ # # # # ]: 0 : r_m = b - UL_m;
356 [ # # ]: 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
357 [ # # # # : 0 : UL_m = lower_m(x);
# # ]
358 [ # # # # ]: 0 : x = r_m - UL_m;
359 [ # # # # : 0 : x = inverse_diagonal_m(x) * x;
# # # # ]
360 : : }
361 [ # # # # : 0 : UL_m = lower_m(x);
# # ]
362 [ # # # # ]: 0 : r_m = b - UL_m;
363 [ # # ]: 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
364 [ # # # # : 0 : UL_m = upper_m(x);
# # ]
365 [ # # # # ]: 0 : x = r_m - UL_m;
366 [ # # # # : 0 : x = inverse_diagonal_m(x) * x;
# # # # ]
367 : : }
368 : : }
369 : 0 : return x;
370 : 0 : }
371 : :
372 : 0 : void init_fields(Field& b) override {
373 : 0 : layout_type& layout = b.getLayout();
374 : 0 : mesh_type& mesh = b.get_mesh();
375 : :
376 [ # # # # ]: 0 : UL_m = Field(mesh, layout);
377 [ # # # # ]: 0 : r_m = Field(mesh, layout);
378 : 0 : }
379 : :
380 : : protected:
381 : : LowerF lower_m;
382 : : UpperF upper_m;
383 : : InvDiagF inverse_diagonal_m;
384 : : unsigned innerloops_m;
385 : : unsigned outerloops_m;
386 : : Field UL_m;
387 : : Field r_m;
388 : : };
389 : :
390 : : /*!
391 : : * Symmetric successive over-relaxation
392 : : */
393 : : template <typename Field, typename LowerF, typename UpperF, typename InvDiagF, typename DiagF>
394 : : struct ssor_preconditioner : public preconditioner<Field> {
395 : : constexpr static unsigned Dim = Field::dim;
396 : : using mesh_type = typename Field::Mesh_t;
397 : : using layout_type = typename Field::Layout_t;
398 : :
399 : 0 : ssor_preconditioner(LowerF&& lower, UpperF&& upper, InvDiagF&& inverse_diagonal,
400 : : DiagF&& diagonal, unsigned innerloops, unsigned outerloops,
401 : : double omega)
402 : : : preconditioner<Field>("ssor")
403 : 0 : , innerloops_m(innerloops)
404 : 0 : , outerloops_m(outerloops)
405 [ # # # # : 0 : , omega_m(omega) {
# # # # ]
406 : 0 : lower_m = std::move(lower);
407 : 0 : upper_m = std::move(upper);
408 : 0 : inverse_diagonal_m = std::move(inverse_diagonal);
409 : 0 : diagonal_m = std::move(diagonal);
410 : 0 : }
411 : :
412 : 0 : Field operator()(Field& b) override {
413 [ # # # # : 0 : static IpplTimings::TimerRef initTimer = IpplTimings::getTimer("SSOR Init");
# # # # ]
414 : 0 : IpplTimings::startTimer(initTimer);
415 : :
416 : : double D;
417 : :
418 : 0 : layout_type& layout = b.getLayout();
419 : 0 : mesh_type& mesh = b.get_mesh();
420 : :
421 : 0 : Field x(mesh, layout);
422 : :
423 [ # # ]: 0 : x = 0; // Initial guess
424 : :
425 [ # # ]: 0 : IpplTimings::stopTimer(initTimer);
426 : :
427 [ # # # # : 0 : static IpplTimings::TimerRef loopTimer = IpplTimings::getTimer("SSOR loop");
# # # # ]
428 [ # # ]: 0 : IpplTimings::startTimer(loopTimer);
429 : :
430 [ # # ]: 0 : for (unsigned int k = 0; k < outerloops_m; ++k) {
431 [ # # # # : 0 : UL_m = upper_m(x);
# # ]
432 [ # # # # ]: 0 : D = diagonal_m(x);
433 [ # # # # : 0 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * D * x;
# # # # #
# ]
434 : :
435 [ # # ]: 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
436 [ # # # # : 0 : UL_m = lower_m(x);
# # ]
437 [ # # # # : 0 : x = r_m - omega_m * UL_m;
# # ]
438 [ # # # # : 0 : x = inverse_diagonal_m(x) * x;
# # # # ]
439 : : }
440 [ # # # # : 0 : UL_m = lower_m(x);
# # ]
441 [ # # # # ]: 0 : D = diagonal_m(x);
442 [ # # # # : 0 : r_m = omega_m * (b - UL_m) + (1.0 - omega_m) * D * x;
# # # # #
# ]
443 [ # # ]: 0 : for (unsigned int j = 0; j < innerloops_m; ++j) {
444 [ # # # # : 0 : UL_m = upper_m(x);
# # ]
445 [ # # # # : 0 : x = r_m - omega_m * UL_m;
# # ]
446 [ # # # # : 0 : x = inverse_diagonal_m(x) * x;
# # # # ]
447 : : }
448 : : }
449 [ # # ]: 0 : IpplTimings::stopTimer(loopTimer);
450 : 0 : return x;
451 : 0 : }
452 : :
453 : 0 : void init_fields(Field& b) override {
454 : 0 : layout_type& layout = b.getLayout();
455 : 0 : mesh_type& mesh = b.get_mesh();
456 : :
457 [ # # # # ]: 0 : UL_m = Field(mesh, layout);
458 [ # # # # ]: 0 : r_m = Field(mesh, layout);
459 : 0 : }
460 : :
461 : : protected:
462 : : LowerF lower_m;
463 : : UpperF upper_m;
464 : : InvDiagF inverse_diagonal_m;
465 : : DiagF diagonal_m;
466 : : unsigned innerloops_m;
467 : : unsigned outerloops_m;
468 : : double omega_m;
469 : : Field UL_m;
470 : : Field r_m;
471 : : };
472 : :
473 : : /*!
474 : : * Computes the largest Eigenvalue of the Functor f
475 : : * @param f Functor
476 : : * @param x_0 initial guess
477 : : * @param max_iter maximum number of iterations
478 : : * @param tol tolerance
479 : : */
480 : : template <typename Field, typename Functor>
481 : : double powermethod(Functor&& f, Field& x_0, unsigned int max_iter = 5000, double tol = 1e-3) {
482 : : unsigned int i = 0;
483 : : using mesh_type = typename Field::Mesh_t;
484 : : using layout_type = typename Field::Layout_t;
485 : : mesh_type& mesh = x_0.get_mesh();
486 : : layout_type& layout = x_0.getLayout();
487 : : Field x_new(mesh, layout);
488 : : Field x_diff(mesh, layout);
489 : : double error = 1.0;
490 : : double lambda = 1.0;
491 : : while (error > tol && i < max_iter) {
492 : : x_new = f(x_0);
493 : : lambda = norm(x_new);
494 : : x_diff = x_new - lambda * x_0;
495 : : error = norm(x_diff);
496 : : x_new = x_new / lambda;
497 : : x_0 = x_new.deepCopy();
498 : : ++i;
499 : : }
500 : : if (i == max_iter) {
501 : : std::cerr << "Powermethod did not converge, lambda_max : " << lambda
502 : : << ", error : " << error << std::endl;
503 : : }
504 : : return lambda;
505 : : }
506 : :
507 : : /*!
508 : : * Computes the smallest Eigenvalue of the Functor f (f must be symmetric positive definite)
509 : : * @param f Functor
510 : : * @param x_0 initial guess
511 : : * @param lambda_max largest Eigenvalue
512 : : * @param max_iter maximum number of iterations
513 : : * @param tol tolerance
514 : : */
515 : : template <typename Field, typename Functor>
516 : : double adapted_powermethod(Functor&& f, Field& x_0, double lambda_max,
517 : : unsigned int max_iter = 5000, double tol = 1e-3) {
518 : : unsigned int i = 0;
519 : : using mesh_type = typename Field::Mesh_t;
520 : : using layout_type = typename Field::Layout_t;
521 : : mesh_type& mesh = x_0.get_mesh();
522 : : layout_type& layout = x_0.getLayout();
523 : : Field x_new(mesh, layout);
524 : : Field x_diff(mesh, layout);
525 : : double error = 1.0;
526 : : double lambda = 1.0;
527 : : while (error > tol && i < max_iter) {
528 : : x_new = f(x_0);
529 : : x_new = x_new - lambda_max * x_0;
530 : : lambda = -norm(x_new); // We know that lambda < 0;
531 : : x_diff = x_new - lambda * x_0;
532 : : error = norm(x_diff);
533 : : x_new = x_new / -lambda;
534 : : x_0 = x_new.deepCopy();
535 : : ++i;
536 : : }
537 : : lambda = lambda + lambda_max;
538 : : if (i == max_iter) {
539 : : std::cerr << "Powermethod did not converge, lambda_min : " << lambda
540 : : << ", error : " << error << std::endl;
541 : : }
542 : : return lambda;
543 : : }
544 : :
545 : : /*
546 : : // Use the powermethod to compute the eigenvalues if no analytical solution is known
547 : : beta = powermethod(std::move(op_m), x_0);
548 : : // Trick for computing the smallest Eigenvalue of SPD Matrix
549 : : alpha = adapted_powermethod(std::move(op_m), x_0, beta);
550 : : */
551 : :
552 : : } // namespace ippl
553 : :
554 : : #endif // IPPL_PRECONDITIONER_H
|
