Alien  1.3.0
Developer documentation
Loading...
Searching...
No Matches
BiCGStab.h
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7
8#pragma once
9
10#include <ostream>
11#include <vector>
12
13namespace Alien
14{
15
16template <typename AlgebraT>
17class BiCGStab
18{
19 public:
20 // clang-format off
21 typedef AlgebraT AlgebraType;
22 typedef typename AlgebraType::Matrix MatrixType;
23 typedef typename AlgebraType::Vector VectorType;
24 typedef typename MatrixType::ValueType ValueType;
25 typedef typename AlgebraType::FutureType FutureType;
26 // clang-format on
27
28 class Iteration
29 {
30 public:
31 Iteration(AlgebraType& algebra,
32 VectorType const& b,
33 ValueType tol,
34 int max_iter,
35 ITraceMng* trace_mng = nullptr)
36 : m_algebra(algebra)
37 , m_value(0)
38 , m_f_value(m_value)
39 , m_max_iteration(max_iter)
40 , m_tol(tol)
41 , m_iter(0)
42 , m_trace_mng(trace_mng)
43 {
44 m_algebra.dot(b, b, m_f_value);
45 m_nrm2_b = m_f_value.get();
46 //m_nrm2_b = m_algebra.dot(b, b) ;
47 if (m_trace_mng)
48 m_trace_mng->info() << "STOP CRITERIA NORME B = " << m_nrm2_b;
49 m_criteria_value = m_tol * m_tol * m_nrm2_b;
50 m_sqrt_nrm2_b = std::sqrt(m_nrm2_b);
51 m_value = m_criteria_value + 1;
52 if (m_nrm2_b == 0)
53 m_status = true;
54 else
55 m_status = false;
56 }
57
58 virtual ~Iteration()
59 {}
60
61 bool nullRhs() const
62 {
63 return m_nrm2_b == 0.;
64 }
65
66 bool first() const
67 {
68 return m_iter == 0;
69 }
70
71 bool stop(VectorType const& r)
72 {
73 if (m_iter >= m_max_iteration)
74 return true;
75 m_algebra.dot(r, r, m_f_value);
76 m_status = m_f_value.get() < m_criteria_value;
77 return m_status;
78 }
79
80 void operator++()
81 {
82 if (m_trace_mng)
83 m_trace_mng->info() << "iteration (" << m_iter << ") criteria = " << getValue();
84 ++m_iter;
85 }
86
87 ValueType getValue() const
88 {
89 if (m_sqrt_nrm2_b == 0)
90 return 0.;
91 else
92 return std::sqrt(m_value) / m_sqrt_nrm2_b;
93 }
94
95 int operator()() const
96 {
97 return m_iter;
98 }
99
100 bool getStatus() const
101 {
102 return m_status;
103 }
104
105 private:
106 // clang-format off
107 AlgebraType& m_algebra;
108 int m_max_iteration = 0;
109 ValueType m_tol = 0.;
110 int m_iter = 0;
111 ValueType m_value = 0.;
112 FutureType m_f_value;
113 ValueType m_criteria_value = 0.;
114 ValueType m_value_init = 0.;
115 ValueType m_nrm2_b = 0.;
116 ValueType m_sqrt_nrm2_b = 0.;
117 bool m_status = false;
118 ITraceMng* m_trace_mng = nullptr;
119 // clang-format on
120 };
121
122 BiCGStab(AlgebraType& algebra, ITraceMng* trace_mng = nullptr)
123 : m_algebra(algebra)
124 , m_trace_mng(trace_mng)
125 {}
126
127 virtual ~BiCGStab()
128 {}
129
130 void setOutputLevel(int level)
131 {
132 m_output_level = level;
133 }
134
135 template <typename PrecondT, typename iterT>
136 int solve(PrecondT& precond, iterT& iter, MatrixType const& A,
137 VectorType const& b, VectorType& x)
138 {
139 if (iter.nullRhs())
140 return 0;
141 ValueType rho(0), rho1(0), alpha(0), beta(0), omega(0);
142 VectorType p, phat, s, shat, t, v, r, r0;
143
144 m_algebra.allocate(AlgebraType::resource(A), p, phat, s, shat, t, v, r, r0);
145
146 // SEQ0
147 // r = b - A * x;
148 m_algebra.copy(b, r);
149 m_algebra.mult(A, x, r0);
150 m_algebra.axpy(-1., r0, r);
151
152 // rtilde = r
153 m_algebra.copy(r, r0);
154 m_algebra.copy(r, p);
155 rho1 = m_algebra.dot(r, r0);
156 if (m_output_level > 2)
157 _print(0, "Seq 0", "rho1", rho1);
158
159 /*
160 phat = solve(M, p);
161 v = A * phat;
162 gamma = dot(r0, v);
163 alpha = rho_1 / gamma;
164 s = r - alpha * v;
165 */
166 // SEQ1
167 m_algebra.exec(precond, p, phat);
168 m_algebra.mult(A, phat, v);
169 alpha = m_algebra.dot(v, r0);
170 if (alpha == 0)
171 throw typename AlgebraType::NullValueException("alpha");
172 alpha = rho1 / alpha;
173 m_algebra.copy(r, s);
174 m_algebra.axpy(-alpha, v, s);
175 if (m_output_level > 2)
176 _print(0, "Seq 1", "alpha", alpha);
177
178 if (iter.stop(s)) {
179 ++iter;
180 m_algebra.axpy(alpha, phat, x);
181 m_algebra.free(p, phat, s, shat, t, v, r, r0);
182 return 0;
183 }
184
185 // SEQ 2
186 m_algebra.exec(precond, s, shat);
187 m_algebra.mult(A, shat, t);
188 omega = m_algebra.dot(t, s);
189 beta = m_algebra.dot(t, t);
190
191 if (beta == 0) {
192 if (iter.stop(r)) {
193 ++iter;
194 m_algebra.axpy(alpha, phat, x);
195 m_algebra.free(p, phat, s, shat, t, v, r, r0);
196 return 0;
197 }
198 else
199 throw typename AlgebraType::NullValueException("beta");
200 }
201 omega = omega / beta;
202 if (m_output_level > 2)
203 _print(iter(), "Seq 2", "beta", beta, "alpha", alpha, "rho1", rho1);
204
205 // SEQ 3
206 m_algebra.axpy(omega, shat, x);
207 m_algebra.axpy(alpha, phat, x);
208 m_algebra.copy(s, r);
209 m_algebra.axpy(-omega, t, r);
210
211 rho = rho1;
212 ++iter;
213 if (m_output_level > 2)
214 _print(iter(), "Seq 3", "beta", beta, "alpha", alpha, "rho1", rho1);
215
216 while (!iter.stop(r)) {
217 //SEQ4
218 rho1 = m_algebra.dot(r, r0);
219 beta = (rho1 / rho) * (alpha / omega);
220 m_algebra.axpy(-omega, v, p);
221 m_algebra.scal(beta, p);
222 m_algebra.axpy(1., r, p);
223 if (m_output_level > 2)
224 _print(iter(), "Seq 4", "beta", beta, "alpha", alpha, "rho1", rho1);
225
226 if (rho == 0)
227 throw typename AlgebraType::NullValueException("rho");
228
229 //m_algebra.process (seq1);
230 m_algebra.exec(precond, p, phat);
231 m_algebra.mult(A, phat, v);
232 alpha = m_algebra.dot(v, r0);
233 if (alpha == 0)
234 throw typename AlgebraType::NullValueException("alpha");
235 else
236 alpha = rho1 / alpha;
237
238 m_algebra.copy(r, s);
239 m_algebra.axpy(-alpha, v, s);
240
241 if (m_output_level > 2)
242 _print(iter(), "Seq 1", "alpha", alpha);
243 if (iter.stop(s)) {
244 m_algebra.axpy(alpha, phat, x);
245 m_algebra.free(p, phat, s, shat, t, v, r, r0);
246 return 0;
247 }
248
249 //m_algebra.process (seq2);
250 m_algebra.exec(precond, s, shat);
251 m_algebra.mult(A, shat, t);
252 omega = m_algebra.dot(t, s);
253 beta = m_algebra.dot(t, t);
254
255 if (m_output_level > 2)
256 _print(iter(), "Seq 2", "beta", beta, "alpha", alpha, "rho1", rho1, "omega", omega);
257 if (beta == 0) {
258 if (iter.stop(s)) {
259 m_algebra.axpy(alpha, phat, x);
260 m_algebra.free(p, phat, s, shat, t, v, r, r0);
261 return 0;
262 }
263 throw typename AlgebraType::NullValueException("beta");
264 }
265 else
266 omega = omega / beta;
267
268 //m_algebra.process (seq3);
269 m_algebra.axpy(omega, shat, x);
270 m_algebra.axpy(alpha, phat, x);
271 m_algebra.copy(s, r);
272 m_algebra.axpy(-omega, t, r);
273
274 rho = rho1;
275
276 ++iter;
277 if (m_output_level > 2)
278 _print(iter(), "end loop", "beta", beta, "alpha", alpha, "rho", rho);
279 }
280
281 m_algebra.free(p, phat, s, shat, t, v, r, r0);
282
283 return 0;
284 }
285
286 template <typename PrecondT, typename iterT>
287 int solve2(PrecondT& precond, iterT& iter, MatrixType const& A,
288 VectorType const& b, VectorType& x)
289 {
290 if (iter.nullRhs())
291 return 0;
292 // clang-format off
293 ValueType rho (0), rho1 (0), alpha (0), beta (0), gamma (0), omega (0);
294 FutureType frho(rho), frho1(rho1), falpha(alpha), fbeta(beta), fgamma(gamma), fomega(omega) ;
295 VectorType p, phat, s, shat, t, v, r, r0;
296 // clang-format on
297
298 m_algebra.allocate(AlgebraType::resource(A), p, phat, s, shat, t, v, r, r0);
299
300 // SEQ0
301 // r = b - A * x;
302 m_algebra.copy(b, r);
303 m_algebra.mult(A, x, r0);
304 m_algebra.axpy(-1., r0, r);
305
306 // rtilde = r
307 m_algebra.copy(r, r0);
308 m_algebra.copy(r, p);
309 m_algebra.dot(r, r0, frho1);
310 if (m_output_level > 2)
311 _print(0, "Seq 0", "rho1", frho1.get());
312
313 /*
314 phat = solve(M, p);
315 v = A * phat;
316 gamma = dot(r0, v);
317 alpha = rho_1 / gamma;
318 s = r - alpha * v;
319 */
320 // SEQ1
321 m_algebra.exec(precond, p, phat);
322 m_algebra.mult(A, phat, v);
323 m_algebra.dot(v, r0, falpha);
324 if (falpha.get() == 0)
325 throw typename AlgebraType::NullValueException("alpha");
326 alpha = frho1.get() / alpha;
327
328 m_algebra.copy(r, s);
329 m_algebra.axpy(-alpha, v, s);
330 if (m_output_level > 2)
331 _print(0, "Seq 1", "alpha", alpha);
332
333 if (iter.stop(s)) {
334 ++iter;
335 m_algebra.axpy(alpha, phat, x);
336 m_algebra.free(p, phat, s, shat, t, v, r, r0);
337 return 0;
338 }
339
340 // SEQ 2
341 m_algebra.exec(precond, s, shat);
342 m_algebra.mult(A, shat, t);
343 m_algebra.dot(t, s, fomega);
344 m_algebra.dot(t, t, fbeta);
345 if (fbeta.get() == 0) {
346 if (iter.stop(r)) {
347 ++iter;
348 m_algebra.axpy(alpha, phat, x);
349 m_algebra.free(p, phat, s, shat, t, v, r, r0);
350 return 0;
351 }
352 else
353 throw typename AlgebraType::NullValueException("beta");
354 }
355 omega = fomega.get() / beta;
356 if (m_output_level > 2)
357 _print(iter(), "Seq 2", "beta", beta, "alpha", alpha, "rho1", rho1);
358
359 // SEQ 3
360 m_algebra.axpy(omega, shat, x);
361 m_algebra.axpy(alpha, phat, x);
362 m_algebra.copy(s, r);
363 m_algebra.axpy(-omega, t, r);
364
365 rho = rho1;
366 ++iter;
367 if (m_output_level > 2)
368 _print(iter(), "Seq 3", "beta", beta, "alpha", alpha, "rho1", rho1);
369
370 while (!iter.stop(r)) {
371 //SEQ4
372 /*
373 beta = (rho_1 / rho_2) * (alpha / omega);
374 p = r + beta * (p - omega * v);
375 */
376 m_algebra.dot(r, r0, frho1);
377 beta = (frho1.get() / rho) * (alpha / omega);
378 m_algebra.axpy(-omega, v, p);
379 m_algebra.scal(beta, p);
380 m_algebra.axpy(1., r, p);
381 if (m_output_level > 2)
382 _print(iter(), "Seq 4", "beta", beta, "alpha", alpha, "rho1", rho1);
383
384 if (rho == 0)
385 throw typename AlgebraType::NullValueException("rho");
386
387 //m_algebra.process (seq1);
388 m_algebra.exec(precond, p, phat);
389 m_algebra.mult(A, phat, v);
390 m_algebra.dot(v, r0, falpha);
391 if (falpha.get() == 0)
392 throw typename AlgebraType::NullValueException("alpha");
393 else
394 alpha = rho1 / alpha;
395
396 m_algebra.copy(r, s);
397 m_algebra.axpy(-alpha, v, s);
398 if (m_output_level > 2)
399 _print(iter(), "Seq 1", "alpha", alpha);
400
401 if (iter.stop(s)) {
402 m_algebra.axpy(alpha, phat, x);
403 m_algebra.free(p, phat, s, shat, t, v, r, r0);
404 return 0;
405 }
406
407 //m_algebra.process (seq2);
408 m_algebra.exec(precond, s, shat);
409 m_algebra.mult(A, shat, t);
410 m_algebra.dot(t, s, fomega);
411 m_algebra.dot(t, t, fbeta);
412 if (m_output_level > 2)
413 _print(iter(), "Seq 2", "beta", beta, "alpha", alpha, "rho1", rho1, "omega", omega);
414 if (fbeta.get() == 0) {
415 if (iter.stop(s)) {
416 m_algebra.axpy(alpha, phat, x);
417 m_algebra.free(p, phat, s, shat, t, v, r, r0);
418 return 0;
419 }
420 throw typename AlgebraType::NullValueException("beta");
421 }
422 else
423 omega = fomega.get() / beta;
424
425 //m_algebra.process (seq3);
426 m_algebra.axpy(omega, shat, x);
427 m_algebra.axpy(alpha, phat, x);
428 m_algebra.copy(s, r);
429 m_algebra.axpy(-omega, t, r);
430
431 rho = rho1;
432
433 ++iter;
434 if (m_output_level > 2)
435 _print(iter(), "end loop", "beta", beta, "alpha", alpha, "rho", rho);
436 }
437
438 m_algebra.free(p, phat, s, shat, t, v, r, r0);
439
440 return 0;
441 }
442
443 private:
444 void
445 _print(int iter, std::string const& msg, std::string const& label0,
446 ValueType value0)
447 {
448 if (m_trace_mng) {
449 m_trace_mng->info() << msg;
450 m_trace_mng->info() << "Iterate: " << iter << " " << label0 << " "
451 << value0;
452 }
453 }
454
455 void
456 _print(int iter, std::string const& msg, std::string const& label0,
457 ValueType value0, std::string const& label1, ValueType value1)
458 {
459 if (m_trace_mng) {
460 _print(iter, msg, label0, value0);
461 m_trace_mng->info() << "Iterate: " << iter << " " << label1 << " "
462 << value1;
463 }
464 }
465
466 void
467 _print(int iter, std::string const& msg, std::string const& label0,
468 ValueType value0, std::string const& label1, ValueType value1,
469 std::string const& label2, ValueType value2)
470 {
471 if (m_trace_mng) {
472 _print(iter, msg, label0, value0, label1, value1);
473 m_trace_mng->info() << "Iterate: " << iter << " " << label2 << " "
474 << value2;
475 }
476 }
477
478 void
479 _print(int iter, std::string const& msg, std::string const& label0,
480 ValueType value0, std::string const& label1, ValueType value1,
481 std::string const& label2, ValueType value2,
482 std::string const& label3, ValueType value3)
483 {
484 if (m_trace_mng) {
485 _print(iter, msg, label0, value0, label1, value1, label2, value2);
486 m_trace_mng->info() << "Iterate: " << iter << " " << label3 << " "
487 << value3;
488 }
489 }
490
491 AlgebraType& m_algebra;
492 ITraceMng* m_trace_mng = nullptr;
493 int m_output_level = 0;
494};
495} // namespace Alien
Alien
-- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature --
Definition BackEnd.h:17