Arcane  4.2.3.0
Developer documentation
Loading...
Searching...
No Matches
CPRPreconditioner.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/* CPRPreconditioner.h (C) 2000-2026 */
9/* */
10/* Two stage preconditioner of the Constrained Pressure Residual type. */
11/*---------------------------------------------------------------------------*/
12#ifndef ARCCORE_ALINA_CPRPRECONDITIONER_H
13#define ARCCORE_ALINA_CPRPRECONDITIONER_H
14/*---------------------------------------------------------------------------*/
15/*---------------------------------------------------------------------------*/
16/*
17 * This file is based on the work on AMGCL library (version march 2026)
18 * which can be found at https://github.com/ddemidov/amgcl.
19 *
20 * Copyright (c) 2012-2022 Denis Demidov <dennis.demidov@gmail.com>
21 * SPDX-License-Identifier: MIT
22 */
23/*---------------------------------------------------------------------------*/
24/*---------------------------------------------------------------------------*/
25
26#include <vector>
27#include <memory>
28#include <cassert>
29
30#include "arccore/alina/BuiltinBackend.h"
31#include "arccore/alina/AlinaUtils.h"
32
33/*---------------------------------------------------------------------------*/
34/*---------------------------------------------------------------------------*/
35
36namespace Arcane::Alina::preconditioner
37{
38
39/*---------------------------------------------------------------------------*/
40/*---------------------------------------------------------------------------*/
44template <class PPrecond, class SPrecond>
45class CPRPreconditioner
46{
47 static_assert(math::static_rows<typename PPrecond::backend_type::value_type>::value == 1,
48 "Pressure backend should have scalar value type!");
49
50 static_assert(backend::backends_compatible<
51 typename PPrecond::backend_type,
52 typename SPrecond::backend_type>::value,
53 "Backends for pressure and flow preconditioners should coincide!");
54
55 public:
56
57 typedef typename SPrecond::backend_type backend_type;
58 typedef typename PPrecond::backend_type backend_type_p;
59
60 typedef typename backend_type::value_type value_type;
61 typedef typename backend_type::matrix matrix;
62 typedef typename backend_type::vector vector;
63 typedef typename backend_type_p::value_type value_type_p;
64 typedef typename backend_type_p::matrix matrix_p;
65 typedef typename backend_type_p::vector vector_p;
66
67 typedef typename backend_type::params backend_params;
68
69 typedef typename BuiltinBackend<value_type>::matrix build_matrix;
70 typedef typename BuiltinBackend<value_type_p>::matrix build_matrix_p;
71
72 typedef typename math::scalar_of<value_type>::type scalar_type;
73
74 struct params
75 {
76 typedef typename PPrecond::params pprecond_params;
77 typedef typename SPrecond::params sprecond_params;
78
79 pprecond_params pprecond;
80 sprecond_params sprecond;
81
82 Int32 block_size = math::static_rows<value_type>::value == 1 ? 2 : math::static_rows<value_type>::value;
83
84 Int32 active_rows = 0;
85
86 params() = default;
87
88 params(const PropertyTree& p)
89 : ARCCORE_ALINA_PARAMS_IMPORT_CHILD(p, pprecond)
90 , ARCCORE_ALINA_PARAMS_IMPORT_CHILD(p, sprecond)
91 , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, block_size)
92 , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, active_rows)
93 {
94 p.check_params({ "pprecond", "sprecond", "block_size", "active_rows" });
95 }
96
97 void get(PropertyTree& p, const std::string& path = "") const
98 {
99 ARCCORE_ALINA_PARAMS_EXPORT_CHILD(p, path, pprecond);
100 ARCCORE_ALINA_PARAMS_EXPORT_CHILD(p, path, sprecond);
101 ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, block_size);
102 ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, active_rows);
103 }
104 };
105
106 template <class Matrix>
107 CPRPreconditioner(const AlinaContext& context,
108 const Matrix& K,
109 const params& prm = params(),
110 const backend_params& bprm = backend_params())
111 : prm(prm)
112 , n(backend::nbRow(K))
113 {
114 init(context, std::make_shared<build_matrix>(K), bprm,
115 std::integral_constant<bool, math::static_rows<value_type>::value == 1>());
116 }
117
118 CPRPreconditioner(const AlinaContext& context,
119 std::shared_ptr<build_matrix> K,
120 const params& prm = params(),
121 const backend_params& bprm = backend_params())
122 : prm(prm)
123 , n(backend::nbRow(*K))
124 {
125 init(context, K, bprm,
126 std::integral_constant<bool, math::static_rows<value_type>::value == 1>());
127 }
128
129 template <class Vec1, class Vec2>
130 void apply(const Vec1& rhs, Vec2&& x) const
131 {
132 const auto one = math::identity<scalar_type>();
133 const auto zero = math::zero<scalar_type>();
134
135 ARCCORE_ALINA_TIC("sprecond");
136 S->apply(rhs, x);
137 ARCCORE_ALINA_TOC("sprecond");
138 backend::residual(rhs, S->system_matrix(), x, *rs);
139
140 backend::spmv(one, *Fpp, *rs, zero, *rp);
141 ARCCORE_ALINA_TIC("pprecond");
142 P->apply(*rp, *xp);
143 ARCCORE_ALINA_TOC("pprecond");
144
145 backend::spmv(one, *Scatter, *xp, one, x);
146 }
147
148 std::shared_ptr<matrix> system_matrix_ptr() const
149 {
150 return S->system_matrix_ptr();
151 }
152
153 const matrix& system_matrix() const
154 {
155 return S->system_matrix();
156 }
157
158 template <class Matrix>
159 void partial_update(const Matrix& K,
160 bool update_transfer_ops = true,
161 const backend_params& bprm = backend_params())
162 {
163 auto K_ptr = std::make_shared<build_matrix>(K);
164 // Update global preconditioner
165 S = std::make_shared<SPrecond>(K_ptr, prm.sprecond, bprm);
166 if (update_transfer_ops) {
167 // Update transfer operator Fpp
168 update_transfer(
169 K_ptr,
170 bprm,
171 std::integral_constant<bool, math::static_rows<value_type>::value == 1>());
172 }
173 }
174
175 params prm;
176
177 private:
178
179 size_t n, np;
180
181 std::shared_ptr<PPrecond> P;
182 std::shared_ptr<SPrecond> S;
183
184 std::shared_ptr<matrix_p> Fpp, Scatter;
185 std::shared_ptr<vector> rs;
186 std::shared_ptr<vector_p> rp, xp;
187
188 // Returns pressure transfer operator fpp and (optionally)
189 // partially constructed pressure system matrix App.
190 std::tuple<std::shared_ptr<build_matrix_p>, std::shared_ptr<build_matrix_p>>
191 first_scalar_pass(std::shared_ptr<build_matrix> K, bool get_app = true)
192 {
193 typedef typename backend::row_iterator<build_matrix>::type row_iterator;
194 const int B = prm.block_size;
195 const ptrdiff_t N = (prm.active_rows ? prm.active_rows : n);
196
197 np = N / B;
198
199 auto fpp = std::make_shared<build_matrix_p>();
200 fpp->set_size(np, N);
201 fpp->set_nonzeros(N);
202 fpp->ptr[0] = 0;
203
204 std::shared_ptr<build_matrix_p> App;
205
206 if (get_app) {
207 App = std::make_shared<build_matrix>();
208 App->set_size(np, np, true);
209 }
210
211 // Get the pressure matrix nonzero pattern,
212 // extract and invert block diagonals.
213 arccoreParallelFor(0, np, ForLoopRunInfo{}, [&](Int32 begin, Int32 size) {
214 std::vector<row_iterator> k;
215 k.reserve(B);
216 multi_array<scalar_type, 2> v(B, B);
217
218 for (ptrdiff_t ip = begin; ip < (begin + size); ++ip) {
219 ptrdiff_t ik = ip * B;
220 bool done = true;
221 ptrdiff_t cur_col = 0;
222
223 k.clear();
224 for (int i = 0; i < B; ++i) {
225 k.push_back(backend::row_begin(*K, ik + i));
226
227 if (k.back() && k.back().col() < N) {
228 ptrdiff_t col = k.back().col() / B;
229 if (done) {
230 cur_col = col;
231 done = false;
232 }
233 else {
234 cur_col = std::min(cur_col, col);
235 }
236 }
237
238 fpp->col[ik + i] = ik + i;
239 }
240 fpp->ptr[ip + 1] = ik + B;
241
242 while (!done) {
243 if (get_app)
244 ++App->ptr[ip + 1];
245
246 ptrdiff_t end = (cur_col + 1) * B;
247
248 if (cur_col == ip) {
249 // This is diagonal block.
250 // Capture its (transposed) value,
251 // invert it and put the relevant row into fpp.
252 for (int i = 0; i < B; ++i)
253 for (int j = 0; j < B; ++j)
254 v(i, j) = 0;
255
256 for (int i = 0; i < B; ++i)
257 for (; k[i] && k[i].col() < end; ++k[i])
258 v(k[i].col() % B, i) = k[i].value();
259
260 invert(v.data(), &fpp->val[ik]);
261
262 if (!get_app)
263 break;
264 }
265 else {
266 // This is off-diagonal block.
267 // Just skip it.
268 for (int i = 0; i < B; ++i)
269 while (k[i] && k[i].col() < end)
270 ++k[i];
271 }
272
273 // Get next column number.
274 done = true;
275 for (int i = 0; i < B; ++i) {
276 if (k[i] && k[i].col() < N) {
277 ptrdiff_t col = k[i].col() / B;
278 if (done) {
279 cur_col = col;
280 done = false;
281 }
282 else {
283 cur_col = std::min(cur_col, col);
284 }
285 }
286 }
287 }
288 }
289 });
290
291 if (get_app)
292 App->set_nonzeros(App->scan_row_sizes());
293
294 return std::make_tuple(fpp, App);
295 }
296
297 // The system matrix has scalar values
298 void init(const AlinaContext& context, std::shared_ptr<build_matrix> K,
299 const backend_params bprm, std::true_type)
300 {
301 typedef typename backend::row_iterator<build_matrix>::type row_iterator;
302 const int B = prm.block_size;
303 const ptrdiff_t N = (prm.active_rows ? prm.active_rows : n);
304
305 np = N / B;
306
307 std::shared_ptr<build_matrix_p> fpp, App;
308 std::tie(fpp, App) = first_scalar_pass(K);
309
310 auto scatter = std::make_shared<build_matrix_p>();
311 scatter->set_size(n, np);
312 scatter->set_nonzeros(np);
313 scatter->ptr[0] = 0;
314
315 arccoreParallelFor(0, np, ForLoopRunInfo{}, [&](Int32 begin, Int32 size) {
316 std::vector<row_iterator> k;
317 k.reserve(B);
318
319 for (ptrdiff_t ip = begin; ip < (begin + size); ++ip) {
320 ptrdiff_t ik = ip * B;
321 ptrdiff_t head = App->ptr[ip];
322 bool done = true;
323 ptrdiff_t cur_col;
324
325 value_type_p* d = &fpp->val[ik];
326
327 k.clear();
328 for (int i = 0; i < B; ++i) {
329 k.push_back(backend::row_begin(*K, ik + i));
330
331 if (k.back() && k.back().col() < N) {
332 ptrdiff_t col = k.back().col() / B;
333 if (done) {
334 cur_col = col;
335 done = false;
336 }
337 else {
338 cur_col = std::min(cur_col, col);
339 }
340 }
341 }
342
343 while (!done) {
344 ptrdiff_t end = (cur_col + 1) * B;
345 value_type_p app = 0;
346
347 for (int i = 0; i < B; ++i) {
348 for (; k[i] && k[i].col() < end; ++k[i]) {
349 if (k[i].col() % B == 0) {
350 app += d[i] * k[i].value();
351 }
352 }
353 }
354
355 App->col[head] = cur_col;
356 App->val[head] = app;
357 ++head;
358
359 // Get next column number.
360 done = true;
361 for (int i = 0; i < B; ++i) {
362 if (k[i] && k[i].col() < N) {
363 ptrdiff_t col = k[i].col() / B;
364 if (done) {
365 cur_col = col;
366 done = false;
367 }
368 else {
369 cur_col = std::min(cur_col, col);
370 }
371 }
372 }
373 }
374
375 scatter->col[ip] = ip;
376 scatter->val[ip] = math::identity<value_type_p>();
377
378 ptrdiff_t nnz = ip;
379 for (int i = 0; i < B; ++i) {
380 if (i == 0)
381 ++nnz;
382 scatter->ptr[ik + i + 1] = nnz;
383 }
384 }
385 });
386
387 for (size_t i = N; i < n; ++i)
388 scatter->ptr[i + 1] = scatter->ptr[i];
389
390 ARCCORE_ALINA_TIC("pprecond");
391 P = std::make_shared<PPrecond>(context, App, prm.pprecond, bprm);
392 ARCCORE_ALINA_TOC("pprecond");
393 ARCCORE_ALINA_TIC("sprecond");
394 S = std::make_shared<SPrecond>(context, K, prm.sprecond, bprm);
395 ARCCORE_ALINA_TOC("sprecond");
396
397 Fpp = backend_type_p::copy_matrix(fpp, bprm);
398 Scatter = backend_type_p::copy_matrix(scatter, bprm);
399
400 rp = backend_type_p::create_vector(np, bprm);
401 xp = backend_type_p::create_vector(np, bprm);
402 rs = backend_type::create_vector(n, bprm);
403 }
404
405 void update_transfer(std::shared_ptr<build_matrix> K, const backend_params bprm, std::true_type)
406 {
407 auto fpp = std::get<0>(first_scalar_pass(K, /*get_app*/ false));
408 Fpp = backend_type_p::copy_matrix(fpp, bprm);
409 }
410
411 // The system matrix has block values
412 void init(const AlinaContext& context, std::shared_ptr<build_matrix> K,
413 const backend_params bprm, std::false_type)
414 {
415 const int B = math::static_rows<value_type>::value;
416 const ptrdiff_t N = (prm.active_rows ? prm.active_rows : n);
417
418 np = N;
419
420 auto fpp = std::make_shared<build_matrix_p>();
421 fpp->set_size(np, np * B);
422 fpp->set_nonzeros(np * B);
423 fpp->ptr[0] = 0;
424
425 auto scatter = std::make_shared<build_matrix_p>();
426 scatter->set_size(np * B, np);
427 scatter->set_nonzeros(np);
428 scatter->ptr[0] = 0;
429
430 auto App = std::make_shared<build_matrix_p>();
431 App->set_size(np, np, true);
432 App->set_nonzeros(K->nbNonZero());
433 App->ptr[0] = 0;
434
435 arccoreParallelFor(0, np, ForLoopRunInfo{}, [&](Int32 begin, Int32 size) {
436 for (ptrdiff_t i = begin; i < (begin + size); ++i) {
437 ptrdiff_t ik = i * B;
438 for (int k = 0; k < B; ++k, ++ik) {
439 fpp->col[ik] = ik;
440 scatter->ptr[ik + 1] = i + 1;
441 }
442
443 fpp->ptr[i + 1] = ik;
444 scatter->col[i] = i;
445 scatter->val[i] = math::identity<value_type_p>();
446
447 ptrdiff_t row_beg = K->ptr[i];
448 ptrdiff_t row_end = K->ptr[i + 1];
449 App->ptr[i + 1] = row_end;
450
451 // Extract and invert block diagonals
452 value_type_p* d = &fpp->val[i * B];
453 for (ptrdiff_t j = row_beg; j < row_end; ++j) {
454 if (K->col[j] == i) {
455 value_type v = math::adjoint(K->val[j]);
456 invert(v.data(), d);
457 break;
458 }
459 }
460
461 for (ptrdiff_t j = row_beg; j < row_end; ++j) {
462 value_type_p app = 0;
463 for (int k = 0; k < B; ++k)
464 app += d[k] * K->val[j](k, 0);
465
466 App->col[j] = K->col[j];
467 App->val[j] = app;
468 }
469 }
470 });
471
472 ARCCORE_ALINA_TIC("pprecond");
473 P = std::make_shared<PPrecond>(context, App, prm.pprecond, bprm);
474 ARCCORE_ALINA_TOC("pprecond");
475 ARCCORE_ALINA_TIC("sprecond");
476 S = std::make_shared<SPrecond>(context, K, prm.sprecond, bprm);
477 ARCCORE_ALINA_TOC("sprecond");
478
479 Fpp = backend_type_p::copy_matrix(fpp, bprm);
480 Scatter = backend_type_p::copy_matrix(scatter, bprm);
481
482 rp = backend_type_p::create_vector(np, bprm);
483 xp = backend_type_p::create_vector(np, bprm);
484 rs = backend_type::create_vector(n, bprm);
485 }
486
487 void update_transfer(std::shared_ptr<build_matrix> K, const backend_params bprm, std::false_type)
488 {
489 const int B = math::static_rows<value_type>::value;
490 const ptrdiff_t N = (prm.active_rows ? prm.active_rows : n);
491
492 np = N;
493
494 auto fpp = std::make_shared<build_matrix_p>();
495 fpp->set_size(np, np * B);
496 fpp->set_nonzeros(np * B);
497 fpp->ptr[0] = 0;
498
499 arccoreParallelFor(0, np, ForLoopRunInfo{}, [&](Int32 begin, Int32 size) {
500 for (ptrdiff_t i = begin; i < (begin + size); ++i) {
501 ptrdiff_t ik = i * B;
502 for (int k = 0; k < B; ++k, ++ik) {
503 fpp->col[ik] = ik;
504 }
505
506 fpp->ptr[i + 1] = ik;
507
508 ptrdiff_t row_beg = K->ptr[i];
509 ptrdiff_t row_end = K->ptr[i + 1];
510
511 // Extract and invert block diagonals
512 value_type_p* d = &fpp->val[i * B];
513 for (ptrdiff_t j = row_beg; j < row_end; ++j) {
514 if (K->col[j] == i) {
515 value_type v = math::adjoint(K->val[j]);
516 invert(v.data(), d);
517 break;
518 }
519 }
520 }
521 });
522 Fpp = backend_type_p::copy_matrix(fpp, bprm);
523 }
524
525 // Inverts dense matrix A;
526 // Returns the first column of the inverted matrix.
527 void invert(scalar_type* A, value_type_p* y)
528 {
529 const int B = math::static_rows<value_type>::value == 1 ? prm.block_size : math::static_rows<value_type>::value;
530
531 // Perform LU-factorization of A in-place
532 for (int k = 0; k < B; ++k) {
533 scalar_type d = A[k * B + k];
534 assert(!math::is_zero(d));
535 for (int i = k + 1; i < B; ++i) {
536 A[i * B + k] /= d;
537 for (int j = k + 1; j < B; ++j)
538 A[i * B + j] -= A[i * B + k] * A[k * B + j];
539 }
540 }
541
542 // Invert unit vector in-place.
543 // Lower triangular solve:
544 for (int i = 0; i < B; ++i) {
545 value_type_p b = static_cast<value_type_p>(i == 0);
546 for (int j = 0; j < i; ++j)
547 b -= A[i * B + j] * y[j];
548 y[i] = b;
549 }
550
551 // Upper triangular solve:
552 for (int i = B; i-- > 0;) {
553 for (int j = i + 1; j < B; ++j)
554 y[i] -= A[i * B + j] * y[j];
555 y[i] /= A[i * B + i];
556 }
557 }
558
559 friend std::ostream& operator<<(std::ostream& os, const CPRPreconditioner& p)
560 {
561 os << "CPR (two-stage preconditioner)\n"
562 "### Pressure preconditioner:\n"
563 << *p.P << "\n"
564 "### Global preconditioner:\n"
565 << *p.S << std::endl;
566 return os;
567 }
568};
569
570/*---------------------------------------------------------------------------*/
571/*---------------------------------------------------------------------------*/
572
573} // namespace Arcane::Alina::preconditioner
574
575/*---------------------------------------------------------------------------*/
576/*---------------------------------------------------------------------------*/
577
578#endif
Arcane::Alina::PropertyTree
Class to store parameters as a hierarchical key/value tree.
Definition AlinaUtils.h:113
Arcane::Alina::preconditioner::CPRPreconditioner
Two stage preconditioner of the Constrained Pressure Residual type.
Definition CPRPreconditioner.h:46
Arcane::Matrix
Matrix class, to be used by user.
Definition matrix/Matrix.h:36
Arcane::arccoreParallelFor
void arccoreParallelFor(const ComplexForLoopRanges< RankValue, IndexType_ > &loop_ranges, const ForLoopRunInfo &run_info, const LambdaType &lambda_function, const ReducerArgs &... reducer_args)
Applies the lambda function lambda_function concurrently over the iteration interval given by loop_ra...
Definition ParallelFor.h:86
Arcane::Int32
std::int32_t Int32
Signed integer type of 32 bits.
Definition ArccoreGlobal.h:233
Arcane::Alina::backend::backends_compatible
Metafunction that checks if two backends are compatible.
Definition BackendInterface.h:73
Arcane::Alina::math::static_rows
Number of rows for statically sized matrix types.
Definition ValueTypeInterface.h:68