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-2016 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
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 /* Specific matrix for hybrid discretization X+C (C for cells) and X=V or F
110 
111  This matrix corresponds to 4 blocks :
112  (0,0) --> square MSR matrix of size #X
113  (1,0) --> rectangular matrix based on a cs_connect_index_t (c2x)
114  (1,1) --> diagonal matrix of size number of cells
115  (0,1) --> rectangular matrix sharing the same pattern as block (1,0)
116  Transposed of (1,0) is the matrix is symmetric
117 
118  The resulting matrix is square matrix of size #C + #X */
119 
120 typedef struct {
121 
122  int flag; // Symmetric, sorted, shared...
123 
124  cs_lnum_t n_x; // n_vertices or n_faces
126  cs_lnum_t n_rows; // n_x + n_cells
127 
128  const cs_connect_index_t *c2x; /* shared with a cs_cdo_connect_t structure
129  which is owner. Enable an easy acces to
130  blocks (1,0) and (0,1) */
131 
132  cs_sla_matrix_t *xx_block; // (0,0) block = MSR-type matrix
133  double *cc_diag; // (1,1) block = diagonal matrix of size n_cells
134  double *cx_vals; // (1,0) block = rectangular based on c2x index
135  double *xc_vals; /* (0,1) block = rectangular based on c2x index
136  Only allocated if matrix is not symmetric */
137 
139 
140 /*============================================================================
141  * Public function prototypes for SLA matrices
142  *============================================================================*/
143 
144 /*----------------------------------------------------------------------------*/
156 /*----------------------------------------------------------------------------*/
157 
160  cs_lnum_t n_cols,
161  int stride,
163  bool sym);
164 
165 /*----------------------------------------------------------------------------*/
177 /*----------------------------------------------------------------------------*/
178 
182  int stride);
183 
184 /*----------------------------------------------------------------------------*/
198 /*----------------------------------------------------------------------------*/
199 
202  bool is_symmetric,
203  bool sorted_idx,
204  int stride);
205 
206 /*----------------------------------------------------------------------------*/
215 /*----------------------------------------------------------------------------*/
216 
219  bool shared);
220 
221 /*----------------------------------------------------------------------------*/
229 /*----------------------------------------------------------------------------*/
230 
233 
234 /*----------------------------------------------------------------------------*/
242 /*----------------------------------------------------------------------------*/
243 
246 
247 /*----------------------------------------------------------------------------*/
256 /*----------------------------------------------------------------------------*/
257 
258 void
260  double threshold,
261  int verbosity);
262 
263 /*----------------------------------------------------------------------------*/
274 /*----------------------------------------------------------------------------*/
275 
276 void
277 cs_sla_matrix_clean(int verbosity,
278  double threshold,
279  cs_sla_matrix_t *m);
280 
281 /*----------------------------------------------------------------------------*/
289 /*----------------------------------------------------------------------------*/
290 
291 size_t
293 
294 /*----------------------------------------------------------------------------*/
300 /*----------------------------------------------------------------------------*/
301 
302 void
304 
305 /*----------------------------------------------------------------------------*/
312 /*----------------------------------------------------------------------------*/
313 
314 void
316  double *p_diag[]);
317 
318 /*----------------------------------------------------------------------------*/
324 /*----------------------------------------------------------------------------*/
325 
326 void
328 
329 /*----------------------------------------------------------------------------*/
336 /*----------------------------------------------------------------------------*/
337 
338 void
340 
341 /*----------------------------------------------------------------------------*/
347 /*----------------------------------------------------------------------------*/
348 
349 void
351 
352 /*----------------------------------------------------------------------------*/
358 /*----------------------------------------------------------------------------*/
359 
360 void
362 
363 /*----------------------------------------------------------------------------*/
369 /*----------------------------------------------------------------------------*/
370 
371 void
373 
374 /*----------------------------------------------------------------------------*/
384 /*----------------------------------------------------------------------------*/
385 
386 void
388  cs_sla_matrix_t *ass,
389  bool only_diag);
390 
391 /*----------------------------------------------------------------------------*/
399 /*----------------------------------------------------------------------------*/
400 
401 void
403  cs_sla_matrix_t *ass);
404 
405 /*----------------------------------------------------------------------------*/
420 /*----------------------------------------------------------------------------*/
421 
423 cs_sla_matrix_pack(cs_lnum_t n_final_rows,
424  cs_lnum_t n_final_cols,
425  const cs_sla_matrix_t *init,
426  const cs_lnum_t *row_z2i_ids,
427  const cs_lnum_t *col_i2z_ids,
428  bool keep_sym);
429 
430 /*----------------------------------------------------------------------------*/
441 /*----------------------------------------------------------------------------*/
442 
444 cs_sla_matrix_add(double alpha,
445  const cs_sla_matrix_t *a,
446  double beta,
447  const cs_sla_matrix_t *b);
448 
449 /*----------------------------------------------------------------------------*/
460 /*----------------------------------------------------------------------------*/
461 
462 void
464  const double v[],
465  double *inout[],
466  bool reset);
467 
468 /*----------------------------------------------------------------------------*/
480 /*----------------------------------------------------------------------------*/
481 
482 void
483 cs_sla_amxby(double alpha,
484  const cs_sla_matrix_t *m,
485  const double x[],
486  double beta,
487  const double y[],
488  double *inout[]);
489 
490 /*----------------------------------------------------------------------------*/
499 /*----------------------------------------------------------------------------*/
500 
503  const cs_sla_matrix_t *b);
504 
505 /*----------------------------------------------------------------------------*/
516 /*----------------------------------------------------------------------------*/
517 
520  const double D[],
521  const cs_sla_matrix_t *A,
522  cs_lnum_t *w);
523 
524 /*----------------------------------------------------------------------------*/
537 /*----------------------------------------------------------------------------*/
538 
541  const cs_sla_matrix_t *a,
542  double beta,
543  const cs_sla_matrix_t *bt,
544  const cs_sla_matrix_t *b);
545 
546 /*----------------------------------------------------------------------------*/
563 /*----------------------------------------------------------------------------*/
564 
565 void
567  const cs_sla_matrix_t *B,
568  const cs_sla_matrix_t *C,
569  const cs_sla_matrix_t *D,
570  const double X[],
571  const double Y[],
572  double *F[],
573  double *G[],
574  bool reset);
575 
576 /*----------------------------------------------------------------------------*/
586 /*----------------------------------------------------------------------------*/
587 
588 void
589 cs_sla_bwrite(const char *name,
590  const cs_sla_matrix_t *m,
591  const double *rhs,
592  const double *sol);
593 
594 /*----------------------------------------------------------------------------*/
604 /*----------------------------------------------------------------------------*/
605 
606 void
607 cs_sla_bread(const char *name,
608  cs_sla_matrix_t **p_mat,
609  double *p_rhs[],
610  double *p_sol[]);
611 
612 /*----------------------------------------------------------------------------*/
620 /*----------------------------------------------------------------------------*/
621 
622 void
623 cs_sla_matrix_summary(const char *name,
624  FILE *f,
625  cs_sla_matrix_t *m);
626 
627 /*----------------------------------------------------------------------------*/
635 /*----------------------------------------------------------------------------*/
636 
637 void
638 cs_sla_matrix_dump(const char *name,
639  FILE *f,
640  const cs_sla_matrix_t *m);
641 
642 /*----------------------------------------------------------------------------*/
651 /*----------------------------------------------------------------------------*/
652 
653 void
654 cs_sla_system_dump(const char *name,
655  FILE *f,
656  const cs_sla_matrix_t *m,
657  const double *rhs);
658 
659 /*----------------------------------------------------------------------------*/
673 /*----------------------------------------------------------------------------*/
674 
677  cs_lnum_t n_cells,
678  bool bktrans,
679  bool bk00sym,
680  const cs_connect_index_t *x2x,
681  const cs_connect_index_t *c2x);
682 
683 /*----------------------------------------------------------------------------*/
691 /*----------------------------------------------------------------------------*/
692 
695 
696 /*----------------------------------------------------------------------------*/
709 /*----------------------------------------------------------------------------*/
710 
711 void
713  const double vx[],
714  const double vc[],
715  double *iox[],
716  double *ioc[],
717  bool reset);
718 
719 /*----------------------------------------------------------------------------*/
729 /*----------------------------------------------------------------------------*/
730 
731 void
733  cs_sla_hmatrix_t *ass);
734 
735 /*----------------------------------------------------------------------------*/
736 
738 
739 #endif /* __CS_SLA_H__ */
int stencil_min
Definition: cs_sla.h:76
double * cc_diag
Definition: cs_sla.h:133
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.c:3555
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.c:1489
Definition: cs_sla.h:65
void cs_sla_hmatvec(const cs_sla_hmatrix_t *hm, const double vx[], const double vc[], double *iox[], double *ioc[], bool reset)
Compute a matrix vector product. If inout is not allocated, allocation is done inside this function I...
Definition: cs_sla.c:3771
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.c:3352
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.c:2902
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.c:3272
Definition: cs_cdo_toolbox.h:60
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.c:2965
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.c:2707
Definition: cs_sla.h:120
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.c:2823
void cs_sla_matrix_clean_zeros(cs_sla_matrix_t *m, double threshold, int verbosity)
Reset to 0 all entries below a given threshold Only available for CSR and MSR matrices with stride = ...
Definition: cs_sla.c:1732
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.c:2561
size_t nnz
Definition: cs_sla.h:81
double * xc_vals
Definition: cs_sla.h:135
#define BEGIN_C_DECLS
Definition: cs_defs.h:448
cs_sla_matrix_t * cs_sla_matrix_transpose(const cs_sla_matrix_t *mat)
Transpose a cs_sla_matrix_t structure.
Definition: cs_sla.c:1572
int stride
Definition: cs_sla.h:92
Definition: cs_field_pointer.h:82
void cs_sla_matrix_clean(int verbosity, double threshold, cs_sla_matrix_t *m)
Set to zero entries in a cs_sla_matrix_t structure if the ratio |a(i,j)| < eps * max|a(i,j)| is below a given threshold. Be careful when using this function since one can loose the symmetry Only available for matrices with stride = 1.
Definition: cs_sla.c:1774
cs_lnum_t n_x
Definition: cs_sla.h:124
double precision, dimension(ncharm), save beta
Definition: cpincl.f90:99
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.c:2023
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.c:2340
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.c:3456
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.c:1430
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.c:2155
cs_lnum_t n_cells
Definition: cs_sla.h:125
short int * sgn
Definition: cs_sla.h:99
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.c:3193
cs_sla_matrix_t * cs_sla_matrix_free(cs_sla_matrix_t *m)
Free a cs_sla_matrix_t structure.
Definition: cs_sla.c:1678
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.c:2789
const cs_connect_index_t * c2x
Definition: cs_sla.h:128
Definition: cs_convection_diffusion.h:47
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.c:1330
Definition: cs_convection_diffusion.h:47
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.c:2084
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.c:1987
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.c:2260
int flag
Definition: cs_sla.h:122
cs_lnum_t n_rows
Definition: cs_sla.h:126
#define END_C_DECLS
Definition: cs_defs.h:449
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.c:2219
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.c:1890
cs_sla_matrix_t * xx_block
Definition: cs_sla.h:132
cs_sla_hmatrix_t * cs_sla_hmatrix_free(cs_sla_hmatrix_t *hm)
Free a cs_sla_hmatrix_t structure.
Definition: cs_sla.c:3735
void cs_sla_assemble_hmat_sym(const cs_locmat_t *loc, cs_sla_hmatrix_t *ass)
Assemble a hybrid matrix from local contributions –> We assume that the local matrices are symmetric...
Definition: cs_sla.c:3868
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.c:1823
cs_lnum_t * didx
Definition: cs_sla.h:102
cs_lnum_t * col_id
Definition: cs_sla.h:97
double * cx_vals
Definition: cs_sla.h:134
cs_sla_hmatrix_t * cs_sla_hmatrix_create(cs_lnum_t n_x, cs_lnum_t n_cells, bool bktrans, bool bk00sym, const cs_connect_index_t *x2x, const cs_connect_index_t *c2x)
Create a cs_sla_hmatrix_t structure This is a square matrix of size n_x+n_cells (stride = 1 up to now...
Definition: cs_sla.c:3668
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.c:1849
Definition: cs_cdo_toolbox.h:71
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.c:1255
double precision, save b
Definition: cs_fuel_incl.f90:146