48 #define CS_SDM_BY_BLOCK (1 << 0) 49 #define CS_SDM_SYMMETRIC (1 << 1) 50 #define CS_SDM_SHARED_VAL (1 << 2) 154 if (x == NULL || y == NULL)
157 for (
int i = 0; i < n; i++)
177 cs_sdm_scalvect(
int n,
182 if (x == NULL || y == NULL)
185 for (
int i = 0; i < n; i++)
203 cs_sdm_add_scalvect(
int n,
208 if (x == NULL || y == NULL)
211 for (
int i = 0; i < n; i++)
225 cs_sdm_symm_ur(cs_sdm_t *mat)
228 for (
int i = 1; i < mat->n_rows; i++) {
229 cs_real_t *m_i = mat->val + i*mat->n_rows;
230 for (
int j = 0; j < i; j++)
231 m_i[j] = mat->val[j*mat->n_rows + i];
309 int n_max_blocks_by_col,
310 const int max_row_block_sizes[],
311 const int max_col_block_sizes[]);
327 int n_max_blocks_by_col);
344 cs_sdm_map_array(
int n_max_rows,
349 assert(array != NULL && m != NULL);
355 m->block_desc = NULL;
391 memset(mat->val, 0, n_rows*n_cols*
sizeof(
cs_real_t));
405 cs_sdm_square_init(
int n_rows,
410 mat->n_rows = mat->n_cols =
n_rows;
411 memset(mat->val, 0, n_rows*n_rows*
sizeof(
cs_real_t));
431 const int row_block_sizes[],
432 const int col_block_sizes[]);
476 cs_sdm_copy(cs_sdm_t *recv,
477 const cs_sdm_t *send)
480 assert(recv->n_max_rows >= send->n_rows);
481 assert(recv->n_max_cols >= send->n_cols);
483 recv->flag = send->flag;
484 recv->n_rows = send->n_rows;
485 recv->n_cols = send->n_cols;
488 memcpy(recv->val, send->val,
sizeof(
cs_real_t)*send->n_rows*send->n_cols);
517 static inline cs_sdm_t *
518 cs_sdm_get_block(
const cs_sdm_t *
const m,
528 return m->block_desc->blocks
529 + row_block_id*m->block_desc->n_col_blocks + col_block_id;
543 cs_sdm_get_col(
const cs_sdm_t *m,
548 assert(m != NULL && col_vals != NULL);
549 assert(col_id < m->n_cols);
552 for(
int i = 0; i < m->n_rows; i++, _col += m->n_cols)
572 cs_sdm_copy_block(
const cs_sdm_t *m,
573 const short int r_id,
574 const short int c_id,
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);
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);
607 for (
int i = 0; i < m->n_rows*m->n_cols; i++)
608 m->val[i] *= scaling;
623 cs_sdm_transpose_and_update(
const cs_sdm_t *m,
626 assert(m != NULL && mt != NULL);
627 assert(m->n_rows == mt->n_cols && m->n_cols == mt->n_rows);
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];
665 cs_sdm_33_matvec(
const cs_sdm_t *mat,
670 assert(mat != NULL && vec != NULL && mv != NULL);
671 assert(mat->n_rows == 3 && mat->n_cols == 3);
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];
693 cs_sdm_multiply_r1c3_rowrow(
const cs_sdm_t *a,
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);
703 c->val[0] += a->val[0]*b->val[0] + a->val[1]*b->val[1] + a->val[2]*b->val[2];
862 const cs_sdm_t *add);
877 const cs_sdm_t *add);
907 const cs_sdm_t *add);
922 const cs_sdm_t *add);
1270 const cs_sdm_t *mat);
1304 const cs_sdm_t *mat);
cs_flag_t flag
Definition: cs_sdm.h:76
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:1428
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:1571
cs_sdm_t * cs_sdm_free(cs_sdm_t *mat)
Free a cs_sdm_t structure.
Definition: cs_sdm.c:331
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:1946
cs_sdm_t * cs_sdm_create_transpose(cs_sdm_t *mat)
Define a new matrix which is its transpose.
Definition: cs_sdm.c:196
double precision, dimension(ncharm), save alpha
Definition: cpincl.f90:99
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:1149
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:1765
int n_max_cols
Definition: cs_sdm.h:83
cs_sdm_t * cs_sdm_create(cs_flag_t flag, int n_max_rows, int n_max_cols)
Allocate and initialize a cs_sdm_t structure Most generic function to create a cs_sdm_t structure...
Definition: cs_sdm.c:138
cs_sdm_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
void cs_sdm_block_square_asymm(cs_sdm_t *mat)
Set the given block matrix into its anti-symmetric part.
Definition: cs_sdm.c:1261
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
#define BEGIN_C_DECLS
Definition: cs_defs.h:510
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
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:2247
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:1390
void cs_sdm_block_33_to_xyz(const cs_sdm_t *mb33, cs_sdm_t *mbxyz)
Convert a matrix with each entry is a 3x3 block into a matrix with a block for each component x...
Definition: cs_sdm.c:457
int n_col_blocks
Definition: cs_sdm.h:63
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:907
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:1006
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:1887
int n_max_blocks_by_row
Definition: cs_sdm.h:60
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:1194
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:1122
double cs_real_t
Floating-point value.
Definition: cs_defs.h:322
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:2198
#define CS_SDM_BY_BLOCK
Definition: cs_sdm.h:48
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:772
#define CS_SDM_SHARED_VAL
Definition: cs_sdm.h:50
int n_cols
Definition: cs_sdm.h:84
double precision, save a
Definition: cs_fuel_incl.f90:146
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:506
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:1618
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:2317
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_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:2141
int n_row_blocks
Definition: cs_sdm.h:61
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:1675
cs_sdm_block_t * block_desc
Definition: cs_sdm.h:89
int n_max_blocks_by_col
Definition: cs_sdm.h:62
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:1983
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:2056
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:316
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:1047
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
void cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a, const cs_sdm_t *b, cs_sdm_t *c)
Compute a row-row matrix product of a and b. It is basically equal to the classical a*b^T...
Definition: cs_sdm.c:668
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:716
#define END_C_DECLS
Definition: cs_defs.h:511
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:622
unsigned short int cs_flag_t
Definition: cs_defs.h:324
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:1328
void cs_sdm_add(cs_sdm_t *mat, const cs_sdm_t *add)
Add two small dense matrices: loc += add.
Definition: cs_sdm.c:1099
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. mv has been previously allocated. mv is updated inside this function. Don't forget to initialize mv if needed.
Definition: cs_sdm.c:940
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:1523
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:865
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
double precision, dimension(:,:,:), allocatable nc
Definition: atimbr.f90:110
cs_real_t * val
Definition: cs_sdm.h:86
void() cs_sdm_matvec_t(const cs_sdm_t *mat, const cs_real_t *vec, cs_real_t *mv)
Generic prototype for computing a dense matrix-vector product mv has been previously allocated...
Definition: cs_sdm.h:125
cs_sdm_t * cs_sdm_block33_create(int n_max_blocks_by_row, int n_max_blocks_by_col)
Allocate and initialize a cs_sdm_t structure by block when the block size is constant and equal to 3...
Definition: cs_sdm.c:291
cs_sdm_t * blocks
Definition: cs_sdm.h:68
void cs_sdm_square_asymm(cs_sdm_t *mat)
Set the given matrix into its anti-symmetric part.
Definition: cs_sdm.c:1225
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:1915
int n_max_rows
Definition: cs_sdm.h:79
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:579
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:829
void cs_sdm_simple_dump(const cs_sdm_t *mat)
Dump a small dense matrix.
Definition: cs_sdm.c:2112
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:970
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:1487
int n_rows
Definition: cs_sdm.h:80
double precision, save b
Definition: cs_fuel_incl.f90:146