Program listing for file numerics/src/tools/NumericsMatrix.h#
Return to documentation for this file
1#ifndef NumericsMatrix_H
2#define NumericsMatrix_H
3
4
5
6#include <assert.h>
7#include <stdio.h>
8#include <stdlib.h>
9#include "CSparseMatrix.h"
10#include "NumericsFwd.h"
11#include "NumericsDataVersion.h"
12#include "NumericsSparseMatrix.h"
13#include "SiconosConfig.h"
14#include "NM_MPI.h"
15#ifndef __cplusplus
16#include <stdbool.h>
17#endif
18
19#ifdef WITH_OPENSSL
20#include <openssl/sha.h>
21#endif
22
23
24typedef struct
25{
26 size_t iWorkSize;
27 void *iWork;
28 size_t sizeof_elt ;
29 size_t dWorkSize;
30 double *dWork;
31 bool isLUfactorized;
32 bool isCholeskyfactorized;
33 bool isLDLTfactorized;
34 bool isInversed;
35#ifdef SICONOS_HAS_MPI
36 MPI_Comm mpi_comm;
37#endif
38#ifdef WITH_OPENSSL
39 unsigned int values_sha1_count;
40 unsigned char values_sha1[SHA_DIGEST_LENGTH];
41#endif
42} NumericsMatrixInternalData;
43
44
45typedef enum NumericsMatrix_types {
46 NM_DENSE,
47 NM_SPARSE_BLOCK,
48 NM_SPARSE,
49 NM_UNKNOWN,
50} NM_types;
51
52
53struct NumericsMatrix
54{
55 NM_types storageType;
56 int size0;
57 int size1;
58 double* matrix0;
59 SparseBlockStructuredMatrix* matrix1;
60 NumericsSparseMatrix* matrix2;
61
62 NumericsMatrixInternalData* internalData;
63
64 NumericsDataVersion version;
65
66 NumericsMatrix* destructible;
67};
68
69typedef struct
70{
71 int size0;
72 int size1;
73 NumericsMatrix* D1;
74 NumericsMatrix* D2;
75 NumericsMatrix* A;
76} BalancingMatrices;
77
78
79typedef NumericsMatrix RawNumericsMatrix;
80
81
82typedef enum {
83 NM_NONE,
84 NM_KEEP_FACTORS,
85 NM_PRESERVE
86} NM_gesv_opts;
87
88#if defined(__cplusplus) && !defined(BUILD_AS_CPP)
89extern "C"
90{
91#endif
92
93
94
95
96
97 RawNumericsMatrix* NM_new(void);
98 RawNumericsMatrix* NM_eye(int size);
99 RawNumericsMatrix* NM_scalar(int size, double s);
100
101
102 RawNumericsMatrix* NM_create(NM_types storageType, int size0, int size1);
103
104
105 RawNumericsMatrix* NM_create_from_data(int storageType, int size0, int size1, void* data);
106 RawNumericsMatrix* NM_create_from_filename(const char *filename);
107 RawNumericsMatrix* NM_create_from_file(FILE *file);
108
109
110 void NM_version_copy(const NumericsMatrix* const A, NumericsMatrix* B);
111
112
113 void NM_copy(const NumericsMatrix* const A, NumericsMatrix* B);
114
115
116 void NM_copy_to_sparse(NumericsMatrix* A, NumericsMatrix* B, double threshold);
117
118
119 RawNumericsMatrix* NM_duplicate(NumericsMatrix* mat);
120
121
122
123 NumericsSparseMatrix* numericsSparseMatrix(NumericsMatrix* A);
124
125
126 CSparseMatrix* NM_triplet(NumericsMatrix* A);
127
128
129 CSparseMatrix* NM_half_triplet(NumericsMatrix* A);
130
131
132 CSparseMatrix* NM_csc(NumericsMatrix *A);
133
134
135 CSparseMatrix* NM_csc_trans(NumericsMatrix* A);
136
137
138 CSparseMatrix* NM_csr(NumericsMatrix *A);
139
140
141 void NM_fill(NumericsMatrix* M, NM_types storageType, int size0, int size1, void* data);
142
143
144 RawNumericsMatrix* NM_new_SBM(int size0, int size1, SparseBlockStructuredMatrix* m1);
145
146
147 RawNumericsMatrix* NM_transpose(NumericsMatrix * A);
148
149
150 void NM_null(NumericsMatrix* A);
151
152
153 bool NM_destructible(const NumericsMatrix* A);
154
155
156 RawNumericsMatrix* NM_preserve(NumericsMatrix* A);
157
158
159 RawNumericsMatrix* NM_unpreserve(NumericsMatrix* A);
160
161
162 bool NM_LU_factorized(const NumericsMatrix* const A);
163
164
165 bool NM_Cholesky_factorized(const NumericsMatrix* const A);
166
167
168 bool NM_LDLT_factorized(const NumericsMatrix* const A);
169
170
171 void NM_set_LU_factorized(NumericsMatrix* A, bool flag);
172 void NM_set_Cholesky_factorized(NumericsMatrix* A, bool flag);
173 void NM_set_LDLT_factorized(NumericsMatrix* A, bool flag);
174
175
176 void NM_update_size(NumericsMatrix* A);
177
178
179 void NM_csc_alloc(NumericsMatrix* A, CS_INT nzmax);
180
181
182 void NM_csc_empty_alloc(NumericsMatrix* A, CS_INT nzmax);
183
184
185 void NM_triplet_alloc(NumericsMatrix* A, CS_INT nzmax);
186
187
188
189
190 void NM_clear(NumericsMatrix* m);
191
192 NumericsMatrix * NM_free(NumericsMatrix* m);
193
194
195 void NM_clear_not_dense(NumericsMatrix* m);
196 NumericsMatrix * NM_free_not_dense(NumericsMatrix* m);
197
198 void NM_clear_not_SBM(NumericsMatrix* m);
199 NumericsMatrix * NM_free_not_SBM(NumericsMatrix* m);
200
201
202
203 void NM_clear_other_storages(NumericsMatrix* M, NM_types storageType);
204
205
206 void NM_zentry(NumericsMatrix* M, int i, int j, double val, double threshold);
207
208
209 void NM_entry(NumericsMatrix* M, int i, int j, double val);
210
211
212 double NM_get_value(const NumericsMatrix* const M, int i, int j);
213
214
215 bool NM_equal(NumericsMatrix* A, NumericsMatrix* B);
216
217
218 bool NM_compare(NumericsMatrix* A, NumericsMatrix* B, double tol);
219
220
221 size_t NM_nnz(const NumericsMatrix* M);
222
223
224 void NM_extract_diag_block(NumericsMatrix* M, int block_row_nb, size_t start_row,
225 int size, double **Block);
226
227
228
229 void NM_extract_diag_block3(NumericsMatrix* M, int block_row_nb, double **Block);
230
231
232
233 void NM_extract_diag_block2(NumericsMatrix* M, int block_row_nb, double **Block);
234
235
236 void NM_extract_diag_block5(NumericsMatrix* M, int block_row_nb, double **Block);
237
238
239 void NM_copy_diag_block3(NumericsMatrix* M, int block_row_nb, double **Block);
240
241
242
243 void NM_insert(NumericsMatrix* A, const NumericsMatrix* const B,
244 const unsigned int start_i, const unsigned int start_j);
245
246
247
248
249 void NM_prod_mv_3x3(int sizeX, int sizeY, NumericsMatrix* A,
250 double* const x, double* y);
251
252
253 void NM_row_prod(int sizeX, int sizeY, int currentRowNumber, NumericsMatrix* A, const double* const x, double* y, int init);
254
255
256 void NM_row_prod_no_diag(size_t sizeX, size_t sizeY, int block_start, size_t row_start, NumericsMatrix* A, double* x, double* y, double* xsave, bool init);
257
258
259 void NM_row_prod_no_diag3(size_t sizeX, int block_start, size_t row_start, NumericsMatrix* A, double* x, double* y, bool init);
260
261
262 void NM_row_prod_no_diag2(size_t sizeX, int block_start, size_t row_start, NumericsMatrix* A, double* x, double* y, bool init);
263
264
265 void NM_row_prod_no_diag1x1(size_t sizeX, int block_start, size_t row_start, NumericsMatrix* A, double* x, double* y, bool init);
266
267
268 void NM_gemv(const double alpha, NumericsMatrix* A, const double *x,
269 const double beta,
270 double *y);
271
272
273 void NM_gemm(const double alpha, NumericsMatrix* A, NumericsMatrix* B,
274 const double beta, NumericsMatrix *C);
275
276
277 RawNumericsMatrix * NM_multiply(NumericsMatrix* A, NumericsMatrix* B);
278
279
280 void NM_tgemv(const double alpha, NumericsMatrix* A, const double *x,
281 const double beta,
282 double *y);
283
284
285
286
287 void NM_dense_to_sparse(NumericsMatrix* A, NumericsMatrix* B, double threshold);
288
289
290 int NM_to_dense(NumericsMatrix* A, NumericsMatrix* B);
291
292
293 void NM_dense_display_matlab(double * m, int nRow, int nCol, int lDim);
294
295
296 void NM_dense_display(double * m, int nRow, int nCol, int lDim);
297
298
299 void NM_vector_display(double * m, int nRow);
300
301
302
303 void NM_display(const NumericsMatrix* const M);
304
305
306 void NM_display_storageType(const NumericsMatrix* const M);
307
308
309
310 void NM_display_row_by_row(const NumericsMatrix* const m);
311
312
313
314
315
316
317 void NM_write_in_filename(const NumericsMatrix* const M, const char *filename);
318
319
320 void NM_read_in_filename(NumericsMatrix* const M, const char *filename);
321
322
323
324 void NM_write_in_file(const NumericsMatrix* const M, FILE* file);
325
326
327 void NM_read_in_file(NumericsMatrix* const M, FILE *file);
328
329
330 RawNumericsMatrix* NM_new_from_file(FILE *file);
331 RawNumericsMatrix* NM_new_from_filename(const char * filename);
332
333
334 void NM_write_in_file_scilab(const NumericsMatrix* const M, FILE* file);
335
336
337 void NM_write_in_file_python(const NumericsMatrix* const M, FILE* file);
338
339
340 void NM_read_in_file_scilab(NumericsMatrix* const M, FILE *file);
341
342
343 void NM_clearDense(NumericsMatrix* A);
344
345
346 void NM_clearSparseBlock(NumericsMatrix* A);
347
348
349 void NM_clearSparse(NumericsMatrix* A);
350
351
352 void NM_clearTriplet(NumericsMatrix* A);
353
354
355 void NM_clearHalfTriplet(NumericsMatrix* A);
356
357
358 void NM_clearCSC(NumericsMatrix* A);
359
360
361 void NM_clearCSCTranspose(NumericsMatrix* A);
362
363
364 void NM_clearCSR(NumericsMatrix* A);
365
366
367 void NM_clearSparseStorage(NumericsMatrix *A);
368
369
370
371
372
373
374 int NM_LU_factorize(NumericsMatrix* A);
375 int NM_Cholesky_factorize(NumericsMatrix* A);
376 int NM_LDLT_factorize(NumericsMatrix* A);
377
378
379 int NM_LU_solve(NumericsMatrix* A, double *b, unsigned int nrhs);
380 int NM_LU_solve_matrix_rhs(NumericsMatrix* Ao, NumericsMatrix* B);
381 int NM_LU_refine(NumericsMatrix* A, double *x, double tol, int max_iter, double *residu);
382 int NM_Cholesky_solve(NumericsMatrix* A, double *b, unsigned int nrhs);
383 int NM_Cholesky_solve_matrix_rhs(NumericsMatrix* Ao, NumericsMatrix* B);
384 int NM_LDLT_solve(NumericsMatrix* A, double *b, unsigned int nrhs);
385 int NM_LDLT_refine(NumericsMatrix* Ao, double *x , double *b, unsigned int nrhs, double tol, int maxitref, int job );
386
387
388 int NM_gesv_expert(NumericsMatrix* A, double *b, unsigned keep);
389 int NM_posv_expert(NumericsMatrix* A, double *b, unsigned keep);
390
391 int NM_gesv_expert_multiple_rhs(NumericsMatrix* A, double *b, unsigned int n_rhs, unsigned keep);
392
393 int NM_Linear_solver_finalize(NumericsMatrix* Ao);
394
395
396 NumericsMatrix* NM_LU_inv(NumericsMatrix* A);
397
398
399 int NM_inverse_diagonal_block_matrix_in_place(NumericsMatrix* A);
400
401
402 NumericsMatrix * NM_inverse_diagonal_block_matrix(NumericsMatrix* A, unsigned int block_number, unsigned int * blocksizes);
403
404
405 static inline int NM_gesv(NumericsMatrix* A, double *b, bool preserve)
406 {
407 return NM_gesv_expert(A, b, preserve ? NM_PRESERVE : NM_NONE);
408 }
409
410
411 NumericsMatrix* NM_gesv_inv(NumericsMatrix* A);
412
413
414
415
416
417 void NM_setSparseSolver(NumericsMatrix* A, NSM_linear_solver solver_id);
418
419
420 NumericsMatrixInternalData* NM_internalData(NumericsMatrix* A);
421
422
423 void NM_internalData_new(NumericsMatrix* M);
424
425
426 void NM_internalData_copy(const NumericsMatrix* const A, NumericsMatrix* B );
427
428
429 void* NM_iWork(NumericsMatrix *A, size_t size, size_t sizeof_elt);
430
431
432 double* NM_dWork(NumericsMatrix *A, int size);
433
434
435 void NM_add_to_diag3(NumericsMatrix* M, double alpha);
436
437
438 void NM_add_to_diag5(NumericsMatrix* M, double alpha);
439
440
441 RawNumericsMatrix * NM_add(double alpha, NumericsMatrix* A, double beta, NumericsMatrix* B);
442
443
444 void NM_scal(double alpha, NumericsMatrix* A);
445
446
447 static inline void NM_assert(NM_types type, NumericsMatrix* M)
448 {
449#ifndef NDEBUG
450 assert(M && "NM_assert :: the matrix is NULL");
451 assert(M->storageType == type && "NM_assert :: the matrix has the wrong type");
452 switch(type)
453 {
454 case NM_DENSE:
455 assert(M->matrix0);
456 break;
457 case NM_SPARSE_BLOCK:
458 assert(M->matrix1);
459 break;
460 case NM_SPARSE:
461 assert(M->matrix2);
462 break;
463 default:
464 assert(0 && "NM_assert :: unknown storageType");
465 }
466#endif
467 }
468
469
470 int NM_check(const NumericsMatrix* const A);
471
472
473 double NM_norm_1(NumericsMatrix* const A);
474
475
476 double NM_norm_inf(NumericsMatrix* const A);
477
478 int NM_is_symmetric(NumericsMatrix* A);
479 double NM_symmetry_discrepancy(NumericsMatrix* A);
480
481
482
483 static inline NumericsMatrix* NM_convert(NumericsMatrix* A)
484 {
485 return A;
486 }
487
488
489 double NM_iterated_power_method(NumericsMatrix* A, double tol, int itermax);
490
491
492 int NM_max_by_columns(NumericsMatrix *A, double * max);
493
494
495 int NM_max_by_rows(NumericsMatrix *A, double * max);
496
497
498 int NM_max_abs_by_columns(NumericsMatrix *A, double * max);
499
500
501 int NM_max_abs_by_rows(NumericsMatrix *A, double * max);
502
503
504 int NM_compute_balancing_matrices(NumericsMatrix* A, double tol, int itermax, BalancingMatrices * B);
505
506
507 BalancingMatrices * NM_BalancingMatrices_new(NumericsMatrix* A);
508
509
510
511 BalancingMatrices * NM_BalancingMatrices_free(BalancingMatrices* A);
512
513
514 void NM_reset_version(NumericsMatrix* M, NM_types id);
515
516
517 void NM_reset_versions(NumericsMatrix* M);
518
519
520 void NM_version_sync(NumericsMatrix* M);
521
522
523 int NM_isnan(NumericsMatrix* M);
524
525#ifdef WITH_OPENSSL
526
527 void NM_compute_values_sha1(NumericsMatrix* A, unsigned char * digest);
528
529
530 unsigned char* NM_values_sha1(NumericsMatrix* A);
531
532
533 void NM_set_values_sha1(NumericsMatrix* A);
534
535
536 void NM_clear_values_sha1(NumericsMatrix* A);
537
538
539 bool NM_check_values_sha1(NumericsMatrix* A);
540
541
542 bool NM_equal_values_sha1(NumericsMatrix* A, NumericsMatrix* B);
543
544#endif
545
546#if defined(__cplusplus) && !defined(BUILD_AS_CPP)
547}
548#endif
549
550#endif