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