7.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: Iterative solvers
6  *============================================================================*/
7 
8 /*
9  This file is part of Code_Saturne, a general-purpose CFD tool.
10 
11  Copyright (C) 1998-2021 EDF S.A.
12 
13  This program is free software; you can redistribute it and/or modify it under
14  the terms of the GNU General Public License as published by the Free Software
15  Foundation; either version 2 of the License, or (at your option) any later
16  version.
17 
18  This program is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21  details.
22 
23  You should have received a copy of the GNU General Public License along with
24  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25  Street, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27 
28 /*----------------------------------------------------------------------------*/
29 
30 #include "cs_defs.h"
31 
32 /*----------------------------------------------------------------------------
33  * Standard C library headers
34  *----------------------------------------------------------------------------*/
35 
36 #include <stdarg.h>
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <string.h>
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 8
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_residue; /* last initial residue 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 update_stats; /* do stats need to be updated ? */
159  bool ignore_convergence; /* ignore convergence for some
160  solvers used as preconditioners */
161 
162  int n_max_iter; /* maximum number of iterations */
163  int restart_interval; /* maximum number of iterations
164  before restarting the algorithm
165  (only applicable for GMRES or GCR
166  algorithm up to now) */
167 
168  cs_sles_it_solve_t *solve; /* pointer to solve function */
169 
170  cs_sles_pc_t *pc; /* pointer to possibly shared
171  preconditioner object */
172  cs_sles_pc_t *_pc; /* pointer to owned
173  preconditioner object */
174 
175  /* Performance data */
176 
177  unsigned n_setups; /* Number of times system setup */
178  unsigned n_solves; /* Number of times system solved */
179 
180  unsigned n_iterations_last; /* Number of iterations for last
181  system resolution */
182  unsigned n_iterations_min; /* Minimum number ot iterations
183  in system resolution history */
184  unsigned n_iterations_max; /* Maximum number ot iterations
185  in system resolution history */
186  unsigned long long n_iterations_tot; /* Total accumulated number of
187  iterations */
188 
189  cs_timer_counter_t t_setup; /* Total setup */
190  cs_timer_counter_t t_solve; /* Total time used */
191 
192  /* Plot info */
193 
194  int plot_time_stamp; /* Plot time stamp */
195  cs_time_plot_t *plot; /* Pointer to plot structure,
196  which may be owned or shared */
197  cs_time_plot_t *_plot; /* Pointer to own plot structure */
198 
199  /* Communicator used for reduction operations
200  (if left at NULL, main communicator will be used) */
201 
202 # if defined(HAVE_MPI)
203  MPI_Comm comm;
204  MPI_Comm caller_comm;
205  int caller_n_ranks;
206 # endif
207 
208  /* Solver setup */
209 
210  const struct _cs_sles_it_t *shared; /* pointer to context sharing some
211  setup and preconditioner data,
212  or NULL */
213 
214  cs_sles_it_add_t *add_data; /* additional data */
215 
216  cs_sles_it_setup_t *setup_data; /* setup data */
217 
218  /* Alternative solvers (fallback or heuristics) */
219 
220  cs_sles_convergence_state_t fallback_cvg; /* threshold for fallback
221  convergence */
222  cs_sles_it_t *fallback; /* fallback solver */
223 
224 };
225 
226 /* Convergence testing and tracking */
227 /*----------------------------------*/
228 
229 struct _cs_sles_it_convergence_t {
230 
231  const char *name; /* Pointer to name string */
232 
233  int verbosity; /* Verbosity level */
234 
235  unsigned n_iterations; /* Current number of iterations */
236  unsigned n_iterations_max; /* Maximum number of iterations */
237 
238  double precision; /* Precision limit */
239  double r_norm; /* Residue normalization */
240  double residue; /* Current residue */
241 
242 };
243 
244 /*============================================================================
245  * Inline static function definitions
246  *============================================================================*/
247 
248 /*----------------------------------------------------------------------------
249  * Compute dot product, summing result over all ranks.
250  *
251  * parameters:
252  * c <-- pointer to solver context info
253  * x <-- first vector in s = x.y
254  * y <-- second vector in s = x.y
255  *
256  * returns:
257  * result of s = x.y
258  *----------------------------------------------------------------------------*/
259 
260 inline static double
261 _dot_product(const cs_sles_it_t *c,
262  const cs_real_t *x,
263  const cs_real_t *y)
264 {
265  double s = cs_dot(c->setup_data->n_rows, x, y);
266 
267 #if defined(HAVE_MPI)
268 
269  if (c->comm != MPI_COMM_NULL) {
270  double _sum;
271  MPI_Allreduce(&s, &_sum, 1, MPI_DOUBLE, MPI_SUM, c->comm);
272  s = _sum;
273  }
274 
275 #endif /* defined(HAVE_MPI) */
276 
277  return s;
278 }
279 
280 /*----------------------------------------------------------------------------
281  * Compute dot product x.x, summing result over all ranks.
282  *
283  * parameters:
284  * c <-- pointer to solver context info
285  * x <-- vector in s = x.x
286  *
287  * returns:
288  * result of s = x.x
289  *----------------------------------------------------------------------------*/
290 
291 inline static double
292 _dot_product_xx(const cs_sles_it_t *c,
293  const cs_real_t *x)
294 {
295  double s;
296 
297  s = cs_dot_xx(c->setup_data->n_rows, x);
298 
299 #if defined(HAVE_MPI)
300 
301  if (c->comm != MPI_COMM_NULL) {
302  double _sum;
303  MPI_Allreduce(&s, &_sum, 1, MPI_DOUBLE, MPI_SUM, c->comm);
304  s = _sum;
305  }
306 
307 #endif /* defined(HAVE_MPI) */
308 
309  return s;
310 }
311 
312 /*----------------------------------------------------------------------------
313  * Compute 2 dot products x.x and x.y, summing result over all ranks.
314  *
315  * parameters:
316  * c <-- pointer to solver context info
317  * x <-- vector in s1 = x.x and s2 = x.y
318  * y <-- vector in s2 = x.y
319  * s1 --> result of s1 = x.x
320  * s2 --> result of s2 = x.y
321  *----------------------------------------------------------------------------*/
322 
323 inline static void
324 _dot_products_xx_xy(const cs_sles_it_t *c,
325  const cs_real_t *x,
326  const cs_real_t *y,
327  double *s1,
328  double *s2)
329 {
330  double s[2];
331 
332  cs_dot_xx_xy(c->setup_data->n_rows, x, y, s, s+1);
333 
334 #if defined(HAVE_MPI)
335 
336  if (c->comm != MPI_COMM_NULL) {
337  double _sum[2];
338  MPI_Allreduce(s, _sum, 2, MPI_DOUBLE, MPI_SUM, c->comm);
339  s[0] = _sum[0];
340  s[1] = _sum[1];
341  }
342 
343 #endif /* defined(HAVE_MPI) */
344 
345  *s1 = s[0];
346  *s2 = s[1];
347 }
348 
349 /*----------------------------------------------------------------------------
350  * Compute 2 dot products x.x and x.y, summing result over all ranks.
351  *
352  * parameters:
353  * c <-- pointer to solver context info
354  * x <-- vector in s1 = x.y
355  * y <-- vector in s1 = x.y and s2 = y.z
356  * z <-- vector in s2 = y.z
357  * s1 --> result of s1 = x.y
358  * s2 --> result of s2 = y.z
359  *----------------------------------------------------------------------------*/
360 
361 inline static void
362 _dot_products_xy_yz(const cs_sles_it_t *c,
363  const cs_real_t *x,
364  const cs_real_t *y,
365  const cs_real_t *z,
366  double *s1,
367  double *s2)
368 {
369  double s[2];
370 
371  cs_dot_xy_yz(c->setup_data->n_rows, x, y, z, s, s+1);
372 
373 #if defined(HAVE_MPI)
374 
375  if (c->comm != MPI_COMM_NULL) {
376  double _sum[2];
377  MPI_Allreduce(s, _sum, 2, MPI_DOUBLE, MPI_SUM, c->comm);
378  s[0] = _sum[0];
379  s[1] = _sum[1];
380  }
381 
382 #endif /* defined(HAVE_MPI) */
383 
384  *s1 = s[0];
385  *s2 = s[1];
386 }
387 
388 /*----------------------------------------------------------------------------
389  * Compute 3 dot products, summing result over all ranks.
390  *
391  * parameters:
392  * c <-- pointer to solver context info
393  * x <-- first vector
394  * y <-- second vector
395  * z <-- third vector
396  * s1 --> result of s1 = x.x
397  * s2 --> result of s2 = x.y
398  * s3 --> result of s3 = y.z
399  *----------------------------------------------------------------------------*/
400 
401 inline static void
402 _dot_products_xx_xy_yz(const cs_sles_it_t *c,
403  const cs_real_t *x,
404  const cs_real_t *y,
405  const cs_real_t *z,
406  double *s1,
407  double *s2,
408  double *s3)
409 {
410  double s[3];
411 
412  cs_dot_xx_xy_yz(c->setup_data->n_rows, x, y, z, s, s+1, s+2);
413 
414 #if defined(HAVE_MPI)
415 
416  if (c->comm != MPI_COMM_NULL) {
417  double _sum[3];
418 
419  MPI_Allreduce(s, _sum, 3, MPI_DOUBLE, MPI_SUM, c->comm);
420  s[0] = _sum[0];
421  s[1] = _sum[1];
422  s[2] = _sum[2];
423  }
424 
425 #endif /* defined(HAVE_MPI) */
426 
427  *s1 = s[0];
428  *s2 = s[1];
429  *s3 = s[2];
430 }
431 
432 /*----------------------------------------------------------------------------
433  * Compute 5 dot products, summing result over all ranks.
434  *
435  * parameters:
436  * c <-- pointer to solver context info
437  * x <-- first vector
438  * y <-- second vector
439  * z <-- third vector
440  * xx --> result of x.x
441  * yy --> result of y.y
442  * xy --> result of x.y
443  * xz --> result of x.z
444  * yz --> result of y.z
445  *----------------------------------------------------------------------------*/
446 
447 inline static void
448 _dot_products_xx_yy_xy_xz_yz(const cs_sles_it_t *c,
449  const cs_real_t *x,
450  const cs_real_t *y,
451  const cs_real_t *z,
452  double *xx,
453  double *yy,
454  double *xy,
455  double *xz,
456  double *yz)
457 {
458  double s[5];
459 
460  cs_dot_xx_yy_xy_xz_yz(c->setup_data->n_rows, x, y, z, s, s+1, s+2, s+3, s+4);
461 
462 #if defined(HAVE_MPI)
463 
464  if (c->comm != MPI_COMM_NULL) {
465  double _sum[5];
466  MPI_Allreduce(s, _sum, 5, MPI_DOUBLE, MPI_SUM, c->comm);
467  memcpy(s, _sum, 5*sizeof(double));
468  }
469 
470 #endif /* defined(HAVE_MPI) */
471 
472  *xx = s[0];
473  *yy = s[1];
474  *xy = s[2];
475  *xz = s[3];
476  *yz = s[4];
477 }
478 
479 /*----------------------------------------------------------------------------
480  * Compute 4 dot products, summing result over all ranks.
481  *
482  * parameters:
483  * c <-- pointer to solver context info
484  * v <-- first vector
485  * r <-- second vector
486  * w <-- third vector
487  * q <-- fourth vector
488  * s1 --> result of s1 = v.r
489  * s2 --> result of s2 = v.w
490  * s3 --> result of s3 = v.q
491  * s4 --> result of s4 = r.r
492  *----------------------------------------------------------------------------*/
493 
494 inline static void
495 _dot_products_vr_vw_vq_rr(const cs_sles_it_t *c,
496  const cs_real_t *v,
497  const cs_real_t *r,
498  const cs_real_t *w,
499  const cs_real_t *q,
500  double *s1,
501  double *s2,
502  double *s3,
503  double *s4)
504 {
505  double s[4];
506 
507  /* Use two separate call as cs_blas.c does not yet hav matching call */
508 
509  cs_dot_xy_yz(c->setup_data->n_rows, w, v, q, s+1, s+2);
510  cs_dot_xx_xy(c->setup_data->n_rows, r, v, s+3, s);
511 
512 #if defined(HAVE_MPI)
513 
514  if (c->comm != MPI_COMM_NULL) {
515  double _sum[4];
516  MPI_Allreduce(s, _sum, 4, MPI_DOUBLE, MPI_SUM, c->comm);
517  memcpy(s, _sum, 4*sizeof(double));
518  }
519 
520 #endif /* defined(HAVE_MPI) */
521 
522  *s1 = s[0];
523  *s2 = s[1];
524  *s3 = s[2];
525  *s4 = s[3];
526 }
527 
528 /*----------------------------------------------------------------------------
529  * Block Jacobi utilities.
530  * Compute forward and backward to solve an LU 3*3 system.
531  *
532  * parameters:
533  * mat <-- 3*3*dim matrix
534  * x --> solution
535  * b --> 1st part of RHS (c - b)
536  * c --> 2nd part of RHS (c - b)
537  *----------------------------------------------------------------------------*/
538 
539 inline static void
540 _fw_and_bw_lu33(const cs_real_t mat[],
541  cs_real_t x[restrict],
542  const cs_real_t b[restrict],
543  const cs_real_t c[restrict])
544 {
545  cs_real_t aux[3];
546 
547  aux[0] = (c[0] - b[0]);
548  aux[1] = (c[1] - b[1]) - aux[0]*mat[3];
549  aux[2] = (c[2] - b[2]) - aux[0]*mat[6] - aux[1]*mat[7];
550 
551  x[2] = aux[2]/mat[8];
552  x[1] = (aux[1] - mat[5]*x[2])/mat[4];
553  x[0] = (aux[0] - mat[1]*x[1] - mat[2]*x[2])/mat[0];
554 }
555 
556 /*----------------------------------------------------------------------------
557  * Block Jacobi 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 --> 1st part of RHS (c - b)
565  * c --> 2nd part of RHS (c - b)
566  *----------------------------------------------------------------------------*/
567 
568 inline static void
569 _fw_and_bw_lu(const cs_real_t mat[],
570  int db_size,
571  cs_real_t x[restrict],
572  const cs_real_t b[restrict],
573  const cs_real_t c[restrict])
574 {
575  assert(db_size <= DB_SIZE_MAX);
576  cs_real_t aux[DB_SIZE_MAX];
577 
578  /* forward */
579  for (int ii = 0; ii < db_size; ii++) {
580  aux[ii] = (c[ii] - b[ii]);
581  for (int jj = 0; jj < ii; jj++) {
582  aux[ii] -= aux[jj]*mat[ii*db_size + jj];
583  }
584  }
585 
586  /* backward */
587  for (int ii = db_size - 1; ii >= 0; ii-=1) {
588  x[ii] = aux[ii];
589  for (int jj = db_size - 1; jj > ii; jj-=1) {
590  x[ii] -= x[jj]*mat[ii*db_size + jj];
591  }
592  x[ii] /= mat[ii*(db_size + 1)];
593  }
594 }
595 
596 /*----------------------------------------------------------------------------
597  * Block Gauss-Seidel utilities.
598  * Compute forward and backward to solve an LU P*P system.
599  *
600  * parameters:
601  * mat <-- P*P*dim matrix
602  * db_size <-- matrix size
603  * x --> solution
604  * b <-> RHS in, work array
605  *----------------------------------------------------------------------------*/
606 
607 inline static void
608 _fw_and_bw_lu_gs(const cs_real_t mat[],
609  int db_size,
610  cs_real_t x[restrict],
611  const cs_real_t b[restrict])
612 {
613  assert(db_size <= DB_SIZE_MAX);
614 
615  /* forward */
616  for (int ii = 0; ii < db_size; ii++) {
617  x[ii] = b[ii];
618  for (int jj = 0; jj < ii; jj++)
619  x[ii] -= x[jj]*mat[ii*db_size + jj];
620  }
621 
622  /* backward */
623  for (int ii = db_size - 1; ii >= 0; ii--) {
624  for (int jj = db_size - 1; jj > ii; jj--)
625  x[ii] -= x[jj]*mat[ii*db_size + jj];
626  x[ii] /= mat[ii*(db_size + 1)];
627  }
628 }
629 
630 /*============================================================================
631  * Public function definitions
632  *============================================================================*/
633 
634 /*----------------------------------------------------------------------------*/
648 /*----------------------------------------------------------------------------*/
649 
650 void
651 cs_sles_it_convergence_init(cs_sles_it_convergence_t *convergence,
652  const char *solver_name,
653  int verbosity,
654  unsigned n_iter_max,
655  double precision,
656  double r_norm,
657  double *residue);
658 
659 /*----------------------------------------------------------------------------
660  * Setup context for iterative linear solver.
661  *
662  * This function is common to most solvers
663  *
664  * parameters:
665  * c <-> pointer to solver context info
666  * name <-- pointer to system name
667  * a <-- matrix
668  * verbosity <-- verbosity level
669  * diag_block_size <-- diagonal block size
670  * block_nn_inverse <-- if diagonal block size is 3 or 6, compute inverse of
671  * block if true, inverse of block diagonal otherwise
672  *----------------------------------------------------------------------------*/
673 
674 void
675 cs_sles_it_setup_priv(cs_sles_it_t *c,
676  const char *name,
677  const cs_matrix_t *a,
678  int verbosity,
679  int diag_block_size,
680  bool block_nn_inverse);
681 
684 /*----------------------------------------------------------------------------*/
685 
687 
688 #endif /* __CS_SLES_IT_PRIV_H__ */
#define restrict
Definition: cs_defs.h:142
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_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
cs_sles_it_type_t
Definition: cs_sles_it.h:55
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
#define BEGIN_C_DECLS
Definition: cs_defs.h:510
struct _cs_sles_pc_t cs_sles_pc_t
Definition: cs_sles_pc.h:66
struct _cs_sles_it_t cs_sles_it_t
Definition: cs_sles_it.h:86
double cs_real_t
Floating-point value.
Definition: cs_defs.h:322
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
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:143
struct _cs_matrix_t cs_matrix_t
Definition: cs_matrix.h:93
struct _cs_time_plot_t cs_time_plot_t
Definition: cs_time_plot.h:48
cs_sles_convergence_state_t
Convergence status indicator.
Definition: cs_sles.h:56
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
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
struct _cs_sles_it_convergence_t cs_sles_it_convergence_t
Definition: cs_sles_it.h:90
double precision, save a
Definition: cs_fuel_incl.f90:146
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
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:316
#define END_C_DECLS
Definition: cs_defs.h:511
Definition: cs_timer.h:55
double precision, save b
Definition: cs_fuel_incl.f90:146