My Project
programmer's documentation
Loading...
Searching...
No Matches
cs_sdm.h
Go to the documentation of this file.
1#ifndef __CS_SDM_H__
2#define __CS_SDM_H__
3
4/*============================================================================
5 * Set of operations to handle Small Dense Matrices (SDM)
6 *============================================================================*/
7
8/*
9 This file is part of Code_Saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-2019 EDF S.A.
12
13 This program is free software; you can redistribute it and/or modify it under
14 the terms of the GNU General Public License as published by the Free Software
15 Foundation; either version 2 of the License, or (at your option) any later
16 version.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21 details.
22
23 You should have received a copy of the GNU General Public License along with
24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25 Street, Fifth Floor, Boston, MA 02110-1301, USA.
26*/
27
28/*----------------------------------------------------------------------------
29 * Standard C library headers
30 *----------------------------------------------------------------------------*/
31
32#include <stdio.h>
33#include <string.h>
34#include <math.h>
35
36/*----------------------------------------------------------------------------
37 * Local headers
38 *----------------------------------------------------------------------------*/
39
40/*----------------------------------------------------------------------------*/
41
43
44/*============================================================================
45 * Macro definitions
46 *============================================================================*/
47
48#define CS_SDM_BY_BLOCK (1 << 0) /* Matrix is defined by block */
49#define CS_SDM_SYMMETRIC (1 << 1) /* Matrix is symmetric by construction */
50#define CS_SDM_SHARED_VAL (1 << 2) /* Matrix is not owner of its values */
51
52/*============================================================================
53 * Type definitions
54 *============================================================================*/
55
56typedef struct _cs_sdm_t cs_sdm_t;
57
58typedef struct {
59
64
65 /* Allocated to n_max_blocks_by_row*n_max_blocks_by_col
66 Cast in cs_sdm_t where values are shared with the master cs_sdm_t struct.
67 */
68 cs_sdm_t *blocks;
69
71
72/* Structure enabling the repeated usage of Small Dense Matrices (SDM) during
73 the building of the linear system by a cellwise process */
74struct _cs_sdm_t {
75
76 cs_flag_t flag; /* Metadata */
77
78 /* Row-related members */
79 int n_max_rows; // max number of entities by primal cells
80 int n_rows; // current number of entities
81
82 /* Column-related members. Not useful if the matrix is square */
83 int n_max_cols; // max number of entries in a column
84 int n_cols; // current number of columns
85
86 cs_real_t *val; // small dense matrix (size: n_max_rows*n_max_cols)
87
88 /* Structure describing the matrix in terms of blocks */
90
91};
92
93/*============================================================================
94 * Prototypes for pointer of functions
95 *============================================================================*/
96
97/*----------------------------------------------------------------------------*/
106/*----------------------------------------------------------------------------*/
107
108typedef void
109(cs_sdm_product_t) (const cs_sdm_t *a,
110 const cs_sdm_t *b,
111 cs_sdm_t *c);
112
113/*----------------------------------------------------------------------------*/
122/*----------------------------------------------------------------------------*/
123
124typedef void
125(cs_sdm_matvec_t) (const cs_sdm_t *mat,
126 const cs_real_t *vec,
127 cs_real_t *mv);
128
129/*============================================================================
130 * Public function prototypes
131 *============================================================================*/
132
133/*----------------------------------------------------------------------------*/
145/*----------------------------------------------------------------------------*/
146
147static inline cs_real_t
149 const cs_real_t x[],
150 const cs_real_t y[])
151{
152 cs_real_t dp = 0;
153
154 if (x == NULL || y == NULL)
155 return dp;
156
157 for (int i = 0; i < n; i++)
158 dp += x[i]*y[i];
159
160 return dp;
161}
162
163/*----------------------------------------------------------------------------*/
174/*----------------------------------------------------------------------------*/
175
176static inline void
178 const cs_real_t s,
179 const cs_real_t x[],
180 cs_real_t y[])
181{
182 if (x == NULL || y == NULL)
183 return;
184
185 for (int i = 0; i < n; i++)
186 y[i] = s * x[i];
187}
188
189/*----------------------------------------------------------------------------*/
200/*----------------------------------------------------------------------------*/
201
202static inline void
204 const cs_real_t s,
205 const cs_real_t x[],
206 cs_real_t y[])
207{
208 if (x == NULL || y == NULL)
209 return;
210
211 for (int i = 0; i < n; i++)
212 y[i] += s * x[i];
213}
214
215/*----------------------------------------------------------------------------*/
226/*----------------------------------------------------------------------------*/
227
228cs_sdm_t *
230 int n_max_rows,
231 int n_max_cols);
232
233/*----------------------------------------------------------------------------*/
242/*----------------------------------------------------------------------------*/
243
244cs_sdm_t *
245cs_sdm_square_create(int n_max_rows);
246
247/*----------------------------------------------------------------------------*/
256/*----------------------------------------------------------------------------*/
257
258cs_sdm_t *
259cs_sdm_create_copy(const cs_sdm_t *m);
260
261/*----------------------------------------------------------------------------*/
269/*----------------------------------------------------------------------------*/
270
271cs_sdm_t *
272cs_sdm_create_transpose(cs_sdm_t *mat);
273
274/*----------------------------------------------------------------------------*/
285/*----------------------------------------------------------------------------*/
286
287cs_sdm_t *
288cs_sdm_block_create(int n_max_blocks_by_row,
289 int n_max_blocks_by_col,
290 const int max_row_block_sizes[],
291 const int max_col_block_sizes[]);
292
293/*----------------------------------------------------------------------------*/
303/*----------------------------------------------------------------------------*/
304
305cs_sdm_t *
306cs_sdm_block33_create(int n_max_blocks_by_row,
307 int n_max_blocks_by_col);
308
309/*----------------------------------------------------------------------------*/
321/*----------------------------------------------------------------------------*/
322
323static inline void
324cs_sdm_map_array(int n_max_rows,
325 int n_max_cols,
326 cs_sdm_t *m,
327 cs_real_t *array)
328{
329 assert(array != NULL && m != NULL); /* Sanity check */
330
331 m->flag = CS_SDM_SHARED_VAL;
332 m->n_rows = m->n_max_rows = n_max_rows;
333 m->n_cols = m->n_max_cols = n_max_cols;
334 m->val = array;
335 m->block_desc = NULL;
336}
337
338/*----------------------------------------------------------------------------*/
346/*----------------------------------------------------------------------------*/
347
348cs_sdm_t *
349cs_sdm_free(cs_sdm_t *mat);
350
351/*----------------------------------------------------------------------------*/
360/*----------------------------------------------------------------------------*/
361
362static inline void
363cs_sdm_init(int n_rows,
364 int n_cols,
365 cs_sdm_t *mat)
366{
367 assert(mat != NULL);
368
369 mat->n_rows = n_rows;
370 mat->n_cols = n_cols;
371 memset(mat->val, 0, n_rows*n_cols*sizeof(cs_real_t));
372}
373
374/*----------------------------------------------------------------------------*/
382/*----------------------------------------------------------------------------*/
383
384static inline void
386 cs_sdm_t *mat)
387{
388 assert(mat != NULL);
389
390 mat->n_rows = mat->n_cols = n_rows; /* square matrix */
391 memset(mat->val, 0, n_rows*n_rows*sizeof(cs_real_t));
392}
393
394/*----------------------------------------------------------------------------*/
405/*----------------------------------------------------------------------------*/
406
407void
408cs_sdm_block_init(cs_sdm_t *m,
409 int n_row_blocks,
410 int n_col_blocks,
411 const int row_block_sizes[],
412 const int col_block_sizes[]);
413
414/*----------------------------------------------------------------------------*/
423/*----------------------------------------------------------------------------*/
424
425void
426cs_sdm_block33_init(cs_sdm_t *m,
427 int n_row_blocks,
428 int n_col_blocks);
429
430/*----------------------------------------------------------------------------*/
438/*----------------------------------------------------------------------------*/
439
440static inline void
441cs_sdm_copy(cs_sdm_t *recv,
442 const cs_sdm_t *send)
443{
444 /* Sanity check */
445 assert(recv->n_max_rows >= send->n_rows);
446 assert(recv->n_max_cols >= send->n_cols);
447
448 recv->flag = send->flag;
449 recv->n_rows = send->n_rows;
450 recv->n_cols = send->n_cols;
451
452 /* Copy values */
453 memcpy(recv->val, send->val, sizeof(cs_real_t)*send->n_rows*send->n_cols);
454}
455
456/*----------------------------------------------------------------------------*/
465/*----------------------------------------------------------------------------*/
466
467cs_sdm_t *
468cs_sdm_block_create_copy(const cs_sdm_t *mref);
469
470/*----------------------------------------------------------------------------*/
480/*----------------------------------------------------------------------------*/
481
482static inline cs_sdm_t *
483cs_sdm_get_block(const cs_sdm_t *const m,
484 int row_block_id,
485 int col_block_id)
486{
487 /* Sanity checks */
488 assert(m != NULL);
489 assert(m->flag & CS_SDM_BY_BLOCK && m->block_desc != NULL);
490 assert(col_block_id < m->block_desc->n_col_blocks);
491 assert(row_block_id < m->block_desc->n_row_blocks);
492
493 return m->block_desc->blocks
494 + row_block_id*m->block_desc->n_col_blocks + col_block_id;
495}
496
497/*----------------------------------------------------------------------------*/
505/*----------------------------------------------------------------------------*/
506
507static inline void
508cs_sdm_get_col(const cs_sdm_t *m,
509 int col_id,
510 cs_real_t *col_vals)
511{
512 /* Sanity checks */
513 assert(m != NULL && col_vals != NULL);
514 assert(col_id < m->n_cols);
515
516 const cs_real_t *_col = m->val + col_id;
517 for(int i = 0; i < m->n_rows; i++, _col += m->n_cols)
518 col_vals[i] = *_col;
519}
520
521/*----------------------------------------------------------------------------*/
534/*----------------------------------------------------------------------------*/
535
536static inline void
537cs_sdm_copy_block(const cs_sdm_t *m,
538 const short int r_id,
539 const short int c_id,
540 const short int nr,
541 const short int nc,
542 cs_sdm_t *b)
543{
544 /* Sanity checks */
545 assert(m != NULL && b != NULL);
546 assert(r_id >= 0 && c_id >= 0);
547 assert((r_id + nr) <= m->n_rows);
548 assert((c_id + nc) <= m->n_cols);
549 assert(nr == b->n_rows && nc == b->n_cols);
550
551 const cs_real_t *_start = m->val + c_id + r_id*m->n_cols;
552 for (short int i = 0; i < nr; i++, _start += m->n_cols)
553 memcpy(b->val + i*nc, _start, sizeof(cs_real_t)*nc);
554}
555
556/*----------------------------------------------------------------------------*/
565/*----------------------------------------------------------------------------*/
566
567static inline void
569 cs_sdm_t *mt)
570{
571 assert(m != NULL && mt != NULL);
572 assert(m->n_rows == mt->n_cols && m->n_cols == mt->n_rows);
573
574 for (short int i = 0; i < m->n_rows; i++) {
575 const cs_real_t *m_i = m->val + i*m->n_cols;
576 for (short int j = 0; j < m->n_cols; j++)
577 mt->val[j*mt->n_cols + i] += m_i[j];
578 }
579}
580
581/*----------------------------------------------------------------------------*/
591/*----------------------------------------------------------------------------*/
592
593void
594cs_sdm_multiply(const cs_sdm_t *a,
595 const cs_sdm_t *b,
596 cs_sdm_t *c);
597
598/*----------------------------------------------------------------------------*/
610/*----------------------------------------------------------------------------*/
611
612static inline void
614 const cs_sdm_t *b,
615 cs_sdm_t *c)
616{
617 /* Sanity check */
618 assert(a != NULL && b != NULL && c != NULL);
619 assert(a->n_cols == 3 && b->n_cols == 3 &&
620 a->n_rows == 1 && c->n_rows == 1 &&
621 c->n_cols == 1 && b->n_rows == 1);
622
623 c->val[0] += a->val[0]*b->val[0] + a->val[1]*b->val[1] + a->val[2]*b->val[2];
624}
625
626/*----------------------------------------------------------------------------*/
638/*----------------------------------------------------------------------------*/
639
640void
641cs_sdm_multiply_rowrow(const cs_sdm_t *a,
642 const cs_sdm_t *b,
643 cs_sdm_t *c);
644
645/*----------------------------------------------------------------------------*/
658/*----------------------------------------------------------------------------*/
659
660void
661cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a,
662 const cs_sdm_t *b,
663 cs_sdm_t *c);
664
665/*----------------------------------------------------------------------------*/
677/*----------------------------------------------------------------------------*/
678
679void
680cs_sdm_block_multiply_rowrow(const cs_sdm_t *a,
681 const cs_sdm_t *b,
682 cs_sdm_t *c);
683
684/*----------------------------------------------------------------------------*/
697/*----------------------------------------------------------------------------*/
698
699void
700cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a,
701 const cs_sdm_t *b,
702 cs_sdm_t *c);
703
704/*----------------------------------------------------------------------------*/
713/*----------------------------------------------------------------------------*/
714
715void
716cs_sdm_square_matvec(const cs_sdm_t *mat,
717 const cs_real_t *vec,
718 cs_real_t *mv);
719
720/*----------------------------------------------------------------------------*/
729/*----------------------------------------------------------------------------*/
730
731void
732cs_sdm_matvec(const cs_sdm_t *mat,
733 const cs_real_t *vec,
734 cs_real_t *mv);
735
736/*----------------------------------------------------------------------------*/
746/*----------------------------------------------------------------------------*/
747
748void
749cs_sdm_update_matvec(const cs_sdm_t *mat,
750 const cs_real_t *vec,
751 cs_real_t *mv);
752
753/*----------------------------------------------------------------------------*/
764/*----------------------------------------------------------------------------*/
765
766void
767cs_sdm_matvec_transposed(const cs_sdm_t *mat,
768 const cs_real_t *vec,
769 cs_real_t *mv);
770
771/*----------------------------------------------------------------------------*/
778/*----------------------------------------------------------------------------*/
779
780void
781cs_sdm_block_add(cs_sdm_t *mat,
782 const cs_sdm_t *add);
783
784/*----------------------------------------------------------------------------*/
792/*----------------------------------------------------------------------------*/
793
794void
795cs_sdm_block_add_mult(cs_sdm_t *mat,
796 cs_real_t mult_coef,
797 const cs_sdm_t *add);
798
799/*----------------------------------------------------------------------------*/
809/*----------------------------------------------------------------------------*/
810
811void
812cs_sdm_block_matvec(const cs_sdm_t *mat,
813 const cs_real_t *vec,
814 cs_real_t *mv);
815
816/*----------------------------------------------------------------------------*/
823/*----------------------------------------------------------------------------*/
824
825void
826cs_sdm_add(cs_sdm_t *mat,
827 const cs_sdm_t *add);
828
829/*----------------------------------------------------------------------------*/
837/*----------------------------------------------------------------------------*/
838
839void
840cs_sdm_add_mult(cs_sdm_t *mat,
841 cs_real_t alpha,
842 const cs_sdm_t *add);
843
844/*----------------------------------------------------------------------------*/
852/*----------------------------------------------------------------------------*/
853
854void
855cs_sdm_square_add_transpose(cs_sdm_t *mat,
856 cs_sdm_t *tr);
857
858/*----------------------------------------------------------------------------*/
865/*----------------------------------------------------------------------------*/
866
867void
868cs_sdm_square_2symm(cs_sdm_t *mat);
869
870/*----------------------------------------------------------------------------*/
876/*----------------------------------------------------------------------------*/
877
878void
879cs_sdm_square_asymm(cs_sdm_t *mat);
880
881/*----------------------------------------------------------------------------*/
887/*----------------------------------------------------------------------------*/
888
889void
890cs_sdm_block_square_asymm(cs_sdm_t *mat);
891
892/*----------------------------------------------------------------------------*/
909/*----------------------------------------------------------------------------*/
910
911void
913 cs_real_t Qt[9],
914 cs_real_t R[6]);
915
916/*----------------------------------------------------------------------------*/
927/*----------------------------------------------------------------------------*/
928
929void
930cs_sdm_33_lu_compute(const cs_sdm_t *m,
931 cs_real_t facto[9]);
932
933/*----------------------------------------------------------------------------*/
946/*----------------------------------------------------------------------------*/
947
948void
949cs_sdm_lu_compute(const cs_sdm_t *m,
950 cs_real_t facto[]);
951
952/*----------------------------------------------------------------------------*/
966/*----------------------------------------------------------------------------*/
967
968void
969cs_sdm_33_lu_solve(const cs_real_t facto[9],
970 const cs_real_t rhs[3],
971 cs_real_t sol[3]);
972
973/*----------------------------------------------------------------------------*/
988/*----------------------------------------------------------------------------*/
989
990void
992 const cs_real_t facto[],
993 const cs_real_t *rhs,
994 cs_real_t *sol);
995
996/*----------------------------------------------------------------------------*/
1010/*----------------------------------------------------------------------------*/
1011
1012void
1013cs_sdm_33_ldlt_compute(const cs_sdm_t *m,
1014 cs_real_t facto[6]);
1015
1016/*----------------------------------------------------------------------------*/
1030/*----------------------------------------------------------------------------*/
1031
1032void
1033cs_sdm_44_ldlt_compute(const cs_sdm_t *m,
1034 cs_real_t facto[10]);
1035
1036/*----------------------------------------------------------------------------*/
1050/*----------------------------------------------------------------------------*/
1051
1052void
1053cs_sdm_66_ldlt_compute(const cs_sdm_t *m,
1054 cs_real_t facto[21]);
1055
1056/*----------------------------------------------------------------------------*/
1071/*----------------------------------------------------------------------------*/
1072
1073void
1074cs_sdm_ldlt_compute(const cs_sdm_t *m,
1075 cs_real_t *facto,
1076 cs_real_t *dkk);
1077
1078/*----------------------------------------------------------------------------*/
1088/*----------------------------------------------------------------------------*/
1089
1090void
1091cs_sdm_33_ldlt_solve(const cs_real_t facto[6],
1092 const cs_real_t rhs[3],
1093 cs_real_t sol[3]);
1094
1095/*----------------------------------------------------------------------------*/
1105/*----------------------------------------------------------------------------*/
1106
1107void
1108cs_sdm_44_ldlt_solve(const cs_real_t facto[10],
1109 const cs_real_t rhs[4],
1110 cs_real_t x[4]);
1111
1112/*----------------------------------------------------------------------------*/
1122/*----------------------------------------------------------------------------*/
1123
1124void
1126 const cs_real_t b[6],
1127 cs_real_t x[6]);
1128
1129/*----------------------------------------------------------------------------*/
1141/*----------------------------------------------------------------------------*/
1142
1143void
1144cs_sdm_ldlt_solve(int n_rows,
1145 const cs_real_t *facto,
1146 const cs_real_t *rhs,
1147 cs_real_t *sol);
1148
1149/*----------------------------------------------------------------------------*/
1159/*----------------------------------------------------------------------------*/
1160
1161double
1162cs_sdm_test_symmetry(const cs_sdm_t *mat);
1163
1164/*----------------------------------------------------------------------------*/
1170/*----------------------------------------------------------------------------*/
1171
1172void
1173cs_sdm_simple_dump(const cs_sdm_t *mat);
1174
1175/*----------------------------------------------------------------------------*/
1184/*----------------------------------------------------------------------------*/
1185
1186void
1187cs_sdm_dump(cs_lnum_t parent_id,
1188 const cs_lnum_t *row_ids,
1189 const cs_lnum_t *col_ids,
1190 const cs_sdm_t *mat);
1191
1192/*----------------------------------------------------------------------------*/
1205/*----------------------------------------------------------------------------*/
1206
1207void
1208cs_sdm_fprintf(FILE *fp,
1209 const char *fname,
1210 cs_real_t thd,
1211 const cs_sdm_t *m);
1212
1213/*----------------------------------------------------------------------------*/
1220/*----------------------------------------------------------------------------*/
1221
1222void
1224 const cs_sdm_t *mat);
1225
1226/*----------------------------------------------------------------------------*/
1239/*----------------------------------------------------------------------------*/
1240
1241void
1242cs_sdm_block_fprintf(FILE *fp,
1243 const char *fname,
1244 cs_real_t thd,
1245 const cs_sdm_t *m);
1246
1247/*----------------------------------------------------------------------------*/
1248
1250
1251#endif /* __CS_SDM_H__ */
#define BEGIN_C_DECLS
Definition cs_defs.h:467
double cs_real_t
Floating-point value.
Definition cs_defs.h:302
#define END_C_DECLS
Definition cs_defs.h:468
int cs_lnum_t
local mesh entity id
Definition cs_defs.h:298
unsigned short int cs_flag_t
Definition cs_defs.h:304
void cs_sdm_44_ldlt_compute(const cs_sdm_t *m, cs_real_t facto[10])
LDL^T: Modified Cholesky decomposition of a 4x4 SPD matrix. For more reference, see for instance http...
Definition cs_sdm.c:1570
void cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T....
Definition cs_sdm.c:724
void cs_sdm_ldlt_compute(const cs_sdm_t *m, cs_real_t *facto, cs_real_t *dkk)
LDL^T: Modified Cholesky decomposition of a SPD matrix. For more reference, see for instance http://m...
Definition cs_sdm.c:1717
void cs_sdm_block_add_mult(cs_sdm_t *mat, cs_real_t mult_coef, const cs_sdm_t *add)
Add two matrices defined by block: loc += mult_coef * add.
Definition cs_sdm.c:958
void cs_sdm_44_ldlt_solve(const cs_real_t facto[10], const cs_real_t rhs[4], cs_real_t x[4])
Solve a 4x4 matrix with a modified Cholesky decomposition (L.D.L^T) The solution should be already al...
Definition cs_sdm.c:1867
void cs_sdm_33_ldlt_compute(const cs_sdm_t *m, cs_real_t facto[6])
LDL^T: Modified Cholesky decomposition of a 3x3 SPD matrix. For more reference, see for instance http...
Definition cs_sdm.c:1523
static void cs_sdm_scalvect(int n, const cs_real_t s, const cs_real_t x[], cs_real_t y[])
Multiply a small vector by a scalar coefficient: y = s x For very small array sizes (3,...
Definition cs_sdm.h:177
cs_sdm_t * cs_sdm_block33_create(int n_max_blocks_by_row, int n_max_blocks_by_col)
Allocate and initialize a cs_sdm_t structure by block when the block size is constant and equal to 3.
Definition cs_sdm.c:291
void cs_sdm_66_ldlt_solve(const cs_real_t f[21], const cs_real_t b[6], cs_real_t x[6])
Solve a 6x6 matrix with a modified Cholesky decomposition (L.D.L^T) The solution should be already al...
Definition cs_sdm.c:1898
cs_sdm_t * cs_sdm_create(cs_flag_t flag, int n_max_rows, int n_max_cols)
Allocate and initialize a cs_sdm_t structure Most generic function to create a cs_sdm_t structure.
Definition cs_sdm.c:138
static void cs_sdm_transpose_and_update(const cs_sdm_t *m, cs_sdm_t *mt)
transpose and copy a matrix into another one already shaped sub-matrix starting from (r_id,...
Definition cs_sdm.h:568
void cs_sdm_block_init(cs_sdm_t *m, int n_row_blocks, int n_col_blocks, const int row_block_sizes[], const int col_block_sizes[])
Initialize the pattern of cs_sdm_t structure defined by block The matrix should have been allocated b...
Definition cs_sdm.c:363
static cs_sdm_t * cs_sdm_get_block(const cs_sdm_t *const m, int row_block_id, int col_block_id)
Get a specific block in a cs_sdm_t structure defined by block.
Definition cs_sdm.h:483
void cs_sdm_lu_compute(const cs_sdm_t *m, cs_real_t facto[])
LU factorization of a small dense matrix. Small means that the number m->n_rows is less than 100 for ...
Definition cs_sdm.c:1380
cs_sdm_t * cs_sdm_block_create_copy(const cs_sdm_t *mref)
Allocate and initialize a cs_sdm_t structure w.r.t. to a given matrix.
Definition cs_sdm.c:458
cs_sdm_t * cs_sdm_free(cs_sdm_t *mat)
Free a cs_sdm_t structure.
Definition cs_sdm.c:331
void cs_sdm_block_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix defined by block mv has been previousl...
Definition cs_sdm.c:999
static void cs_sdm_map_array(int n_max_rows, int n_max_cols, cs_sdm_t *m, cs_real_t *array)
Map an array into a predefined cs_sdm_t structure. This array is shared and the lifecycle of this arr...
Definition cs_sdm.h:324
void() cs_sdm_matvec_t(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Generic prototype for computing a dense matrix-vector product mv has been previously allocated.
Definition cs_sdm.h:125
void cs_sdm_update_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix mv has been previously allocated and i...
Definition cs_sdm.c:859
void cs_sdm_block_fprintf(FILE *fp, const char *fname, cs_real_t thd, const cs_sdm_t *m)
Print a cs_sdm_t structure which is defined by block Print into the file f if given otherwise open a ...
Definition cs_sdm.c:2264
void cs_sdm_33_sym_qr_compute(const cs_real_t m[9], cs_real_t Qt[9], cs_real_t R[6])
Decompose a matrix into the matrix product Q.R Case of a 3x3 symmetric matrix.
Definition cs_sdm.c:1280
void() cs_sdm_product_t(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Generic prototype for computing a local dense matrix-product c = a*b where c has been previously allo...
Definition cs_sdm.h:109
void cs_sdm_66_ldlt_compute(const cs_sdm_t *m, cs_real_t facto[21])
LDL^T: Modified Cholesky decomposition of a 6x6 SPD matrix. For more reference, see for instance http...
Definition cs_sdm.c:1627
void cs_sdm_add(cs_sdm_t *mat, const cs_sdm_t *add)
Add two small dense matrices: loc += add.
Definition cs_sdm.c:1051
static cs_real_t cs_sdm_dot(int n, const cs_real_t x[], const cs_real_t y[])
Basic dot product for a small vector For very small array sizes (3, 4, 6) prefer functions in cs_math...
Definition cs_sdm.h:148
static void cs_sdm_copy_block(const cs_sdm_t *m, const short int r_id, const short int c_id, const short int nr, const short int nc, cs_sdm_t *b)
copy a block of a matrix into a sub-matrix starting from (r_id, c_id) with a size of nr rows and nc c...
Definition cs_sdm.h:537
double cs_sdm_test_symmetry(const cs_sdm_t *mat)
Test if a matrix is symmetric. Return 0. if the extradiagonal differences are lower thann the machine...
Definition cs_sdm.c:2008
void cs_sdm_lu_solve(cs_lnum_t n_rows, const cs_real_t facto[], const cs_real_t *rhs, cs_real_t *sol)
Solve a system A.sol = rhs using a LU factorization of A (a small dense matrix).
Definition cs_sdm.c:1475
void cs_sdm_dump(cs_lnum_t parent_id, const cs_lnum_t *row_ids, const cs_lnum_t *col_ids, const cs_sdm_t *mat)
Dump a small dense matrix.
Definition cs_sdm.c:2093
cs_sdm_t * cs_sdm_square_create(int n_max_rows)
Allocate and initialize a cs_sdm_t structure Case of a square matrix.
Definition cs_sdm.c:157
void cs_sdm_square_asymm(cs_sdm_t *mat)
Set the given matrix into its anti-symmetric part.
Definition cs_sdm.c:1177
#define CS_SDM_SHARED_VAL
Definition cs_sdm.h:50
void cs_sdm_multiply(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a local dense matrix-product c = a*b c has been previously allocated.
Definition cs_sdm.c:531
cs_sdm_t * cs_sdm_create_transpose(cs_sdm_t *mat)
Define a new matrix which is its transpose.
Definition cs_sdm.c:196
void cs_sdm_block33_init(cs_sdm_t *m, int n_row_blocks, int n_col_blocks)
Initialize the pattern of cs_sdm_t structure defined by 3x3 block The matrix should have been allocat...
Definition cs_sdm.c:419
static void cs_sdm_add_scalvect(int n, const cs_real_t s, const cs_real_t x[], cs_real_t y[])
Multiply a small vector by a scalar coefficient: y += s x For very small array sizes (3,...
Definition cs_sdm.h:203
void cs_sdm_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix mv has been previously allocated.
Definition cs_sdm.c:817
void cs_sdm_ldlt_solve(int n_rows, const cs_real_t *facto, const cs_real_t *rhs, cs_real_t *sol)
Solve a SPD matrix with a L.D.L^T (Modified Cholesky decomposition) The solution should be already al...
Definition cs_sdm.c:1935
static void cs_sdm_get_col(const cs_sdm_t *m, int col_id, cs_real_t *col_vals)
Get a copy of a column in a preallocated vector.
Definition cs_sdm.h:508
void cs_sdm_block_multiply_rowrow(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T....
Definition cs_sdm.c:668
void cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T....
Definition cs_sdm.c:620
void cs_sdm_33_lu_compute(const cs_sdm_t *m, cs_real_t facto[9])
LU factorization of a small dense 3x3 matrix.
Definition cs_sdm.c:1342
void cs_sdm_matvec_transposed(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a rectangular matrix which is transposed....
Definition cs_sdm.c:892
void cs_sdm_fprintf(FILE *fp, const char *fname, cs_real_t thd, const cs_sdm_t *m)
Print a cs_sdm_t structure not defined by block Print into the file f if given otherwise open a new f...
Definition cs_sdm.c:2145
void cs_sdm_block_dump(cs_lnum_t parent_id, const cs_sdm_t *mat)
Dump a small dense matrix defined by blocks.
Definition cs_sdm.c:2194
void cs_sdm_multiply_rowrow(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T....
Definition cs_sdm.c:574
static void cs_sdm_square_init(int n_rows, cs_sdm_t *mat)
Initialize a cs_sdm_t structure Case of a square matrix.
Definition cs_sdm.h:385
void cs_sdm_add_mult(cs_sdm_t *mat, cs_real_t alpha, const cs_sdm_t *add)
Add two small dense matrices: loc += alpha*add.
Definition cs_sdm.c:1074
static void cs_sdm_multiply_r1c3_rowrow(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T....
Definition cs_sdm.h:613
void cs_sdm_simple_dump(const cs_sdm_t *mat)
Dump a small dense matrix.
Definition cs_sdm.c:2064
#define CS_SDM_BY_BLOCK
Definition cs_sdm.h:48
static void cs_sdm_copy(cs_sdm_t *recv, const cs_sdm_t *send)
Copy a cs_sdm_t structure into another cs_sdm_t structure which has been already allocated.
Definition cs_sdm.h:441
void cs_sdm_block_square_asymm(cs_sdm_t *mat)
Set the given block matrix into its anti-symmetric part.
Definition cs_sdm.c:1213
void cs_sdm_square_add_transpose(cs_sdm_t *mat, cs_sdm_t *tr)
Define a new matrix by adding the given matrix with its transpose. Keep the transposed matrix for a f...
Definition cs_sdm.c:1101
void cs_sdm_33_lu_solve(const cs_real_t facto[9], const cs_real_t rhs[3], cs_real_t sol[3])
Solve a system A.sol = rhs using a LU factorization of A (a small 3x3 dense matrix).
Definition cs_sdm.c:1439
cs_sdm_t * cs_sdm_block_create(int n_max_blocks_by_row, int n_max_blocks_by_col, const int max_row_block_sizes[], const int max_col_block_sizes[])
Allocate and initialize a cs_sdm_t structure.
Definition cs_sdm.c:231
cs_sdm_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
static void cs_sdm_init(int n_rows, int n_cols, cs_sdm_t *mat)
Initialize a cs_sdm_t structure Case of a square matrix.
Definition cs_sdm.h:363
void cs_sdm_33_ldlt_solve(const cs_real_t facto[6], const cs_real_t rhs[3], cs_real_t sol[3])
Solve a 3x3 matrix with a modified Cholesky decomposition (L.D.L^T) The solution should be already al...
Definition cs_sdm.c:1839
void cs_sdm_square_matvec(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Compute a dense matrix-vector product for a small square matrix mv has been previously allocated.
Definition cs_sdm.c:781
void cs_sdm_block_add(cs_sdm_t *mat, const cs_sdm_t *add)
Add two matrices defined by block: loc += add.
Definition cs_sdm.c:922
void cs_sdm_square_2symm(cs_sdm_t *mat)
Set the given matrix to two times its symmetric part mat --> mat + mat_tr = 2*symm(mat)
Definition cs_sdm.c:1146
Definition cs_sdm.h:74
int n_rows
Definition cs_sdm.h:80
int n_cols
Definition cs_sdm.h:84
int n_max_cols
Definition cs_sdm.h:83
cs_flag_t flag
Definition cs_sdm.h:76
cs_sdm_block_t * block_desc
Definition cs_sdm.h:89
cs_real_t * val
Definition cs_sdm.h:86
int n_max_rows
Definition cs_sdm.h:79
Definition cs_sdm.h:58
int n_max_blocks_by_row
Definition cs_sdm.h:60
cs_sdm_t * blocks
Definition cs_sdm.h:68
int n_row_blocks
Definition cs_sdm.h:61
int n_max_blocks_by_col
Definition cs_sdm.h:62
int n_col_blocks
Definition cs_sdm.h:63