9.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-2026 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 <cassert>
33#include <cstring>
34
35/*----------------------------------------------------------------------------
36 * Local headers
37 *----------------------------------------------------------------------------*/
38
39
40/*============================================================================
41 * Macro definitions
42 *============================================================================*/
43
44#define CS_SDM_BY_BLOCK (1 << 0) /* Matrix is defined by block */
45#define CS_SDM_SYMMETRIC (1 << 1) /* Matrix is symmetric by construction */
46#define CS_SDM_SHARED_VAL (1 << 2) /* Matrix is not owner of its values */
47
48/*============================================================================
49 * Type definitions
50 *============================================================================*/
51
52typedef struct _cs_sdm_t cs_sdm_t;
53
54typedef struct {
55
56 int n_max_blocks_by_row;
57 int n_row_blocks;
58 int n_max_blocks_by_col;
59 int n_col_blocks;
60
61 /* Allocated to n_max_blocks_by_row*n_max_blocks_by_col
62 Cast in cs_sdm_t where values are shared with the master cs_sdm_t struct.
63 */
64 cs_sdm_t *blocks;
65
66} cs_sdm_block_t;
67
68/* Structure enabling the repeated usage of Small Dense Matrices (SDM) during
69 the building of the linear system by a cellwise process */
70struct _cs_sdm_t {
71
72 cs_flag_t flag; /* Metadata */
73
74 /* Row-related members */
75 int n_max_rows; // max number of entries in a row
76 int n_rows; // current number of entities
77
78 /* Column-related members. Not useful if the matrix is square */
79 int n_max_cols; // max number of entries in a column
80 int n_cols; // current number of columns
81
82 cs_real_t *val; // small dense matrix (size: n_max_rows*n_max_cols)
83 // storage is row-major
84
85 /* Structure describing the matrix in terms of blocks */
86 cs_sdm_block_t *block_desc;
87
88 /*----------------------------------------------------------------------------*/
93 /*----------------------------------------------------------------------------*/
94 ~_cs_sdm_t();
95
96 /*----------------------------------------------------------------------------*/
103 /*----------------------------------------------------------------------------*/
104 _cs_sdm_t(const int n_max_rows_);
105
106 /*----------------------------------------------------------------------------*/
115 /*----------------------------------------------------------------------------*/
116 _cs_sdm_t(const cs_flag_t flag_,
117 const int n_max_rows_,
118 const int n_max_cols_);
119
120 /*----------------------------------------------------------------------------*/
127 /*----------------------------------------------------------------------------*/
128 _cs_sdm_t(const _cs_sdm_t &other);
129
130 /*----------------------------------------------------------------------------*/
135 /*----------------------------------------------------------------------------*/
136
137 inline int
138 size() const
139 {
140 return n_rows * n_cols;
141 };
142
143 /*----------------------------------------------------------------------------*/
151 /*----------------------------------------------------------------------------*/
152
153 inline cs_real_t *
154 row(const int row_i) const
155 {
156 assert(0 <= row_i && row_i < n_rows);
157 return val + row_i * n_cols;
158 }
159
160 /* Deleted for the moment */
161 _cs_sdm_t &
162 operator=(const _cs_sdm_t &) = delete;
163
164 /*----------------------------------------------------------------------------*/
173 /*----------------------------------------------------------------------------*/
174 cs_real_t &
175 operator()(const int i, const int j);
176
177 /*----------------------------------------------------------------------------*/
186 /*----------------------------------------------------------------------------*/
187 const cs_real_t &
188 operator()(const int i, const int j) const;
189
190 /*----------------------------------------------------------------------------*/
198 /*----------------------------------------------------------------------------*/
199
200 void
201 init(const int nrows, const int ncols);
202
203 /*----------------------------------------------------------------------------*/
210 /*----------------------------------------------------------------------------*/
211
212 void
213 init(const int nrows);
214
215 /*----------------------------------------------------------------------------*/
220 /*----------------------------------------------------------------------------*/
221
222 inline cs_real_t *
223 data() const
224 {
225 return val;
226 };
227
228 /*----------------------------------------------------------------------------*/
234 /*----------------------------------------------------------------------------*/
235
236 _cs_sdm_t
237 transpose() const;
238
239 /*----------------------------------------------------------------------------*/
245 /*----------------------------------------------------------------------------*/
246
247 void
248 symmetrize_ur();
249
250 /*----------------------------------------------------------------------------*/
259 /*----------------------------------------------------------------------------*/
260
261 inline cs_sdm_t *
262 get_block(const int row_block_id, const int col_block_id) const
263 {
264 /* Sanity checks */
265 assert(flag & CS_SDM_BY_BLOCK && block_desc != nullptr);
266 assert(col_block_id < block_desc->n_col_blocks);
267 assert(row_block_id < block_desc->n_row_blocks);
268
269 return block_desc->blocks + row_block_id * block_desc->n_col_blocks +
270 col_block_id;
271 };
272
273 /*----------------------------------------------------------------------------*/
279 /*----------------------------------------------------------------------------*/
280
281 _cs_sdm_t &
282 operator*=(const cs_real_t &scaling);
283
284 /*----------------------------------------------------------------------------*/
290 /*----------------------------------------------------------------------------*/
291
292 _cs_sdm_t &
293 operator+=(const _cs_sdm_t &add);
294
295 /*----------------------------------------------------------------------------*/
303 /*----------------------------------------------------------------------------*/
304
305 void
306 dot(const cs_real_t vec[], cs_real_t mv[]) const;
307
308 /*----------------------------------------------------------------------------*/
312 /*----------------------------------------------------------------------------*/
313
314 void
315 dump() const;
316
317 /*----------------------------------------------------------------------------*/
323 /*----------------------------------------------------------------------------*/
324
325 void
326 dump(cs_lnum_t parent_id) const;
327
328 /*----------------------------------------------------------------------------*/
337 /*----------------------------------------------------------------------------*/
338
339 void
340 dump(cs_lnum_t parent_id,
341 const cs_lnum_t *row_ids,
342 const cs_lnum_t *col_ids) const;
343
344 /*----------------------------------------------------------------------------*/
356 /*----------------------------------------------------------------------------*/
357
358 void
359 dump(FILE *fp, const char *fname, cs_real_t thd) const;
360
361 /*----------------------------------------------------------------------------*/
372 /*----------------------------------------------------------------------------*/
373
374 void
375 map_array(int n_max_rows_, int n_max_cols_, cs_real_t *array);
376};
377
378namespace cs {
379using sdm_t = cs_sdm_t;
380}
381
382/*============================================================================
383 * Public function prototypes
384 *============================================================================*/
385
386/*----------------------------------------------------------------------------*/
398/*----------------------------------------------------------------------------*/
399
400static inline cs_real_t
401cs_sdm_dot(int n,
402 const cs_real_t x[],
403 const cs_real_t y[])
404{
405 cs_real_t dp = 0;
406
407 if (x == NULL || y == NULL)
408 return dp;
409
410 for (int i = 0; i < n; i++)
411 dp += x[i]*y[i];
412
413 return dp;
414}
415
416/*----------------------------------------------------------------------------*/
427/*----------------------------------------------------------------------------*/
428
429static inline void
430cs_sdm_scalvect(int n,
431 const cs_real_t s,
432 const cs_real_t x[],
433 cs_real_t y[])
434{
435 if (x == NULL || y == NULL)
436 return;
437
438 for (int i = 0; i < n; i++)
439 y[i] = s * x[i];
440}
441
442/*----------------------------------------------------------------------------*/
453/*----------------------------------------------------------------------------*/
454
455static inline void
456cs_sdm_add_scalvect(int n,
457 const cs_real_t s,
458 const cs_real_t x[],
459 cs_real_t y[])
460{
461 if (x == NULL || y == NULL)
462 return;
463
464 for (int i = 0; i < n; i++)
465 y[i] += s * x[i];
466}
467
468/*----------------------------------------------------------------------------*/
479/*----------------------------------------------------------------------------*/
480
481cs_sdm_t *
482cs_sdm_create(cs_flag_t flag,
483 int n_max_rows,
484 int n_max_cols);
485
486/*----------------------------------------------------------------------------*/
495/*----------------------------------------------------------------------------*/
496
497cs_sdm_t *
498cs_sdm_square_create(int n_max_rows);
499
500/*----------------------------------------------------------------------------*/
509/*----------------------------------------------------------------------------*/
510
511cs_sdm_t *
512cs_sdm_create_copy(const cs_sdm_t *m);
513
514/*----------------------------------------------------------------------------*/
522/*----------------------------------------------------------------------------*/
523
524cs_sdm_t *
525cs_sdm_create_transpose(cs_sdm_t *mat);
526
527/*----------------------------------------------------------------------------*/
538/*----------------------------------------------------------------------------*/
539
540cs_sdm_t *
541cs_sdm_block_create(int n_max_blocks_by_row,
542 int n_max_blocks_by_col,
543 const int max_row_block_sizes[],
544 const int max_col_block_sizes[]);
545
546/*----------------------------------------------------------------------------*/
556/*----------------------------------------------------------------------------*/
557
558cs_sdm_t *
559cs_sdm_block33_create(int n_max_blocks_by_row,
560 int n_max_blocks_by_col);
561
562/*----------------------------------------------------------------------------*/
570/*----------------------------------------------------------------------------*/
571
572cs_sdm_t *
573cs_sdm_free(cs_sdm_t *mat);
574
575/*----------------------------------------------------------------------------*/
586/*----------------------------------------------------------------------------*/
587
588void
589cs_sdm_block_init(cs_sdm_t *m,
590 int n_row_blocks,
591 int n_col_blocks,
592 const int row_block_sizes[],
593 const int col_block_sizes[]);
594
595/*----------------------------------------------------------------------------*/
604/*----------------------------------------------------------------------------*/
605
606void
607cs_sdm_block33_init(cs_sdm_t *m,
608 int n_row_blocks,
609 int n_col_blocks);
610
611
612/*----------------------------------------------------------------------------*/
620/*----------------------------------------------------------------------------*/
621
622void
623cs_sdm_block_33_to_xyz(const cs_sdm_t *mb33,
624 cs_sdm_t *mbxyz);
625
626/*----------------------------------------------------------------------------*/
634/*----------------------------------------------------------------------------*/
635
636static inline void
637cs_sdm_copy(cs_sdm_t *recv,
638 const cs_sdm_t *send)
639{
640 /* Sanity check */
641 assert(recv->n_max_rows >= send->n_rows);
642 assert(recv->n_max_cols >= send->n_cols);
643
644 recv->flag = send->flag;
645 recv->n_rows = send->n_rows;
646 recv->n_cols = send->n_cols;
647
648 /* Copy values */
649 std::memcpy(recv->val,
650 send->val,
651 sizeof(cs_real_t) * send->n_rows * send->n_cols);
652}
653
654/*----------------------------------------------------------------------------*/
663/*----------------------------------------------------------------------------*/
664
665cs_sdm_t *
666cs_sdm_block_create_copy(const cs_sdm_t *mref);
667
668/*----------------------------------------------------------------------------*/
680/*----------------------------------------------------------------------------*/
681
682static inline void
683cs_sdm_copy_block(const cs_sdm_t *m,
684 const short int r_id,
685 const short int c_id,
686 const short int nr,
687 const short int nc,
688 cs_sdm_t *b)
689{
690 /* Sanity checks */
691 assert(m != NULL && b != NULL);
692 assert(r_id >= 0 && c_id >= 0);
693 assert((r_id + nr) <= m->n_rows);
694 assert((c_id + nc) <= m->n_cols);
695 assert(nr == b->n_rows && nc == b->n_cols);
696
697 const cs_real_t *_start = m->val + c_id + r_id*m->n_cols;
698 for (short int i = 0; i < nr; i++, _start += m->n_cols)
699 std::memcpy(b->val + i * nc, _start, sizeof(cs_real_t) * nc);
700}
701
702/*----------------------------------------------------------------------------*/
711/*----------------------------------------------------------------------------*/
712
713static inline void
714cs_sdm_transpose_and_update(const cs_sdm_t *m,
715 cs_sdm_t *mt)
716{
717 assert(m != NULL && mt != NULL);
718 assert(m->n_rows == mt->n_cols && m->n_cols == mt->n_rows);
719
720 for (short int i = 0; i < m->n_rows; i++) {
721 const cs_real_t *m_i = m->val + i*m->n_cols;
722 for (short int j = 0; j < m->n_cols; j++)
723 mt->val[j*mt->n_cols + i] += m_i[j];
724 }
725}
726
727/*----------------------------------------------------------------------------*/
737/*----------------------------------------------------------------------------*/
738
739void
740cs_sdm_multiply(const cs_sdm_t *a,
741 const cs_sdm_t *b,
742 cs_sdm_t *c);
743
744/*----------------------------------------------------------------------------*/
756/*----------------------------------------------------------------------------*/
757
758static inline void
759cs_sdm_multiply_r1c3_rowrow(const cs_sdm_t *a,
760 const cs_sdm_t *b,
761 cs_sdm_t *c)
762{
763 /* Sanity check */
764 assert(a != NULL && b != NULL && c != NULL);
765 assert(a->n_cols == 3 && b->n_cols == 3 &&
766 a->n_rows == 1 && c->n_rows == 1 &&
767 c->n_cols == 1 && b->n_rows == 1);
768
769 c->val[0] += a->val[0]*b->val[0] + a->val[1]*b->val[1] + a->val[2]*b->val[2];
770}
771
772/*----------------------------------------------------------------------------*/
784/*----------------------------------------------------------------------------*/
785
786void
787cs_sdm_multiply_rowrow(const cs_sdm_t *a,
788 const cs_sdm_t *b,
789 cs_sdm_t *c);
790
791/*----------------------------------------------------------------------------*/
804/*----------------------------------------------------------------------------*/
805
806void
807cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a,
808 const cs_sdm_t *b,
809 cs_sdm_t *c);
810
811/*----------------------------------------------------------------------------*/
823/*----------------------------------------------------------------------------*/
824
825void
826cs_sdm_block_multiply_rowrow(const cs_sdm_t *a,
827 const cs_sdm_t *b,
828 cs_sdm_t *c);
829
830/*----------------------------------------------------------------------------*/
843/*----------------------------------------------------------------------------*/
844
845void
846cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a,
847 const cs_sdm_t *b,
848 cs_sdm_t *c);
849
850/*----------------------------------------------------------------------------*/
860/*----------------------------------------------------------------------------*/
861
862void
863cs_sdm_update_matvec(const cs_sdm_t *mat,
864 const cs_real_t *vec,
865 cs_real_t *mv);
866
867/*----------------------------------------------------------------------------*/
878/*----------------------------------------------------------------------------*/
879
880void
881cs_sdm_matvec_transposed(const cs_sdm_t *mat,
882 const cs_real_t *vec,
883 cs_real_t *mv);
884
885/*----------------------------------------------------------------------------*/
892/*----------------------------------------------------------------------------*/
893
894void
895cs_sdm_block_add(cs_sdm_t *mat,
896 const cs_sdm_t *add);
897
898/*----------------------------------------------------------------------------*/
906/*----------------------------------------------------------------------------*/
907
908void
909cs_sdm_block_add_mult(cs_sdm_t *mat,
910 cs_real_t mult_coef,
911 const cs_sdm_t *add);
912
913/*----------------------------------------------------------------------------*/
923/*----------------------------------------------------------------------------*/
924
925void
926cs_sdm_block_matvec(const cs_sdm_t *mat,
927 const cs_real_t *vec,
928 cs_real_t *mv);
929
930/*----------------------------------------------------------------------------*/
938/*----------------------------------------------------------------------------*/
939
940void
941cs_sdm_add_mult(cs_sdm_t *mat,
942 cs_real_t alpha,
943 const cs_sdm_t *add);
944
945/*----------------------------------------------------------------------------*/
957/*----------------------------------------------------------------------------*/
958
959void cs_sdm_add_block(cs_sdm_t *mat,
960 const short int r_id,
961 const short int c_id,
962 const short int nr,
963 const short int nc,
964 const cs_sdm_t *add);
965
966/*----------------------------------------------------------------------------*/
976/*----------------------------------------------------------------------------*/
977
978void cs_sdm_add_block_topleft(cs_sdm_t *mat,
979 const short int nr,
980 const short int nc,
981 const cs_sdm_t *add);
982
983/*----------------------------------------------------------------------------*/
991/*----------------------------------------------------------------------------*/
992
993void
994cs_sdm_square_add_transpose(cs_sdm_t *mat,
995 cs_sdm_t *tr);
996
997/*----------------------------------------------------------------------------*/
1004/*----------------------------------------------------------------------------*/
1005
1006void
1007cs_sdm_square_2symm(cs_sdm_t *mat);
1008
1009/*----------------------------------------------------------------------------*/
1015/*----------------------------------------------------------------------------*/
1016
1017void
1018cs_sdm_square_asymm(cs_sdm_t *mat);
1019
1020/*----------------------------------------------------------------------------*/
1026/*----------------------------------------------------------------------------*/
1027
1028void
1029cs_sdm_block_square_asymm(cs_sdm_t *mat);
1030
1031/*----------------------------------------------------------------------------*/
1048/*----------------------------------------------------------------------------*/
1049
1050void
1051cs_sdm_33_sym_qr_compute(const cs_real_t m[9],
1052 cs_real_t Qt[9],
1053 cs_real_t R[6]);
1054
1055/*----------------------------------------------------------------------------*/
1066/*----------------------------------------------------------------------------*/
1067
1068void
1069cs_sdm_33_lu_compute(const cs_sdm_t *m,
1070 cs_real_t facto[9]);
1071
1072/*----------------------------------------------------------------------------*/
1085/*----------------------------------------------------------------------------*/
1086
1087void
1088cs_sdm_lu_compute(const cs_sdm_t *m,
1089 cs_real_t facto[]);
1090
1091/*----------------------------------------------------------------------------*/
1105/*----------------------------------------------------------------------------*/
1106
1107void
1108cs_sdm_33_lu_solve(const cs_real_t facto[9],
1109 const cs_real_t rhs[3],
1110 cs_real_t sol[3]);
1111
1112/*----------------------------------------------------------------------------*/
1127/*----------------------------------------------------------------------------*/
1128
1129void
1130cs_sdm_lu_solve(cs_lnum_t n_rows,
1131 const cs_real_t facto[],
1132 const cs_real_t *rhs,
1133 cs_real_t *sol);
1134
1135/*----------------------------------------------------------------------------*/
1149/*----------------------------------------------------------------------------*/
1150
1151void
1152cs_sdm_33_ldlt_compute(const cs_sdm_t *m,
1153 cs_real_t facto[6]);
1154
1155/*----------------------------------------------------------------------------*/
1169/*----------------------------------------------------------------------------*/
1170
1171void
1172cs_sdm_44_ldlt_compute(const cs_sdm_t *m,
1173 cs_real_t facto[10]);
1174
1175/*----------------------------------------------------------------------------*/
1189/*----------------------------------------------------------------------------*/
1190
1191void
1192cs_sdm_66_ldlt_compute(const cs_sdm_t *m,
1193 cs_real_t facto[21]);
1194
1195/*----------------------------------------------------------------------------*/
1210/*----------------------------------------------------------------------------*/
1211
1212void
1213cs_sdm_ldlt_compute(const cs_sdm_t *m,
1214 cs_real_t *facto,
1215 cs_real_t *dkk);
1216
1217/*----------------------------------------------------------------------------*/
1227/*----------------------------------------------------------------------------*/
1228
1229void
1230cs_sdm_33_ldlt_solve(const cs_real_t facto[6],
1231 const cs_real_t rhs[3],
1232 cs_real_t sol[3]);
1233
1234/*----------------------------------------------------------------------------*/
1244/*----------------------------------------------------------------------------*/
1245
1246void
1247cs_sdm_44_ldlt_solve(const cs_real_t facto[10],
1248 const cs_real_t rhs[4],
1249 cs_real_t x[4]);
1250
1251/*----------------------------------------------------------------------------*/
1261/*----------------------------------------------------------------------------*/
1262
1263void
1264cs_sdm_66_ldlt_solve(const cs_real_t f[21],
1265 const cs_real_t b[6],
1266 cs_real_t x[6]);
1267
1268/*----------------------------------------------------------------------------*/
1280/*----------------------------------------------------------------------------*/
1281
1282void
1283cs_sdm_ldlt_solve(int n_rows,
1284 const cs_real_t *facto,
1285 const cs_real_t *rhs,
1286 cs_real_t *sol);
1287
1288/*----------------------------------------------------------------------------*/
1298/*----------------------------------------------------------------------------*/
1299
1300double
1301cs_sdm_test_symmetry(const cs_sdm_t *mat);
1302
1303/*----------------------------------------------------------------------------*/
1310/*----------------------------------------------------------------------------*/
1311
1312void
1313cs_sdm_block_dump(cs_lnum_t parent_id,
1314 const cs_sdm_t *mat);
1315
1316/*----------------------------------------------------------------------------*/
1329/*----------------------------------------------------------------------------*/
1330
1331void
1332cs_sdm_block_fprintf(FILE *fp,
1333 const char *fname,
1334 cs_real_t thd,
1335 const cs_sdm_t *m);
1336
1337#endif /* CS_SDM_H */
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
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.cpp:1546
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:766
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:1700
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:911
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:1850
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:1497
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:430
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:288
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:1880
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:133
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:714
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:363
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:1353
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:506
cs_sdm_free
cs_sdm_t * cs_sdm_free(cs_sdm_t *mat)
Free a cs_sdm_t structure.
Definition: cs_sdm.cpp:331
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:948
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:819
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:2113
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:1254
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:1064
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:1606
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:401
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:683
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:1987
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:1450
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:152
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:1155
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:191
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:420
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:456
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:1916
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:1315
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:849
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:2043
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_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:995
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:759
CS_SDM_BY_BLOCK
#define CS_SDM_BY_BLOCK
Definition: cs_sdm.h:44
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:637
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:1189
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:1083
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:1025
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:1413
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:225
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:169
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:1823
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:878
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:1126
cs
Definition: cs_algorithm.h:51
cs::sdm_t
cs_sdm_t sdm_t
Definition: cs_sdm.h:379
_cs_sdm_t
Definition: cs_sdm.h:70
_cs_sdm_t::n_rows
int n_rows
Definition: cs_sdm.h:76
_cs_sdm_t::_cs_sdm_t
_cs_sdm_t(const int n_max_rows_)
Constructor for a square matrix.
Definition: cs_sdm.cpp:2189
_cs_sdm_t::n_cols
int n_cols
Definition: cs_sdm.h:80
_cs_sdm_t::operator()
cs_real_t & operator()(const int i, const int j)
Modifiy the value (i,j) of the matric.
Definition: cs_sdm.cpp:2253
_cs_sdm_t::row
cs_real_t * row(const int row_i) const
Return a pointer on row.
Definition: cs_sdm.h:154
_cs_sdm_t::map_array
void map_array(int n_max_rows_, int n_max_cols_, 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.cpp:2572
_cs_sdm_t::n_max_cols
int n_max_cols
Definition: cs_sdm.h:79
_cs_sdm_t::operator+=
_cs_sdm_t & operator+=(const _cs_sdm_t &add)
Add two small dense matrices: loc += add.
Definition: cs_sdm.cpp:2386
_cs_sdm_t::symmetrize_ur
void symmetrize_ur()
Set the lower left according to the upper right part in order to get a symmetric matrix.
Definition: cs_sdm.cpp:2348
_cs_sdm_t::get_block
cs_sdm_t * get_block(const int row_block_id, const int col_block_id) const
Get a specific block in a cs_sdm_t structure defined by block.
Definition: cs_sdm.h:262
_cs_sdm_t::init
void init(const int nrows, const int ncols)
Initialize matrix.
Definition: cs_sdm.cpp:2291
_cs_sdm_t::flag
cs_flag_t flag
Definition: cs_sdm.h:72
_cs_sdm_t::dot
void dot(const cs_real_t vec[], cs_real_t mv[]) const
Compute a dense matrix-vector product for a small square matrix mv has been previously allocated.
Definition: cs_sdm.cpp:2410
_cs_sdm_t::operator=
_cs_sdm_t & operator=(const _cs_sdm_t &)=delete
_cs_sdm_t::block_desc
cs_sdm_block_t * block_desc
Definition: cs_sdm.h:86
_cs_sdm_t::val
cs_real_t * val
Definition: cs_sdm.h:82
_cs_sdm_t::n_max_rows
int n_max_rows
Definition: cs_sdm.h:75
_cs_sdm_t::data
cs_real_t * data() const
Return a pointer on the values.
Definition: cs_sdm.h:223
_cs_sdm_t::dump
void dump() const
Dump a small dense matrix.
Definition: cs_sdm.cpp:2438
_cs_sdm_t::size
int size() const
Return the size of the matrix.
Definition: cs_sdm.h:138
_cs_sdm_t::operator*=
_cs_sdm_t & operator*=(const cs_real_t &scaling)
Multiply a matrix with the scaling factor given as parameter.
Definition: cs_sdm.cpp:2367
_cs_sdm_t::~_cs_sdm_t
~_cs_sdm_t()
Destructor.
Definition: cs_sdm.cpp:2169
_cs_sdm_t::transpose
_cs_sdm_t transpose() const
Transpose a matrix.
Definition: cs_sdm.cpp:2326
cs_sdm_block_t
Definition: cs_sdm.h:54
cs_sdm_block_t::n_max_blocks_by_row
int n_max_blocks_by_row
Definition: cs_sdm.h:56
cs_sdm_block_t::blocks
cs_sdm_t * blocks
Definition: cs_sdm.h:64
cs_sdm_block_t::n_row_blocks
int n_row_blocks
Definition: cs_sdm.h:57
cs_sdm_block_t::n_max_blocks_by_col
int n_max_blocks_by_col
Definition: cs_sdm.h:58
cs_sdm_block_t::n_col_blocks
int n_col_blocks
Definition: cs_sdm.h:59