programmer's documentation
cs_sdm.h
Go to the documentation of this file.
1 #ifndef __CS_SDM_H__
2 #define __CS_SDM_H__
3 
4 /*============================================================================
5  * Set of operations to handle Small Dense Matrices (SDM)
6  *============================================================================*/
7 
8 /*
9  This file is part of Code_Saturne, a general-purpose CFD tool.
10 
11  Copyright (C) 1998-2018 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  * Standard C library headers
30  *----------------------------------------------------------------------------*/
31 
32 #include <stdio.h>
33 #include <string.h>
34 #include <math.h>
35 
36 /*----------------------------------------------------------------------------
37  * Local headers
38  *----------------------------------------------------------------------------*/
39 
40 /*----------------------------------------------------------------------------*/
41 
43 
44 /*============================================================================
45  * Macro definitions
46  *============================================================================*/
47 
48 #define CS_SDM_BY_BLOCK (1 << 0) /* Matrix is defined by block */
49 #define CS_SDM_SYMMETRIC (1 << 1) /* Matrix is symmetric by construction */
50 #define CS_SDM_SHARED_VAL (1 << 2) /* Matrix is not owner of its values */
51 
52 /*============================================================================
53  * Type definitions
54  *============================================================================*/
55 
56 typedef struct _cs_sdm_t cs_sdm_t;
57 
58 typedef struct {
59 
64 
65  /* Allocated to n_max_blocks_by_row*n_max_blocks_by_col
66  Cast in cs_sdm_t where values are shared with the master cs_sdm_t struct.
67  */
68  cs_sdm_t *blocks;
69 
71 
72 /* Structure enabling the repeated usage of Small Dense Matrices (SDM) during
73  the building of the linear system by a cellwise process */
74 struct _cs_sdm_t {
75 
76  cs_flag_t flag; /* Metadata */
77 
78  /* Row-related members */
79  int n_max_rows; // max number of entities by primal cells
80  int n_rows; // current number of entities
81 
82  /* Column-related members. Not useful if the matrix is square */
83  int n_max_cols; // max number of entries in a column
84  int n_cols; // current number of columns
85 
86  cs_real_t *val; // small dense matrix (size: n_max_rows*n_max_cols)
87 
88  /* Structure describing the matrix in terms of blocks */
90 
91 };
92 
93 /*============================================================================
94  * Prototypes for pointer of functions
95  *============================================================================*/
96 
97 /*----------------------------------------------------------------------------*/
106 /*----------------------------------------------------------------------------*/
107 
108 typedef void
109 (cs_sdm_product_t) (const cs_sdm_t *a,
110  const cs_sdm_t *b,
111  cs_sdm_t *c);
112 
113 /*----------------------------------------------------------------------------*/
122 /*----------------------------------------------------------------------------*/
123 
124 typedef void
125 (cs_sdm_matvec_t) (const cs_sdm_t *mat,
126  const cs_real_t *vec,
127  cs_real_t *mv);
128 
129 /*============================================================================
130  * Public function prototypes
131  *============================================================================*/
132 
133 /*----------------------------------------------------------------------------*/
145 /*----------------------------------------------------------------------------*/
146 
147 static inline cs_real_t
149  const cs_real_t x[],
150  const cs_real_t y[])
151 {
152  cs_real_t dp = 0;
153 
154  if (x == NULL || y == NULL)
155  return dp;
156 
157  for (int i = 0; i < n; i++)
158  dp += x[i]*y[i];
159 
160  return dp;
161 }
162 
163 /*----------------------------------------------------------------------------*/
174 /*----------------------------------------------------------------------------*/
175 
176 static inline void
178  const cs_real_t s,
179  const cs_real_t x[],
180  cs_real_t y[])
181 {
182  if (x == NULL || y == NULL)
183  return;
184 
185  for (int i = 0; i < n; i++)
186  y[i] = s * x[i];
187 }
188 
189 /*----------------------------------------------------------------------------*/
200 /*----------------------------------------------------------------------------*/
201 
202 static inline void
204  const cs_real_t s,
205  const cs_real_t x[],
206  cs_real_t y[])
207 {
208  if (x == NULL || y == NULL)
209  return;
210 
211  for (int i = 0; i < n; i++)
212  y[i] += s * x[i];
213 }
214 
215 /*----------------------------------------------------------------------------*/
226 /*----------------------------------------------------------------------------*/
227 
228 cs_sdm_t *
230  int n_max_rows,
231  int n_max_cols);
232 
233 /*----------------------------------------------------------------------------*/
242 /*----------------------------------------------------------------------------*/
243 
244 cs_sdm_t *
245 cs_sdm_square_create(int n_max_rows);
246 
247 /*----------------------------------------------------------------------------*/
256 /*----------------------------------------------------------------------------*/
257 
258 cs_sdm_t *
259 cs_sdm_create_copy(const cs_sdm_t *m);
260 
261 /*----------------------------------------------------------------------------*/
270 /*----------------------------------------------------------------------------*/
271 
272 cs_sdm_t *
273 cs_sdm_create_transpose(cs_sdm_t *mat);
274 
275 /*----------------------------------------------------------------------------*/
286 /*----------------------------------------------------------------------------*/
287 
288 cs_sdm_t *
289 cs_sdm_block_create(int n_max_blocks_by_row,
290  int n_max_blocks_by_col,
291  const short int max_row_block_sizes[],
292  const short int max_col_block_sizes[]);
293 
294 /*----------------------------------------------------------------------------*/
306 /*----------------------------------------------------------------------------*/
307 
308 static inline void
309 cs_sdm_map_array(int n_max_rows,
310  int n_max_cols,
311  cs_sdm_t *m,
312  cs_real_t *array)
313 {
314  assert(array != NULL && m != NULL); /* Sanity check */
315 
316  m->flag = CS_SDM_SHARED_VAL;
317  m->n_rows = m->n_max_rows = n_max_rows;
318  m->n_cols = m->n_max_cols = n_max_cols;
319  m->val = array;
320  m->block_desc = NULL;
321 }
322 
323 /*----------------------------------------------------------------------------*/
331 /*----------------------------------------------------------------------------*/
332 
333 cs_sdm_t *
334 cs_sdm_free(cs_sdm_t *mat);
335 
336 /*----------------------------------------------------------------------------*/
345 /*----------------------------------------------------------------------------*/
346 
347 static inline void
349  int n_cols,
350  cs_sdm_t *mat)
351 {
352  assert(mat != NULL);
353 
354  mat->n_rows = n_rows;
355  mat->n_cols = n_cols;
356  memset(mat->val, 0, n_rows*n_cols*sizeof(cs_real_t));
357 }
358 
359 /*----------------------------------------------------------------------------*/
367 /*----------------------------------------------------------------------------*/
368 
369 static inline void
371  cs_sdm_t *mat)
372 {
373  assert(mat != NULL);
374 
375  mat->n_rows = mat->n_cols = n_rows; /* square matrix */
376  memset(mat->val, 0, n_rows*n_rows*sizeof(cs_real_t));
377 }
378 
379 /*----------------------------------------------------------------------------*/
390 /*----------------------------------------------------------------------------*/
391 
392 void
393 cs_sdm_block_init(cs_sdm_t *m,
394  int n_row_blocks,
395  int n_col_blocks,
396  const short int row_block_sizes[],
397  const short int col_block_sizes[]);
398 
399 /*----------------------------------------------------------------------------*/
407 /*----------------------------------------------------------------------------*/
408 
409 static inline void
410 cs_sdm_copy(cs_sdm_t *recv,
411  const cs_sdm_t *send)
412 {
413  /* Sanity check */
414  assert(recv->n_max_rows >= send->n_rows);
415  assert(recv->n_max_cols >= send->n_cols);
416 
417  recv->flag = send->flag;
418  recv->n_rows = send->n_rows;
419  recv->n_cols = send->n_cols;
420 
421  /* Copy values */
422  memcpy(recv->val, send->val, sizeof(cs_real_t)*send->n_rows*send->n_cols);
423 }
424 
425 /*----------------------------------------------------------------------------*/
434 /*----------------------------------------------------------------------------*/
435 
436 cs_sdm_t *
437 cs_sdm_block_create_copy(const cs_sdm_t *mref);
438 
439 /*----------------------------------------------------------------------------*/
449 /*----------------------------------------------------------------------------*/
450 
451 static inline cs_sdm_t *
452 cs_sdm_get_block(const cs_sdm_t *m,
453  int row_block_id,
454  int col_block_id)
455 {
456  /* Sanity checks */
457  assert(m != NULL);
458  assert(m->flag & CS_SDM_BY_BLOCK && m->block_desc != NULL);
459  assert(col_block_id < m->block_desc->n_col_blocks);
460  assert(row_block_id < m->block_desc->n_row_blocks);
461 
462  const cs_sdm_block_t *bd = m->block_desc;
463 
464  return bd->blocks + row_block_id*bd->n_col_blocks + col_block_id;
465 }
466 
467 /*----------------------------------------------------------------------------*/
475 /*----------------------------------------------------------------------------*/
476 
477 static inline void
478 cs_sdm_get_col(const cs_sdm_t *m,
479  int col_id,
480  cs_real_t *col_vals)
481 {
482  /* Sanity checks */
483  assert(m != NULL && col_vals != NULL);
484  assert(col_id < m->n_cols);
485 
486  const cs_real_t *_col = m->val + col_id;
487  for(int i = 0; i < m->n_rows; i++, _col += m->n_cols)
488  col_vals[i] = *_col;
489 }
490 
491 /*----------------------------------------------------------------------------*/
504 /*----------------------------------------------------------------------------*/
505 
506 static inline void
507 cs_sdm_copy_block(const cs_sdm_t *m,
508  const short int r_id,
509  const short int c_id,
510  const short int nr,
511  const short int nc,
512  cs_sdm_t *b)
513 {
514  /* Sanity checks */
515  assert(m != NULL && b != NULL);
516  assert(r_id >= 0 && c_id >= 0);
517  assert((r_id + nr) <= m->n_rows);
518  assert((c_id + nc) <= m->n_cols);
519  assert(nr == b->n_rows && nc == b->n_cols);
520 
521  const cs_real_t *_start = m->val + c_id + r_id*m->n_cols;
522  for (short int i = 0; i < nr; i++, _start += m->n_cols)
523  memcpy(b->val + i*nc, _start, sizeof(cs_real_t)*nc);
524 }
525 
526 /*----------------------------------------------------------------------------*/
535 /*----------------------------------------------------------------------------*/
536 
537 static inline void
538 cs_sdm_transpose_and_update(const cs_sdm_t *m,
539  cs_sdm_t *mt)
540 {
541  assert(m != NULL && mt != NULL);
542  assert(m->n_rows == mt->n_cols && m->n_cols == mt->n_rows);
543 
544  for (short int i = 0; i < m->n_rows; i++) {
545  const cs_real_t *m_i = m->val + i*m->n_cols;
546  for (short int j = 0; j < m->n_cols; j++)
547  mt->val[j*mt->n_cols + i] += m_i[j];
548  }
549 }
550 
551 /*----------------------------------------------------------------------------*/
561 /*----------------------------------------------------------------------------*/
562 
563 void
564 cs_sdm_multiply(const cs_sdm_t *a,
565  const cs_sdm_t *b,
566  cs_sdm_t *c);
567 
568 /*----------------------------------------------------------------------------*/
580 /*----------------------------------------------------------------------------*/
581 
582 static inline void
583 cs_sdm_multiply_r1c3_rowrow(const cs_sdm_t *a,
584  const cs_sdm_t *b,
585  cs_sdm_t *c)
586 {
587  /* Sanity check */
588  assert(a != NULL && b != NULL && c != NULL);
589  assert(a->n_cols == 3 && b->n_cols == 3 &&
590  a->n_rows == 1 && c->n_rows == 1 &&
591  c->n_cols == 1 && b->n_rows == 1);
592 
593  c->val[0] += a->val[0]*b->val[0] + a->val[1]*b->val[1] + a->val[2]*b->val[2];
594 }
595 
596 /*----------------------------------------------------------------------------*/
608 /*----------------------------------------------------------------------------*/
609 
610 void
611 cs_sdm_multiply_rowrow(const cs_sdm_t *a,
612  const cs_sdm_t *b,
613  cs_sdm_t *c);
614 
615 /*----------------------------------------------------------------------------*/
628 /*----------------------------------------------------------------------------*/
629 
630 void
631 cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a,
632  const cs_sdm_t *b,
633  cs_sdm_t *c);
634 
635 /*----------------------------------------------------------------------------*/
647 /*----------------------------------------------------------------------------*/
648 
649 void
650 cs_sdm_block_multiply_rowrow(const cs_sdm_t *a,
651  const cs_sdm_t *b,
652  cs_sdm_t *c);
653 
654 /*----------------------------------------------------------------------------*/
667 /*----------------------------------------------------------------------------*/
668 
669 void
670 cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a,
671  const cs_sdm_t *b,
672  cs_sdm_t *c);
673 
674 /*----------------------------------------------------------------------------*/
683 /*----------------------------------------------------------------------------*/
684 
685 void
686 cs_sdm_square_matvec(const cs_sdm_t *mat,
687  const cs_real_t *vec,
688  cs_real_t *mv);
689 
690 /*----------------------------------------------------------------------------*/
699 /*----------------------------------------------------------------------------*/
700 
701 void
702 cs_sdm_matvec(const cs_sdm_t *mat,
703  const cs_real_t *vec,
704  cs_real_t *mv);
705 
706 /*----------------------------------------------------------------------------*/
716 /*----------------------------------------------------------------------------*/
717 
718 void
719 cs_sdm_update_matvec(const cs_sdm_t *mat,
720  const cs_real_t *vec,
721  cs_real_t *mv);
722 
723 /*----------------------------------------------------------------------------*/
734 /*----------------------------------------------------------------------------*/
735 
736 void
737 cs_sdm_matvec_transposed(const cs_sdm_t *mat,
738  const cs_real_t *vec,
739  cs_real_t *mv);
740 
741 /*----------------------------------------------------------------------------*/
748 /*----------------------------------------------------------------------------*/
749 
750 void
751 cs_sdm_block_add(cs_sdm_t *mat,
752  const cs_sdm_t *add);
753 
754 /*----------------------------------------------------------------------------*/
762 /*----------------------------------------------------------------------------*/
763 
764 void
765 cs_sdm_block_add_mult(cs_sdm_t *mat,
766  cs_real_t mult_coef,
767  const cs_sdm_t *add);
768 
769 /*----------------------------------------------------------------------------*/
779 /*----------------------------------------------------------------------------*/
780 
781 void
782 cs_sdm_block_matvec(const cs_sdm_t *mat,
783  const cs_real_t *vec,
784  cs_real_t *mv);
785 
786 /*----------------------------------------------------------------------------*/
793 /*----------------------------------------------------------------------------*/
794 
795 void
796 cs_sdm_add(cs_sdm_t *mat,
797  const cs_sdm_t *add);
798 
799 /*----------------------------------------------------------------------------*/
807 /*----------------------------------------------------------------------------*/
808 
809 void
810 cs_sdm_add_mult(cs_sdm_t *mat,
812  const cs_sdm_t *add);
813 
814 /*----------------------------------------------------------------------------*/
822 /*----------------------------------------------------------------------------*/
823 
824 void
825 cs_sdm_square_add_transpose(cs_sdm_t *mat,
826  cs_sdm_t *tr);
827 
828 /*----------------------------------------------------------------------------*/
834 /*----------------------------------------------------------------------------*/
835 
836 void
837 cs_sdm_square_asymm(cs_sdm_t *mat);
838 
839 /*----------------------------------------------------------------------------*/
845 /*----------------------------------------------------------------------------*/
846 
847 void
848 cs_sdm_block_square_asymm(cs_sdm_t *mat);
849 
850 /*----------------------------------------------------------------------------*/
867 /*----------------------------------------------------------------------------*/
868 
869 void
871  cs_real_t Qt[9],
872  cs_real_t R[6]);
873 
874 /*----------------------------------------------------------------------------*/
888 /*----------------------------------------------------------------------------*/
889 
890 void
891 cs_sdm_33_ldlt_compute(const cs_sdm_t *m,
892  cs_real_t facto[6]);
893 
894 /*----------------------------------------------------------------------------*/
908 /*----------------------------------------------------------------------------*/
909 
910 void
911 cs_sdm_44_ldlt_compute(const cs_sdm_t *m,
912  cs_real_t facto[10]);
913 
914 /*----------------------------------------------------------------------------*/
928 /*----------------------------------------------------------------------------*/
929 
930 void
931 cs_sdm_66_ldlt_compute(const cs_sdm_t *m,
932  cs_real_t facto[21]);
933 
934 /*----------------------------------------------------------------------------*/
949 /*----------------------------------------------------------------------------*/
950 
951 void
952 cs_sdm_ldlt_compute(const cs_sdm_t *m,
953  cs_real_t *facto,
954  cs_real_t *dkk);
955 
956 /*----------------------------------------------------------------------------*/
966 /*----------------------------------------------------------------------------*/
967 
968 void
969 cs_sdm_33_ldlt_solve(const cs_real_t facto[6],
970  const cs_real_t rhs[3],
971  cs_real_t sol[3]);
972 
973 /*----------------------------------------------------------------------------*/
983 /*----------------------------------------------------------------------------*/
984 
985 void
986 cs_sdm_44_ldlt_solve(const cs_real_t facto[10],
987  const cs_real_t rhs[4],
988  cs_real_t x[4]);
989 
990 /*----------------------------------------------------------------------------*/
1000 /*----------------------------------------------------------------------------*/
1001 
1002 void
1003 cs_sdm_66_ldlt_solve(const cs_real_t f[21],
1004  const cs_real_t b[6],
1005  cs_real_t x[6]);
1006 
1007 /*----------------------------------------------------------------------------*/
1019 /*----------------------------------------------------------------------------*/
1020 
1021 void
1023  const cs_real_t *facto,
1024  const cs_real_t *rhs,
1025  cs_real_t *sol);
1026 
1027 /*----------------------------------------------------------------------------*/
1037 /*----------------------------------------------------------------------------*/
1038 
1039 double
1040 cs_sdm_test_symmetry(const cs_sdm_t *mat);
1041 
1042 /*----------------------------------------------------------------------------*/
1048 /*----------------------------------------------------------------------------*/
1049 
1050 void
1051 cs_sdm_simple_dump(const cs_sdm_t *mat);
1052 
1053 /*----------------------------------------------------------------------------*/
1062 /*----------------------------------------------------------------------------*/
1063 
1064 void
1065 cs_sdm_dump(cs_lnum_t parent_id,
1066  const cs_lnum_t *row_ids,
1067  const cs_lnum_t *col_ids,
1068  const cs_sdm_t *mat);
1069 
1070 /*----------------------------------------------------------------------------*/
1083 /*----------------------------------------------------------------------------*/
1084 
1085 void
1086 cs_sdm_fprintf(FILE *fp,
1087  const char *fname,
1088  cs_real_t thd,
1089  const cs_sdm_t *m);
1090 
1091 /*----------------------------------------------------------------------------*/
1098 /*----------------------------------------------------------------------------*/
1099 
1100 void
1101 cs_sdm_block_dump(cs_lnum_t parent_id,
1102  const cs_sdm_t *mat);
1103 
1104 /*----------------------------------------------------------------------------*/
1117 /*----------------------------------------------------------------------------*/
1118 
1119 void
1120 cs_sdm_block_fprintf(FILE *fp,
1121  const char *fname,
1122  cs_real_t thd,
1123  const cs_sdm_t *m);
1124 
1125 /*----------------------------------------------------------------------------*/
1126 
1128 
1129 #endif /* __CS_SDM_H__ */
cs_flag_t flag
Definition: cs_sdm.h:76
void cs_sdm_33_ldlt_compute(const cs_sdm_t *m, cs_real_t facto[6])
LDL^T: Modified Cholesky decomposition of a 3x3 SPD matrix. For more reference, see for instance http...
Definition: cs_sdm.c:1230
cs_sdm_t * cs_sdm_free(cs_sdm_t *mat)
Free a cs_sdm_t structure.
Definition: cs_sdm.c:290
void cs_sdm_66_ldlt_solve(const cs_real_t f[21], const cs_real_t b[6], cs_real_t x[6])
Solve a 6x6 matrix with a modified Cholesky decomposition (L.D.L^T) The solution should be already al...
Definition: cs_sdm.c:1605
cs_sdm_t * cs_sdm_create_transpose(cs_sdm_t *mat)
Define a new matrix by adding the given matrix with its transpose. Keep the transposed matrix for a f...
Definition: cs_sdm.c:197
void cs_sdm_square_add_transpose(cs_sdm_t *mat, cs_sdm_t *tr)
Define a new matrix by adding the given matrix with its transpose. Keep the transposed matrix for a f...
Definition: cs_sdm.c:1019
void cs_sdm_ldlt_compute(const cs_sdm_t *m, cs_real_t *facto, cs_real_t *dkk)
LDL^T: Modified Cholesky decomposition of a SPD matrix. For more reference, see for instance http://m...
Definition: cs_sdm.c:1424
int n_max_cols
Definition: cs_sdm.h:83
cs_sdm_t * cs_sdm_create(cs_flag_t flag, int n_max_rows, int n_max_cols)
Allocate and initialize a cs_sdm_t structure Most generic function to create a cs_sdm_t structure...
Definition: cs_sdm.c:138
static void cs_sdm_init(int n_rows, int n_cols, cs_sdm_t *mat)
Initialize a cs_sdm_t structure Case of a square matrix.
Definition: cs_sdm.h:348
cs_sdm_t * cs_sdm_create_copy(const cs_sdm_t *m)
Allocate a cs_sdm_t structure and initialized it with the copy of the matrix m in input...
Definition: cs_sdm.c:174
void cs_sdm_block_square_asymm(cs_sdm_t *mat)
Set the given block matrix into its anti-symmetric part.
Definition: cs_sdm.c:1098
#define BEGIN_C_DECLS
Definition: cs_defs.h:462
Definition: cs_field_pointer.h:83
void cs_sdm_block_dump(cs_lnum_t parent_id, const cs_sdm_t *mat)
Dump a small dense matrix defined by blocks.
Definition: cs_sdm.c:1901
int n_col_blocks
Definition: cs_sdm.h:63
void cs_sdm_update_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix mv has been previously allocated and i...
Definition: cs_sdm.c:778
void cs_sdm_block_add_mult(cs_sdm_t *mat, cs_real_t mult_coef, const cs_sdm_t *add)
Add two matrices defined by block: loc += mult_coef * add.
Definition: cs_sdm.c:877
void cs_sdm_33_ldlt_solve(const cs_real_t facto[6], const cs_real_t rhs[3], cs_real_t sol[3])
Solve a 3x3 matrix with a modified Cholesky decomposition (L.D.L^T) The solution should be already al...
Definition: cs_sdm.c:1546
int n_max_blocks_by_row
Definition: cs_sdm.h:60
static void cs_sdm_copy(cs_sdm_t *recv, const cs_sdm_t *send)
Copy a cs_sdm_t structure into another cs_sdm_t structure which has been already allocated.
Definition: cs_sdm.h:410
void cs_sdm_add_mult(cs_sdm_t *mat, cs_real_t alpha, const cs_sdm_t *add)
Add two small dense matrices: loc += alpha*add.
Definition: cs_sdm.c:992
double cs_real_t
Floating-point value.
Definition: cs_defs.h:297
void cs_sdm_fprintf(FILE *fp, const char *fname, cs_real_t thd, const cs_sdm_t *m)
Print a cs_sdm_t structure not defined by block Print into the file f if given otherwise open a new f...
Definition: cs_sdm.c:1852
static void cs_sdm_map_array(int n_max_rows, int n_max_cols, cs_sdm_t *m, cs_real_t *array)
Map an array into a predefined cs_sdm_t structure. This array is shared and the lifecycle of this arr...
Definition: cs_sdm.h:309
#define CS_SDM_BY_BLOCK
Definition: cs_sdm.h:48
void cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T...
Definition: cs_sdm.c:643
static void cs_sdm_get_col(const cs_sdm_t *m, int col_id, cs_real_t *col_vals)
Get a copy of a column in a preallocated vector.
Definition: cs_sdm.h:478
#define CS_SDM_SHARED_VAL
Definition: cs_sdm.h:50
static cs_sdm_t * cs_sdm_get_block(const cs_sdm_t *m, int row_block_id, int col_block_id)
Get a specific block in a cs_sdm_t structure defined by block.
Definition: cs_sdm.h:452
int n_cols
Definition: cs_sdm.h:84
double precision, save a
Definition: cs_fuel_incl.f90:146
cs_sdm_t * cs_sdm_block_create_copy(const cs_sdm_t *mref)
Allocate and initialize a cs_sdm_t structure w.r.t. to a given matrix.
Definition: cs_sdm.c:378
void cs_sdm_44_ldlt_compute(const cs_sdm_t *m, cs_real_t facto[10])
LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http...
Definition: cs_sdm.c:1277
static cs_real_t cs_sdm_dot(int n, const cs_real_t x[], const cs_real_t y[])
Basic dot product for a small vector For very small array sizes (3, 4, 6) prefer functions in cs_math...
Definition: cs_sdm.h:148
void cs_sdm_block_fprintf(FILE *fp, const char *fname, cs_real_t thd, const cs_sdm_t *m)
Print a cs_sdm_t structure which is defined by block Print into the file f if given otherwise open a ...
Definition: cs_sdm.c:1971
cs_sdm_t * cs_sdm_square_create(int n_max_rows)
Allocate and initialize a cs_sdm_t structure Case of a square matrix.
Definition: cs_sdm.c:157
void cs_sdm_dump(cs_lnum_t parent_id, const cs_lnum_t *row_ids, const cs_lnum_t *col_ids, const cs_sdm_t *mat)
Dump a small dense matrix.
Definition: cs_sdm.c:1800
void cs_sdm_block_init(cs_sdm_t *m, int n_row_blocks, int n_col_blocks, const short int row_block_sizes[], const short int col_block_sizes[])
Initialize the pattern of cs_sdm_t structure defined by block The matrix should have been allocated b...
Definition: cs_sdm.c:322
int n_row_blocks
Definition: cs_sdm.h:61
void cs_sdm_66_ldlt_compute(const cs_sdm_t *m, cs_real_t facto[21])
LDL^T: Modified Cholesky decomposition of a 6x6 SPD matrix. For more reference, see for instance http...
Definition: cs_sdm.c:1334
cs_sdm_block_t * block_desc
Definition: cs_sdm.h:89
cs_sdm_t * cs_sdm_block_create(int n_max_blocks_by_row, int n_max_blocks_by_col, const short int max_row_block_sizes[], const short int max_col_block_sizes[])
Allocate and initialize a cs_sdm_t structure.
Definition: cs_sdm.c:232
int n_max_blocks_by_col
Definition: cs_sdm.h:62
void cs_sdm_ldlt_solve(int n_rows, const cs_real_t *facto, const cs_real_t *rhs, cs_real_t *sol)
Solve a SPD matrix with a L.D.L^T (Modified Cholesky decomposition) The solution should be already al...
Definition: cs_sdm.c:1642
static void cs_sdm_transpose_and_update(const cs_sdm_t *m, cs_sdm_t *mt)
transpose and copy a matrix into another one already shaped sub-matrix starting from (r_id...
Definition: cs_sdm.h:538
double cs_sdm_test_symmetry(const cs_sdm_t *mat)
Test if a matrix is symmetric. Return 0. if the extradiagonal differences are lower thann the machine...
Definition: cs_sdm.c:1715
Definition: cs_sdm.h:74
Definition: cs_sdm.h:58
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:293
void cs_sdm_block_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix defined by block mv has been previousl...
Definition: cs_sdm.c:917
void cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T...
Definition: cs_sdm.c:540
void cs_sdm_block_multiply_rowrow(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T...
Definition: cs_sdm.c:588
#define END_C_DECLS
Definition: cs_defs.h:463
void cs_sdm_multiply_rowrow(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T...
Definition: cs_sdm.c:494
unsigned short int cs_flag_t
Definition: cs_defs.h:299
void cs_sdm_33_sym_qr_compute(const cs_real_t m[9], cs_real_t Qt[9], cs_real_t R[6])
Decompose a matrix into the matrix product Q.R Case of a 3x3 symmetric matrix.
Definition: cs_sdm.c:1165
void cs_sdm_add(cs_sdm_t *mat, const cs_sdm_t *add)
Add two small dense matrices: loc += add.
Definition: cs_sdm.c:969
static void cs_sdm_scalvect(int n, const cs_real_t s, const cs_real_t x[], cs_real_t y[])
Multiply a small vector by a scalar coefficient: y = s x For very small array sizes (3...
Definition: cs_sdm.h:177
void cs_sdm_matvec_transposed(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix which is transposed. mv has been previously allocated. mv is updated inside this function. Don&#39;t forget to initialize mv if needed.
Definition: cs_sdm.c:811
void cs_sdm_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix mv has been previously allocated...
Definition: cs_sdm.c:736
void() cs_sdm_product_t(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Generic prototype for computing a local dense matrix-product c = a*b where c has been previously allo...
Definition: cs_sdm.h:109
double precision, dimension(:,:,:), allocatable nc
Definition: atimbr.f90:110
cs_real_t * val
Definition: cs_sdm.h:86
void() cs_sdm_matvec_t(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Generic prototype for computing a dense matrix-vector product mv has been previously allocated...
Definition: cs_sdm.h:125
static void cs_sdm_add_scalvect(int n, const cs_real_t s, const cs_real_t x[], cs_real_t y[])
Multiply a small vector by a scalar coefficient: y += s x For very small array sizes (3...
Definition: cs_sdm.h:203
static void cs_sdm_multiply_r1c3_rowrow(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T...
Definition: cs_sdm.h:583
cs_sdm_t * blocks
Definition: cs_sdm.h:68
void cs_sdm_square_asymm(cs_sdm_t *mat)
Set the given matrix into its anti-symmetric part.
Definition: cs_sdm.c:1062
static void cs_sdm_copy_block(const cs_sdm_t *m, const short int r_id, const short int c_id, const short int nr, const short int nc, cs_sdm_t *b)
copy a block of a matrix into a sub-matrix starting from (r_id, c_id) with a size of nr rows and nc c...
Definition: cs_sdm.h:507
void cs_sdm_44_ldlt_solve(const cs_real_t facto[10], const cs_real_t rhs[4], cs_real_t x[4])
Solve a 4x4 matrix with a modified Cholesky decomposition (L.D.L^T) The solution should be already al...
Definition: cs_sdm.c:1574
int n_max_rows
Definition: cs_sdm.h:79
void cs_sdm_multiply(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a local dense matrix-product c = a*b c has been previously allocated.
Definition: cs_sdm.c:451
void cs_sdm_square_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a small square matrix mv has been previously allocated...
Definition: cs_sdm.c:700
void cs_sdm_simple_dump(const cs_sdm_t *mat)
Dump a small dense matrix.
Definition: cs_sdm.c:1771
void cs_sdm_block_add(cs_sdm_t *mat, const cs_sdm_t *add)
Add two matrices defined by block: loc += add.
Definition: cs_sdm.c:841
int n_rows
Definition: cs_sdm.h:80
double precision, save b
Definition: cs_fuel_incl.f90:146
static void cs_sdm_square_init(int n_rows, cs_sdm_t *mat)
Initialize a cs_sdm_t structure Case of a square matrix.
Definition: cs_sdm.h:370