8.3
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
43BEGIN_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
57typedef struct _cs_sdm_t cs_sdm_t;
58
59typedef 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 */
75struct _cs_sdm_t {
76
77 cs_flag_t flag; /* Metadata */
78
79 /* Row-related members */
80 int n_max_rows; // max number of entries in a row
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 // storage is row-major
89
90 /* Structure describing the matrix in terms of blocks */
91 cs_sdm_block_t *block_desc;
92
93};
94
95/*============================================================================
96 * Prototypes for pointer of functions
97 *============================================================================*/
98
99/*----------------------------------------------------------------------------*/
108/*----------------------------------------------------------------------------*/
109
110typedef void
111(cs_sdm_product_t)(const cs_sdm_t *a,
112 const cs_sdm_t *b,
113 cs_sdm_t *c);
114
115/*----------------------------------------------------------------------------*/
124/*----------------------------------------------------------------------------*/
125
126typedef void
127(cs_sdm_matvec_t)(const cs_sdm_t *mat,
128 const cs_real_t *vec,
129 cs_real_t *mv);
130
131/*============================================================================
132 * Public function prototypes
133 *============================================================================*/
134
135/*----------------------------------------------------------------------------*/
147/*----------------------------------------------------------------------------*/
148
149static inline cs_real_t
150cs_sdm_dot(int n,
151 const cs_real_t x[],
152 const cs_real_t y[])
153{
154 cs_real_t dp = 0;
155
156 if (x == NULL || y == NULL)
157 return dp;
158
159 for (int i = 0; i < n; i++)
160 dp += x[i]*y[i];
161
162 return dp;
163}
164
165/*----------------------------------------------------------------------------*/
176/*----------------------------------------------------------------------------*/
177
178static inline void
179cs_sdm_scalvect(int n,
180 const cs_real_t s,
181 const cs_real_t x[],
182 cs_real_t y[])
183{
184 if (x == NULL || y == NULL)
185 return;
186
187 for (int i = 0; i < n; i++)
188 y[i] = s * x[i];
189}
190
191/*----------------------------------------------------------------------------*/
202/*----------------------------------------------------------------------------*/
203
204static inline void
205cs_sdm_add_scalvect(int n,
206 const cs_real_t s,
207 const cs_real_t x[],
208 cs_real_t y[])
209{
210 if (x == NULL || y == NULL)
211 return;
212
213 for (int i = 0; i < n; i++)
214 y[i] += s * x[i];
215}
216
217/*----------------------------------------------------------------------------*/
224/*----------------------------------------------------------------------------*/
225
226static inline void
227cs_sdm_symm_ur(cs_sdm_t *mat)
228{
229 assert(mat != NULL);
230 for (int i = 1; i < mat->n_rows; i++) {
231 cs_real_t *m_i = mat->val + i*mat->n_rows;
232 for (int j = 0; j < i; j++)
233 m_i[j] = mat->val[j*mat->n_rows + i];
234 }
235}
236
237/*----------------------------------------------------------------------------*/
248/*----------------------------------------------------------------------------*/
249
250cs_sdm_t *
251cs_sdm_create(cs_flag_t flag,
252 int n_max_rows,
253 int n_max_cols);
254
255/*----------------------------------------------------------------------------*/
264/*----------------------------------------------------------------------------*/
265
266cs_sdm_t *
267cs_sdm_square_create(int n_max_rows);
268
269/*----------------------------------------------------------------------------*/
278/*----------------------------------------------------------------------------*/
279
280cs_sdm_t *
281cs_sdm_create_copy(const cs_sdm_t *m);
282
283/*----------------------------------------------------------------------------*/
291/*----------------------------------------------------------------------------*/
292
293cs_sdm_t *
294cs_sdm_create_transpose(cs_sdm_t *mat);
295
296/*----------------------------------------------------------------------------*/
307/*----------------------------------------------------------------------------*/
308
309cs_sdm_t *
310cs_sdm_block_create(int n_max_blocks_by_row,
311 int n_max_blocks_by_col,
312 const int max_row_block_sizes[],
313 const int max_col_block_sizes[]);
314
315/*----------------------------------------------------------------------------*/
325/*----------------------------------------------------------------------------*/
326
327cs_sdm_t *
328cs_sdm_block33_create(int n_max_blocks_by_row,
329 int n_max_blocks_by_col);
330
331/*----------------------------------------------------------------------------*/
343/*----------------------------------------------------------------------------*/
344
345static inline void
346cs_sdm_map_array(int n_max_rows,
347 int n_max_cols,
348 cs_sdm_t *m,
349 cs_real_t *array)
350{
351 assert(array != NULL && m != NULL); /* Sanity check */
352
353 m->flag = CS_SDM_SHARED_VAL;
354 m->n_rows = m->n_max_rows = n_max_rows;
355 m->n_cols = m->n_max_cols = n_max_cols;
356 m->val = array;
357 m->block_desc = NULL;
358}
359
360/*----------------------------------------------------------------------------*/
368/*----------------------------------------------------------------------------*/
369
370cs_sdm_t *
371cs_sdm_free(cs_sdm_t *mat);
372
373/*----------------------------------------------------------------------------*/
382/*----------------------------------------------------------------------------*/
383
384static inline void
385cs_sdm_init(int n_rows,
386 int n_cols,
387 cs_sdm_t *mat)
388{
389 assert(mat != NULL);
390 assert(n_rows <= mat->n_max_rows);
391 assert(n_cols <= mat->n_max_cols);
392
393 mat->n_rows = n_rows;
394 mat->n_cols = n_cols;
395 memset(mat->val, 0, n_rows*n_cols*sizeof(cs_real_t));
396}
397
398/*----------------------------------------------------------------------------*/
406/*----------------------------------------------------------------------------*/
407
408static inline void
409cs_sdm_square_init(int n_rows,
410 cs_sdm_t *mat)
411{
412 cs_sdm_init(n_rows, n_rows, mat);
413}
414
415/*----------------------------------------------------------------------------*/
426/*----------------------------------------------------------------------------*/
427
428void
429cs_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
446void
447cs_sdm_block33_init(cs_sdm_t *m,
448 int n_row_blocks,
449 int n_col_blocks);
450
451
452/*----------------------------------------------------------------------------*/
460/*----------------------------------------------------------------------------*/
461
462void
463cs_sdm_block_33_to_xyz(const cs_sdm_t *mb33,
464 cs_sdm_t *mbxyz);
465
466/*----------------------------------------------------------------------------*/
474/*----------------------------------------------------------------------------*/
475
476static inline void
477cs_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
503cs_sdm_t *
504cs_sdm_block_create_copy(const cs_sdm_t *mref);
505
506/*----------------------------------------------------------------------------*/
516/*----------------------------------------------------------------------------*/
517
518static inline cs_sdm_t *
519cs_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
543static inline void
544cs_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
571static inline void
572cs_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
600static inline void
601cs_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
622static inline void
623cs_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
648void
649cs_sdm_multiply(const cs_sdm_t *a,
650 const cs_sdm_t *b,
651 cs_sdm_t *c);
652
653/*----------------------------------------------------------------------------*/
662/*----------------------------------------------------------------------------*/
663
664static inline void
665cs_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
692static inline void
693cs_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
720void
721cs_sdm_multiply_rowrow(const cs_sdm_t *a,
722 const cs_sdm_t *b,
723 cs_sdm_t *c);
724
725/*----------------------------------------------------------------------------*/
738/*----------------------------------------------------------------------------*/
739
740void
741cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a,
742 const cs_sdm_t *b,
743 cs_sdm_t *c);
744
745/*----------------------------------------------------------------------------*/
757/*----------------------------------------------------------------------------*/
758
759void
760cs_sdm_block_multiply_rowrow(const cs_sdm_t *a,
761 const cs_sdm_t *b,
762 cs_sdm_t *c);
763
764/*----------------------------------------------------------------------------*/
777/*----------------------------------------------------------------------------*/
778
779void
780cs_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
795void
796cs_sdm_square_matvec(const cs_sdm_t *mat,
797 const cs_real_t *vec,
798 cs_real_t *mv);
799
800/*----------------------------------------------------------------------------*/
809/*----------------------------------------------------------------------------*/
810
811void
812cs_sdm_matvec(const cs_sdm_t *mat,
813 const cs_real_t *vec,
814 cs_real_t *mv);
815
816/*----------------------------------------------------------------------------*/
826/*----------------------------------------------------------------------------*/
827
828void
829cs_sdm_update_matvec(const cs_sdm_t *mat,
830 const cs_real_t *vec,
831 cs_real_t *mv);
832
833/*----------------------------------------------------------------------------*/
844/*----------------------------------------------------------------------------*/
845
846void
847cs_sdm_matvec_transposed(const cs_sdm_t *mat,
848 const cs_real_t *vec,
849 cs_real_t *mv);
850
851/*----------------------------------------------------------------------------*/
858/*----------------------------------------------------------------------------*/
859
860void
861cs_sdm_block_add(cs_sdm_t *mat,
862 const cs_sdm_t *add);
863
864/*----------------------------------------------------------------------------*/
872/*----------------------------------------------------------------------------*/
873
874void
875cs_sdm_block_add_mult(cs_sdm_t *mat,
876 cs_real_t mult_coef,
877 const cs_sdm_t *add);
878
879/*----------------------------------------------------------------------------*/
889/*----------------------------------------------------------------------------*/
890
891void
892cs_sdm_block_matvec(const cs_sdm_t *mat,
893 const cs_real_t *vec,
894 cs_real_t *mv);
895
896/*----------------------------------------------------------------------------*/
903/*----------------------------------------------------------------------------*/
904
905void
906cs_sdm_add(cs_sdm_t *mat,
907 const cs_sdm_t *add);
908
909/*----------------------------------------------------------------------------*/
917/*----------------------------------------------------------------------------*/
918
919void
920cs_sdm_add_mult(cs_sdm_t *mat,
921 cs_real_t alpha,
922 const cs_sdm_t *add);
923
924/*----------------------------------------------------------------------------*/
936/*----------------------------------------------------------------------------*/
937
938void cs_sdm_add_block(cs_sdm_t *mat,
939 const short int r_id,
940 const short int c_id,
941 const short int nr,
942 const short int nc,
943 const cs_sdm_t *add);
944
945/*----------------------------------------------------------------------------*/
955/*----------------------------------------------------------------------------*/
956
957void cs_sdm_add_block_topleft(cs_sdm_t *mat,
958 const short int nr,
959 const short int nc,
960 const cs_sdm_t *add);
961
962/*----------------------------------------------------------------------------*/
970/*----------------------------------------------------------------------------*/
971
972void
973cs_sdm_square_add_transpose(cs_sdm_t *mat,
974 cs_sdm_t *tr);
975
976/*----------------------------------------------------------------------------*/
983/*----------------------------------------------------------------------------*/
984
985void
986cs_sdm_square_2symm(cs_sdm_t *mat);
987
988/*----------------------------------------------------------------------------*/
994/*----------------------------------------------------------------------------*/
995
996void
997cs_sdm_square_asymm(cs_sdm_t *mat);
998
999/*----------------------------------------------------------------------------*/
1005/*----------------------------------------------------------------------------*/
1006
1007void
1008cs_sdm_block_square_asymm(cs_sdm_t *mat);
1009
1010/*----------------------------------------------------------------------------*/
1027/*----------------------------------------------------------------------------*/
1028
1029void
1030cs_sdm_33_sym_qr_compute(const cs_real_t m[9],
1031 cs_real_t Qt[9],
1032 cs_real_t R[6]);
1033
1034/*----------------------------------------------------------------------------*/
1045/*----------------------------------------------------------------------------*/
1046
1047void
1048cs_sdm_33_lu_compute(const cs_sdm_t *m,
1049 cs_real_t facto[9]);
1050
1051/*----------------------------------------------------------------------------*/
1064/*----------------------------------------------------------------------------*/
1065
1066void
1067cs_sdm_lu_compute(const cs_sdm_t *m,
1068 cs_real_t facto[]);
1069
1070/*----------------------------------------------------------------------------*/
1084/*----------------------------------------------------------------------------*/
1085
1086void
1087cs_sdm_33_lu_solve(const cs_real_t facto[9],
1088 const cs_real_t rhs[3],
1089 cs_real_t sol[3]);
1090
1091/*----------------------------------------------------------------------------*/
1106/*----------------------------------------------------------------------------*/
1107
1108void
1109cs_sdm_lu_solve(cs_lnum_t n_rows,
1110 const cs_real_t facto[],
1111 const cs_real_t *rhs,
1112 cs_real_t *sol);
1113
1114/*----------------------------------------------------------------------------*/
1128/*----------------------------------------------------------------------------*/
1129
1130void
1131cs_sdm_33_ldlt_compute(const cs_sdm_t *m,
1132 cs_real_t facto[6]);
1133
1134/*----------------------------------------------------------------------------*/
1148/*----------------------------------------------------------------------------*/
1149
1150void
1151cs_sdm_44_ldlt_compute(const cs_sdm_t *m,
1152 cs_real_t facto[10]);
1153
1154/*----------------------------------------------------------------------------*/
1168/*----------------------------------------------------------------------------*/
1169
1170void
1171cs_sdm_66_ldlt_compute(const cs_sdm_t *m,
1172 cs_real_t facto[21]);
1173
1174/*----------------------------------------------------------------------------*/
1189/*----------------------------------------------------------------------------*/
1190
1191void
1192cs_sdm_ldlt_compute(const cs_sdm_t *m,
1193 cs_real_t *facto,
1194 cs_real_t *dkk);
1195
1196/*----------------------------------------------------------------------------*/
1206/*----------------------------------------------------------------------------*/
1207
1208void
1209cs_sdm_33_ldlt_solve(const cs_real_t facto[6],
1210 const cs_real_t rhs[3],
1211 cs_real_t sol[3]);
1212
1213/*----------------------------------------------------------------------------*/
1223/*----------------------------------------------------------------------------*/
1224
1225void
1226cs_sdm_44_ldlt_solve(const cs_real_t facto[10],
1227 const cs_real_t rhs[4],
1228 cs_real_t x[4]);
1229
1230/*----------------------------------------------------------------------------*/
1240/*----------------------------------------------------------------------------*/
1241
1242void
1243cs_sdm_66_ldlt_solve(const cs_real_t f[21],
1244 const cs_real_t b[6],
1245 cs_real_t x[6]);
1246
1247/*----------------------------------------------------------------------------*/
1259/*----------------------------------------------------------------------------*/
1260
1261void
1262cs_sdm_ldlt_solve(int n_rows,
1263 const cs_real_t *facto,
1264 const cs_real_t *rhs,
1265 cs_real_t *sol);
1266
1267/*----------------------------------------------------------------------------*/
1277/*----------------------------------------------------------------------------*/
1278
1279double
1280cs_sdm_test_symmetry(const cs_sdm_t *mat);
1281
1282/*----------------------------------------------------------------------------*/
1288/*----------------------------------------------------------------------------*/
1289
1290void
1291cs_sdm_simple_dump(const cs_sdm_t *mat);
1292
1293/*----------------------------------------------------------------------------*/
1302/*----------------------------------------------------------------------------*/
1303
1304void
1305cs_sdm_dump(cs_lnum_t parent_id,
1306 const cs_lnum_t *row_ids,
1307 const cs_lnum_t *col_ids,
1308 const cs_sdm_t *mat);
1309
1310/*----------------------------------------------------------------------------*/
1323/*----------------------------------------------------------------------------*/
1324
1325void
1326cs_sdm_fprintf(FILE *fp,
1327 const char *fname,
1328 cs_real_t thd,
1329 const cs_sdm_t *m);
1330
1331/*----------------------------------------------------------------------------*/
1338/*----------------------------------------------------------------------------*/
1339
1340void
1341cs_sdm_block_dump(cs_lnum_t parent_id,
1342 const cs_sdm_t *mat);
1343
1344/*----------------------------------------------------------------------------*/
1357/*----------------------------------------------------------------------------*/
1358
1359void
1360cs_sdm_block_fprintf(FILE *fp,
1361 const char *fname,
1362 cs_real_t thd,
1363 const cs_sdm_t *m);
1364
1365/*----------------------------------------------------------------------------*/
1366
1367END_C_DECLS
1368
1369#endif /* __CS_SDM_H__ */
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:542
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:342
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:543
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:335
cs_flag_t
unsigned short int cs_flag_t
Definition: cs_defs.h:344
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.cpp:1662
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.cpp:768
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.cpp:1816
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.cpp:458
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.cpp: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.cpp:1966
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.cpp:1613
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:179
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.cpp:291
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.cpp:1996
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.cpp:137
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.cpp:364
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_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.cpp:1468
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.cpp:507
cs_sdm_free
cs_sdm_t * cs_sdm_free(cs_sdm_t *mat)
Free a cs_sdm_t structure.
Definition: cs_sdm.cpp:332
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.cpp: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:346
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:127
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.cpp: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.cpp:2364
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.cpp:1369
cs_sdm_add_block_topleft
void cs_sdm_add_block_topleft(cs_sdm_t *mat, const short int nr, const short int nc, const cs_sdm_t *add)
Add a block of a matrix into a sub-matrix starting from (0, 0) with a size of nr rows and nc cols.
Definition: cs_sdm.cpp:1178
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:111
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.cpp:1722
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.cpp: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:150
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.cpp:2103
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.cpp:1565
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.cpp:2188
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.cpp:156
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.cpp:1269
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.cpp:583
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.cpp:195
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.cpp: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:205
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.cpp: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.cpp:2032
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.cpp: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.cpp: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.cpp:1430
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.cpp: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.cpp:2244
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.cpp:2293
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:227
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.cpp: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:409
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.cpp: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.cpp:2159
CS_SDM_BY_BLOCK
#define CS_SDM_BY_BLOCK
Definition: cs_sdm.h:49
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_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.cpp:1303
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.cpp:1197
cs_sdm_add_block
void cs_sdm_add_block(cs_sdm_t *mat, const short int r_id, const short int c_id, const short int nr, const short int nc, const cs_sdm_t *add)
Add a block of a matrix into a sub-matrix starting from (r_id, c_id) with a size of nr rows and nc co...
Definition: cs_sdm.cpp:1139
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.cpp:1528
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.cpp:229
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.cpp:173
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:385
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.cpp:1939
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.cpp: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.cpp: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.cpp:1240
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:91
_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