Program listing for file numerics/src/tools/NumericsMatrix.h

Program listing for file numerics/src/tools/NumericsMatrix.h#

  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