8.3
general documentation
cs_sles_it_priv.h
Go to the documentation of this file.
1#ifndef __CS_SLES_IT_PRIV_H__
2#define __CS_SLES_IT_PRIV_H__
3
4/*============================================================================
5 * Sparse Linear Equation Solvers: private elements.
6 *
7 * These elements are shared between iterative solvers and smoother
8 * both for host and device implementations, but are not accessible to
9 * calling code.
10 *============================================================================*/
11
12/*
13 This file is part of code_saturne, a general-purpose CFD tool.
14
15 Copyright (C) 1998-2024 EDF S.A.
16
17 This program is free software; you can redistribute it and/or modify it under
18 the terms of the GNU General Public License as published by the Free Software
19 Foundation; either version 2 of the License, or (at your option) any later
20 version.
21
22 This program is distributed in the hope that it will be useful, but WITHOUT
23 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
24 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
25 details.
26
27 You should have received a copy of the GNU General Public License along with
28 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
29 Street, Fifth Floor, Boston, MA 02110-1301, USA.
30*/
31
32/*----------------------------------------------------------------------------*/
33
34#include "cs_defs.h"
35
36/*----------------------------------------------------------------------------
37 * Standard C library headers
38 *----------------------------------------------------------------------------*/
39
40#include <assert.h>
41#include <math.h>
42
43#if defined(HAVE_MPI)
44#include <mpi.h>
45#endif
46
47/*----------------------------------------------------------------------------
48 * Local headers
49 *----------------------------------------------------------------------------*/
50
51#include "bft_mem.h"
52#include "bft_error.h"
53#include "bft_printf.h"
54
55#include "cs_base.h"
56#include "cs_blas.h"
57#include "cs_file.h"
58#include "cs_log.h"
59#include "cs_halo.h"
60#include "cs_mesh.h"
61#include "cs_matrix.h"
62#include "cs_matrix_default.h"
63#include "cs_matrix_util.h"
64#include "cs_post.h"
65#include "cs_timer.h"
66#include "cs_time_plot.h"
67
68/*----------------------------------------------------------------------------
69 * Header for the current file
70 *----------------------------------------------------------------------------*/
71
72#include "cs_sles.h"
73#include "cs_sles_it.h"
74#include "cs_sles_pc.h"
75
76/*----------------------------------------------------------------------------*/
77
79
82/*=============================================================================
83 * Local Macro Definitions
84 *============================================================================*/
85
86#if !defined(HUGE_VAL)
87#define HUGE_VAL 1.E+12
88#endif
89
90#define DB_SIZE_MAX 9
91
92/*=============================================================================
93 * Local Structure Definitions
94 *============================================================================*/
95
96/*----------------------------------------------------------------------------
97 * Function pointer for actual resolution of a linear system.
98 *
99 * parameters:
100 * c <-- pointer to solver context info
101 * a <-- linear equation matrix
102 * convergence <-- convergence information structure
103 * rhs <-- right hand side
104 * vx_ini <-- initial system solution
105 * (vx if nonzero, nullptr if zero)
106 * vx <-> system solution
107 * aux_size <-- number of elements in aux_vectors (in bytes)
108 * aux_vectors --- optional working area (allocation otherwise)
109 *
110 * returns:
111 * convergence status
112 *----------------------------------------------------------------------------*/
113
115(cs_sles_it_solve_t) (cs_sles_it_t *c,
116 const cs_matrix_t *a,
117 cs_lnum_t diag_block_size,
118 cs_sles_it_convergence_t *convergence,
119 const cs_real_t *rhs,
120 cs_real_t *vx_ini,
121 cs_real_t *vx,
122 size_t aux_size,
123 void *aux_vectors);
124
125/* Solver setup data */
126/*-------------------*/
127
128typedef struct _cs_sles_it_setup_t {
129
130 double initial_residual; /* last initial residual value */
131
132 cs_lnum_t n_rows; /* number of associated rows */
133
134 const cs_real_t *ad_inv; /* pointer to diagonal inverse */
135 cs_real_t *_ad_inv; /* private pointer to
136 diagonal inverse */
137
138 void *pc_context; /* preconditioner context */
139 cs_sles_pc_apply_t *pc_apply; /* preconditioner apply */
140
141} cs_sles_it_setup_t;
142
143/* Solver additional data */
144/*------------------------*/
145
146typedef struct _cs_sles_it_add_t {
147
148 cs_lnum_t *order; /* ordering */
149
150} cs_sles_it_add_t;
151
152/* Basic per linear system options and logging */
153/*---------------------------------------------*/
154
155struct _cs_sles_it_t {
156
157 /* Base settings */
158
159 cs_sles_it_type_t type; /* Solver type */
160
161 bool on_device; /* SpMV on device ? */
162
163 bool update_stats; /* do stats need to be updated ? */
164 bool ignore_convergence; /* ignore convergence for some
165 solvers used as preconditioners */
166
167 int n_max_iter; /* maximum number of iterations */
168 int restart_interval; /* maximum number of iterations
169 before restarting the algorithm
170 (only applicable for GMRES or GCR
171 algorithm up to now) */
172
173 cs_sles_it_solve_t *solve; /* pointer to solve function */
174
175 cs_sles_pc_t *pc; /* pointer to possibly shared
176 preconditioner object */
177 cs_sles_pc_t *_pc; /* pointer to owned
178 preconditioner object */
179
180 /* Performance data */
181
182 unsigned n_setups; /* Number of times system setup */
183 unsigned n_solves; /* Number of times system solved */
184
185 unsigned n_iterations_last; /* Number of iterations for last
186 system resolution */
187 unsigned n_iterations_min; /* Minimum number ot iterations
188 in system resolution history */
189 unsigned n_iterations_max; /* Maximum number ot iterations
190 in system resolution history */
191 unsigned long long n_iterations_tot; /* Total accumulated number of
192 iterations */
193
194 cs_timer_counter_t t_setup; /* Total setup */
195 cs_timer_counter_t t_solve; /* Total time used */
196
197 /* Plot info */
198
199 int plot_time_stamp; /* Plot time stamp */
200 cs_time_plot_t *plot; /* Pointer to plot structure,
201 which may be owned or shared */
202 cs_time_plot_t *_plot; /* Pointer to own plot structure */
203
204 /* Communicator used for reduction operations
205 (if left at NULL, main communicator will be used) */
206
207# if defined(HAVE_MPI)
208 MPI_Comm comm;
209 MPI_Comm caller_comm;
210 int caller_n_ranks;
211# endif
212
213 /* Solver setup */
214
215 const struct _cs_sles_it_t *shared; /* pointer to context sharing some
216 setup and preconditioner data,
217 or NULL */
218
219 cs_sles_it_add_t *add_data; /* additional data */
220
221 cs_sles_it_setup_t *setup_data; /* setup data */
222
223 /* Alternative solvers (fallback or heuristics) */
224
225 cs_sles_convergence_state_t fallback_cvg; /* threshold for fallback
226 convergence */
227 int fallback_n_max_iter; /* number of maximum iteration
228 for fallback solver */
229
230 cs_sles_it_t *fallback; /* fallback solver */
231
232
233};
234
235/* Convergence testing and tracking */
236/*----------------------------------*/
237
238struct _cs_sles_it_convergence_t {
239
240 const char *name; /* Pointer to name string */
241
242 int verbosity; /* Verbosity level */
243
244 unsigned n_iterations; /* Current number of iterations */
245 unsigned n_iterations_max; /* Maximum number of iterations */
246
247 double precision; /* Precision limit */
248 double r_norm; /* Residual normalization */
249 double residual; /* Current residual */
250
251};
252
253/*============================================================================
254 * Inline static function definitions
255 *============================================================================*/
256
257/*----------------------------------------------------------------------------
258 * Compute dot product, summing result over all ranks.
259 *
260 * parameters:
261 * c <-- pointer to solver context info
262 * x <-- first vector in s = x.y
263 * y <-- second vector in s = x.y
264 *
265 * returns:
266 * result of s = x.y
267 *----------------------------------------------------------------------------*/
268
269inline static double
270_dot_product(const cs_sles_it_t *c,
271 const cs_real_t *x,
272 const cs_real_t *y)
273{
274 double s = cs_dot(c->setup_data->n_rows, x, y);
275
276#if defined(HAVE_MPI)
277
278 if (c->comm != MPI_COMM_NULL) {
279 double _sum;
280 MPI_Allreduce(&s, &_sum, 1, MPI_DOUBLE, MPI_SUM, c->comm);
281 s = _sum;
282 }
283
284#endif /* defined(HAVE_MPI) */
285
286 return s;
287}
288
289/*----------------------------------------------------------------------------
290 * Compute dot product x.x, summing result over all ranks.
291 *
292 * parameters:
293 * c <-- pointer to solver context info
294 * x <-- vector in s = x.x
295 *
296 * returns:
297 * result of s = x.x
298 *----------------------------------------------------------------------------*/
299
300inline static double
301_dot_product_xx(const cs_sles_it_t *c,
302 const cs_real_t *x)
303{
304 double s;
305
306 s = cs_dot_xx(c->setup_data->n_rows, x);
307
308#if defined(HAVE_MPI)
309
310 if (c->comm != MPI_COMM_NULL) {
311 double _sum;
312 MPI_Allreduce(&s, &_sum, 1, MPI_DOUBLE, MPI_SUM, c->comm);
313 s = _sum;
314 }
315
316#endif /* defined(HAVE_MPI) */
317
318 return s;
319}
320
321/*----------------------------------------------------------------------------
322 * Compute 2 dot products x.x and x.y, summing result over all ranks.
323 *
324 * parameters:
325 * c <-- pointer to solver context info
326 * x <-- vector in s1 = x.x and s2 = x.y
327 * y <-- vector in s2 = x.y
328 * s1 --> result of s1 = x.x
329 * s2 --> result of s2 = x.y
330 *----------------------------------------------------------------------------*/
331
332inline static void
333_dot_products_xx_xy(const cs_sles_it_t *c,
334 const cs_real_t *x,
335 const cs_real_t *y,
336 double *s1,
337 double *s2)
338{
339 double s[2];
340
341 cs_dot_xx_xy(c->setup_data->n_rows, x, y, s, s+1);
342
343#if defined(HAVE_MPI)
344
345 if (c->comm != MPI_COMM_NULL) {
346 double _sum[2];
347 MPI_Allreduce(s, _sum, 2, MPI_DOUBLE, MPI_SUM, c->comm);
348 s[0] = _sum[0];
349 s[1] = _sum[1];
350 }
351
352#endif /* defined(HAVE_MPI) */
353
354 *s1 = s[0];
355 *s2 = s[1];
356}
357
358/*----------------------------------------------------------------------------
359 * Compute 2 dot products x.x and x.y, summing result over all ranks.
360 *
361 * parameters:
362 * c <-- pointer to solver context info
363 * x <-- vector in s1 = x.y
364 * y <-- vector in s1 = x.y and s2 = y.z
365 * z <-- vector in s2 = y.z
366 * s1 --> result of s1 = x.y
367 * s2 --> result of s2 = y.z
368 *----------------------------------------------------------------------------*/
369
370inline static void
371_dot_products_xy_yz(const cs_sles_it_t *c,
372 const cs_real_t *x,
373 const cs_real_t *y,
374 const cs_real_t *z,
375 double *s1,
376 double *s2)
377{
378 double s[2];
379
380 cs_dot_xy_yz(c->setup_data->n_rows, x, y, z, s, s+1);
381
382#if defined(HAVE_MPI)
383
384 if (c->comm != MPI_COMM_NULL) {
385 double _sum[2];
386 MPI_Allreduce(s, _sum, 2, MPI_DOUBLE, MPI_SUM, c->comm);
387 s[0] = _sum[0];
388 s[1] = _sum[1];
389 }
390
391#endif /* defined(HAVE_MPI) */
392
393 *s1 = s[0];
394 *s2 = s[1];
395}
396
397/*----------------------------------------------------------------------------
398 * Compute 3 dot products, summing result over all ranks.
399 *
400 * parameters:
401 * c <-- pointer to solver context info
402 * x <-- first vector
403 * y <-- second vector
404 * z <-- third vector
405 * s1 --> result of s1 = x.x
406 * s2 --> result of s2 = x.y
407 * s3 --> result of s3 = y.z
408 *----------------------------------------------------------------------------*/
409
410inline static void
411_dot_products_xx_xy_yz(const cs_sles_it_t *c,
412 const cs_real_t *x,
413 const cs_real_t *y,
414 const cs_real_t *z,
415 double *s1,
416 double *s2,
417 double *s3)
418{
419 double s[3];
420
421 cs_dot_xx_xy_yz(c->setup_data->n_rows, x, y, z, s, s+1, s+2);
422
423#if defined(HAVE_MPI)
424
425 if (c->comm != MPI_COMM_NULL) {
426 double _sum[3];
427
428 MPI_Allreduce(s, _sum, 3, MPI_DOUBLE, MPI_SUM, c->comm);
429 s[0] = _sum[0];
430 s[1] = _sum[1];
431 s[2] = _sum[2];
432 }
433
434#endif /* defined(HAVE_MPI) */
435
436 *s1 = s[0];
437 *s2 = s[1];
438 *s3 = s[2];
439}
440
441/*----------------------------------------------------------------------------
442 * Compute 5 dot products, summing result over all ranks.
443 *
444 * parameters:
445 * c <-- pointer to solver context info
446 * x <-- first vector
447 * y <-- second vector
448 * z <-- third vector
449 * xx --> result of x.x
450 * yy --> result of y.y
451 * xy --> result of x.y
452 * xz --> result of x.z
453 * yz --> result of y.z
454 *----------------------------------------------------------------------------*/
455
456inline static void
457_dot_products_xx_yy_xy_xz_yz(const cs_sles_it_t *c,
458 const cs_real_t *x,
459 const cs_real_t *y,
460 const cs_real_t *z,
461 double *restrict xx,
462 double *restrict yy,
463 double *restrict xy,
464 double *restrict xz,
465 double *restrict yz)
466{
467 double s[5];
468
469 cs_dot_xx_yy_xy_xz_yz(c->setup_data->n_rows, x, y, z, s, s+1, s+2, s+3, s+4);
470
471#if defined(HAVE_MPI)
472
473 if (c->comm != MPI_COMM_NULL) {
474 double _sum[5];
475 MPI_Allreduce(s, _sum, 5, MPI_DOUBLE, MPI_SUM, c->comm);
476 memcpy(s, _sum, 5*sizeof(double));
477 }
478
479#endif /* defined(HAVE_MPI) */
480
481 *xx = s[0];
482 *yy = s[1];
483 *xy = s[2];
484 *xz = s[3];
485 *yz = s[4];
486}
487
488/*----------------------------------------------------------------------------
489 * Block Jacobi utilities.
490 * Compute forward and backward to solve an LU 3*3 system.
491 *
492 * parameters:
493 * mat <-- 3*3*dim matrix
494 * x --> solution
495 * b --> 1st part of RHS (c - b)
496 * c --> 2nd part of RHS (c - b)
497 *----------------------------------------------------------------------------*/
498
499inline static void
500_fw_and_bw_lu33(const cs_real_t mat[],
502 const cs_real_t *restrict b,
503 const cs_real_t *restrict c)
504{
505 cs_real_t aux[3];
506
507 aux[0] = (c[0] - b[0]);
508 aux[1] = (c[1] - b[1]) - aux[0]*mat[3];
509 aux[2] = (c[2] - b[2]) - aux[0]*mat[6] - aux[1]*mat[7];
510
511 x[2] = aux[2]/mat[8];
512 x[1] = (aux[1] - mat[5]*x[2])/mat[4];
513 x[0] = (aux[0] - mat[1]*x[1] - mat[2]*x[2])/mat[0];
514}
515
516/*----------------------------------------------------------------------------
517 * Block Jacobi utilities.
518 * Compute forward and backward to solve an LU P*P system.
519 *
520 * parameters:
521 * mat <-- P*P*dim matrix
522 * db_size <-- matrix size
523 * x --> solution
524 * b --> 1st part of RHS (c - b)
525 * c --> 2nd part of RHS (c - b)
526 *----------------------------------------------------------------------------*/
527
528inline static void
529_fw_and_bw_lu(const cs_real_t mat[],
530 int db_size,
532 const cs_real_t *restrict b,
533 const cs_real_t *restrict c)
534{
535 assert(db_size <= DB_SIZE_MAX);
536 cs_real_t aux[DB_SIZE_MAX];
537
538 /* forward */
539 for (int ii = 0; ii < db_size; ii++) {
540 aux[ii] = (c[ii] - b[ii]);
541 for (int jj = 0; jj < ii; jj++) {
542 aux[ii] -= aux[jj]*mat[ii*db_size + jj];
543 }
544 }
545
546 /* backward */
547 for (int ii = db_size - 1; ii >= 0; ii-=1) {
548 x[ii] = aux[ii];
549 for (int jj = db_size - 1; jj > ii; jj-=1) {
550 x[ii] -= x[jj]*mat[ii*db_size + jj];
551 }
552 x[ii] /= mat[ii*(db_size + 1)];
553 }
554}
555
556/*----------------------------------------------------------------------------
557 * Block Gauss-Seidel utilities.
558 * Compute forward and backward to solve an LU P*P system.
559 *
560 * parameters:
561 * mat <-- P*P*dim matrix
562 * db_size <-- matrix size
563 * x --> solution
564 * b <-> RHS in, work array
565 *----------------------------------------------------------------------------*/
566
567inline static void
568_fw_and_bw_lu_gs(const cs_real_t mat[],
569 int db_size,
571 const cs_real_t *restrict b)
572{
573 assert(db_size <= DB_SIZE_MAX);
574
575 /* forward */
576 for (int ii = 0; ii < db_size; ii++) {
577 x[ii] = b[ii];
578 for (int jj = 0; jj < ii; jj++)
579 x[ii] -= x[jj]*mat[ii*db_size + jj];
580 }
581
582 /* backward */
583 for (int ii = db_size - 1; ii >= 0; ii--) {
584 for (int jj = db_size - 1; jj > ii; jj--)
585 x[ii] -= x[jj]*mat[ii*db_size + jj];
586 x[ii] /= mat[ii*(db_size + 1)];
587 }
588}
589
590/*============================================================================
591 * Public function definitions
592 *============================================================================*/
593
594/*----------------------------------------------------------------------------*/
608/*----------------------------------------------------------------------------*/
609
610void
611cs_sles_it_convergence_init(cs_sles_it_convergence_t *convergence,
612 const char *solver_name,
613 int verbosity,
614 unsigned n_iter_max,
615 double precision,
616 double r_norm,
617 double *residual);
618
619/*----------------------------------------------------------------------------
620 * Setup context for iterative linear solver.
621 *
622 * This function is common to most solvers
623 *
624 * parameters:
625 * c <-> pointer to solver context info
626 * name <-- pointer to system name
627 * a <-- matrix
628 * verbosity <-- verbosity level
629 * diag_block_size <-- diagonal block size
630 * block_nn_inverse <-- if diagonal block size is 3 or 6, compute inverse of
631 * block if true, inverse of block diagonal otherwise
632 *----------------------------------------------------------------------------*/
633
634void
635cs_sles_it_setup_priv(cs_sles_it_t *c,
636 const char *name,
637 const cs_matrix_t *a,
638 int verbosity,
639 int diag_block_size,
640 bool block_nn_inverse);
641
644/*----------------------------------------------------------------------------*/
645
647
648#endif /* __CS_SLES_IT_PRIV_H__ */
double cs_dot(cs_lnum_t n, const cs_real_t *x, const cs_real_t *y)
Return the dot product of 2 vectors: x.y.
Definition: cs_blas.cpp:1495
void cs_dot_xx_yy_xy_xz_yz(cs_lnum_t n, const cs_real_t *restrict x, const cs_real_t *restrict y, const cs_real_t *restrict z, double *xx, double *yy, double *xy, double *xz, double *yz)
Return 5 dot products of 3 vectors: x.x, y.y, x.y, x.z, and y.z.
Definition: cs_blas.cpp:1684
double cs_dot_xx(cs_lnum_t n, const cs_real_t *x)
Return dot products of a vector with itself: x.x.
Definition: cs_blas.cpp:1516
void cs_dot_xy_yz(cs_lnum_t n, const cs_real_t *restrict x, const cs_real_t *restrict y, const cs_real_t *restrict z, double *xy, double *yz)
Return 2 dot products of 3 vectors: x.y, and y.z.
Definition: cs_blas.cpp:1623
void cs_dot_xx_xy_yz(cs_lnum_t n, const cs_real_t *restrict x, const cs_real_t *restrict y, const cs_real_t *restrict z, double *xx, double *xy, double *yz)
Return 3 dot products of 3 vectors: x.x, x.y, and y.z.
Definition: cs_blas.cpp:1652
void cs_dot_xx_xy(cs_lnum_t n, const cs_real_t *restrict x, const cs_real_t *restrict y, double *xx, double *xy)
Return 2 dot products of 2 vectors: x.x, and x.y.
Definition: cs_blas.cpp:1596
#define restrict
Definition: cs_defs.h:145
#define BEGIN_C_DECLS
Definition: cs_defs.h:542
double cs_real_t
Floating-point value.
Definition: cs_defs.h:342
#define END_C_DECLS
Definition: cs_defs.h:543
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:335
struct _cs_matrix_t cs_matrix_t
Definition: cs_matrix.h:110
cs_sles_convergence_state_t
Definition: cs_sles.h:56
struct _cs_sles_it_t cs_sles_it_t
Definition: cs_sles_it.h:86
struct _cs_sles_it_convergence_t cs_sles_it_convergence_t
Definition: cs_sles_it.h:90
cs_sles_it_type_t
Definition: cs_sles_it.h:55
cs_sles_pc_state_t() cs_sles_pc_apply_t(void *context, const cs_real_t *x_in, cs_real_t *x_out)
Function pointer for application of a preconditioner.
Definition: cs_sles_pc.h:145
struct _cs_sles_pc_t cs_sles_pc_t
Definition: cs_sles_pc.h:66
struct _cs_time_plot_t cs_time_plot_t
Definition: cs_time_plot.h:48
Definition: cs_timer.h:55