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 #define CS_SLA_MATRIX_INFO (1 << 3) /* 8: info struct. is set */
58 
59 /*============================================================================
60  * Type definitions for matrices
61  *============================================================================*/
62 
63 typedef enum {
64 
70 
72 
73 typedef struct {
74 
75  /* Stencil */
78  double stencil_mean;
79 
80  /* Fill-in */
81  size_t nnz;
82  double fillin;
83 
85 
86 typedef struct {
87 
88  cs_sla_matrix_type_t type; /* Type of matrix storage */
89  cs_sla_matrix_info_t info; /* Information on how matrix is filled */
90  int flag; /* Symmetric, sorted, shared... */
91 
92  int stride; /* Number of entries in "val" for each couple A(i,j) */
93  int n_rows;
94  int n_cols;
95 
98 
99  short int *sgn; /* Only for DEC type: -1 or 1 according to orientation */
100  double *val; /* Only for CSR and MSR type */
101 
102  cs_lnum_t *didx; /* Diagonal index: used to flag diagonal index and to
103  separate lower and upper part */
104  double *diag; /* Used in MSR type but can also be used for a
105  given preconditioning technique */
106 
108 
109 /*============================================================================
110  * Type definitions for iterative solvers
111  *============================================================================*/
112 
113 typedef struct {
114 
115  int code; // Convergence code
116  int iter; // Current iteration
117  double residual; // Current residual norm computed
118 
120 
121 /*============================================================================
122  * Public function prototypes for SLA matrices
123  *============================================================================*/
124 
125 /*----------------------------------------------------------------------------*/
137 /*----------------------------------------------------------------------------*/
138 
141  cs_lnum_t n_cols,
142  int stride,
144  bool sym);
145 
146 /*----------------------------------------------------------------------------*/
158 /*----------------------------------------------------------------------------*/
159 
163  int stride);
164 
165 /*----------------------------------------------------------------------------*/
179 /*----------------------------------------------------------------------------*/
180 
183  bool is_symmetric,
184  bool sorted_idx,
185  int stride);
186 
187 /*----------------------------------------------------------------------------*/
196 /*----------------------------------------------------------------------------*/
197 
200  bool shared);
201 
202 /*----------------------------------------------------------------------------*/
210 /*----------------------------------------------------------------------------*/
211 
214 
215 /*----------------------------------------------------------------------------*/
223 /*----------------------------------------------------------------------------*/
224 
227 
228 /*----------------------------------------------------------------------------*/
235 /*----------------------------------------------------------------------------*/
236 
237 void
239 
240 /*----------------------------------------------------------------------------*/
250 /*----------------------------------------------------------------------------*/
251 
252 void
254  double eps);
255 
256 /*----------------------------------------------------------------------------*/
264 /*----------------------------------------------------------------------------*/
265 
266 size_t
268 
269 /*----------------------------------------------------------------------------*/
275 /*----------------------------------------------------------------------------*/
276 
277 void
279 
280 /*----------------------------------------------------------------------------*/
287 /*----------------------------------------------------------------------------*/
288 
289 void
291  double *p_diag[]);
292 
293 /*----------------------------------------------------------------------------*/
299 /*----------------------------------------------------------------------------*/
300 
301 void
303 
304 /*----------------------------------------------------------------------------*/
311 /*----------------------------------------------------------------------------*/
312 
313 void
315 
316 /*----------------------------------------------------------------------------*/
322 /*----------------------------------------------------------------------------*/
323 
324 void
326 
327 /*----------------------------------------------------------------------------*/
333 /*----------------------------------------------------------------------------*/
334 
335 void
337 
338 /*----------------------------------------------------------------------------*/
344 /*----------------------------------------------------------------------------*/
345 
346 void
348 
349 /*----------------------------------------------------------------------------*/
359 /*----------------------------------------------------------------------------*/
360 
361 void
363  cs_sla_matrix_t *ass,
364  bool only_diag);
365 
366 /*----------------------------------------------------------------------------*/
374 /*----------------------------------------------------------------------------*/
375 
376 void
378  cs_sla_matrix_t *ass);
379 
380 /*----------------------------------------------------------------------------*/
395 /*----------------------------------------------------------------------------*/
396 
398 cs_sla_matrix_pack(cs_lnum_t n_final_rows,
399  cs_lnum_t n_final_cols,
400  const cs_sla_matrix_t *init,
401  const cs_lnum_t *row_z2i_ids,
402  const cs_lnum_t *col_i2z_ids,
403  bool keep_sym);
404 
405 /*----------------------------------------------------------------------------*/
416 /*----------------------------------------------------------------------------*/
417 
419 cs_sla_matrix_add(double alpha,
420  const cs_sla_matrix_t *a,
421  double beta,
422  const cs_sla_matrix_t *b);
423 
424 /*----------------------------------------------------------------------------*/
435 /*----------------------------------------------------------------------------*/
436 
437 void
439  const double v[],
440  double *inout[],
441  bool reset);
442 
443 /*----------------------------------------------------------------------------*/
455 /*----------------------------------------------------------------------------*/
456 
457 void
458 cs_sla_amxby(double alpha,
459  const cs_sla_matrix_t *m,
460  const double x[],
461  double beta,
462  const double y[],
463  double *inout[]);
464 
465 /*----------------------------------------------------------------------------*/
474 /*----------------------------------------------------------------------------*/
475 
478  const cs_sla_matrix_t *b);
479 
480 /*----------------------------------------------------------------------------*/
491 /*----------------------------------------------------------------------------*/
492 
495  const double D[],
496  const cs_sla_matrix_t *A,
497  cs_lnum_t *w);
498 
499 /*----------------------------------------------------------------------------*/
512 /*----------------------------------------------------------------------------*/
513 
516  const cs_sla_matrix_t *a,
517  double beta,
518  const cs_sla_matrix_t *bt,
519  const cs_sla_matrix_t *b);
520 
521 /*----------------------------------------------------------------------------*/
538 /*----------------------------------------------------------------------------*/
539 
540 void
542  const cs_sla_matrix_t *B,
543  const cs_sla_matrix_t *C,
544  const cs_sla_matrix_t *D,
545  const double X[],
546  const double Y[],
547  double *F[],
548  double *G[],
549  bool reset);
550 
551 /*----------------------------------------------------------------------------*/
561 /*----------------------------------------------------------------------------*/
562 
563 void
564 cs_sla_bwrite(const char *name,
565  const cs_sla_matrix_t *m,
566  const double *rhs,
567  const double *sol);
568 
569 /*----------------------------------------------------------------------------*/
579 /*----------------------------------------------------------------------------*/
580 
581 void
582 cs_sla_bread(const char *name,
583  cs_sla_matrix_t **p_mat,
584  double *p_rhs[],
585  double *p_sol[]);
586 
587 /*----------------------------------------------------------------------------*/
595 /*----------------------------------------------------------------------------*/
596 
597 void
598 cs_sla_matrix_summary(const char *name,
599  FILE *f,
600  cs_sla_matrix_t *m);
601 
602 /*----------------------------------------------------------------------------*/
610 /*----------------------------------------------------------------------------*/
611 
612 void
613 cs_sla_matrix_dump(const char *name,
614  FILE *f,
615  const cs_sla_matrix_t *m);
616 
617 /*----------------------------------------------------------------------------*/
626 /*----------------------------------------------------------------------------*/
627 
628 void
629 cs_sla_system_dump(const char *name,
630  FILE *f,
631  const cs_sla_matrix_t *m,
632  const double *rhs);
633 
634 /*----------------------------------------------------------------------------*/
635 
637 
638 #endif /* __CS_SLA_H__ */
int stencil_min
Definition: cs_sla.h:76
int iter
Definition: cs_sla.h:116
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:3567
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_matrix.c:1486
Definition: cs_sla.h:65
void cs_sla_matrix_summary(const char *name, FILE *f, cs_sla_matrix_t *m)
Synthesis of a cs_sla_matrix_t structure.
Definition: cs_sla_matrix.c:3374
cs_sla_matrix_t * cs_sla_multiply_AtDA(const cs_sla_matrix_t *At, const double D[], const cs_sla_matrix_t *A, cs_lnum_t *w)
Compute the product C = At * Diag * A.
Definition: cs_sla_matrix.c:2924
double * diag
Definition: cs_sla.h:104
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:3294
Definition: cs_cdo_toolbox.h:82
double residual
Definition: cs_sla.h:117
double * val
Definition: cs_sla.h:100
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 where alpha and beta are scalar...
Definition: cs_sla_matrix.c:2987
void cs_sla_matvec(const cs_sla_matrix_t *m, const double v[], double *inout[], bool reset)
Compute a matrix vector product. If inout is not allocated, allocation is done inside this function I...
Definition: cs_sla_matrix.c:2730
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:2845
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:1796
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:2584
size_t nnz
Definition: cs_sla.h:81
#define BEGIN_C_DECLS
Definition: cs_defs.h:429
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:1569
int stride
Definition: cs_sla.h:92
Definition: cs_field_pointer.h:82
double precision, dimension(ncharm), save beta
Definition: cpincl.f90:97
Definition: cs_sla.h:113
void cs_sla_matrix_set_info(cs_sla_matrix_t *m)
Compute general information related to a cs_sla_matrix_t structure and store it into a cs_sla_matrix_...
Definition: cs_sla_matrix.c:2071
cs_sla_matrix_info_t info
Definition: cs_sla.h:89
void cs_sla_assemble_msr(const cs_locmat_t *loc, cs_sla_matrix_t *ass)
Assemble a MSR matrix from local contributions –> We assume that the assembled matrix has its column...
Definition: cs_sla_matrix.c:2375
cs_lnum_t * idx
Definition: cs_sla.h:96
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:3468
cs_sla_matrix_t * cs_sla_matrix_create_msr_from_index(const cs_connect_index_t *connect_idx, bool is_symmetric, bool sorted_idx, int stride)
Create a cs_sla_matrix_t structure with MSR type from an existing connectivity index. Be aware of removing the diagonal entry in the connectivity index before the call to this routine.
Definition: cs_sla_matrix.c:1429
double stencil_mean
Definition: cs_sla.h:78
void cs_sla_matrix_csr2msr(cs_sla_matrix_t *a)
Change matrix representation from CSR to MSR.
Definition: cs_sla_matrix.c:2203
short int * sgn
Definition: cs_sla.h:99
int code
Definition: cs_sla.h:115
Definition: cs_sla.h:86
double precision, dimension(:,:,:), allocatable v
Definition: atimbr.f90:114
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:3215
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:1675
void cs_sla_amxby(double alpha, const cs_sla_matrix_t *m, const double x[], double beta, const double y[], double *inout[])
Compute an array resulting from alpha * M(x) + beta * y If inout is not allocated, allocation is done inside this function.
Definition: cs_sla_matrix.c:2811
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
cs_sla_matrix_t * cs_sla_matrix_create_from_ref(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. idx, didx and col_id are shared with ref...
Definition: cs_sla_matrix.c:1332
void cs_sla_matrix_rmzeros(cs_sla_matrix_t *m)
Remove entries with zero values Only available for CSR and MSR matrices with stride = 1...
Definition: cs_sla_matrix.c:1726
int n_cols
Definition: cs_sla.h:94
Definition: cs_sla.h:73
int n_rows
Definition: cs_sla.h:93
Definition: cs_sla.h:67
void cs_sla_matrix_msr2csr(cs_sla_matrix_t *a)
Change matrix representation from MSR to CSR.
Definition: cs_sla_matrix.c:2132
cs_sla_matrix_type_t type
Definition: cs_sla.h:88
Definition: cs_sla.h:69
double fillin
Definition: cs_sla.h:82
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:2035
void cs_sla_assemble_msr_sym(const cs_locmat_t *loc, cs_sla_matrix_t *ass, bool only_diag)
Assemble a MSR matrix from local contributions –> We assume that the local matrices are symmetric –...
Definition: cs_sla_matrix.c:2308
#define END_C_DECLS
Definition: cs_defs.h:430
cs_sla_matrix_type_t
Definition: cs_sla.h:63
void cs_sla_matrix_share2own(cs_sla_matrix_t *a)
Allocate its own pattern if shared.
Definition: cs_sla_matrix.c:2267
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:66
int flag
Definition: cs_sla.h:90
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:1938
Definition: cs_field_pointer.h:70
Definition: cs_sla.h:68
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:1871
cs_lnum_t * didx
Definition: cs_sla.h:102
cs_lnum_t * col_id
Definition: cs_sla.h:97
int stencil_max
Definition: cs_sla.h:77
void cs_sla_matrix_diag_idx(cs_sla_matrix_t *m)
Build diagonal index.
Definition: cs_sla_matrix.c:1897
Definition: cs_cdo_toolbox.h:72
cs_sla_matrix_t * cs_sla_matrix_create(cs_lnum_t n_rows, cs_lnum_t n_cols, int stride, cs_sla_matrix_type_t type, bool sym)
Create a cs_sla_matrix_t structure.
Definition: cs_sla_matrix.c:1257
double precision, save b
Definition: cs_fuel_incl.f90:146