8.2
general 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-2024 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 <assert.h>
33 #include <stdio.h>
34 #include <string.h>
35 #include <math.h>
36 
37 /*----------------------------------------------------------------------------
38  * Local headers
39  *----------------------------------------------------------------------------*/
40 
41 /*----------------------------------------------------------------------------*/
42 
43 BEGIN_C_DECLS
44 
45 /*============================================================================
46  * Macro definitions
47  *============================================================================*/
48 
49 #define CS_SDM_BY_BLOCK (1 << 0) /* Matrix is defined by block */
50 #define CS_SDM_SYMMETRIC (1 << 1) /* Matrix is symmetric by construction */
51 #define CS_SDM_SHARED_VAL (1 << 2) /* Matrix is not owner of its values */
52 
53 /*============================================================================
54  * Type definitions
55  *============================================================================*/
56 
57 typedef struct _cs_sdm_t cs_sdm_t;
58 
59 typedef struct {
60 
61  int n_max_blocks_by_row;
62  int n_row_blocks;
63  int n_max_blocks_by_col;
64  int n_col_blocks;
65 
66  /* Allocated to n_max_blocks_by_row*n_max_blocks_by_col
67  Cast in cs_sdm_t where values are shared with the master cs_sdm_t struct.
68  */
69  cs_sdm_t *blocks;
70 
71 } cs_sdm_block_t;
72 
73 /* Structure enabling the repeated usage of Small Dense Matrices (SDM) during
74  the building of the linear system by a cellwise process */
75 struct _cs_sdm_t {
76 
77  cs_flag_t flag; /* Metadata */
78 
79  /* Row-related members */
80  int n_max_rows; // max number of entities by primal cells
81  int n_rows; // current number of entities
82 
83  /* Column-related members. Not useful if the matrix is square */
84  int n_max_cols; // max number of entries in a column
85  int n_cols; // current number of columns
86 
87  cs_real_t *val; // small dense matrix (size: n_max_rows*n_max_cols)
88 
89  /* Structure describing the matrix in terms of blocks */
90  cs_sdm_block_t *block_desc;
91 
92 };
93 
94 /*============================================================================
95  * Prototypes for pointer of functions
96  *============================================================================*/
97 
98 /*----------------------------------------------------------------------------*/
107 /*----------------------------------------------------------------------------*/
108 
109 typedef void
110 (cs_sdm_product_t)(const cs_sdm_t *a,
111  const cs_sdm_t *b,
112  cs_sdm_t *c);
113 
114 /*----------------------------------------------------------------------------*/
123 /*----------------------------------------------------------------------------*/
124 
125 typedef void
126 (cs_sdm_matvec_t)(const cs_sdm_t *mat,
127  const cs_real_t *vec,
128  cs_real_t *mv);
129 
130 /*============================================================================
131  * Public function prototypes
132  *============================================================================*/
133 
134 /*----------------------------------------------------------------------------*/
146 /*----------------------------------------------------------------------------*/
147 
148 static inline cs_real_t
149 cs_sdm_dot(int n,
150  const cs_real_t x[],
151  const cs_real_t y[])
152 {
153  cs_real_t dp = 0;
154 
155  if (x == NULL || y == NULL)
156  return dp;
157 
158  for (int i = 0; i < n; i++)
159  dp += x[i]*y[i];
160 
161  return dp;
162 }
163 
164 /*----------------------------------------------------------------------------*/
175 /*----------------------------------------------------------------------------*/
176 
177 static inline void
178 cs_sdm_scalvect(int n,
179  const cs_real_t s,
180  const cs_real_t x[],
181  cs_real_t y[])
182 {
183  if (x == NULL || y == NULL)
184  return;
185 
186  for (int i = 0; i < n; i++)
187  y[i] = s * x[i];
188 }
189 
190 /*----------------------------------------------------------------------------*/
201 /*----------------------------------------------------------------------------*/
202 
203 static inline void
204 cs_sdm_add_scalvect(int n,
205  const cs_real_t s,
206  const cs_real_t x[],
207  cs_real_t y[])
208 {
209  if (x == NULL || y == NULL)
210  return;
211 
212  for (int i = 0; i < n; i++)
213  y[i] += s * x[i];
214 }
215 
216 /*----------------------------------------------------------------------------*/
223 /*----------------------------------------------------------------------------*/
224 
225 static inline void
226 cs_sdm_symm_ur(cs_sdm_t *mat)
227 {
228  assert(mat != NULL);
229  for (int i = 1; i < mat->n_rows; i++) {
230  cs_real_t *m_i = mat->val + i*mat->n_rows;
231  for (int j = 0; j < i; j++)
232  m_i[j] = mat->val[j*mat->n_rows + i];
233  }
234 }
235 
236 /*----------------------------------------------------------------------------*/
247 /*----------------------------------------------------------------------------*/
248 
249 cs_sdm_t *
250 cs_sdm_create(cs_flag_t flag,
251  int n_max_rows,
252  int n_max_cols);
253 
254 /*----------------------------------------------------------------------------*/
263 /*----------------------------------------------------------------------------*/
264 
265 cs_sdm_t *
266 cs_sdm_square_create(int n_max_rows);
267 
268 /*----------------------------------------------------------------------------*/
277 /*----------------------------------------------------------------------------*/
278 
279 cs_sdm_t *
280 cs_sdm_create_copy(const cs_sdm_t *m);
281 
282 /*----------------------------------------------------------------------------*/
290 /*----------------------------------------------------------------------------*/
291 
292 cs_sdm_t *
293 cs_sdm_create_transpose(cs_sdm_t *mat);
294 
295 /*----------------------------------------------------------------------------*/
306 /*----------------------------------------------------------------------------*/
307 
308 cs_sdm_t *
309 cs_sdm_block_create(int n_max_blocks_by_row,
310  int n_max_blocks_by_col,
311  const int max_row_block_sizes[],
312  const int max_col_block_sizes[]);
313 
314 /*----------------------------------------------------------------------------*/
324 /*----------------------------------------------------------------------------*/
325 
326 cs_sdm_t *
327 cs_sdm_block33_create(int n_max_blocks_by_row,
328  int n_max_blocks_by_col);
329 
330 /*----------------------------------------------------------------------------*/
342 /*----------------------------------------------------------------------------*/
343 
344 static inline void
345 cs_sdm_map_array(int n_max_rows,
346  int n_max_cols,
347  cs_sdm_t *m,
348  cs_real_t *array)
349 {
350  assert(array != NULL && m != NULL); /* Sanity check */
351 
352  m->flag = CS_SDM_SHARED_VAL;
353  m->n_rows = m->n_max_rows = n_max_rows;
354  m->n_cols = m->n_max_cols = n_max_cols;
355  m->val = array;
356  m->block_desc = NULL;
357 }
358 
359 /*----------------------------------------------------------------------------*/
367 /*----------------------------------------------------------------------------*/
368 
369 cs_sdm_t *
370 cs_sdm_free(cs_sdm_t *mat);
371 
372 /*----------------------------------------------------------------------------*/
381 /*----------------------------------------------------------------------------*/
382 
383 static inline void
384 cs_sdm_init(int n_rows,
385  int n_cols,
386  cs_sdm_t *mat)
387 {
388  assert(mat != NULL);
389 
390  mat->n_rows = n_rows;
391  mat->n_cols = n_cols;
392  memset(mat->val, 0, n_rows*n_cols*sizeof(cs_real_t));
393 }
394 
395 /*----------------------------------------------------------------------------*/
403 /*----------------------------------------------------------------------------*/
404 
405 static inline void
406 cs_sdm_square_init(int n_rows,
407  cs_sdm_t *mat)
408 {
409  assert(mat != NULL);
410 
411  mat->n_rows = mat->n_cols = n_rows; /* square matrix */
412  memset(mat->val, 0, n_rows*n_rows*sizeof(cs_real_t));
413 }
414 
415 /*----------------------------------------------------------------------------*/
426 /*----------------------------------------------------------------------------*/
427 
428 void
429 cs_sdm_block_init(cs_sdm_t *m,
430  int n_row_blocks,
431  int n_col_blocks,
432  const int row_block_sizes[],
433  const int col_block_sizes[]);
434 
435 /*----------------------------------------------------------------------------*/
444 /*----------------------------------------------------------------------------*/
445 
446 void
447 cs_sdm_block33_init(cs_sdm_t *m,
448  int n_row_blocks,
449  int n_col_blocks);
450 
451 
452 /*----------------------------------------------------------------------------*/
460 /*----------------------------------------------------------------------------*/
461 
462 void
463 cs_sdm_block_33_to_xyz(const cs_sdm_t *mb33,
464  cs_sdm_t *mbxyz);
465 
466 /*----------------------------------------------------------------------------*/
474 /*----------------------------------------------------------------------------*/
475 
476 static inline void
477 cs_sdm_copy(cs_sdm_t *recv,
478  const cs_sdm_t *send)
479 {
480  /* Sanity check */
481  assert(recv->n_max_rows >= send->n_rows);
482  assert(recv->n_max_cols >= send->n_cols);
483 
484  recv->flag = send->flag;
485  recv->n_rows = send->n_rows;
486  recv->n_cols = send->n_cols;
487 
488  /* Copy values */
489  memcpy(recv->val, send->val, sizeof(cs_real_t)*send->n_rows*send->n_cols);
490 }
491 
492 /*----------------------------------------------------------------------------*/
501 /*----------------------------------------------------------------------------*/
502 
503 cs_sdm_t *
504 cs_sdm_block_create_copy(const cs_sdm_t *mref);
505 
506 /*----------------------------------------------------------------------------*/
516 /*----------------------------------------------------------------------------*/
517 
518 static inline cs_sdm_t *
519 cs_sdm_get_block(const cs_sdm_t *const m,
520  int row_block_id,
521  int col_block_id)
522 {
523  /* Sanity checks */
524  assert(m != NULL);
525  assert(m->flag & CS_SDM_BY_BLOCK && m->block_desc != NULL);
526  assert(col_block_id < m->block_desc->n_col_blocks);
527  assert(row_block_id < m->block_desc->n_row_blocks);
528 
529  return m->block_desc->blocks
530  + row_block_id*m->block_desc->n_col_blocks + col_block_id;
531 }
532 
533 /*----------------------------------------------------------------------------*/
541 /*----------------------------------------------------------------------------*/
542 
543 static inline void
544 cs_sdm_get_col(const cs_sdm_t *m,
545  int col_id,
546  cs_real_t *col_vals)
547 {
548  /* Sanity checks */
549  assert(m != NULL && col_vals != NULL);
550  assert(col_id < m->n_cols);
551 
552  const cs_real_t *_col = m->val + col_id;
553  for(int i = 0; i < m->n_rows; i++, _col += m->n_cols)
554  col_vals[i] = *_col;
555 }
556 
557 /*----------------------------------------------------------------------------*/
569 /*----------------------------------------------------------------------------*/
570 
571 static inline void
572 cs_sdm_copy_block(const cs_sdm_t *m,
573  const short int r_id,
574  const short int c_id,
575  const short int nr,
576  const short int nc,
577  cs_sdm_t *b)
578 {
579  /* Sanity checks */
580  assert(m != NULL && b != NULL);
581  assert(r_id >= 0 && c_id >= 0);
582  assert((r_id + nr) <= m->n_rows);
583  assert((c_id + nc) <= m->n_cols);
584  assert(nr == b->n_rows && nc == b->n_cols);
585 
586  const cs_real_t *_start = m->val + c_id + r_id*m->n_cols;
587  for (short int i = 0; i < nr; i++, _start += m->n_cols)
588  memcpy(b->val + i*nc, _start, sizeof(cs_real_t)*nc);
589 }
590 
591 /*----------------------------------------------------------------------------*/
598 /*----------------------------------------------------------------------------*/
599 
600 static inline void
601 cs_sdm_scale(cs_real_t scaling,
602  cs_sdm_t *m)
603 {
604  /* Sanity checks */
605  assert(m != NULL);
606 
607  for (int i = 0; i < m->n_rows*m->n_cols; i++)
608  m->val[i] *= scaling;
609 }
610 
611 /*----------------------------------------------------------------------------*/
620 /*----------------------------------------------------------------------------*/
621 
622 static inline void
623 cs_sdm_transpose_and_update(const cs_sdm_t *m,
624  cs_sdm_t *mt)
625 {
626  assert(m != NULL && mt != NULL);
627  assert(m->n_rows == mt->n_cols && m->n_cols == mt->n_rows);
628 
629  for (short int i = 0; i < m->n_rows; i++) {
630  const cs_real_t *m_i = m->val + i*m->n_cols;
631  for (short int j = 0; j < m->n_cols; j++)
632  mt->val[j*mt->n_cols + i] += m_i[j];
633  }
634 }
635 
636 /*----------------------------------------------------------------------------*/
646 /*----------------------------------------------------------------------------*/
647 
648 void
649 cs_sdm_multiply(const cs_sdm_t *a,
650  const cs_sdm_t *b,
651  cs_sdm_t *c);
652 
653 /*----------------------------------------------------------------------------*/
662 /*----------------------------------------------------------------------------*/
663 
664 static inline void
665 cs_sdm_33_matvec(const cs_sdm_t *mat,
666  const cs_real_t *vec,
667  cs_real_t *mv)
668 {
669  /* Sanity checks */
670  assert(mat != NULL && vec != NULL && mv != NULL);
671  assert(mat->n_rows == 3 && mat->n_cols == 3);
672 
673  mv[0] = vec[0]*mat->val[0] + vec[1]*mat->val[1] + vec[2]*mat->val[2];
674  mv[1] = vec[0]*mat->val[3] + vec[1]*mat->val[4] + vec[2]*mat->val[5];
675  mv[2] = vec[0]*mat->val[6] + vec[1]*mat->val[7] + vec[2]*mat->val[8];
676 }
677 
678 /*----------------------------------------------------------------------------*/
690 /*----------------------------------------------------------------------------*/
691 
692 static inline void
693 cs_sdm_multiply_r1c3_rowrow(const cs_sdm_t *a,
694  const cs_sdm_t *b,
695  cs_sdm_t *c)
696 {
697  /* Sanity check */
698  assert(a != NULL && b != NULL && c != NULL);
699  assert(a->n_cols == 3 && b->n_cols == 3 &&
700  a->n_rows == 1 && c->n_rows == 1 &&
701  c->n_cols == 1 && b->n_rows == 1);
702 
703  c->val[0] += a->val[0]*b->val[0] + a->val[1]*b->val[1] + a->val[2]*b->val[2];
704 }
705 
706 /*----------------------------------------------------------------------------*/
718 /*----------------------------------------------------------------------------*/
719 
720 void
721 cs_sdm_multiply_rowrow(const cs_sdm_t *a,
722  const cs_sdm_t *b,
723  cs_sdm_t *c);
724 
725 /*----------------------------------------------------------------------------*/
738 /*----------------------------------------------------------------------------*/
739 
740 void
741 cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a,
742  const cs_sdm_t *b,
743  cs_sdm_t *c);
744 
745 /*----------------------------------------------------------------------------*/
757 /*----------------------------------------------------------------------------*/
758 
759 void
760 cs_sdm_block_multiply_rowrow(const cs_sdm_t *a,
761  const cs_sdm_t *b,
762  cs_sdm_t *c);
763 
764 /*----------------------------------------------------------------------------*/
777 /*----------------------------------------------------------------------------*/
778 
779 void
780 cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a,
781  const cs_sdm_t *b,
782  cs_sdm_t *c);
783 
784 /*----------------------------------------------------------------------------*/
793 /*----------------------------------------------------------------------------*/
794 
795 void
796 cs_sdm_square_matvec(const cs_sdm_t *mat,
797  const cs_real_t *vec,
798  cs_real_t *mv);
799 
800 /*----------------------------------------------------------------------------*/
809 /*----------------------------------------------------------------------------*/
810 
811 void
812 cs_sdm_matvec(const cs_sdm_t *mat,
813  const cs_real_t *vec,
814  cs_real_t *mv);
815 
816 /*----------------------------------------------------------------------------*/
826 /*----------------------------------------------------------------------------*/
827 
828 void
829 cs_sdm_update_matvec(const cs_sdm_t *mat,
830  const cs_real_t *vec,
831  cs_real_t *mv);
832 
833 /*----------------------------------------------------------------------------*/
844 /*----------------------------------------------------------------------------*/
845 
846 void
847 cs_sdm_matvec_transposed(const cs_sdm_t *mat,
848  const cs_real_t *vec,
849  cs_real_t *mv);
850 
851 /*----------------------------------------------------------------------------*/
858 /*----------------------------------------------------------------------------*/
859 
860 void
861 cs_sdm_block_add(cs_sdm_t *mat,
862  const cs_sdm_t *add);
863 
864 /*----------------------------------------------------------------------------*/
872 /*----------------------------------------------------------------------------*/
873 
874 void
875 cs_sdm_block_add_mult(cs_sdm_t *mat,
876  cs_real_t mult_coef,
877  const cs_sdm_t *add);
878 
879 /*----------------------------------------------------------------------------*/
889 /*----------------------------------------------------------------------------*/
890 
891 void
892 cs_sdm_block_matvec(const cs_sdm_t *mat,
893  const cs_real_t *vec,
894  cs_real_t *mv);
895 
896 /*----------------------------------------------------------------------------*/
903 /*----------------------------------------------------------------------------*/
904 
905 void
906 cs_sdm_add(cs_sdm_t *mat,
907  const cs_sdm_t *add);
908 
909 /*----------------------------------------------------------------------------*/
917 /*----------------------------------------------------------------------------*/
918 
919 void
920 cs_sdm_add_mult(cs_sdm_t *mat,
921  cs_real_t alpha,
922  const cs_sdm_t *add);
923 
924 /*----------------------------------------------------------------------------*/
932 /*----------------------------------------------------------------------------*/
933 
934 void
935 cs_sdm_square_add_transpose(cs_sdm_t *mat,
936  cs_sdm_t *tr);
937 
938 /*----------------------------------------------------------------------------*/
945 /*----------------------------------------------------------------------------*/
946 
947 void
948 cs_sdm_square_2symm(cs_sdm_t *mat);
949 
950 /*----------------------------------------------------------------------------*/
956 /*----------------------------------------------------------------------------*/
957 
958 void
959 cs_sdm_square_asymm(cs_sdm_t *mat);
960 
961 /*----------------------------------------------------------------------------*/
967 /*----------------------------------------------------------------------------*/
968 
969 void
970 cs_sdm_block_square_asymm(cs_sdm_t *mat);
971 
972 /*----------------------------------------------------------------------------*/
989 /*----------------------------------------------------------------------------*/
990 
991 void
992 cs_sdm_33_sym_qr_compute(const cs_real_t m[9],
993  cs_real_t Qt[9],
994  cs_real_t R[6]);
995 
996 /*----------------------------------------------------------------------------*/
1007 /*----------------------------------------------------------------------------*/
1008 
1009 void
1010 cs_sdm_33_lu_compute(const cs_sdm_t *m,
1011  cs_real_t facto[9]);
1012 
1013 /*----------------------------------------------------------------------------*/
1026 /*----------------------------------------------------------------------------*/
1027 
1028 void
1029 cs_sdm_lu_compute(const cs_sdm_t *m,
1030  cs_real_t facto[]);
1031 
1032 /*----------------------------------------------------------------------------*/
1046 /*----------------------------------------------------------------------------*/
1047 
1048 void
1049 cs_sdm_33_lu_solve(const cs_real_t facto[9],
1050  const cs_real_t rhs[3],
1051  cs_real_t sol[3]);
1052 
1053 /*----------------------------------------------------------------------------*/
1068 /*----------------------------------------------------------------------------*/
1069 
1070 void
1071 cs_sdm_lu_solve(cs_lnum_t n_rows,
1072  const cs_real_t facto[],
1073  const cs_real_t *rhs,
1074  cs_real_t *sol);
1075 
1076 /*----------------------------------------------------------------------------*/
1090 /*----------------------------------------------------------------------------*/
1091 
1092 void
1093 cs_sdm_33_ldlt_compute(const cs_sdm_t *m,
1094  cs_real_t facto[6]);
1095 
1096 /*----------------------------------------------------------------------------*/
1110 /*----------------------------------------------------------------------------*/
1111 
1112 void
1113 cs_sdm_44_ldlt_compute(const cs_sdm_t *m,
1114  cs_real_t facto[10]);
1115 
1116 /*----------------------------------------------------------------------------*/
1130 /*----------------------------------------------------------------------------*/
1131 
1132 void
1133 cs_sdm_66_ldlt_compute(const cs_sdm_t *m,
1134  cs_real_t facto[21]);
1135 
1136 /*----------------------------------------------------------------------------*/
1151 /*----------------------------------------------------------------------------*/
1152 
1153 void
1154 cs_sdm_ldlt_compute(const cs_sdm_t *m,
1155  cs_real_t *facto,
1156  cs_real_t *dkk);
1157 
1158 /*----------------------------------------------------------------------------*/
1168 /*----------------------------------------------------------------------------*/
1169 
1170 void
1171 cs_sdm_33_ldlt_solve(const cs_real_t facto[6],
1172  const cs_real_t rhs[3],
1173  cs_real_t sol[3]);
1174 
1175 /*----------------------------------------------------------------------------*/
1185 /*----------------------------------------------------------------------------*/
1186 
1187 void
1188 cs_sdm_44_ldlt_solve(const cs_real_t facto[10],
1189  const cs_real_t rhs[4],
1190  cs_real_t x[4]);
1191 
1192 /*----------------------------------------------------------------------------*/
1202 /*----------------------------------------------------------------------------*/
1203 
1204 void
1205 cs_sdm_66_ldlt_solve(const cs_real_t f[21],
1206  const cs_real_t b[6],
1207  cs_real_t x[6]);
1208 
1209 /*----------------------------------------------------------------------------*/
1221 /*----------------------------------------------------------------------------*/
1222 
1223 void
1224 cs_sdm_ldlt_solve(int n_rows,
1225  const cs_real_t *facto,
1226  const cs_real_t *rhs,
1227  cs_real_t *sol);
1228 
1229 /*----------------------------------------------------------------------------*/
1239 /*----------------------------------------------------------------------------*/
1240 
1241 double
1242 cs_sdm_test_symmetry(const cs_sdm_t *mat);
1243 
1244 /*----------------------------------------------------------------------------*/
1250 /*----------------------------------------------------------------------------*/
1251 
1252 void
1253 cs_sdm_simple_dump(const cs_sdm_t *mat);
1254 
1255 /*----------------------------------------------------------------------------*/
1264 /*----------------------------------------------------------------------------*/
1265 
1266 void
1267 cs_sdm_dump(cs_lnum_t parent_id,
1268  const cs_lnum_t *row_ids,
1269  const cs_lnum_t *col_ids,
1270  const cs_sdm_t *mat);
1271 
1272 /*----------------------------------------------------------------------------*/
1285 /*----------------------------------------------------------------------------*/
1286 
1287 void
1288 cs_sdm_fprintf(FILE *fp,
1289  const char *fname,
1290  cs_real_t thd,
1291  const cs_sdm_t *m);
1292 
1293 /*----------------------------------------------------------------------------*/
1300 /*----------------------------------------------------------------------------*/
1301 
1302 void
1303 cs_sdm_block_dump(cs_lnum_t parent_id,
1304  const cs_sdm_t *mat);
1305 
1306 /*----------------------------------------------------------------------------*/
1319 /*----------------------------------------------------------------------------*/
1320 
1321 void
1322 cs_sdm_block_fprintf(FILE *fp,
1323  const char *fname,
1324  cs_real_t thd,
1325  const cs_sdm_t *m);
1326 
1327 /*----------------------------------------------------------------------------*/
1328 
1329 END_C_DECLS
1330 
1331 #endif /* __CS_SDM_H__ */
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:528
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:529
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
cs_flag_t
unsigned short int cs_flag_t
Definition: cs_defs.h:334
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:1600
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:768
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:292
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:1754
cs_sdm_block_33_to_xyz
void cs_sdm_block_33_to_xyz(const cs_sdm_t *mb33, cs_sdm_t *mbxyz)
Convert a matrix with each entry is a 3x3 block into a matrix with a block for each component x,...
Definition: cs_sdm.c:458
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:230
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_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:997
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:1904
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:1551
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:178
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:1934
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:623
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:365
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:1406
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_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:1037
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:345
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:126
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:901
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:2302
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:1307
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:110
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:1660
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_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:1087
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:149
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:572
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:2041
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:1503
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:2126
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:1207
CS_SDM_SHARED_VAL
#define CS_SDM_SHARED_VAL
Definition: cs_sdm.h:51
cs_sdm_scale
static void cs_sdm_scale(cs_real_t scaling, cs_sdm_t *m)
Multiply a matrix with the scaling factor given as parameter.
Definition: cs_sdm.h:601
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:583
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:421
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:204
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:859
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:1970
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:544
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:714
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:668
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:1368
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:933
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:2182
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:519
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:2231
cs_sdm_symm_ur
static void cs_sdm_symm_ur(cs_sdm_t *mat)
Set the lower left according to the upper right part in order to get a symmetric matrix.
Definition: cs_sdm.h:226
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:624
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:406
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:1109
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:693
cs_sdm_simple_dump
void cs_sdm_simple_dump(const cs_sdm_t *mat)
Dump a small dense matrix.
Definition: cs_sdm.c:2097
CS_SDM_BY_BLOCK
#define CS_SDM_BY_BLOCK
Definition: cs_sdm.h:49
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:507
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:477
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_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:1241
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:1135
cs_sdm_free
cs_sdm_t * cs_sdm_free(cs_sdm_t *mat)
Free a cs_sdm_t structure.
Definition: cs_sdm.c:333
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:1466
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:384
cs_sdm_33_matvec
static void cs_sdm_33_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense 3x3 matrix-vector product mv has been previously allocated.
Definition: cs_sdm.h:665
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:1877
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:823
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:962
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:1178
atimbr::nc
double precision, dimension(:,:,:), allocatable nc
Definition: atimbr.f90:109
_cs_sdm_t
Definition: cs_sdm.h:75
_cs_sdm_t::n_rows
int n_rows
Definition: cs_sdm.h:81
_cs_sdm_t::n_cols
int n_cols
Definition: cs_sdm.h:85
_cs_sdm_t::n_max_cols
int n_max_cols
Definition: cs_sdm.h:84
_cs_sdm_t::flag
cs_flag_t flag
Definition: cs_sdm.h:77
_cs_sdm_t::block_desc
cs_sdm_block_t * block_desc
Definition: cs_sdm.h:90
_cs_sdm_t::val
cs_real_t * val
Definition: cs_sdm.h:87
_cs_sdm_t::n_max_rows
int n_max_rows
Definition: cs_sdm.h:80
cs_sdm_block_t
Definition: cs_sdm.h:59
cs_sdm_block_t::n_max_blocks_by_row
int n_max_blocks_by_row
Definition: cs_sdm.h:61
cs_sdm_block_t::blocks
cs_sdm_t * blocks
Definition: cs_sdm.h:69
cs_sdm_block_t::n_row_blocks
int n_row_blocks
Definition: cs_sdm.h:62
cs_sdm_block_t::n_max_blocks_by_col
int n_max_blocks_by_col
Definition: cs_sdm.h:63
cs_sdm_block_t::n_col_blocks
int n_col_blocks
Definition: cs_sdm.h:64