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