programmer's documentation
cs_sles_it.h
Go to the documentation of this file.
1 #ifndef __CS_SLES_IT_H__
2 #define __CS_SLES_IT_H__
3 
4 /*============================================================================
5  * Sparse Linear Equation Solvers
6  *============================================================================*/
7 
8 /*
9  This file is part of Code_Saturne, a general-purpose CFD tool.
10 
11  Copyright (C) 1998-2015 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 /*----------------------------------------------------------------------------
31  * Local headers
32  *----------------------------------------------------------------------------*/
33 
34 #include "cs_base.h"
35 #include "cs_halo_perio.h"
36 #include "cs_matrix.h"
37 #include "cs_time_plot.h"
38 #include "cs_sles.h"
39 
40 /*----------------------------------------------------------------------------*/
41 
43 
44 /*============================================================================
45  * Macro definitions
46  *============================================================================*/
47 
48 /*============================================================================
49  * Type definitions
50  *============================================================================*/
51 
52 /*----------------------------------------------------------------------------
53  * Solver types
54  *----------------------------------------------------------------------------*/
55 
56 typedef enum {
57 
58  CS_SLES_PCG, /* Preconditionned conjugate gradient */
59  CS_SLES_JACOBI, /* Jacobi */
60  CS_SLES_BICGSTAB, /* Bi-conjugate gradient stabilized */
61  CS_SLES_BICGSTAB2, /* Bi-conjugate gradient stabilized - 2*/
62  CS_SLES_GMRES, /* Generalized minimal residual */
63  CS_SLES_P_GAUSS_SEIDEL, /* Process-local Gauss-Seidel */
64  CS_SLES_PCR3, /* 3-layer conjugate residual */
65  CS_SLES_N_IT_TYPES /* Number of resolution algorithms */
66 
68 
69 /* Iterative linear solver context (opaque) */
70 
71 typedef struct _cs_sles_it_t cs_sles_it_t;
72 
73 /*============================================================================
74  * Global variables
75  *============================================================================*/
76 
77 /* Short names for matrix types */
78 
79 extern const char *cs_sles_it_type_name[];
80 
81 /*=============================================================================
82  * Public function prototypes
83  *============================================================================*/
84 
85 /*----------------------------------------------------------------------------
86  * Define and associate an iterative sparse linear system solver
87  * for a given field or equation name.
88  *
89  * If this system did not previously exist, it is added to the list of
90  * "known" systems. Otherwise, its definition is replaced by the one
91  * defined here.
92  *
93  * This is a utility function: if finer control is needed, see
94  * cs_sles_define() and cs_sles_it_create().
95  *
96  * Note that this function returns a pointer directly to the iterative solver
97  * management structure. This may be used to set further options,
98  * for example using cs_sles_set_plot_options(). If needed, cs_sles_find()
99  * may be used to obtain a pointer to the matching cs_sles_t container.
100  *
101  * parameters:
102  * f_id <-- associated field id, or < 0
103  * name <-- associated name if f_id < 0, or NULL
104  * solver_type <-- type of solver (PCG, Jacobi, ...)
105  * poly_degree <-- preconditioning polynomial degree
106  * (0: diagonal; -1: non-preconditioned)
107  * n_max_iter <-- maximum number of iterations
108  *
109  * returns:
110  * pointer to newly created iterative solver info object.
111  *----------------------------------------------------------------------------*/
112 
113 cs_sles_it_t *
114 cs_sles_it_define(int f_id,
115  const char *name,
116  cs_sles_it_type_t solver_type,
117  int poly_degree,
118  int n_max_iter);
119 
120 /*----------------------------------------------------------------------------
121  * Create iterative sparse linear system solver info and context.
122  *
123  * parameters:
124  * solver_type <-- type of solver (PCG, Jacobi, ...)
125  * poly_degree <-- preconditioning polynomial degree
126  * (0: diagonal; -1: non-preconditioned)
127  * n_max_iter <-- maximum number of iterations
128  * update_stats <-- automatic solver statistics indicator
129  *
130  * returns:
131  * pointer to newly created solver info object.
132  *----------------------------------------------------------------------------*/
133 
134 cs_sles_it_t *
136  int poly_degree,
137  int n_max_iter,
138  bool update_stats);
139 
140 /*----------------------------------------------------------------------------
141  * Destroy iterative sparse linear system solver info and context.
142  *
143  * parameters:
144  * context <-> pointer to iterative sparse linear solver info
145  * (actual type: cs_sles_it_t **)
146  *----------------------------------------------------------------------------*/
147 
148 void
149 cs_sles_it_destroy(void **context);
150 
151 /*----------------------------------------------------------------------------
152  * Create iterative sparse linear system solver info and context
153  * based on existing info and context.
154  *
155  * parameters:
156  * context <-- pointer to reference info and context
157  * (actual type: cs_sles_it_t *)
158  *
159  * returns:
160  * pointer to newly created solver info object
161  * (actual type: cs_sles_it_t *)
162  *----------------------------------------------------------------------------*/
163 
164 void *
165 cs_sles_it_copy(const void *context);
166 
167 /*----------------------------------------------------------------------------
168  * Setup iterative sparse linear equation solver.
169  *
170  * parameters:
171  * context <-> pointer to iterative sparse linear solver info
172  * (actual type: cs_sles_it_t *)
173  * name <-- pointer to system name
174  * a <-- associated matrix
175  * verbosity <-- verbosity level
176  *----------------------------------------------------------------------------*/
177 
178 void
179 cs_sles_it_setup(void *context,
180  const char *name,
181  const cs_matrix_t *a,
182  int verbosity);
183 
184 /*----------------------------------------------------------------------------
185  * Call iterative sparse linear equation solver.
186  *
187  * parameters:
188  * context <-> pointer to iterative sparse linear solver info
189  * (actual type: cs_sles_it_t *)
190  * name <-- pointer to system name
191  * a <-- matrix
192  * verbosity <-- verbosity level
193  * rotation_mode <-- halo update option for rotational periodicity
194  * precision <-- solver precision
195  * r_norm <-- residue normalization
196  * n_iter --> number of iterations
197  * residue --> residue
198  * rhs <-- right hand side
199  * vx <-> system solution
200  * aux_size <-- number of elements in aux_vectors (in bytes)
201  * aux_vectors --- optional working area (internal allocation if NULL)
202  *
203  * returns:
204  * convergence state
205  *----------------------------------------------------------------------------*/
206 
208 cs_sles_it_solve(void *context,
209  const char *name,
210  const cs_matrix_t *a,
211  int verbosity,
212  cs_halo_rotation_t rotation_mode,
213  double precision,
214  double r_norm,
215  int *n_iter,
216  double *residue,
217  const cs_real_t *rhs,
218  cs_real_t *vx,
219  size_t aux_size,
220  void *aux_vectors);
221 
222 /*----------------------------------------------------------------------------
223  * Free iterative sparse linear equation solver setup context.
224  *
225  * This function frees resolution-related data, such as
226  * buffers and preconditioning but does not free the whole context,
227  * as info used for logging (especially performance data) is maintained.
228 
229  * parameters:
230  * context <-> pointer to iterative sparse linear solver info
231  * (actual type: cs_sles_it_t *)
232  *----------------------------------------------------------------------------*/
233 
234 void
235 cs_sles_it_free(void *context);
236 
237 /*----------------------------------------------------------------------------
238  * Log sparse linear equation solver info.
239  *
240  * parameters:
241  * context <-> pointer to iterative sparse linear solver info
242  * (actual type: cs_sles_it_t *)
243  * log_type <-- log type
244  *----------------------------------------------------------------------------*/
245 
246 void
247 cs_sles_it_log(const void *context,
248  cs_log_t log_type);
249 
250 /*----------------------------------------------------------------------------
251  * Return iterative solver type.
252  *
253  * parameters:
254  * context <-- pointer to iterative solver info and context
255  *
256  * returns:
257  * selected solver type
258  *----------------------------------------------------------------------------*/
259 
261 cs_sles_it_get_type(const cs_sles_it_t *context);
262 
263 /*----------------------------------------------------------------------------
264  * Return the initial residue for the previous solve operation with a solver.
265  *
266  * This is useful for convergence tests when this solver is used as
267  * a preconditioning smoother.
268  *
269  * This operation is only valid between calls to cs_sles_it_setup()
270  * (or cs_sles_it_solve()) and cs_sles_it_free().
271  * It returns -1 otherwise.
272  *
273  * parameters:
274  * context <-- pointer to iterative solver info and context
275  *
276  * returns:
277  * initial residue from last call to \ref cs_sles_solve with this solver
278  *----------------------------------------------------------------------------*/
279 
280 double
282 
283 /*----------------------------------------------------------------------------
284  * Associate a similar info and context object with which some setup
285  * data may be shared.
286  *
287  * This is especially useful for sharing preconditioning data between
288  * similar solver contexts (for example ascending and descending multigrid
289  * smoothers based on the same matrix).
290  *
291  * For preconditioning data to be effectively shared, cs_sles_it_setup()
292  * (or cs_sles_it_solve()) must be called on "shareable" before being
293  * called on "context" (without cs_sles_it_free() being called in between,
294  * of course).
295  *
296  * It is the caller's responsibility to ensure the context is not used
297  * for a cs_sles_it_setup() or cs_sles_it_solve() operation after the
298  * shareable object has been destroyed (normally by cs_sles_it_destroy()).
299  *
300  * parameters:
301  * context <-> pointer to iterative sparse linear system solver info
302  * shareable <-- pointer to iterative solver info and context
303  * whose context may be shared
304  *----------------------------------------------------------------------------*/
305 
306 void
308  const cs_sles_it_t *shareable);
309 
310 #if defined(HAVE_MPI)
311 
312 /*----------------------------------------------------------------------------
313  * Set MPI communicator for dot products.
314  *
315  * parameters:
316  * context <-> pointer to iterative sparse linear system solver info
317  * comm <-- MPI communicator
318  *----------------------------------------------------------------------------*/
319 
320 void
321 cs_sles_it_set_mpi_reduce_comm(cs_sles_it_t *context,
322  MPI_Comm comm);
323 
324 #endif /* defined(HAVE_MPI) */
325 
326 /*----------------------------------------------------------------------------
327  * Assign ordering to iterative solver.
328  *
329  * The solver context takes ownership of the order array (i.e. it will
330  * handle its later deallocation).
331  *
332  * This is useful only for Block Gauss-Seidel.
333  *
334  * parameters:
335  * context <-> pointer to iterative solver info and context
336  * order <-> pointer to ordering array
337  *----------------------------------------------------------------------------*/
338 
339 void
341  cs_lnum_t **order);
342 
343 /*----------------------------------------------------------------------------
344  * Assign symmetric option to iterative solver.
345  *
346  * This is useful only for Process-local Gauss-Seidel.
347  *
348  * parameters:
349  * context <-> pointer to iterative solver info and context
350  * symmetric <-- symmetric if true
351  *----------------------------------------------------------------------------*/
352 
353 void
355  bool symmetric);
356 
357 /*----------------------------------------------------------------------------
358  * Query mean number of rows under which Conjugate Gradient algorithm
359  * uses the single-reduction variant.
360  *
361  * The single-reduction variant requires only one parallel sum per
362  * iteration (instead of 2), at the cost of additional vector operations,
363  * so it tends to be more expensive when the number of matrix rows per
364  * MPI rank is high, then becomes cheaper when the MPI latency cost becomes
365  * more significant.
366  *
367  * This option is ignored for non-parallel runs, so 0 is returned.
368  *
369  * return:
370  * mean number of rows per active rank under which the
371  * single-reduction variant will be used
372  *----------------------------------------------------------------------------*/
373 
374 cs_lnum_t
376 
377 /*----------------------------------------------------------------------------
378  * Set mean number of rows under which Conjugate Gradient algorithm
379  * should use the single-reduction variant.
380  *
381  * The single-reduction variant requires only one parallel sum per
382  * iteration (instead of 2), at the cost of additional vector operations,
383  * so it tends to be more expensive when the number of matrix rows per
384  * MPI rank is high, then becomes cheaper when the MPI latency cost becomes
385  * more significant.
386  *
387  * This option is ignored for non-parallel runs.
388  *
389  * parameters:
390  * threshold <-- mean number of rows per active rank under which the
391  * single-reduction variant will be used
392  *----------------------------------------------------------------------------*/
393 
394 void
396 
397 /*----------------------------------------------------------------------------
398  * Log the current global settings relative to parallelism.
399  *----------------------------------------------------------------------------*/
400 
401 void
403 
404 /*----------------------------------------------------------------------------
405  * Error handler for iterative sparse linear equation solver.
406  *
407  * In case of divergence or breakdown, this error handler outputs
408  * postprocessing data to assist debugging, then aborts the run.
409  * It does nothing in case the maximum iteration count is reached.
410  *
411  * parameters:
412  * context <-> pointer to iterative sparse linear system solver info
413  * (actual type: cs_sles_it_t *)
414  * state <-- convergence state
415  * name <-- pointer to name of linear system
416  * a <-- matrix
417  * rotation_mode <-- halo update option for rotational periodicity
418  * rhs <-- right hand side
419  * vx <-> system solution
420  */
421 /*----------------------------------------------------------------------------*/
422 
423 void
424 cs_sles_it_error_post_and_abort(void *context,
426  const char *name,
427  const cs_matrix_t *a,
428  cs_halo_rotation_t rotation_mode,
429  const cs_real_t *rhs,
430  cs_real_t *vx);
431 
432 /*----------------------------------------------------------------------------
433  * Set plotting options for an iterative sparse linear equation solver.
434  *
435  * parameters:
436  * context <-> pointer to iterative solver info and context
437  * base_name <-- base plot name to activate, NULL otherwise
438  * use_iteration <-- if true, use iteration as time stamp
439  * otherwise, use wall clock time
440  *----------------------------------------------------------------------------*/
441 
442 void
444  const char *base_name,
445  bool use_iteration);
446 
447 /*----------------------------------------------------------------------------
448  * Assign existing time plot to iterative sparse linear equation solver.
449  *
450  * This is useful mainly when a time plot has a longer lifecycle than
451  * the linear solver context, such as inside a multigrid solver.
452  *
453  * parameters:
454  * context <-> pointer to iterative solver info and context
455  * time_plot <-- pointer to time plot structure
456  * time_stamp <-- associated time stamp
457  *----------------------------------------------------------------------------*/
458 
459 void
461  cs_time_plot_t *time_plot,
462  int time_stamp);
463 
464 /*----------------------------------------------------------------------------*/
465 
467 
468 #endif /* __CS_SLES_IT_H__ */
void cs_sles_it_setup(void *context, const char *name, const cs_matrix_t *a, int verbosity)
Setup iterative sparse linear equation solver.
Definition: cs_sles_it.c:3967
cs_halo_rotation_t
Definition: cs_halo.h:59
cs_sles_it_type_t
Iterative solver types.
Definition: cs_sles_it.h:56
Definition: cs_sles_it.h:62
void cs_sles_it_error_post_and_abort(void *context, cs_sles_convergence_state_t state, const char *name, const cs_matrix_t *a, cs_halo_rotation_t rotation_mode, const cs_real_t *rhs, cs_real_t *vx)
Error handler for iterative sparse linear equation solver.
Definition: cs_sles_it.c:4558
double cs_sles_it_get_last_initial_residue(const cs_sles_it_t *context)
Return the initial residue for the previous solve operation with a solver.
Definition: cs_sles_it.c:4332
#define BEGIN_C_DECLS
Definition: cs_defs.h:419
void cs_sles_it_set_symmetric(cs_sles_it_t *context, bool symmetric)
Assign symmetric option to iterative solver.
Definition: cs_sles_it.c:4458
cs_lnum_t cs_sles_it_get_pcg_single_reduction(void)
Query mean number of rows under which Conjugate Gradient algorithm uses the single-reduction variant...
Definition: cs_sles_it.c:4484
struct _cs_sles_it_t cs_sles_it_t
Definition: cs_sles_it.h:71
struct _cs_matrix_t cs_matrix_t
Definition: cs_matrix.h:86
void cs_sles_it_destroy(void **context)
Destroy iterative sparse linear system solver info and context.
Definition: cs_sles_it.c:3840
struct _cs_time_plot_t cs_time_plot_t
Definition: cs_time_plot.h:48
void cs_sles_it_free(void *context)
Free iterative sparse linear equation solver setup context.
Definition: cs_sles_it.c:4277
void cs_sles_it_log(const void *context, cs_log_t log_type)
Log sparse linear equation solver info.
Definition: cs_sles_it.c:3901
cs_sles_convergence_state_t
Convergence status indicator.
Definition: cs_sles.h:55
Definition: cs_sles_it.h:65
void cs_sles_it_assign_plot(cs_sles_it_t *context, cs_time_plot_t *time_plot, int time_stamp)
Assign existing time plot to iterative sparse linear equation solver.
Definition: cs_sles_it.c:4648
Definition: cs_sles_it.h:61
Definition: cs_sles_it.h:58
double precision, save a
Definition: cs_fuel_incl.f90:146
Definition: cs_sles_it.h:64
const char * cs_sles_it_type_name[]
cs_sles_convergence_state_t cs_sles_it_solve(void *context, const char *name, const cs_matrix_t *a, int verbosity, cs_halo_rotation_t rotation_mode, double precision, double r_norm, int *n_iter, double *residue, const cs_real_t *rhs, cs_real_t *vx, size_t aux_size, void *aux_vectors)
Call iterative sparse linear equation solver.
Definition: cs_sles_it.c:4008
cs_sles_it_type_t cs_sles_it_get_type(const cs_sles_it_t *context)
Return iterative solver type.
Definition: cs_sles_it.c:4307
cs_log_t
Definition: cs_log.h:47
void cs_sles_it_log_parallel_options(void)
Log the current global settings relative to parallelism.
Definition: cs_sles_it.c:4526
void cs_sles_it_set_plot_options(cs_sles_it_t *context, const char *base_name, bool use_iteration)
Set plotting options for an iterative sparse linear equation solver.
Definition: cs_sles_it.c:4604
void cs_sles_it_set_pcg_single_reduction(cs_lnum_t threshold)
Set mean number of rows under which Conjugate Gradient algorithm should use the single-reduction vari...
Definition: cs_sles_it.c:4512
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:292
Definition: cs_sles_it.h:60
void cs_sles_it_set_shareable(cs_sles_it_t *context, const cs_sles_it_t *shareable)
Associate a similar info and context object with which some setup data may be shared.
Definition: cs_sles_it.c:4366
#define END_C_DECLS
Definition: cs_defs.h:420
double cs_real_t
Definition: cs_defs.h:296
void cs_sles_it_assign_order(cs_sles_it_t *context, cs_lnum_t **order)
Assign ordering to iterative solver.
Definition: cs_sles_it.c:4424
cs_sles_it_t * cs_sles_it_create(cs_sles_it_type_t solver_type, int poly_degree, int n_max_iter, bool update_stats)
Create iterative sparse linear system solver info and context.
Definition: cs_sles_it.c:3783
cs_sles_it_t * cs_sles_it_define(int f_id, const char *name, cs_sles_it_type_t solver_type, int poly_degree, int n_max_iter)
Define and associate an iterative sparse linear system solver for a given field or equation name...
Definition: cs_sles_it.c:3738
Definition: cs_sles_it.h:63
void * cs_sles_it_copy(const void *context)
Create iterative sparse linear system solver info and context based on existing info and context...
Definition: cs_sles_it.c:3872
Definition: cs_sles_it.h:59