7.1
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-2021 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 
56 typedef struct _cs_sdm_t cs_sdm_t;
57 
58 typedef 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 */
74 struct _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 
108 typedef 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 
124 typedef 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 
147 static inline cs_real_t
148 cs_sdm_dot(int n,
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 
176 static inline void
177 cs_sdm_scalvect(int n,
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 
202 static inline void
203 cs_sdm_add_scalvect(int n,
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 /*----------------------------------------------------------------------------*/
222 /*----------------------------------------------------------------------------*/
223 
224 static inline void
225 cs_sdm_symm_ur(cs_sdm_t *mat)
226 {
227  assert(mat != NULL);
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];
232  }
233 }
234 
235 /*----------------------------------------------------------------------------*/
246 /*----------------------------------------------------------------------------*/
247 
248 cs_sdm_t *
250  int n_max_rows,
251  int n_max_cols);
252 
253 /*----------------------------------------------------------------------------*/
262 /*----------------------------------------------------------------------------*/
263 
264 cs_sdm_t *
265 cs_sdm_square_create(int n_max_rows);
266 
267 /*----------------------------------------------------------------------------*/
276 /*----------------------------------------------------------------------------*/
277 
278 cs_sdm_t *
279 cs_sdm_create_copy(const cs_sdm_t *m);
280 
281 /*----------------------------------------------------------------------------*/
289 /*----------------------------------------------------------------------------*/
290 
291 cs_sdm_t *
292 cs_sdm_create_transpose(cs_sdm_t *mat);
293 
294 /*----------------------------------------------------------------------------*/
305 /*----------------------------------------------------------------------------*/
306 
307 cs_sdm_t *
308 cs_sdm_block_create(int n_max_blocks_by_row,
309  int n_max_blocks_by_col,
310  const int max_row_block_sizes[],
311  const int max_col_block_sizes[]);
312 
313 /*----------------------------------------------------------------------------*/
323 /*----------------------------------------------------------------------------*/
324 
325 cs_sdm_t *
326 cs_sdm_block33_create(int n_max_blocks_by_row,
327  int n_max_blocks_by_col);
328 
329 /*----------------------------------------------------------------------------*/
341 /*----------------------------------------------------------------------------*/
342 
343 static inline void
344 cs_sdm_map_array(int n_max_rows,
345  int n_max_cols,
346  cs_sdm_t *m,
347  cs_real_t *array)
348 {
349  assert(array != NULL && m != NULL); /* Sanity check */
350 
351  m->flag = CS_SDM_SHARED_VAL;
352  m->n_rows = m->n_max_rows = n_max_rows;
353  m->n_cols = m->n_max_cols = n_max_cols;
354  m->val = array;
355  m->block_desc = NULL;
356 }
357 
358 /*----------------------------------------------------------------------------*/
366 /*----------------------------------------------------------------------------*/
367 
368 cs_sdm_t *
369 cs_sdm_free(cs_sdm_t *mat);
370 
371 /*----------------------------------------------------------------------------*/
380 /*----------------------------------------------------------------------------*/
381 
382 static inline void
383 cs_sdm_init(int n_rows,
384  int n_cols,
385  cs_sdm_t *mat)
386 {
387  assert(mat != NULL);
388 
389  mat->n_rows = n_rows;
390  mat->n_cols = n_cols;
391  memset(mat->val, 0, n_rows*n_cols*sizeof(cs_real_t));
392 }
393 
394 /*----------------------------------------------------------------------------*/
402 /*----------------------------------------------------------------------------*/
403 
404 static inline void
405 cs_sdm_square_init(int n_rows,
406  cs_sdm_t *mat)
407 {
408  assert(mat != NULL);
409 
410  mat->n_rows = mat->n_cols = n_rows; /* square matrix */
411  memset(mat->val, 0, n_rows*n_rows*sizeof(cs_real_t));
412 }
413 
414 /*----------------------------------------------------------------------------*/
425 /*----------------------------------------------------------------------------*/
426 
427 void
428 cs_sdm_block_init(cs_sdm_t *m,
429  int n_row_blocks,
430  int n_col_blocks,
431  const int row_block_sizes[],
432  const int col_block_sizes[]);
433 
434 /*----------------------------------------------------------------------------*/
443 /*----------------------------------------------------------------------------*/
444 
445 void
446 cs_sdm_block33_init(cs_sdm_t *m,
447  int n_row_blocks,
448  int n_col_blocks);
449 
450 
451 /*----------------------------------------------------------------------------*/
459 /*----------------------------------------------------------------------------*/
460 
461 void
462 cs_sdm_block_33_to_xyz(const cs_sdm_t *mb33,
463  cs_sdm_t *mbxyz);
464 
465 /*----------------------------------------------------------------------------*/
473 /*----------------------------------------------------------------------------*/
474 
475 static inline void
476 cs_sdm_copy(cs_sdm_t *recv,
477  const cs_sdm_t *send)
478 {
479  /* Sanity check */
480  assert(recv->n_max_rows >= send->n_rows);
481  assert(recv->n_max_cols >= send->n_cols);
482 
483  recv->flag = send->flag;
484  recv->n_rows = send->n_rows;
485  recv->n_cols = send->n_cols;
486 
487  /* Copy values */
488  memcpy(recv->val, send->val, sizeof(cs_real_t)*send->n_rows*send->n_cols);
489 }
490 
491 /*----------------------------------------------------------------------------*/
500 /*----------------------------------------------------------------------------*/
501 
502 cs_sdm_t *
503 cs_sdm_block_create_copy(const cs_sdm_t *mref);
504 
505 /*----------------------------------------------------------------------------*/
515 /*----------------------------------------------------------------------------*/
516 
517 static inline cs_sdm_t *
518 cs_sdm_get_block(const cs_sdm_t *const m,
519  int row_block_id,
520  int col_block_id)
521 {
522  /* Sanity checks */
523  assert(m != NULL);
524  assert(m->flag & CS_SDM_BY_BLOCK && m->block_desc != NULL);
525  assert(col_block_id < m->block_desc->n_col_blocks);
526  assert(row_block_id < m->block_desc->n_row_blocks);
527 
528  return m->block_desc->blocks
529  + row_block_id*m->block_desc->n_col_blocks + col_block_id;
530 }
531 
532 /*----------------------------------------------------------------------------*/
540 /*----------------------------------------------------------------------------*/
541 
542 static inline void
543 cs_sdm_get_col(const cs_sdm_t *m,
544  int col_id,
545  cs_real_t *col_vals)
546 {
547  /* Sanity checks */
548  assert(m != NULL && col_vals != NULL);
549  assert(col_id < m->n_cols);
550 
551  const cs_real_t *_col = m->val + col_id;
552  for(int i = 0; i < m->n_rows; i++, _col += m->n_cols)
553  col_vals[i] = *_col;
554 }
555 
556 /*----------------------------------------------------------------------------*/
569 /*----------------------------------------------------------------------------*/
570 
571 static inline void
572 cs_sdm_copy_block(const cs_sdm_t *m,
573  const short int r_id,
574  const short int c_id,
575  const short int nr,
576  const short int nc,
577  cs_sdm_t *b)
578 {
579  /* Sanity checks */
580  assert(m != NULL && b != NULL);
581  assert(r_id >= 0 && c_id >= 0);
582  assert((r_id + nr) <= m->n_rows);
583  assert((c_id + nc) <= m->n_cols);
584  assert(nr == b->n_rows && nc == b->n_cols);
585 
586  const cs_real_t *_start = m->val + c_id + r_id*m->n_cols;
587  for (short int i = 0; i < nr; i++, _start += m->n_cols)
588  memcpy(b->val + i*nc, _start, sizeof(cs_real_t)*nc);
589 }
590 
591 /*----------------------------------------------------------------------------*/
598 /*----------------------------------------------------------------------------*/
599 
600 static inline void
601 cs_sdm_scale(cs_real_t scaling,
602  cs_sdm_t *m)
603 {
604  /* Sanity checks */
605  assert(m != NULL);
606 
607  for (int i = 0; i < m->n_rows*m->n_cols; i++)
608  m->val[i] *= scaling;
609 }
610 
611 /*----------------------------------------------------------------------------*/
620 /*----------------------------------------------------------------------------*/
621 
622 static inline void
623 cs_sdm_transpose_and_update(const cs_sdm_t *m,
624  cs_sdm_t *mt)
625 {
626  assert(m != NULL && mt != NULL);
627  assert(m->n_rows == mt->n_cols && m->n_cols == mt->n_rows);
628 
629  for (short int i = 0; i < m->n_rows; i++) {
630  const cs_real_t *m_i = m->val + i*m->n_cols;
631  for (short int j = 0; j < m->n_cols; j++)
632  mt->val[j*mt->n_cols + i] += m_i[j];
633  }
634 }
635 
636 /*----------------------------------------------------------------------------*/
646 /*----------------------------------------------------------------------------*/
647 
648 void
649 cs_sdm_multiply(const cs_sdm_t *a,
650  const cs_sdm_t *b,
651  cs_sdm_t *c);
652 
653 /*----------------------------------------------------------------------------*/
662 /*----------------------------------------------------------------------------*/
663 
664 static inline void
665 cs_sdm_33_matvec(const cs_sdm_t *mat,
666  const cs_real_t *vec,
667  cs_real_t *mv)
668 {
669  /* Sanity checks */
670  assert(mat != NULL && vec != NULL && mv != NULL);
671  assert(mat->n_rows == 3 && mat->n_cols == 3);
672 
673  mv[0] = vec[0]*mat->val[0] + vec[1]*mat->val[1] + vec[2]*mat->val[2];
674  mv[1] = vec[0]*mat->val[3] + vec[1]*mat->val[4] + vec[2]*mat->val[5];
675  mv[2] = vec[0]*mat->val[6] + vec[1]*mat->val[7] + vec[2]*mat->val[8];
676 }
677 
678 /*----------------------------------------------------------------------------*/
690 /*----------------------------------------------------------------------------*/
691 
692 static inline void
693 cs_sdm_multiply_r1c3_rowrow(const cs_sdm_t *a,
694  const cs_sdm_t *b,
695  cs_sdm_t *c)
696 {
697  /* Sanity check */
698  assert(a != NULL && b != NULL && c != NULL);
699  assert(a->n_cols == 3 && b->n_cols == 3 &&
700  a->n_rows == 1 && c->n_rows == 1 &&
701  c->n_cols == 1 && b->n_rows == 1);
702 
703  c->val[0] += a->val[0]*b->val[0] + a->val[1]*b->val[1] + a->val[2]*b->val[2];
704 }
705 
706 /*----------------------------------------------------------------------------*/
718 /*----------------------------------------------------------------------------*/
719 
720 void
721 cs_sdm_multiply_rowrow(const cs_sdm_t *a,
722  const cs_sdm_t *b,
723  cs_sdm_t *c);
724 
725 /*----------------------------------------------------------------------------*/
738 /*----------------------------------------------------------------------------*/
739 
740 void
741 cs_sdm_multiply_rowrow_sym(const cs_sdm_t *a,
742  const cs_sdm_t *b,
743  cs_sdm_t *c);
744 
745 /*----------------------------------------------------------------------------*/
757 /*----------------------------------------------------------------------------*/
758 
759 void
760 cs_sdm_block_multiply_rowrow(const cs_sdm_t *a,
761  const cs_sdm_t *b,
762  cs_sdm_t *c);
763 
764 /*----------------------------------------------------------------------------*/
777 /*----------------------------------------------------------------------------*/
778 
779 void
780 cs_sdm_block_multiply_rowrow_sym(const cs_sdm_t *a,
781  const cs_sdm_t *b,
782  cs_sdm_t *c);
783 
784 /*----------------------------------------------------------------------------*/
793 /*----------------------------------------------------------------------------*/
794 
795 void
796 cs_sdm_square_matvec(const cs_sdm_t *mat,
797  const cs_real_t *vec,
798  cs_real_t *mv);
799 
800 /*----------------------------------------------------------------------------*/
809 /*----------------------------------------------------------------------------*/
810 
811 void
812 cs_sdm_matvec(const cs_sdm_t *mat,
813  const cs_real_t *vec,
814  cs_real_t *mv);
815 
816 /*----------------------------------------------------------------------------*/
826 /*----------------------------------------------------------------------------*/
827 
828 void
829 cs_sdm_update_matvec(const cs_sdm_t *mat,
830  const cs_real_t *vec,
831  cs_real_t *mv);
832 
833 /*----------------------------------------------------------------------------*/
844 /*----------------------------------------------------------------------------*/
845 
846 void
847 cs_sdm_matvec_transposed(const cs_sdm_t *mat,
848  const cs_real_t *vec,
849  cs_real_t *mv);
850 
851 /*----------------------------------------------------------------------------*/
858 /*----------------------------------------------------------------------------*/
859 
860 void
861 cs_sdm_block_add(cs_sdm_t *mat,
862  const cs_sdm_t *add);
863 
864 /*----------------------------------------------------------------------------*/
872 /*----------------------------------------------------------------------------*/
873 
874 void
875 cs_sdm_block_add_mult(cs_sdm_t *mat,
876  cs_real_t mult_coef,
877  const cs_sdm_t *add);
878 
879 /*----------------------------------------------------------------------------*/
889 /*----------------------------------------------------------------------------*/
890 
891 void
892 cs_sdm_block_matvec(const cs_sdm_t *mat,
893  const cs_real_t *vec,
894  cs_real_t *mv);
895 
896 /*----------------------------------------------------------------------------*/
903 /*----------------------------------------------------------------------------*/
904 
905 void
906 cs_sdm_add(cs_sdm_t *mat,
907  const cs_sdm_t *add);
908 
909 /*----------------------------------------------------------------------------*/
917 /*----------------------------------------------------------------------------*/
918 
919 void
920 cs_sdm_add_mult(cs_sdm_t *mat,
922  const cs_sdm_t *add);
923 
924 /*----------------------------------------------------------------------------*/
932 /*----------------------------------------------------------------------------*/
933 
934 void
935 cs_sdm_square_add_transpose(cs_sdm_t *mat,
936  cs_sdm_t *tr);
937 
938 /*----------------------------------------------------------------------------*/
945 /*----------------------------------------------------------------------------*/
946 
947 void
948 cs_sdm_square_2symm(cs_sdm_t *mat);
949 
950 /*----------------------------------------------------------------------------*/
956 /*----------------------------------------------------------------------------*/
957 
958 void
959 cs_sdm_square_asymm(cs_sdm_t *mat);
960 
961 /*----------------------------------------------------------------------------*/
967 /*----------------------------------------------------------------------------*/
968 
969 void
970 cs_sdm_block_square_asymm(cs_sdm_t *mat);
971 
972 /*----------------------------------------------------------------------------*/
989 /*----------------------------------------------------------------------------*/
990 
991 void
993  cs_real_t Qt[9],
994  cs_real_t R[6]);
995 
996 /*----------------------------------------------------------------------------*/
1007 /*----------------------------------------------------------------------------*/
1008 
1009 void
1010 cs_sdm_33_lu_compute(const cs_sdm_t *m,
1011  cs_real_t facto[9]);
1012 
1013 /*----------------------------------------------------------------------------*/
1026 /*----------------------------------------------------------------------------*/
1027 
1028 void
1029 cs_sdm_lu_compute(const cs_sdm_t *m,
1030  cs_real_t facto[]);
1031 
1032 /*----------------------------------------------------------------------------*/
1046 /*----------------------------------------------------------------------------*/
1047 
1048 void
1049 cs_sdm_33_lu_solve(const cs_real_t facto[9],
1050  const cs_real_t rhs[3],
1051  cs_real_t sol[3]);
1052 
1053 /*----------------------------------------------------------------------------*/
1068 /*----------------------------------------------------------------------------*/
1069 
1070 void
1071 cs_sdm_lu_solve(cs_lnum_t n_rows,
1072  const cs_real_t facto[],
1073  const cs_real_t *rhs,
1074  cs_real_t *sol);
1075 
1076 /*----------------------------------------------------------------------------*/
1090 /*----------------------------------------------------------------------------*/
1091 
1092 void
1093 cs_sdm_33_ldlt_compute(const cs_sdm_t *m,
1094  cs_real_t facto[6]);
1095 
1096 /*----------------------------------------------------------------------------*/
1110 /*----------------------------------------------------------------------------*/
1111 
1112 void
1113 cs_sdm_44_ldlt_compute(const cs_sdm_t *m,
1114  cs_real_t facto[10]);
1115 
1116 /*----------------------------------------------------------------------------*/
1130 /*----------------------------------------------------------------------------*/
1131 
1132 void
1133 cs_sdm_66_ldlt_compute(const cs_sdm_t *m,
1134  cs_real_t facto[21]);
1135 
1136 /*----------------------------------------------------------------------------*/
1151 /*----------------------------------------------------------------------------*/
1152 
1153 void
1154 cs_sdm_ldlt_compute(const cs_sdm_t *m,
1155  cs_real_t *facto,
1156  cs_real_t *dkk);
1157 
1158 /*----------------------------------------------------------------------------*/
1168 /*----------------------------------------------------------------------------*/
1169 
1170 void
1171 cs_sdm_33_ldlt_solve(const cs_real_t facto[6],
1172  const cs_real_t rhs[3],
1173  cs_real_t sol[3]);
1174 
1175 /*----------------------------------------------------------------------------*/
1185 /*----------------------------------------------------------------------------*/
1186 
1187 void
1188 cs_sdm_44_ldlt_solve(const cs_real_t facto[10],
1189  const cs_real_t rhs[4],
1190  cs_real_t x[4]);
1191 
1192 /*----------------------------------------------------------------------------*/
1202 /*----------------------------------------------------------------------------*/
1203 
1204 void
1205 cs_sdm_66_ldlt_solve(const cs_real_t f[21],
1206  const cs_real_t b[6],
1207  cs_real_t x[6]);
1208 
1209 /*----------------------------------------------------------------------------*/
1221 /*----------------------------------------------------------------------------*/
1222 
1223 void
1224 cs_sdm_ldlt_solve(int n_rows,
1225  const cs_real_t *facto,
1226  const cs_real_t *rhs,
1227  cs_real_t *sol);
1228 
1229 /*----------------------------------------------------------------------------*/
1239 /*----------------------------------------------------------------------------*/
1240 
1241 double
1242 cs_sdm_test_symmetry(const cs_sdm_t *mat);
1243 
1244 /*----------------------------------------------------------------------------*/
1250 /*----------------------------------------------------------------------------*/
1251 
1252 void
1253 cs_sdm_simple_dump(const cs_sdm_t *mat);
1254 
1255 /*----------------------------------------------------------------------------*/
1264 /*----------------------------------------------------------------------------*/
1265 
1266 void
1267 cs_sdm_dump(cs_lnum_t parent_id,
1268  const cs_lnum_t *row_ids,
1269  const cs_lnum_t *col_ids,
1270  const cs_sdm_t *mat);
1271 
1272 /*----------------------------------------------------------------------------*/
1285 /*----------------------------------------------------------------------------*/
1286 
1287 void
1288 cs_sdm_fprintf(FILE *fp,
1289  const char *fname,
1290  cs_real_t thd,
1291  const cs_sdm_t *m);
1292 
1293 /*----------------------------------------------------------------------------*/
1300 /*----------------------------------------------------------------------------*/
1301 
1302 void
1303 cs_sdm_block_dump(cs_lnum_t parent_id,
1304  const cs_sdm_t *mat);
1305 
1306 /*----------------------------------------------------------------------------*/
1319 /*----------------------------------------------------------------------------*/
1320 
1321 void
1322 cs_sdm_block_fprintf(FILE *fp,
1323  const char *fname,
1324  cs_real_t thd,
1325  const cs_sdm_t *m);
1326 
1327 /*----------------------------------------------------------------------------*/
1328 
1330 
1331 #endif /* __CS_SDM_H__ */
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
Definition: cs_sdm.h:74
Definition: cs_sdm.h:58
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&#39;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