programmer's documentation
cs_matrix_priv.h
Go to the documentation of this file.
1 #ifndef __CS_MATRIX_PRIV_H__
2 #define __CS_MATRIX_PRIV_H__
3 
4 /*============================================================================
5  * Private types for sparse matrix representation and operations.
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_defs.h"
35 
36 #include "cs_matrix.h"
37 
38 /*----------------------------------------------------------------------------*/
39 
41 
44 /*=============================================================================
45  * Macro definitions
46  *============================================================================*/
47 
48 /*============================================================================
49  * Type definitions
50  *============================================================================*/
51 
52 /* Formats currently supported:
53  *
54  * - Native
55  * - Compressed Sparse Row (CSR)
56  * - Symmetric Compressed Sparse Row (CSR_SYM)
57  */
58 
59 /*----------------------------------------------------------------------------
60  * Function pointer types
61  *----------------------------------------------------------------------------*/
62 
63 typedef void
64 (cs_matrix_set_coeffs_t) (cs_matrix_t *matrix,
65  bool symmetric,
66  bool copy,
67  cs_lnum_t n_edges,
68  const cs_lnum_2_t *restrict edges,
69  const cs_real_t *restrict da,
70  const cs_real_t *restrict xa);
71 
72 typedef void
73 (cs_matrix_release_coeffs_t) (cs_matrix_t *matrix);
74 
75 typedef void
76 (cs_matrix_copy_diagonal_t) (const cs_matrix_t *matrix,
77  cs_real_t *restrict da);
78 
79 typedef void
80 (cs_matrix_vector_product_t) (bool exclude_diag,
81  const cs_matrix_t *matrix,
82  const cs_real_t *restrict x,
83  cs_real_t *restrict y);
84 
85 /*----------------------------------------------------------------------------
86  * Matrix types
87  *----------------------------------------------------------------------------*/
88 
89 /* Native matrix structure representation */
90 /*----------------------------------------*/
91 
92 /* Note: the members of this structure are already available through the top
93  * matrix structure, but are replicated here in case of future removal
94  * from the top structure (which would require computation/assignment of
95  * matrix coefficients in another form) */
96 
97 typedef struct _cs_matrix_struct_native_t {
98 
99  cs_lnum_t n_rows; /* Local number of rows */
100  cs_lnum_t n_cols_ext; /* Local number of columns + ghosts */
101  cs_lnum_t n_edges; /* Local number of graph edges
102  (for extra-diagonal terms) */
103 
104  /* Pointers to shared arrays */
105 
106  const cs_lnum_2_t *edges; /* Edges (symmetric row <-> column)
107  connectivity */
108 
109 } cs_matrix_struct_native_t;
110 
111 /* Native matrix coefficients */
112 /*----------------------------*/
113 
114 typedef struct _cs_matrix_coeff_native_t {
115 
116  bool symmetric; /* Symmetry indicator */
117  int max_db_size; /* Current max allocated diag block size */
118  int max_eb_size; /* Current max allocated extradiag block size */
119 
120  /* Pointers to possibly shared arrays */
121 
122  const cs_real_t *da; /* Diagonal terms */
123  const cs_real_t *xa; /* Extra-diagonal terms */
124 
125  /* Pointers to private arrays (NULL if shared) */
126 
127  cs_real_t *_da; /* Diagonal terms */
128  cs_real_t *_xa; /* Extra-diagonal terms */
129 
130 } cs_matrix_coeff_native_t;
131 
132 /* CSR (Compressed Sparse Row) matrix structure representation */
133 /*-------------------------------------------------------------*/
134 
135 typedef struct _cs_matrix_struct_csr_t {
136 
137  cs_lnum_t n_rows; /* Local number of rows */
138  cs_lnum_t n_cols_ext; /* Local number of columns + ghosts */
139  cs_lnum_t n_cols_max; /* Maximum number of nonzero values
140  on a given row */
141 
142  /* Pointers to structure arrays and info (row_index, col_id) */
143 
144  bool have_diag; /* Has non-zero diagonal */
145  bool direct_assembly; /* True if each value corresponds to
146  a unique face ; false if multiple
147  faces contribute to the same
148  value (i.e. we have split faces) */
149 
150  const cs_lnum_t *row_index; /* Pointer to row index (0 to n-1) */
151  const cs_lnum_t *col_id; /* Pointer to column id (0 to n-1) */
152 
153  cs_lnum_t *_row_index; /* Row index (0 to n-1), if owner */
154  cs_lnum_t *_col_id; /* Column id (0 to n-1), if owner */
155 
156 } cs_matrix_struct_csr_t;
157 
158 /* CSR matrix coefficients representation */
159 /*----------------------------------------*/
160 
161 typedef struct _cs_matrix_coeff_csr_t {
162 
163  int n_prefetch_rows; /* Number of rows at a time for which
164  the x values in y = Ax should be
165  prefetched (0 if no prefetch) */
166 
167  /* Pointers to possibly shared arrays */
168 
169  const cs_real_t *val; /* Matrix coefficients */
170 
171  /* Pointers to private arrays (NULL if shared) */
172 
173  cs_real_t *_val; /* Diagonal matrix coefficients */
174 
175  cs_real_t *x_prefetch; /* Prefetch array for x in y = Ax */
176 
177  /* Pointers to auxiliary arrays used for queries */
178 
179  const cs_real_t *d_val; /* Pointer to diagonal matrix
180  coefficients, if queried */
181  cs_real_t *_d_val; /* Diagonal matrix coefficients,
182  if queried */
183 
184 } cs_matrix_coeff_csr_t;
185 
186 /* CSR_SYM (Symmetric Compressed Sparse Row) matrix structure representation */
187 /*---------------------------------------------------------------------------*/
188 
189 typedef struct _cs_matrix_struct_csr_sym_t {
190 
191  cs_lnum_t n_rows; /* Local number of rows */
192  cs_lnum_t n_cols; /* Local number of columns
193  (> n_rows in case of ghost columns) */
194  cs_lnum_t n_cols_max; /* Maximum number of nonzero values
195  on a given row */
196 
197  /* Pointers to structure arrays and info (row_index, col_id) */
198 
199  bool have_diag; /* Has non-zero diagonal */
200  bool direct_assembly; /* True if each value corresponds to
201  a unique face ; false if multiple
202  faces contribute to the same
203  value (i.e. we have split faces) */
204 
205  cs_lnum_t *row_index; /* Row index (0 to n-1) */
206  cs_lnum_t *col_id; /* Column id (0 to n-1) */
207 
208 } cs_matrix_struct_csr_sym_t;
209 
210 /* symmetric CSR matrix coefficients representation */
211 /*--------------------------------------------------*/
212 
213 typedef struct _cs_matrix_coeff_csr_sym_t {
214 
215  cs_real_t *val; /* Matrix coefficients */
216 
217  /* Pointers to auxiliary arrays used for queries */
218 
219  const cs_real_t *d_val; /* Pointer to diagonal matrix
220  coefficients, if queried */
221  cs_real_t *_d_val; /* Diagonal matrix coefficients,
222  if queried */
223 
224 } cs_matrix_coeff_csr_sym_t;
225 
226 /* MSR matrix coefficients representation */
227 /*----------------------------------------*/
228 
229 typedef struct _cs_matrix_coeff_msr_t {
230 
231  int n_prefetch_rows; /* Number of rows at a time for which
232  the x values in y = Ax should be
233  prefetched (0 if no prefetch) */
234  int max_db_size; /* Current max allocated block size */
235  int max_eb_size; /* Current max allocated extradiag block size */
236 
237  /* Pointers to possibly shared arrays */
238 
239  const cs_real_t *d_val; /* Diagonal matrix coefficients */
240  const cs_real_t *x_val; /* Extra-diagonal matrix coefficients */
241 
242  /* Pointers to private arrays (NULL if shared) */
243 
244  cs_real_t *_d_val; /* Diagonal matrix coefficients */
245  cs_real_t *_x_val; /* Extra-diagonal matrix coefficients */
246 
247  cs_real_t *x_prefetch; /* Prefetch array for x in y = Ax */
248 
249 } cs_matrix_coeff_msr_t;
250 
251 /* Matrix structure (representation-independent part) */
252 /*----------------------------------------------------*/
253 
254 struct _cs_matrix_structure_t {
255 
256  cs_matrix_type_t type; /* Matrix storage and definition type */
257 
258  cs_lnum_t n_rows; /* Local number of rows */
259  cs_lnum_t n_cols_ext; /* Local number of columns + ghosts */
260 
261  void *structure; /* Matrix structure */
262 
263  /* Pointers to shared arrays from mesh structure
264  (face->cell connectivity for coefficient assignment,
265  local->local cell numbering for future info or renumbering,
266  and halo) */
267 
268  const cs_gnum_t *row_num; /* Global row numbers */
269  const cs_halo_t *halo; /* Parallel or periodic halo */
270  const cs_numbering_t *numbering; /* Vectorization or thread-related
271  numbering information */
272 };
273 
274 /* Structure associated with Matrix (representation-independent part) */
275 /*--------------------------------------------------------------------*/
276 
277 struct _cs_matrix_t {
278 
279  cs_matrix_type_t type; /* Matrix storage and definition type */
280 
281  cs_lnum_t n_rows; /* Local number of rows */
282  cs_lnum_t n_cols_ext; /* Local number of columns + ghosts */
283 
284  cs_matrix_fill_type_t fill_type; /* Matrix fill type */
285 
286  int db_size[4]; /* Diag Block size, including padding:
287  0: useful block size
288  1: vector block extents
289  2: matrix line extents
290  3: matrix line*column extents */
291 
292  int eb_size[4]; /* Extradiag block size, including padding:
293  0: useful block size
294  1: vector block extents
295  2: matrix line extents
296  3: matrix line*column extents */
297 
298  /* Pointer to shared structure */
299 
300  const void *structure; /* Matrix structure */
301 
302  /* Pointers to arrays possibly shared from mesh structure
303  (graph edges: face->cell connectivity for coefficient assignment,
304  rows: local->local cell numbering for future info or renumbering,
305  and halo) */
306 
307  const cs_halo_t *halo; /* Parallel or periodic halo */
308  const cs_numbering_t *numbering; /* Vectorization or thread-related
309  numbering information */
310 
311  /* Pointer to shared arrays from coefficient assignment from
312  "native" type. This should be removed in the future, but requires
313  removing the dependency to the native structure in the multigrid
314  code first. */
315 
316  const cs_real_t *xa; /* Extra-diagonal terms */
317 
318  /* Pointer to private data */
319 
320  void *coeffs; /* Matrix coefficients */
321 
322  /* Function pointers */
323 
324  cs_matrix_set_coeffs_t *set_coefficients;
325  cs_matrix_release_coeffs_t *release_coefficients;
326  cs_matrix_copy_diagonal_t *copy_diagonal;
327 
328  /* Function pointer arrays, with CS_MATRIX_N_FILL_TYPES variants:
329  fill_type*2 + exclude_diagonal_flag */
330 
331  cs_matrix_vector_product_t *vector_multiply[CS_MATRIX_N_FILL_TYPES][2];
332 
333  /* Loop length parameter for some SpMv algorithms */
334 
335  int loop_length[CS_MATRIX_N_FILL_TYPES][2];
336 
337 };
338 
339 /* Structure used for tuning variants */
340 /*------------------------------------*/
341 
342 struct _cs_matrix_variant_t {
343 
344  char name[32]; /* Variant name */
345 
346  cs_matrix_type_t type; /* Matrix storage and definition type */
347 
348  /* Loop length parameter for some SpMv algorithms */
349 
350  int loop_length[CS_MATRIX_N_FILL_TYPES][2];
351 
352  /* Function pointer arrays, with variants:
353  fill_type + exclude_diagonal_flag */
354 
355  cs_matrix_vector_product_t *vector_multiply[CS_MATRIX_N_FILL_TYPES][2];
356 
357  /* Measured structure creation cost, or -1 otherwise */
358 
359  double matrix_create_cost;
360 
361  /* Measured assignment costs for each available fill type, or -1 otherwise */
362 
363  double matrix_assign_cost[CS_MATRIX_N_FILL_TYPES];
364 
365  /* Measured operation costs for each available operation, or -1 otherwise
366  fill_type*2 + exclude_diagonal_flag */
367 
368  double matrix_vector_cost[CS_MATRIX_N_FILL_TYPES][2][2];
369 
370 };
371 
374 /*=============================================================================
375  * Semi-private function prototypes
376  *============================================================================*/
377 
378 /*----------------------------------------------------------------------------
379  * Create CSR matrix coefficients.
380  *
381  * returns:
382  * pointer to allocated CSR coefficients structure.
383  *----------------------------------------------------------------------------*/
384 
385 cs_matrix_coeff_csr_t *
387 
388 /*----------------------------------------------------------------------------
389  * Destroy CSR matrix coefficients.
390  *
391  * parameters:
392  * coeff <-> Pointer to CSR matrix coefficients pointer
393  *----------------------------------------------------------------------------*/
394 
395 void
396 cs_matrix_destroy_coeff_csr(cs_matrix_coeff_csr_t **coeff);
397 
398 /*----------------------------------------------------------------------------
399  * Release shared CSR matrix coefficients.
400  *
401  * parameters:
402  * matrix <-- Pointer to matrix structure
403  *----------------------------------------------------------------------------*/
404 
405 void
407 
408 /*----------------------------------------------------------------------------
409  * Copy diagonal of CSR matrix.
410  *
411  * parameters:
412  * matrix <-- Pointer to matrix structure
413  * da --> Diagonal (pre-allocated, size: n_rows)
414  *----------------------------------------------------------------------------*/
415 
416 void
418  cs_real_t *restrict da);
419 
420 /*----------------------------------------------------------------------------
421  * Destroy CSR matrix structure.
422  *
423  * parameters:
424  * matrix <-> Pointer to CSR matrix structure pointer
425  *----------------------------------------------------------------------------*/
426 
427 void
428 cs_matrix_destroy_struct_csr(cs_matrix_struct_csr_t **matrix);
429 
430 /*----------------------------------------------------------------------------
431  * Local matrix.vector product y = A.x with CSR matrix.
432  *
433  * parameters:
434  * exclude_diag <-- exclude diagonal if true
435  * matrix <-- Pointer to matrix structure
436  * x <-- Multipliying vector values
437  * y --> Resulting vector
438  *----------------------------------------------------------------------------*/
439 
440 void
441 cs_matrix_vec_p_l_csr(bool exclude_diag,
442  const cs_matrix_t *matrix,
443  const cs_real_t *restrict x,
444  cs_real_t *restrict y);
445 
446 #if defined (HAVE_MKL)
447 /*----------------------------------------------------------------------------
448  * Local matrix.vector product y = A.x with MSR matrix, using MKL
449  *
450  * parameters:
451  * exclude_diag <-- exclude diagonal if true
452  * matrix <-- Pointer to matrix structure
453  * x <-- Multipliying vector values
454  * y --> Resulting vector
455  *----------------------------------------------------------------------------*/
456 
457 void
458 cs_matrix_vec_p_l_csr_mkl(bool exclude_diag,
459  const cs_matrix_t *matrix,
460  const cs_real_t *restrict x,
461  cs_real_t *restrict y);
462 
463 #endif /* defined (HAVE_MKL) */
464 
465 /*----------------------------------------------------------------------------
466  * Copy diagonal of native or MSR matrix.
467  *
468  * parameters:
469  * matrix <-- Pointer to matrix structure
470  * da --> Diagonal (pre-allocated, size: n_cols)
471  *----------------------------------------------------------------------------*/
472 
473 void
475  cs_real_t *restrict da);
476 
477 /*----------------------------------------------------------------------------
478  * Create MSR matrix coefficients.
479  *
480  * returns:
481  * pointer to allocated MSR coefficients structure.
482  *----------------------------------------------------------------------------*/
483 
484 cs_matrix_coeff_msr_t *
486 
487 /*----------------------------------------------------------------------------
488  * Release shared MSR matrix coefficients.
489  *
490  * parameters:
491  * matrix <-- Pointer to matrix structure
492  *----------------------------------------------------------------------------*/
493 
494 void
496 
497 /*----------------------------------------------------------------------------
498  * Local matrix.vector product y = A.x with MSR matrix.
499  *
500  * parameters:
501  * exclude_diag <-- exclude diagonal if true
502  * matrix <-- Pointer to matrix structure
503  * x <-- Multipliying vector values
504  * y --> Resulting vector
505  *----------------------------------------------------------------------------*/
506 
507 void
508 cs_matrix_vec_p_l_msr(bool exclude_diag,
509  const cs_matrix_t *matrix,
510  const cs_real_t *restrict x,
511  cs_real_t *restrict y);
512 
513 #if defined (HAVE_MKL)
514 /*----------------------------------------------------------------------------
515  * Local matrix.vector product y = A.x with MSR matrix, using MKL
516  *
517  * parameters:
518  * exclude_diag <-- exclude diagonal if true
519  * matrix <-- Pointer to matrix structure
520  * x <-- Multipliying vector values
521  * y --> Resulting vector
522  *----------------------------------------------------------------------------*/
523 
524 void
525 cs_matrix_vec_p_l_msr_mkl(bool exclude_diag,
526  const cs_matrix_t *matrix,
527  const cs_real_t *restrict x,
528  cs_real_t *restrict y);
529 #endif /* defined (HAVE_MKL) */
530 
531 /*----------------------------------------------------------------------------*/
532 
534 
535 #endif /* __CS_MATRIX_PRIV_H__ */
unsigned long cs_gnum_t
global mesh entity number
Definition: cs_defs.h:280
#define restrict
Definition: cs_defs.h:122
void cs_matrix_copy_diagonal_separate(const cs_matrix_t *matrix, cs_real_t *restrict da)
cs_matrix_coeff_csr_t * cs_matrix_create_coeff_csr(void)
#define BEGIN_C_DECLS
Definition: cs_defs.h:419
Definition: cs_matrix.h:76
Definition: cs_halo.h:70
void cs_matrix_copy_diagonal_csr(const cs_matrix_t *matrix, cs_real_t *restrict da)
void cs_matrix_release_coeffs_csr(cs_matrix_t *matrix)
void cs_matrix_release_coeffs_msr(cs_matrix_t *matrix)
struct _cs_matrix_t cs_matrix_t
Definition: cs_matrix.h:86
void cs_matrix_destroy_coeff_csr(cs_matrix_coeff_csr_t **coeff)
cs_matrix_coeff_msr_t * cs_matrix_create_coeff_msr(void)
void cs_matrix_vec_p_l_msr(bool exclude_diag, const cs_matrix_t *matrix, const cs_real_t *restrict x, cs_real_t *restrict y)
void cs_matrix_vec_p_l_csr(bool exclude_diag, const cs_matrix_t *matrix, const cs_real_t *restrict x, cs_real_t *restrict y)
cs_matrix_type_t
Definition: cs_matrix.h:54
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 cs_lnum_2_t[2]
vector of 2 local mesh-entity ids
Definition: cs_defs.h:301
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:292
#define END_C_DECLS
Definition: cs_defs.h:420
double cs_real_t
Definition: cs_defs.h:296
Definition: cs_numbering.h:78
void cs_matrix_destroy_struct_csr(cs_matrix_struct_csr_t **matrix)
cs_matrix_fill_type_t
Definition: cs_matrix.h:66