8.1
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-2023 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 --> system solution
105  * aux_size <-- number of elements in aux_vectors (in bytes)
106  * aux_vectors --- optional working area (allocation otherwise)
107  *
108  * returns:
109  * convergence status
110  *----------------------------------------------------------------------------*/
111 
113 (cs_sles_it_solve_t) (cs_sles_it_t *c,
114  const cs_matrix_t *a,
115  cs_lnum_t diag_block_size,
116  cs_sles_it_convergence_t *convergence,
117  const cs_real_t *rhs,
118  cs_real_t *restrict vx,
119  size_t aux_size,
120  void *aux_vectors);
121 
122 /* Solver setup data */
123 /*-------------------*/
124 
125 typedef struct _cs_sles_it_setup_t {
126 
127  double initial_residual; /* last initial residual value */
128 
129  cs_lnum_t n_rows; /* number of associated rows */
130 
131  const cs_real_t *ad_inv; /* pointer to diagonal inverse */
132  cs_real_t *_ad_inv; /* private pointer to
133  diagonal inverse */
134 
135  void *pc_context; /* preconditioner context */
136  cs_sles_pc_apply_t *pc_apply; /* preconditioner apply */
137 
138 } cs_sles_it_setup_t;
139 
140 /* Solver additional data */
141 /*------------------------*/
142 
143 typedef struct _cs_sles_it_add_t {
144 
145  cs_lnum_t *order; /* ordering */
146 
147 } cs_sles_it_add_t;
148 
149 /* Basic per linear system options and logging */
150 /*---------------------------------------------*/
151 
152 struct _cs_sles_it_t {
153 
154  /* Base settings */
155 
156  cs_sles_it_type_t type; /* Solver type */
157 
158  bool on_device; /* SpMV on device ? */
159 
160  bool update_stats; /* do stats need to be updated ? */
161  bool ignore_convergence; /* ignore convergence for some
162  solvers used as preconditioners */
163 
164  int n_max_iter; /* maximum number of iterations */
165  int restart_interval; /* maximum number of iterations
166  before restarting the algorithm
167  (only applicable for GMRES or GCR
168  algorithm up to now) */
169 
170  cs_sles_it_solve_t *solve; /* pointer to solve function */
171 
172  cs_sles_pc_t *pc; /* pointer to possibly shared
173  preconditioner object */
174  cs_sles_pc_t *_pc; /* pointer to owned
175  preconditioner object */
176 
177  /* Performance data */
178 
179  unsigned n_setups; /* Number of times system setup */
180  unsigned n_solves; /* Number of times system solved */
181 
182  unsigned n_iterations_last; /* Number of iterations for last
183  system resolution */
184  unsigned n_iterations_min; /* Minimum number ot iterations
185  in system resolution history */
186  unsigned n_iterations_max; /* Maximum number ot iterations
187  in system resolution history */
188  unsigned long long n_iterations_tot; /* Total accumulated number of
189  iterations */
190 
191  cs_timer_counter_t t_setup; /* Total setup */
192  cs_timer_counter_t t_solve; /* Total time used */
193 
194  /* Plot info */
195 
196  int plot_time_stamp; /* Plot time stamp */
197  cs_time_plot_t *plot; /* Pointer to plot structure,
198  which may be owned or shared */
199  cs_time_plot_t *_plot; /* Pointer to own plot structure */
200 
201  /* Communicator used for reduction operations
202  (if left at NULL, main communicator will be used) */
203 
204 # if defined(HAVE_MPI)
205  MPI_Comm comm;
206  MPI_Comm caller_comm;
207  int caller_n_ranks;
208 # endif
209 
210  /* Solver setup */
211 
212  const struct _cs_sles_it_t *shared; /* pointer to context sharing some
213  setup and preconditioner data,
214  or NULL */
215 
216  cs_sles_it_add_t *add_data; /* additional data */
217 
218  cs_sles_it_setup_t *setup_data; /* setup data */
219 
220  /* Alternative solvers (fallback or heuristics) */
221 
222  cs_sles_convergence_state_t fallback_cvg; /* threshold for fallback
223  convergence */
224  int fallback_n_max_iter; /* number of maximum iteration
225  for fallback solver */
226 
227  cs_sles_it_t *fallback; /* fallback solver */
228 
229 
230 };
231 
232 /* Convergence testing and tracking */
233 /*----------------------------------*/
234 
235 struct _cs_sles_it_convergence_t {
236 
237  const char *name; /* Pointer to name string */
238 
239  int verbosity; /* Verbosity level */
240 
241  unsigned n_iterations; /* Current number of iterations */
242  unsigned n_iterations_max; /* Maximum number of iterations */
243 
244  double precision; /* Precision limit */
245  double r_norm; /* Residual normalization */
246  double residual; /* Current residual */
247 
248 };
249 
250 /*============================================================================
251  * Inline static function definitions
252  *============================================================================*/
253 
254 /*----------------------------------------------------------------------------
255  * Compute dot product, summing result over all ranks.
256  *
257  * parameters:
258  * c <-- pointer to solver context info
259  * x <-- first vector in s = x.y
260  * y <-- second vector in s = x.y
261  *
262  * returns:
263  * result of s = x.y
264  *----------------------------------------------------------------------------*/
265 
266 inline static double
267 _dot_product(const cs_sles_it_t *c,
268  const cs_real_t *x,
269  const cs_real_t *y)
270 {
271  double s = cs_dot(c->setup_data->n_rows, x, y);
272 
273 #if defined(HAVE_MPI)
274 
275  if (c->comm != MPI_COMM_NULL) {
276  double _sum;
277  MPI_Allreduce(&s, &_sum, 1, MPI_DOUBLE, MPI_SUM, c->comm);
278  s = _sum;
279  }
280 
281 #endif /* defined(HAVE_MPI) */
282 
283  return s;
284 }
285 
286 /*----------------------------------------------------------------------------
287  * Compute dot product x.x, summing result over all ranks.
288  *
289  * parameters:
290  * c <-- pointer to solver context info
291  * x <-- vector in s = x.x
292  *
293  * returns:
294  * result of s = x.x
295  *----------------------------------------------------------------------------*/
296 
297 inline static double
298 _dot_product_xx(const cs_sles_it_t *c,
299  const cs_real_t *x)
300 {
301  double s;
302 
303  s = cs_dot_xx(c->setup_data->n_rows, x);
304 
305 #if defined(HAVE_MPI)
306 
307  if (c->comm != MPI_COMM_NULL) {
308  double _sum;
309  MPI_Allreduce(&s, &_sum, 1, MPI_DOUBLE, MPI_SUM, c->comm);
310  s = _sum;
311  }
312 
313 #endif /* defined(HAVE_MPI) */
314 
315  return s;
316 }
317 
318 /*----------------------------------------------------------------------------
319  * Compute 2 dot products x.x and x.y, summing result over all ranks.
320  *
321  * parameters:
322  * c <-- pointer to solver context info
323  * x <-- vector in s1 = x.x and s2 = x.y
324  * y <-- vector in s2 = x.y
325  * s1 --> result of s1 = x.x
326  * s2 --> result of s2 = x.y
327  *----------------------------------------------------------------------------*/
328 
329 inline static void
330 _dot_products_xx_xy(const cs_sles_it_t *c,
331  const cs_real_t *x,
332  const cs_real_t *y,
333  double *s1,
334  double *s2)
335 {
336  double s[2];
337 
338  cs_dot_xx_xy(c->setup_data->n_rows, x, y, s, s+1);
339 
340 #if defined(HAVE_MPI)
341 
342  if (c->comm != MPI_COMM_NULL) {
343  double _sum[2];
344  MPI_Allreduce(s, _sum, 2, MPI_DOUBLE, MPI_SUM, c->comm);
345  s[0] = _sum[0];
346  s[1] = _sum[1];
347  }
348 
349 #endif /* defined(HAVE_MPI) */
350 
351  *s1 = s[0];
352  *s2 = s[1];
353 }
354 
355 /*----------------------------------------------------------------------------
356  * Compute 2 dot products x.x and x.y, summing result over all ranks.
357  *
358  * parameters:
359  * c <-- pointer to solver context info
360  * x <-- vector in s1 = x.y
361  * y <-- vector in s1 = x.y and s2 = y.z
362  * z <-- vector in s2 = y.z
363  * s1 --> result of s1 = x.y
364  * s2 --> result of s2 = y.z
365  *----------------------------------------------------------------------------*/
366 
367 inline static void
368 _dot_products_xy_yz(const cs_sles_it_t *c,
369  const cs_real_t *x,
370  const cs_real_t *y,
371  const cs_real_t *z,
372  double *s1,
373  double *s2)
374 {
375  double s[2];
376 
377  cs_dot_xy_yz(c->setup_data->n_rows, x, y, z, s, s+1);
378 
379 #if defined(HAVE_MPI)
380 
381  if (c->comm != MPI_COMM_NULL) {
382  double _sum[2];
383  MPI_Allreduce(s, _sum, 2, MPI_DOUBLE, MPI_SUM, c->comm);
384  s[0] = _sum[0];
385  s[1] = _sum[1];
386  }
387 
388 #endif /* defined(HAVE_MPI) */
389 
390  *s1 = s[0];
391  *s2 = s[1];
392 }
393 
394 /*----------------------------------------------------------------------------
395  * Compute 3 dot products, summing result over all ranks.
396  *
397  * parameters:
398  * c <-- pointer to solver context info
399  * x <-- first vector
400  * y <-- second vector
401  * z <-- third vector
402  * s1 --> result of s1 = x.x
403  * s2 --> result of s2 = x.y
404  * s3 --> result of s3 = y.z
405  *----------------------------------------------------------------------------*/
406 
407 inline static void
408 _dot_products_xx_xy_yz(const cs_sles_it_t *c,
409  const cs_real_t *x,
410  const cs_real_t *y,
411  const cs_real_t *z,
412  double *s1,
413  double *s2,
414  double *s3)
415 {
416  double s[3];
417 
418  cs_dot_xx_xy_yz(c->setup_data->n_rows, x, y, z, s, s+1, s+2);
419 
420 #if defined(HAVE_MPI)
421 
422  if (c->comm != MPI_COMM_NULL) {
423  double _sum[3];
424 
425  MPI_Allreduce(s, _sum, 3, MPI_DOUBLE, MPI_SUM, c->comm);
426  s[0] = _sum[0];
427  s[1] = _sum[1];
428  s[2] = _sum[2];
429  }
430 
431 #endif /* defined(HAVE_MPI) */
432 
433  *s1 = s[0];
434  *s2 = s[1];
435  *s3 = s[2];
436 }
437 
438 /*----------------------------------------------------------------------------
439  * Compute 5 dot products, summing result over all ranks.
440  *
441  * parameters:
442  * c <-- pointer to solver context info
443  * x <-- first vector
444  * y <-- second vector
445  * z <-- third vector
446  * xx --> result of x.x
447  * yy --> result of y.y
448  * xy --> result of x.y
449  * xz --> result of x.z
450  * yz --> result of y.z
451  *----------------------------------------------------------------------------*/
452 
453 inline static void
454 _dot_products_xx_yy_xy_xz_yz(const cs_sles_it_t *c,
455  const cs_real_t *x,
456  const cs_real_t *y,
457  const cs_real_t *z,
458  double *xx,
459  double *yy,
460  double *xy,
461  double *xz,
462  double *yz)
463 {
464  double s[5];
465 
466  cs_dot_xx_yy_xy_xz_yz(c->setup_data->n_rows, x, y, z, s, s+1, s+2, s+3, s+4);
467 
468 #if defined(HAVE_MPI)
469 
470  if (c->comm != MPI_COMM_NULL) {
471  double _sum[5];
472  MPI_Allreduce(s, _sum, 5, MPI_DOUBLE, MPI_SUM, c->comm);
473  memcpy(s, _sum, 5*sizeof(double));
474  }
475 
476 #endif /* defined(HAVE_MPI) */
477 
478  *xx = s[0];
479  *yy = s[1];
480  *xy = s[2];
481  *xz = s[3];
482  *yz = s[4];
483 }
484 
485 /*----------------------------------------------------------------------------
486  * Block Jacobi utilities.
487  * Compute forward and backward to solve an LU 3*3 system.
488  *
489  * parameters:
490  * mat <-- 3*3*dim matrix
491  * x --> solution
492  * b --> 1st part of RHS (c - b)
493  * c --> 2nd part of RHS (c - b)
494  *----------------------------------------------------------------------------*/
495 
496 inline static void
497 _fw_and_bw_lu33(const cs_real_t mat[],
498  cs_real_t x[restrict],
499  const cs_real_t b[restrict],
500  const cs_real_t c[restrict])
501 {
502  cs_real_t aux[3];
503 
504  aux[0] = (c[0] - b[0]);
505  aux[1] = (c[1] - b[1]) - aux[0]*mat[3];
506  aux[2] = (c[2] - b[2]) - aux[0]*mat[6] - aux[1]*mat[7];
507 
508  x[2] = aux[2]/mat[8];
509  x[1] = (aux[1] - mat[5]*x[2])/mat[4];
510  x[0] = (aux[0] - mat[1]*x[1] - mat[2]*x[2])/mat[0];
511 }
512 
513 /*----------------------------------------------------------------------------
514  * Block Jacobi utilities.
515  * Compute forward and backward to solve an LU P*P system.
516  *
517  * parameters:
518  * mat <-- P*P*dim matrix
519  * db_size <-- matrix size
520  * x --> solution
521  * b --> 1st part of RHS (c - b)
522  * c --> 2nd part of RHS (c - b)
523  *----------------------------------------------------------------------------*/
524 
525 inline static void
526 _fw_and_bw_lu(const cs_real_t mat[],
527  int db_size,
528  cs_real_t x[restrict],
529  const cs_real_t b[restrict],
530  const cs_real_t c[restrict])
531 {
532  assert(db_size <= DB_SIZE_MAX);
533  cs_real_t aux[DB_SIZE_MAX];
534 
535  /* forward */
536  for (int ii = 0; ii < db_size; ii++) {
537  aux[ii] = (c[ii] - b[ii]);
538  for (int jj = 0; jj < ii; jj++) {
539  aux[ii] -= aux[jj]*mat[ii*db_size + jj];
540  }
541  }
542 
543  /* backward */
544  for (int ii = db_size - 1; ii >= 0; ii-=1) {
545  x[ii] = aux[ii];
546  for (int jj = db_size - 1; jj > ii; jj-=1) {
547  x[ii] -= x[jj]*mat[ii*db_size + jj];
548  }
549  x[ii] /= mat[ii*(db_size + 1)];
550  }
551 }
552 
553 /*----------------------------------------------------------------------------
554  * Block Gauss-Seidel utilities.
555  * Compute forward and backward to solve an LU P*P system.
556  *
557  * parameters:
558  * mat <-- P*P*dim matrix
559  * db_size <-- matrix size
560  * x --> solution
561  * b <-> RHS in, work array
562  *----------------------------------------------------------------------------*/
563 
564 inline static void
565 _fw_and_bw_lu_gs(const cs_real_t mat[],
566  int db_size,
567  cs_real_t x[restrict],
568  const cs_real_t b[restrict])
569 {
570  assert(db_size <= DB_SIZE_MAX);
571 
572  /* forward */
573  for (int ii = 0; ii < db_size; ii++) {
574  x[ii] = b[ii];
575  for (int jj = 0; jj < ii; jj++)
576  x[ii] -= x[jj]*mat[ii*db_size + jj];
577  }
578 
579  /* backward */
580  for (int ii = db_size - 1; ii >= 0; ii--) {
581  for (int jj = db_size - 1; jj > ii; jj--)
582  x[ii] -= x[jj]*mat[ii*db_size + jj];
583  x[ii] /= mat[ii*(db_size + 1)];
584  }
585 }
586 
587 /*============================================================================
588  * Public function definitions
589  *============================================================================*/
590 
591 /*----------------------------------------------------------------------------*/
605 /*----------------------------------------------------------------------------*/
606 
607 void
608 cs_sles_it_convergence_init(cs_sles_it_convergence_t *convergence,
609  const char *solver_name,
610  int verbosity,
611  unsigned n_iter_max,
612  double precision,
613  double r_norm,
614  double *residual);
615 
616 /*----------------------------------------------------------------------------
617  * Setup context for iterative linear solver.
618  *
619  * This function is common to most solvers
620  *
621  * parameters:
622  * c <-> pointer to solver context info
623  * name <-- pointer to system name
624  * a <-- matrix
625  * verbosity <-- verbosity level
626  * diag_block_size <-- diagonal block size
627  * block_nn_inverse <-- if diagonal block size is 3 or 6, compute inverse of
628  * block if true, inverse of block diagonal otherwise
629  *----------------------------------------------------------------------------*/
630 
631 void
632 cs_sles_it_setup_priv(cs_sles_it_t *c,
633  const char *name,
634  const cs_matrix_t *a,
635  int verbosity,
636  int diag_block_size,
637  bool block_nn_inverse);
638 
641 /*----------------------------------------------------------------------------*/
642 
644 
645 #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:1465
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:1654
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:1486
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:1593
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:1622
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:1566
#define restrict
Definition: cs_defs.h:139
#define BEGIN_C_DECLS
Definition: cs_defs.h:514
double cs_real_t
Floating-point value.
Definition: cs_defs.h:319
#define END_C_DECLS
Definition: cs_defs.h:515
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:313
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