8.2
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 *restrict vx_ini,
121  cs_real_t *restrict vx,
122  size_t aux_size,
123  void *aux_vectors);
124 
125 /* Solver setup data */
126 /*-------------------*/
127 
128 typedef 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 
146 typedef 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 
155 struct _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 
238 struct _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 
269 inline 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 
300 inline 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 
332 inline 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 
370 inline 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 
410 inline 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 
456 inline 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 *xx,
462  double *yy,
463  double *xy,
464  double *xz,
465  double *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 
499 inline static void
500 _fw_and_bw_lu33(const cs_real_t mat[],
501  cs_real_t x[restrict],
502  const cs_real_t b[restrict],
503  const cs_real_t c[restrict])
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 
528 inline static void
529 _fw_and_bw_lu(const cs_real_t mat[],
530  int db_size,
531  cs_real_t x[restrict],
532  const cs_real_t b[restrict],
533  const cs_real_t c[restrict])
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 
567 inline static void
568 _fw_and_bw_lu_gs(const cs_real_t mat[],
569  int db_size,
570  cs_real_t x[restrict],
571  const cs_real_t b[restrict])
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 
610 void
611 cs_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 
634 void
635 cs_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.c:1493
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.c:1682
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.c:1514
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.c:1621
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.c:1650
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.c:1594
#define restrict
Definition: cs_defs.h:141
#define BEGIN_C_DECLS
Definition: cs_defs.h:528
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
#define END_C_DECLS
Definition: cs_defs.h:529
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
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