My Project
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-2019 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 /*----------------------------------------------------------------------------*/
269 /*----------------------------------------------------------------------------*/
270 
271 cs_sdm_t *
272 cs_sdm_create_transpose(cs_sdm_t *mat);
273 
274 /*----------------------------------------------------------------------------*/
285 /*----------------------------------------------------------------------------*/
286 
287 cs_sdm_t *
288 cs_sdm_block_create(int n_max_blocks_by_row,
289  int n_max_blocks_by_col,
290  const int max_row_block_sizes[],
291  const int max_col_block_sizes[]);
292 
293 /*----------------------------------------------------------------------------*/
303 /*----------------------------------------------------------------------------*/
304 
305 cs_sdm_t *
306 cs_sdm_block33_create(int n_max_blocks_by_row,
307  int n_max_blocks_by_col);
308 
309 /*----------------------------------------------------------------------------*/
321 /*----------------------------------------------------------------------------*/
322 
323 static inline void
324 cs_sdm_map_array(int n_max_rows,
325  int n_max_cols,
326  cs_sdm_t *m,
327  cs_real_t *array)
328 {
329  assert(array != NULL && m != NULL); /* Sanity check */
330 
331  m->flag = CS_SDM_SHARED_VAL;
332  m->n_rows = m->n_max_rows = n_max_rows;
333  m->n_cols = m->n_max_cols = n_max_cols;
334  m->val = array;
335  m->block_desc = NULL;
336 }
337 
338 /*----------------------------------------------------------------------------*/
346 /*----------------------------------------------------------------------------*/
347 
348 cs_sdm_t *
349 cs_sdm_free(cs_sdm_t *mat);
350 
351 /*----------------------------------------------------------------------------*/
360 /*----------------------------------------------------------------------------*/
361 
362 static inline void
363 cs_sdm_init(int n_rows,
364  int n_cols,
365  cs_sdm_t *mat)
366 {
367  assert(mat != NULL);
368 
369  mat->n_rows = n_rows;
370  mat->n_cols = n_cols;
371  memset(mat->val, 0, n_rows*n_cols*sizeof(cs_real_t));
372 }
373 
374 /*----------------------------------------------------------------------------*/
382 /*----------------------------------------------------------------------------*/
383 
384 static inline void
386  cs_sdm_t *mat)
387 {
388  assert(mat != NULL);
389 
390  mat->n_rows = mat->n_cols = n_rows; /* square matrix */
391  memset(mat->val, 0, n_rows*n_rows*sizeof(cs_real_t));
392 }
393 
394 /*----------------------------------------------------------------------------*/
405 /*----------------------------------------------------------------------------*/
406 
407 void
408 cs_sdm_block_init(cs_sdm_t *m,
409  int n_row_blocks,
410  int n_col_blocks,
411  const int row_block_sizes[],
412  const int col_block_sizes[]);
413 
414 /*----------------------------------------------------------------------------*/
423 /*----------------------------------------------------------------------------*/
424 
425 void
426 cs_sdm_block33_init(cs_sdm_t *m,
427  int n_row_blocks,
428  int n_col_blocks);
429 
430 /*----------------------------------------------------------------------------*/
438 /*----------------------------------------------------------------------------*/
439 
440 static inline void
441 cs_sdm_copy(cs_sdm_t *recv,
442  const cs_sdm_t *send)
443 {
444  /* Sanity check */
445  assert(recv->n_max_rows >= send->n_rows);
446  assert(recv->n_max_cols >= send->n_cols);
447 
448  recv->flag = send->flag;
449  recv->n_rows = send->n_rows;
450  recv->n_cols = send->n_cols;
451 
452  /* Copy values */
453  memcpy(recv->val, send->val, sizeof(cs_real_t)*send->n_rows*send->n_cols);
454 }
455 
456 /*----------------------------------------------------------------------------*/
465 /*----------------------------------------------------------------------------*/
466 
467 cs_sdm_t *
468 cs_sdm_block_create_copy(const cs_sdm_t *mref);
469 
470 /*----------------------------------------------------------------------------*/
480 /*----------------------------------------------------------------------------*/
481 
482 static inline cs_sdm_t *
483 cs_sdm_get_block(const cs_sdm_t *const m,
484  int row_block_id,
485  int col_block_id)
486 {
487  /* Sanity checks */
488  assert(m != NULL);
489  assert(m->flag & CS_SDM_BY_BLOCK && m->block_desc != NULL);
490  assert(col_block_id < m->block_desc->n_col_blocks);
491  assert(row_block_id < m->block_desc->n_row_blocks);
492 
493  return m->block_desc->blocks
494  + row_block_id*m->block_desc->n_col_blocks + col_block_id;
495 }
496 
497 /*----------------------------------------------------------------------------*/
505 /*----------------------------------------------------------------------------*/
506 
507 static inline void
508 cs_sdm_get_col(const cs_sdm_t *m,
509  int col_id,
510  cs_real_t *col_vals)
511 {
512  /* Sanity checks */
513  assert(m != NULL && col_vals != NULL);
514  assert(col_id < m->n_cols);
515 
516  const cs_real_t *_col = m->val + col_id;
517  for(int i = 0; i < m->n_rows; i++, _col += m->n_cols)
518  col_vals[i] = *_col;
519 }
520 
521 /*----------------------------------------------------------------------------*/
534 /*----------------------------------------------------------------------------*/
535 
536 static inline void
537 cs_sdm_copy_block(const cs_sdm_t *m,
538  const short int r_id,
539  const short int c_id,
540  const short int nr,
541  const short int nc,
542  cs_sdm_t *b)
543 {
544  /* Sanity checks */
545  assert(m != NULL && b != NULL);
546  assert(r_id >= 0 && c_id >= 0);
547  assert((r_id + nr) <= m->n_rows);
548  assert((c_id + nc) <= m->n_cols);
549  assert(nr == b->n_rows && nc == b->n_cols);
550 
551  const cs_real_t *_start = m->val + c_id + r_id*m->n_cols;
552  for (short int i = 0; i < nr; i++, _start += m->n_cols)
553  memcpy(b->val + i*nc, _start, sizeof(cs_real_t)*nc);
554 }
555 
556 /*----------------------------------------------------------------------------*/
565 /*----------------------------------------------------------------------------*/
566 
567 static inline void
568 cs_sdm_transpose_and_update(const cs_sdm_t *m,
569  cs_sdm_t *mt)
570 {
571  assert(m != NULL && mt != NULL);
572  assert(m->n_rows == mt->n_cols && m->n_cols == mt->n_rows);
573 
574  for (short int i = 0; i < m->n_rows; i++) {
575  const cs_real_t *m_i = m->val + i*m->n_cols;
576  for (short int j = 0; j < m->n_cols; j++)
577  mt->val[j*mt->n_cols + i] += m_i[j];
578  }
579 }
580 
581 /*----------------------------------------------------------------------------*/
591 /*----------------------------------------------------------------------------*/
592 
593 void
594 cs_sdm_multiply(const cs_sdm_t *a,
595  const cs_sdm_t *b,
596  cs_sdm_t *c);
597 
598 /*----------------------------------------------------------------------------*/
610 /*----------------------------------------------------------------------------*/
611 
612 static inline void
614  const cs_sdm_t *b,
615  cs_sdm_t *c)
616 {
617  /* Sanity check */
618  assert(a != NULL && b != NULL && c != NULL);
619  assert(a->n_cols == 3 && b->n_cols == 3 &&
620  a->n_rows == 1 && c->n_rows == 1 &&
621  c->n_cols == 1 && b->n_rows == 1);
622 
623  c->val[0] += a->val[0]*b->val[0] + a->val[1]*b->val[1] + a->val[2]*b->val[2];
624 }
625 
626 /*----------------------------------------------------------------------------*/
638 /*----------------------------------------------------------------------------*/
639 
640 void
641 cs_sdm_multiply_rowrow(const cs_sdm_t *a,
642  const cs_sdm_t *b,
643  cs_sdm_t *c);
644 
645 /*----------------------------------------------------------------------------*/
658 /*----------------------------------------------------------------------------*/
659 
660 void
661 cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a,
662  const cs_sdm_t *b,
663  cs_sdm_t *c);
664 
665 /*----------------------------------------------------------------------------*/
677 /*----------------------------------------------------------------------------*/
678 
679 void
680 cs_sdm_block_multiply_rowrow(const cs_sdm_t *a,
681  const cs_sdm_t *b,
682  cs_sdm_t *c);
683 
684 /*----------------------------------------------------------------------------*/
697 /*----------------------------------------------------------------------------*/
698 
699 void
700 cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a,
701  const cs_sdm_t *b,
702  cs_sdm_t *c);
703 
704 /*----------------------------------------------------------------------------*/
713 /*----------------------------------------------------------------------------*/
714 
715 void
716 cs_sdm_square_matvec(const cs_sdm_t *mat,
717  const cs_real_t *vec,
718  cs_real_t *mv);
719 
720 /*----------------------------------------------------------------------------*/
729 /*----------------------------------------------------------------------------*/
730 
731 void
732 cs_sdm_matvec(const cs_sdm_t *mat,
733  const cs_real_t *vec,
734  cs_real_t *mv);
735 
736 /*----------------------------------------------------------------------------*/
746 /*----------------------------------------------------------------------------*/
747 
748 void
749 cs_sdm_update_matvec(const cs_sdm_t *mat,
750  const cs_real_t *vec,
751  cs_real_t *mv);
752 
753 /*----------------------------------------------------------------------------*/
764 /*----------------------------------------------------------------------------*/
765 
766 void
767 cs_sdm_matvec_transposed(const cs_sdm_t *mat,
768  const cs_real_t *vec,
769  cs_real_t *mv);
770 
771 /*----------------------------------------------------------------------------*/
778 /*----------------------------------------------------------------------------*/
779 
780 void
781 cs_sdm_block_add(cs_sdm_t *mat,
782  const cs_sdm_t *add);
783 
784 /*----------------------------------------------------------------------------*/
792 /*----------------------------------------------------------------------------*/
793 
794 void
795 cs_sdm_block_add_mult(cs_sdm_t *mat,
796  cs_real_t mult_coef,
797  const cs_sdm_t *add);
798 
799 /*----------------------------------------------------------------------------*/
809 /*----------------------------------------------------------------------------*/
810 
811 void
812 cs_sdm_block_matvec(const cs_sdm_t *mat,
813  const cs_real_t *vec,
814  cs_real_t *mv);
815 
816 /*----------------------------------------------------------------------------*/
823 /*----------------------------------------------------------------------------*/
824 
825 void
826 cs_sdm_add(cs_sdm_t *mat,
827  const cs_sdm_t *add);
828 
829 /*----------------------------------------------------------------------------*/
837 /*----------------------------------------------------------------------------*/
838 
839 void
840 cs_sdm_add_mult(cs_sdm_t *mat,
842  const cs_sdm_t *add);
843 
844 /*----------------------------------------------------------------------------*/
852 /*----------------------------------------------------------------------------*/
853 
854 void
855 cs_sdm_square_add_transpose(cs_sdm_t *mat,
856  cs_sdm_t *tr);
857 
858 /*----------------------------------------------------------------------------*/
865 /*----------------------------------------------------------------------------*/
866 
867 void
868 cs_sdm_square_2symm(cs_sdm_t *mat);
869 
870 /*----------------------------------------------------------------------------*/
876 /*----------------------------------------------------------------------------*/
877 
878 void
879 cs_sdm_square_asymm(cs_sdm_t *mat);
880 
881 /*----------------------------------------------------------------------------*/
887 /*----------------------------------------------------------------------------*/
888 
889 void
890 cs_sdm_block_square_asymm(cs_sdm_t *mat);
891 
892 /*----------------------------------------------------------------------------*/
909 /*----------------------------------------------------------------------------*/
910 
911 void
913  cs_real_t Qt[9],
914  cs_real_t R[6]);
915 
916 /*----------------------------------------------------------------------------*/
927 /*----------------------------------------------------------------------------*/
928 
929 void
930 cs_sdm_33_lu_compute(const cs_sdm_t *m,
931  cs_real_t facto[9]);
932 
933 /*----------------------------------------------------------------------------*/
946 /*----------------------------------------------------------------------------*/
947 
948 void
949 cs_sdm_lu_compute(const cs_sdm_t *m,
950  cs_real_t facto[]);
951 
952 /*----------------------------------------------------------------------------*/
966 /*----------------------------------------------------------------------------*/
967 
968 void
969 cs_sdm_33_lu_solve(const cs_real_t facto[9],
970  const cs_real_t rhs[3],
971  cs_real_t sol[3]);
972 
973 /*----------------------------------------------------------------------------*/
988 /*----------------------------------------------------------------------------*/
989 
990 void
992  const cs_real_t facto[],
993  const cs_real_t *rhs,
994  cs_real_t *sol);
995 
996 /*----------------------------------------------------------------------------*/
1010 /*----------------------------------------------------------------------------*/
1011 
1012 void
1013 cs_sdm_33_ldlt_compute(const cs_sdm_t *m,
1014  cs_real_t facto[6]);
1015 
1016 /*----------------------------------------------------------------------------*/
1030 /*----------------------------------------------------------------------------*/
1031 
1032 void
1033 cs_sdm_44_ldlt_compute(const cs_sdm_t *m,
1034  cs_real_t facto[10]);
1035 
1036 /*----------------------------------------------------------------------------*/
1050 /*----------------------------------------------------------------------------*/
1051 
1052 void
1053 cs_sdm_66_ldlt_compute(const cs_sdm_t *m,
1054  cs_real_t facto[21]);
1055 
1056 /*----------------------------------------------------------------------------*/
1071 /*----------------------------------------------------------------------------*/
1072 
1073 void
1074 cs_sdm_ldlt_compute(const cs_sdm_t *m,
1075  cs_real_t *facto,
1076  cs_real_t *dkk);
1077 
1078 /*----------------------------------------------------------------------------*/
1088 /*----------------------------------------------------------------------------*/
1089 
1090 void
1091 cs_sdm_33_ldlt_solve(const cs_real_t facto[6],
1092  const cs_real_t rhs[3],
1093  cs_real_t sol[3]);
1094 
1095 /*----------------------------------------------------------------------------*/
1105 /*----------------------------------------------------------------------------*/
1106 
1107 void
1108 cs_sdm_44_ldlt_solve(const cs_real_t facto[10],
1109  const cs_real_t rhs[4],
1110  cs_real_t x[4]);
1111 
1112 /*----------------------------------------------------------------------------*/
1122 /*----------------------------------------------------------------------------*/
1123 
1124 void
1125 cs_sdm_66_ldlt_solve(const cs_real_t f[21],
1126  const cs_real_t b[6],
1127  cs_real_t x[6]);
1128 
1129 /*----------------------------------------------------------------------------*/
1141 /*----------------------------------------------------------------------------*/
1142 
1143 void
1144 cs_sdm_ldlt_solve(int n_rows,
1145  const cs_real_t *facto,
1146  const cs_real_t *rhs,
1147  cs_real_t *sol);
1148 
1149 /*----------------------------------------------------------------------------*/
1159 /*----------------------------------------------------------------------------*/
1160 
1161 double
1162 cs_sdm_test_symmetry(const cs_sdm_t *mat);
1163 
1164 /*----------------------------------------------------------------------------*/
1170 /*----------------------------------------------------------------------------*/
1171 
1172 void
1173 cs_sdm_simple_dump(const cs_sdm_t *mat);
1174 
1175 /*----------------------------------------------------------------------------*/
1184 /*----------------------------------------------------------------------------*/
1185 
1186 void
1187 cs_sdm_dump(cs_lnum_t parent_id,
1188  const cs_lnum_t *row_ids,
1189  const cs_lnum_t *col_ids,
1190  const cs_sdm_t *mat);
1191 
1192 /*----------------------------------------------------------------------------*/
1205 /*----------------------------------------------------------------------------*/
1206 
1207 void
1208 cs_sdm_fprintf(FILE *fp,
1209  const char *fname,
1210  cs_real_t thd,
1211  const cs_sdm_t *m);
1212 
1213 /*----------------------------------------------------------------------------*/
1220 /*----------------------------------------------------------------------------*/
1221 
1222 void
1223 cs_sdm_block_dump(cs_lnum_t parent_id,
1224  const cs_sdm_t *mat);
1225 
1226 /*----------------------------------------------------------------------------*/
1239 /*----------------------------------------------------------------------------*/
1240 
1241 void
1242 cs_sdm_block_fprintf(FILE *fp,
1243  const char *fname,
1244  cs_real_t thd,
1245  const cs_sdm_t *m);
1246 
1247 /*----------------------------------------------------------------------------*/
1248 
1250 
1251 #endif /* __CS_SDM_H__ */
cs_sdm_block_add
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:922
cs_sdm_create_copy
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
cs_sdm_block33_create
cs_sdm_t * cs_sdm_block33_create(int n_max_blocks_by_row, int n_max_blocks_by_col)
Allocate and initialize a cs_sdm_t structure by block when the block size is constant and equal to 3.
Definition: cs_sdm.c:291
cs_sdm_block_init
void cs_sdm_block_init(cs_sdm_t *m, int n_row_blocks, int n_col_blocks, const int row_block_sizes[], const 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:363
cs_sdm_block_t::n_max_blocks_by_row
int n_max_blocks_by_row
Definition: cs_sdm.h:60
cs_sdm_square_asymm
void cs_sdm_square_asymm(cs_sdm_t *mat)
Set the given matrix into its anti-symmetric part.
Definition: cs_sdm.c:1177
_cs_sdm_t::block_desc
cs_sdm_block_t * block_desc
Definition: cs_sdm.h:89
cs_sdm_create_transpose
cs_sdm_t * cs_sdm_create_transpose(cs_sdm_t *mat)
Define a new matrix which is its transpose.
Definition: cs_sdm.c:196
cs_sdm_square_matvec
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:781
cs_sdm_multiply
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:531
_cs_sdm_t::val
cs_real_t * val
Definition: cs_sdm.h:86
cs_sdm_block_dump
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:2194
cs_sdm_lu_compute
void cs_sdm_lu_compute(const cs_sdm_t *m, cs_real_t facto[])
LU factorization of a small dense matrix. Small means that the number m->n_rows is less than 100 for ...
Definition: cs_sdm.c:1380
cs_sdm_block_create
cs_sdm_t * cs_sdm_block_create(int n_max_blocks_by_row, int n_max_blocks_by_col, const int max_row_block_sizes[], const int max_col_block_sizes[])
Allocate and initialize a cs_sdm_t structure.
Definition: cs_sdm.c:231
cs_sdm_block_t::blocks
cs_sdm_t * blocks
Definition: cs_sdm.h:68
cs_sdm_create
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
cs_sdm_init
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:363
_cs_sdm_t::n_max_rows
int n_max_rows
Definition: cs_sdm.h:79
cs_fuel_incl::a
double precision, save a
Definition: cs_fuel_incl.f90:146
cs_sdm_block_square_asymm
void cs_sdm_block_square_asymm(cs_sdm_t *mat)
Set the given block matrix into its anti-symmetric part.
Definition: cs_sdm.c:1213
cs_sdm_66_ldlt_solve
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:1898
cs_sdm_update_matvec
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:859
cs_sdm_square_init
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:385
_cs_sdm_t
Definition: cs_sdm.h:74
cs_sdm_square_add_transpose
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:1101
cs_sdm_33_ldlt_compute
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:1523
cs_sdm_33_lu_compute
void cs_sdm_33_lu_compute(const cs_sdm_t *m, cs_real_t facto[9])
LU factorization of a small dense 3x3 matrix.
Definition: cs_sdm.c:1342
cs_fuel_incl::b
double precision, save b
Definition: cs_fuel_incl.f90:146
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:468
cs_sdm_block_add_mult
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:958
cs_sdm_33_ldlt_solve
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:1839
cs_sdm_copy
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:441
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:302
cs_sdm_add_mult
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:1074
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:467
CS_SDM_SHARED_VAL
#define CS_SDM_SHARED_VAL
Definition: cs_sdm.h:50
cs_sdm_map_array
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:324
cs_sdm_block_t
Definition: cs_sdm.h:58
cs_sdm_dump
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:2093
cs_sdm_block_t::n_max_blocks_by_col
int n_max_blocks_by_col
Definition: cs_sdm.h:62
cs_sdm_66_ldlt_compute
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:1627
cs_sdm_block_t::n_col_blocks
int n_col_blocks
Definition: cs_sdm.h:63
cs_sdm_square_2symm
void cs_sdm_square_2symm(cs_sdm_t *mat)
Set the given matrix to two times its symmetric part mat --> mat + mat_tr = 2*symm(mat)
Definition: cs_sdm.c:1146
cs_sdm_block_create_copy
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:458
cs_sdm_44_ldlt_compute
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:1570
cs_sdm_get_block
static cs_sdm_t * cs_sdm_get_block(const cs_sdm_t *const 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:483
cs_sdm_test_symmetry
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:2008
cs_sdm_block_multiply_rowrow_sym
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:724
cs_sdm_fprintf
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:2145
y
void const cs_lnum_t *const const cs_real_t *const const cs_real_t *const const cs_real_t *const const cs_real_t *const y
Definition: cs_wall_functions.h:1147
cs_sdm_get_col
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:508
cs_sdm_block_multiply_rowrow
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:668
cs_sdm_ldlt_solve
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:1935
cs_sdm_block_fprintf
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:2264
CS_SDM_BY_BLOCK
#define CS_SDM_BY_BLOCK
Definition: cs_sdm.h:48
cs_sdm_33_sym_qr_compute
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:1280
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:298
cs_sdm_block33_init
void cs_sdm_block33_init(cs_sdm_t *m, int n_row_blocks, int n_col_blocks)
Initialize the pattern of cs_sdm_t structure defined by 3x3 block The matrix should have been allocat...
Definition: cs_sdm.c:419
cs_sdm_block_t::n_row_blocks
int n_row_blocks
Definition: cs_sdm.h:61
cs_sdm_multiply_r1c3_rowrow
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:613
cs_sdm_dot
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
cs_sdm_copy_block
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:537
cs_sdm_transpose_and_update
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:568
cs_sdm_square_create
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
_cs_sdm_t::n_cols
int n_cols
Definition: cs_sdm.h:84
cs_sdm_block_matvec
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:999
cs_sdm_product_t
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
cs_sdm_multiply_rowrow_sym
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:620
cs_sdm_matvec_t
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
cs_sdm_add_scalvect
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
_cs_sdm_t::flag
cs_flag_t flag
Definition: cs_sdm.h:76
cs_flag_t
unsigned short int cs_flag_t
Definition: cs_defs.h:304
cs_sdm_matvec_transposed
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....
Definition: cs_sdm.c:892
cs_sdm_simple_dump
void cs_sdm_simple_dump(const cs_sdm_t *mat)
Dump a small dense matrix.
Definition: cs_sdm.c:2064
cs_sdm_lu_solve
void cs_sdm_lu_solve(cs_lnum_t n_rows, const cs_real_t facto[], const cs_real_t *rhs, cs_real_t *sol)
Solve a system A.sol = rhs using a LU factorization of A (a small dense matrix).
Definition: cs_sdm.c:1475
cs_sdm_matvec
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:817
cpincl::alpha
double precision, dimension(ncharm), save alpha
Definition: cpincl.f90:99
cs_sdm_multiply_rowrow
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:574
cs_sdm_33_lu_solve
void cs_sdm_33_lu_solve(const cs_real_t facto[9], const cs_real_t rhs[3], cs_real_t sol[3])
Solve a system A.sol = rhs using a LU factorization of A (a small 3x3 dense matrix).
Definition: cs_sdm.c:1439
cs_sdm_free
cs_sdm_t * cs_sdm_free(cs_sdm_t *mat)
Free a cs_sdm_t structure.
Definition: cs_sdm.c:331
_cs_sdm_t::n_rows
int n_rows
Definition: cs_sdm.h:80
cs_sdm_44_ldlt_solve
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:1867
cs_sdm_scalvect
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
cs_sdm_add
void cs_sdm_add(cs_sdm_t *mat, const cs_sdm_t *add)
Add two small dense matrices: loc += add.
Definition: cs_sdm.c:1051
cs_sdm_ldlt_compute
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:1717
_cs_sdm_t::n_max_cols
int n_max_cols
Definition: cs_sdm.h:83
atimbr::nc
double precision, dimension(:,:,:), allocatable nc
Definition: atimbr.f90:110