programmer's documentation
cs_sla.h
Go to the documentation of this file.
1 #ifndef __CS_SLA_H__
2 #define __CS_SLA_H__
3 
4 /*============================================================================
5  * Sparse linear algebra routines
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  * Standard C library headers
32  *----------------------------------------------------------------------------*/
33 
34 #include <stdio.h>
35 
36 /*----------------------------------------------------------------------------
37  * Local headers
38  *----------------------------------------------------------------------------*/
39 
40 #include "cs_base.h"
41 #include "cs_cdo.h"
42 #include "cs_cdo_toolbox.h"
43 #include "cs_param.h"
44 
45 /*----------------------------------------------------------------------------*/
46 
48 
49 /*============================================================================
50  * Macro definitions
51  *============================================================================*/
52 
53 /* Matrix flag */
54 #define CS_SLA_MATRIX_SYM (1 << 0) /* 1: symmetric */
55 #define CS_SLA_MATRIX_SORTED (1 << 1) /* 2: sorted */
56 #define CS_SLA_MATRIX_SHARED (1 << 2) /* 4: share pattern */
57 
58 /*============================================================================
59  * Type definitions for matrices
60  *============================================================================*/
61 
62 typedef enum {
63 
69 
71 
72 typedef struct {
73 
74  /* Stencil */
77  double stencil_mean;
78 
79  /* Fill-in */
80  size_t nnz;
81  double fillin;
82 
84 
85 typedef struct _sla_mat_prop_t cs_sla_mat_prop_t; /* Private definition */
86 
87 typedef struct {
88 
90  cs_sla_mat_prop_t *properties;
91  int flag; /* Symmetric, sorted, shared... */
92 
93  int stride; /* Number of entries in "val" for each couple A(i,j) */
94  int n_rows;
95  int n_cols;
96 
99 
100  short int *sgn; /* Only for DEC type: -1 or 1 according to orientation */
101  double *val; /* Only for CSR and MSR type */
102 
103  cs_lnum_t *didx; /* Diagonal index: used to flag diagonal index and to
104  separate lower and upper part */
105  double *diag; /* Used in MSR type but can also be used for a
106  given preconditioning technique */
107 
109 
110 /*============================================================================
111  * Type definitions for iterative solvers
112  *============================================================================*/
113 
114 typedef struct {
115 
116  int code; // Convergence code
117  int iter; // Current iteration
118  double residual; // Current residual norm computed
119 
121 
122 
123 /*============================================================================
124  * Public function prototypes for SLA matrices
125  *============================================================================*/
126 
127 /*----------------------------------------------------------------------------*/
137 /*----------------------------------------------------------------------------*/
138 
139 void
140 cs_sla_bwrite(const char *name,
141  const cs_sla_matrix_t *m,
142  const double *rhs,
143  const double *sol);
144 
145 /*----------------------------------------------------------------------------*/
155 /*----------------------------------------------------------------------------*/
156 
157 void
158 cs_sla_bread(const char *name,
159  cs_sla_matrix_t **p_mat,
160  double *p_rhs[],
161  double *p_sol[]);
162 
163 /*----------------------------------------------------------------------------*/
170 /*----------------------------------------------------------------------------*/
171 
172 void
174 
175 /*----------------------------------------------------------------------------*/
181 /*----------------------------------------------------------------------------*/
182 
183 void
185 
186 /*----------------------------------------------------------------------------*/
192 /*----------------------------------------------------------------------------*/
193 
194 void
196 
197 /*----------------------------------------------------------------------------*/
212 /*----------------------------------------------------------------------------*/
213 
215 cs_sla_matrix_pack(cs_lnum_t n_final_rows,
216  cs_lnum_t n_final_cols,
217  const cs_sla_matrix_t *init,
218  const cs_lnum_t *row_z2i_ids,
219  const cs_lnum_t *col_i2z_ids,
220  bool keep_sym);
221 
222 /*----------------------------------------------------------------------------*/
228 /*----------------------------------------------------------------------------*/
229 
230 void
232 
233 /*----------------------------------------------------------------------------*/
240 /*----------------------------------------------------------------------------*/
241 
242 void
244  double *p_diag[]);
245 
246 /*----------------------------------------------------------------------------*/
252 /*----------------------------------------------------------------------------*/
253 
254 void
256 
257 /*----------------------------------------------------------------------------
258  * Create a new matrix from a block description (not done for all matrix
259  * combinations, see below).
260  *
261  * |A | B|
262  * M = |--|--|
263  * |C | D|
264  *----------------------------------------------------------------------------*/
265 
268  const cs_sla_matrix_t *B,
269  const cs_sla_matrix_t *C,
270  const cs_sla_matrix_t *D,
271  bool sym);
272 
273 /* Operations on matrices */
274 void cs_sla_matvec(const cs_sla_matrix_t *m,
275  const double v[],
276  double *inout[],
277  bool reset);
278 
279 void cs_sla_amxby(double alpha,
280  const cs_sla_matrix_t *m,
281  const double x[],
282  double beta,
283  const double y[],
284  double *inout[]);
285 
286 /*----------------------------------------------------------------------------*/
303 /*----------------------------------------------------------------------------*/
304 
305 void
307  const cs_sla_matrix_t *B,
308  const cs_sla_matrix_t *C,
309  const cs_sla_matrix_t *D,
310  const double X[],
311  const double Y[],
312  double *F[],
313  double *G[],
314  bool reset);
315 
316 /*----------------------------------------------------------------------------*/
327 /*----------------------------------------------------------------------------*/
328 
330 cs_sla_matrix_add(double alpha,
331  const cs_sla_matrix_t *a,
332  double beta,
333  const cs_sla_matrix_t *b);
334 
335 /*----------------------------------------------------------------------------*/
344 /*----------------------------------------------------------------------------*/
345 
348  const cs_sla_matrix_t *b);
349 
350 /*----------------------------------------------------------------------------*/
363 /*----------------------------------------------------------------------------*/
364 
367  const cs_sla_matrix_t *a,
368  double beta,
369  const cs_sla_matrix_t *bt,
370  const cs_sla_matrix_t *b);
371 
372 /*----------------------------------------------------------------------------*/
383 /*----------------------------------------------------------------------------*/
384 
386  const double D[],
387  const cs_sla_matrix_t *A,
388  int *w);
389 
391  int iter_max,
392  double epsilon,
393  double *lambda);
394 
395 double cs_sla_get_matrix_norm(const cs_sla_matrix_t *m);
396 
397 /*----------------------------------------------------------------------------*/
409 /*----------------------------------------------------------------------------*/
410 
412 cs_sla_matrix_create(int n_rows,
413  int n_cols,
414  int stride,
416  bool sym);
417 
418 /*----------------------------------------------------------------------------*/
428 /*----------------------------------------------------------------------------*/
429 
433  int stride);
434 
435 /*----------------------------------------------------------------------------*/
444 /*----------------------------------------------------------------------------*/
445 
448  bool shared);
449 
450 /*----------------------------------------------------------------------------*/
458 /*----------------------------------------------------------------------------*/
459 
462 
463 /*----------------------------------------------------------------------------*/
471 /*----------------------------------------------------------------------------*/
472 
475 
476 /*----------------------------------------------------------------------------*/
484 /*----------------------------------------------------------------------------*/
485 
486 void
488  double eps);
489 
490 /*----------------------------------------------------------------------------*/
498 /*----------------------------------------------------------------------------*/
499 
500 size_t
502 
503 /*----------------------------------------------------------------------------*/
509 /*----------------------------------------------------------------------------*/
510 
513 
514 /*----------------------------------------------------------------------------*/
522 /*----------------------------------------------------------------------------*/
523 
524 void
525 cs_sla_matrix_resume(const char *name,
526  FILE *f,
527  const cs_sla_matrix_t *m);
528 
529 /*----------------------------------------------------------------------------*/
537 /*----------------------------------------------------------------------------*/
538 
539 void
540 cs_sla_matrix_dump(const char *name,
541  FILE *f,
542  const cs_sla_matrix_t *m);
543 
544 /*----------------------------------------------------------------------------*/
553 /*----------------------------------------------------------------------------*/
554 
555 void
556 cs_sla_system_dump(const char *name,
557  FILE *f,
558  const cs_sla_matrix_t *m,
559  const double *rhs);
560 
561 /*----------------------------------------------------------------------------*/
570 /*----------------------------------------------------------------------------*/
571 
572 void
574  cs_sla_matrix_t *ass);
575 
576 /*----------------------------------------------------------------------------*/
577 
579 
580 #endif /* __CS_SLA_H__ */
int stencil_min
Definition: cs_sla.h:75
int iter
Definition: cs_sla.h:117
void cs_sla_system_dump(const char *name, FILE *f, const cs_sla_matrix_t *m, const double *rhs)
Dump a cs_sla_matrix_t structure and its related right-hand side.
Definition: cs_sla_matrix.c:3464
cs_sla_matrix_t * cs_sla_matrix_copy(const cs_sla_matrix_t *a, bool shared)
Create a new matrix structure from the copy of an existing one.
Definition: cs_sla.h:64
double * diag
Definition: cs_sla.h:105
void cs_sla_bwrite(const char *name, const cs_sla_matrix_t *m, const double *rhs, const double *sol)
Write in binary format a matrix in CSR format, its righ hand side and the solution.
Definition: cs_sla_matrix.c:2686
double residual
Definition: cs_sla.h:118
void cs_sla_assemble_msr(const cs_toolbox_locmat_t *loc, cs_sla_matrix_t *ass)
Assemble a MSR matrix from local contributions –> We assume that the local matrices are symmetric –...
Definition: cs_sla_matrix.c:3568
double * val
Definition: cs_sla.h:101
cs_sla_matrix_t * cs_sla_matrix_combine(double alpha, const cs_sla_matrix_t *a, double beta, const cs_sla_matrix_t *bt, const cs_sla_matrix_t *b)
Specific matrix multiplication. Compute Bt * beta * B + alpha * A alpha and beta are scalar...
Definition: cs_sla_matrix.c:1732
void cs_sla_matvec(const cs_sla_matrix_t *m, const double v[], double *inout[], bool reset)
cs_sla_matrix_t * cs_sla_matrix_create_from_pattern(const cs_sla_matrix_t *ref, cs_sla_matrix_type_t type, int stride)
Create a cs_sla_matrix_t structure from an existing one.
Definition: cs_sla_matrix.c:2838
cs_sla_matrix_t * cs_sla_matrix_multiply(const cs_sla_matrix_t *a, const cs_sla_matrix_t *b)
Main subroutine to multiply two sparse matrices a and b. c= a*b.
Definition: cs_sla_matrix.c:1651
void cs_sla_matrix_clean(cs_sla_matrix_t *m, double eps)
Remove entries in a cs_sla_matrix_t structure below a given threshold. |a(i,j)| < eps * max|a(i...
Definition: cs_sla_matrix.c:3118
cs_sla_matrix_t * cs_sla_matrix_add(double alpha, const cs_sla_matrix_t *a, double beta, const cs_sla_matrix_t *b)
Add two sparse matrices a and b. c = alpha*a + beta*b.
Definition: cs_sla_matrix.c:1507
size_t nnz
Definition: cs_sla.h:80
cs_sla_matrix_t * cs_sla_multiply_AtDA(const cs_sla_matrix_t *At, const double D[], const cs_sla_matrix_t *A, int *w)
Compute the product C = At * Diag * A.
Definition: cs_sla_matrix.c:1879
#define BEGIN_C_DECLS
Definition: cs_defs.h:419
cs_sla_matrix_t * cs_sla_matrix_transpose(const cs_sla_matrix_t *mat)
Transpose a cs_sla_matrix_t structure.
Definition: cs_sla_matrix.c:2956
int stride
Definition: cs_sla.h:93
Definition: cs_field_pointer.h:81
double precision, dimension(ncharm), save beta
Definition: cpincl.f90:97
Definition: cs_sla.h:114
cs_lnum_t * idx
Definition: cs_sla.h:97
void cs_sla_matrix_dump(const char *name, FILE *f, const cs_sla_matrix_t *m)
Dump a cs_sla_matrix_t structure.
Definition: cs_sla_matrix.c:3370
double stencil_mean
Definition: cs_sla.h:77
void cs_sla_matrix_csr2msr(cs_sla_matrix_t *a)
Change matrix representation from CSR to MSR.
Definition: cs_sla_matrix.c:2191
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:108
short int * sgn
Definition: cs_sla.h:100
double cs_sla_get_matrix_norm(const cs_sla_matrix_t *m)
Definition: cs_sla_matrix.c:1277
cs_sla_matrix_t * cs_sla_matrix_create(int n_rows, int n_cols, int stride, cs_sla_matrix_type_t type, bool sym)
Create a cs_sla_matrix_t structure.
int code
Definition: cs_sla.h:116
Definition: cs_sla.h:87
cs_sla_mat_prop_t * properties
Definition: cs_sla.h:90
void cs_sla_bread(const char *name, cs_sla_matrix_t **p_mat, double *p_rhs[], double *p_sol[])
Read from a binary file a matrix in CSR format, its righ hand side and the solution. Matrix must have a stride equal to 1.
Definition: cs_sla_matrix.c:2607
cs_sla_matrix_t * cs_sla_matrix_free(cs_sla_matrix_t *m)
Free a cs_sla_matrix_t structure.
Definition: cs_sla_matrix.c:3062
void cs_sla_amxby(double alpha, const cs_sla_matrix_t *m, const double x[], double beta, const double y[], double *inout[])
Definition: cs_sla_matrix.c:2007
int cs_sla_get_mat_spectral_radius(const cs_sla_matrix_t *matrix, int iter_max, double epsilon, double *lambda)
Definition: cs_sla_matrix.c:1311
cs_sla_matrix_t * cs_sla_matrix_pack(cs_lnum_t n_final_rows, cs_lnum_t n_final_cols, const cs_sla_matrix_t *init, const cs_lnum_t *row_z2i_ids, const cs_lnum_t *col_i2z_ids, bool keep_sym)
Build a new matrix resulting from the extraction of some listed rows and columns. The output is a new...
double precision, save a
Definition: cs_fuel_incl.f90:146
void matrix(const cs_int_t *const iconvp, const cs_int_t *const idiffp, const cs_int_t *const ndircp, const cs_int_t *const isym, const cs_real_t *const thetap, const cs_int_t *const imucpp, const cs_real_t coefbp[], const cs_real_t cofbfp[], const cs_real_t rovsdt[], const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], const cs_real_t xcpp[], cs_real_t da[], cs_real_t xa[])
Definition: cs_matrix_building.c:111
int n_cols
Definition: cs_sla.h:95
Definition: cs_sla.h:72
int n_rows
Definition: cs_sla.h:94
Definition: cs_sla.h:66
void cs_sla_matrix_msr2csr(cs_sla_matrix_t *a)
Change matrix representation from MSR to CSR.
Definition: cs_sla_matrix.c:2120
cs_sla_matrix_t * cs_sla_block2mat(cs_sla_matrix_t *A, const cs_sla_matrix_t *B, const cs_sla_matrix_t *C, const cs_sla_matrix_t *D, bool sym)
Definition: cs_cdo_toolbox.h:72
cs_sla_matrix_info_t cs_sla_matrix_analyse(const cs_sla_matrix_t *m)
Compute general information related to a cs_sla_matrix_t structure.
Definition: cs_sla_matrix.c:3215
cs_sla_matrix_type_t type
Definition: cs_sla.h:89
Definition: cs_sla.h:68
double fillin
Definition: cs_sla.h:81
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:292
void cs_sla_matrix_sort(cs_sla_matrix_t *m)
Sort each row by increasing colomn number.
Definition: cs_sla_matrix.c:2568
#define END_C_DECLS
Definition: cs_defs.h:420
cs_sla_matrix_type_t
Definition: cs_sla.h:62
void cs_sla_matrix_share2own(cs_sla_matrix_t *a)
Allocate its own pattern if shared.
Definition: cs_sla_matrix.c:2255
void cs_sla_matvec_block2(const cs_sla_matrix_t *A, const cs_sla_matrix_t *B, const cs_sla_matrix_t *C, const cs_sla_matrix_t *D, const double X[], const double Y[], double *F[], double *G[], bool reset)
Matrix block 2x2 multiply by a vector.
Definition: cs_sla.h:65
int flag
Definition: cs_sla.h:91
void cs_sla_matrix_get_diag(const cs_sla_matrix_t *m, double *p_diag[])
Get the diagonal entries of a given matrix.
Definition: cs_sla_matrix.c:2471
void cs_sla_matrix_resume(const char *name, FILE *f, const cs_sla_matrix_t *m)
Synthesis of a cs_sla_matrix_t structure.
Definition: cs_sla_matrix.c:3275
Definition: cs_field_pointer.h:70
Definition: cs_sla.h:67
size_t cs_sla_matrix_get_nnz(const cs_sla_matrix_t *m)
Retrieve the number of non-zeros (nnz) elements in a matrix.
Definition: cs_sla_matrix.c:3189
cs_lnum_t * didx
Definition: cs_sla.h:103
cs_lnum_t * col_id
Definition: cs_sla.h:98
int stencil_max
Definition: cs_sla.h:76
void cs_sla_matrix_diag_idx(cs_sla_matrix_t *m)
Build diagonal index.
Definition: cs_sla_matrix.c:2430
double precision, save b
Definition: cs_fuel_incl.f90:146