1 /* Siconos is a program dedicated to modeling, simulation and control 2 * of non smooth dynamical systems. 3 * 4 * Copyright 2021 INRIA. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 #ifndef NumericsMatrix_H 20 #define NumericsMatrix_H 21 22 /*!\file NumericsMatrix.h 23 \brief Structure definition and functions related to matrix storage in Numerics 24 */ 25 26 #include <assert.h> // for assert 27 #include <stdio.h> // for size_t, FILE, NULL 28 #include <stdlib.h> // for malloc 29 #include "CSparseMatrix.h" // for CS_INT, CSparseMatrix 30 #include "NumericsFwd.h" // for NumericsMatrix, NumericsSparseMatrix, Spa... 31 #include "NumericsDataVersion.h" // Versioning 32 #include "NumericsSparseMatrix.h" // for NSM_linear_solver typedef 33 #include "SiconosConfig.h" // for BUILD_AS_CPP, SICONOS_HAS_MP // IWYU pragma: keep 34 #include "NM_MPI.h" 35 #ifndef __cplusplus 36 #include <stdbool.h> // for bool 37 #endif 38 39 #ifdef WITH_OPENSSL 40 #include <openssl/sha.h> 41 #endif 42 43 /** \struct NumericsMatrixInternalData NumericsMatrix.h 44 * Structure for simple workspaces 45 */ 46 typedef struct 47 { 48 size_t iWorkSize; /**< size of iWork */ 49 void *iWork; /**< integer workspace */ 50 size_t sizeof_elt ; /**< sizeof_elt of an element in bytes (result of sizeof for instance)*/ 51 size_t dWorkSize; /**< size of dWork */ 52 double *dWork; /**< double workspace */ 53 bool isLUfactorized; /**< true if the matrix has already been LU-factorized */ 54 bool isCholeskyfactorized; /**< true if the matrix has already been Cholesky factorized */ 55 bool isLDLTfactorized; /**< true if the matrix has already been LDLT factorized */ 56 bool isInversed; /**< true if the matrix contains its inverse (in place inversion) */ 57 #ifdef SICONOS_HAS_MPI 58 MPI_Comm mpi_comm; /**< optional mpi communicator */ 59 #endif 60 #ifdef WITH_OPENSSL 61 unsigned int values_sha1_count; /**< counter for sha1 */ 62 unsigned char values_sha1[SHA_DIGEST_LENGTH]; /**< sha1 hash of 63 * values. Matrices of 64 * differents sizes may have 65 * the same hash */ 66 #endif 67 } NumericsMatrixInternalData; 68 69 /*! Available types of storage for NumericsMatrix */ 70 typedef enum NumericsMatrix_types { 71 NM_DENSE, /**< dense format */ 72 NM_SPARSE_BLOCK, /**< sparse block format */ 73 NM_SPARSE, /**< compressed column format */ 74 NM_UNKNOWN, /** unset. Used in NM_null */ 75 } NM_types; 76 77 /** \struct NumericsMatrix NumericsMatrix.h 78 Interface to different type of matrices in numerics component. 79 80 See NM_* functions for linear algebra operations on dense, sparse block and sparse storage. 81 */ 82 struct NumericsMatrix 83 { 84 NM_types storageType; /**< the type of storage: 85 0: dense (double*), 86 1: SparseBlockStructuredMatrix, 87 2: classical sparse (csc, csr or triplet) via CSparse (from T. Davis)*/ 88 int size0; /**< number of rows */ 89 int size1; /**< number of columns */ 90 double* matrix0; /**< dense storage */ 91 SparseBlockStructuredMatrix* matrix1; /**< sparse block storage */ 92 NumericsSparseMatrix* matrix2; /**< csc, csr or triplet storage */ 93 94 NumericsMatrixInternalData* internalData; /**< internal storage, used for workspace among other things */ 95 96 NumericsDataVersion version; /*< version of dense storage */ 97 98 NumericsMatrix* destructible; /**< pointer on the destructible 99 * matrix, by default points toward 100 * the matrix itself */ 101 }; 102 103 typedef struct 104 { 105 int size0; 106 int size1; 107 NumericsMatrix* D1; 108 NumericsMatrix* D2; 109 NumericsMatrix* A; 110 } BalancingMatrices; 111 112 /*! RawNumericsMatrix is used without conversion in python */ 113 typedef NumericsMatrix RawNumericsMatrix; 114 115 116 typedef enum { 117 NM_NONE, /**< keep nothing */ 118 NM_KEEP_FACTORS, /**< keep all the factorization data (useful to reuse the factorization) */ 119 NM_PRESERVE /**< keep the matrix as-is (useful for the dense case) */ 120 } NM_gesv_opts; 121 122 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 123 extern "C" 124 { 125 #endif 126 /**************************************************/ 127 /** Constructors and destructors ****************/ 128 /**************************************************/ 129 130 /** Creation of an empty NumericsMatrix. 131 * \return a pointer to allocated space 132 */ 133 RawNumericsMatrix* NM_new(void); 134 RawNumericsMatrix* NM_eye(int size); 135 136 /** create a NumericsMatrix and allocate the memory according to the matrix type 137 * \param storageType the type of storage 138 * \param size0 number of rows 139 * \param size1 number of columns 140 * \return a pointer to a NumericsMatrix 141 */ 142 RawNumericsMatrix* NM_create(NM_types storageType, int size0, int size1); 143 144 /** create a NumericsMatrix and possibly set the data 145 * \param storageType the type of storage 146 * \param size0 number of rows 147 * \param size1 number of columns 148 * \param data pointer to the matrix data. If NULL, all matrixX fields are 149 * set to NULL 150 * \return a pointer to a NumericsMatrix 151 */ 152 RawNumericsMatrix* NM_create_from_data(int storageType, int size0, int size1, void* data); 153 RawNumericsMatrix* NM_create_from_filename(const char *filename); 154 RawNumericsMatrix* NM_create_from_file(FILE *file); 155 156 157 /** Copy a NumericsMatrix inside another NumericsMatrix (deep). 158 * Reallocations are performed if B cannot hold a copy of A 159 * \param[in] A a NumericsMatrix 160 * \param[in,out] B a NumericsMatrix 161 */ 162 void NM_copy(const NumericsMatrix* const A, NumericsMatrix* B); 163 164 /** Copy a NumericsMatrix to s sparse one. 165 * Allocation or reallocation are performed on B 166 * \warning It is assumed that B has been properly initialized: its storageType must 167 * be set to NM_SPARSE. 168 * \param[in] A a NumericsMatrix 169 * \param[in,out] B a NumericsMatrix 170 * \param threshold if the original matrix is dense, a threshold can be applied 171 * on the absolute value of the entries 172 */ 173 void NM_copy_to_sparse(const NumericsMatrix* const A, NumericsMatrix* B, double threshold); 174 175 /** create a NumericsMatrix similar to the another one. The structure is the same 176 * \param mat the model matrix 177 * \return a pointer to a NumericsMatrix 178 */ 179 RawNumericsMatrix* NM_duplicate(NumericsMatrix* mat); 180 181 182 /** Creation, if needed, of sparse matrix storage. 183 * \param[in,out] A a NumericsMatrix 184 * \return a pointer on the sparse matrix storage 185 */ 186 NumericsSparseMatrix* numericsSparseMatrix(NumericsMatrix* A); 187 188 /** Creation, if needed, of triplet storage from sparse block storage. 189 * \param[in,out] A a NumericsMatrix initialized with sparsed block storage. 190 * \return the triplet sparse Matrix created in A. 191 */ 192 CSparseMatrix* NM_triplet(NumericsMatrix* A); 193 194 /** Creation, if needed, of half triplet storage from sparse block storage. 195 * \param[in,out] A a NumericsMatrix initialized with sparsed block storage. 196 * \return the triplet sparse Matrix created in A. 197 */ 198 CSparseMatrix* NM_half_triplet(NumericsMatrix* A); 199 200 /** Creation, if needed, of compress column storage of a NumericsMatrix. 201 * \param[in,out] A a NumericsMatrix with sparse block storage initialized 202 * \return the compressed column CSparseMatrix created in A. 203 */ 204 CSparseMatrix* NM_csc(NumericsMatrix *A); 205 206 /** Creation, if needed, of the transposed compress column storage 207 * from compress column storage. 208 * \param[in,out] A a NumericsMatrix with sparse block storage. 209 * \return the transposed compressed column matrix created in A. 210 */ 211 CSparseMatrix* NM_csc_trans(NumericsMatrix* A); 212 213 /** Creation, if needed, of compress row storage of a NumericsMatrix 214 * \warning This rely on the MKL 215 * \param[in,out] A a NumericsMatrix with sparse block storage initialized 216 * \return the compressed row CSparseMatrix created in A. 217 */ 218 CSparseMatrix* NM_csr(NumericsMatrix *A); 219 220 /** fill an existing NumericsMatrix struct 221 * \param[in,out] M the struct to fill 222 * \param storageType the type of storage 223 * \param size0 number of rows 224 * \param size1 number of columns 225 * \param data pointer to the matrix data. If NULL, all matrixX fields are 226 * set to NULL 227 */ 228 void NM_fill(NumericsMatrix* M, NM_types storageType, int size0, int size1, void* data); 229 230 /** new NumericsMatrix with sparse storage from minimal set of data 231 * \param[in] size0 number of rows 232 * \param[in] size1 number of columns 233 * \param[in] m1 the SparseBlockStructuredMatrix 234 * \return a pointer to a NumericsMatrix 235 */ 236 RawNumericsMatrix* NM_new_SBM(int size0, int size1, SparseBlockStructuredMatrix* m1); 237 238 /** new NumericsMatrix equal to the transpose of a given matrix 239 * \param[in] A 240 * \return a pointer to a NumericsMatrix 241 */ 242 RawNumericsMatrix* NM_transpose(NumericsMatrix * A); 243 244 /** set NumericsMatrix fields to NULL 245 * \param A a matrix 246 */ 247 void NM_null(NumericsMatrix* A); 248 249 /** Check if a matrix is destructible. 250 * \param[in] A the NumericsMatrix 251 * \return true if the matrix is destructible */ 252 bool NM_destructible(const NumericsMatrix* A); 253 254 /** Preservation of a matrix before in-place transformations such as 255 * factorizations. 256 * \param[in] A the NumericsMatrix 257 * \return a pointer on the preserved Matrix; 258 */ 259 RawNumericsMatrix* NM_preserve(NumericsMatrix* A); 260 261 /** Set the matrix as destructible, clear the preserved data. 262 * \param[in] A the NumericsMatrix 263 * \return a pointer on the Matrix; 264 */ 265 RawNumericsMatrix* NM_unpreserve(NumericsMatrix* A); 266 267 /** Check for a previous LU factorization. 268 * \param[in] A the NumericsMatrix 269 * \return true if the matrix has been LU factorized. 270 */ 271 bool NM_LU_factorized(const NumericsMatrix* const A); 272 273 /** Check for a previous Cholesky factorization. 274 * \param[in] A the NumericsMatrix 275 * \return true if the matrix has been Cholesky factorized. 276 */ 277 bool NM_Cholesky_factorized(const NumericsMatrix* const A); 278 279 /** Check for a previous LDLT factorization. 280 * \param[in] A the NumericsMatrix 281 * \return true if the matrix has been Cholesky factorized. 282 */ 283 bool NM_LDLT_factorized(const NumericsMatrix* const A); 284 285 /** Set the factorization flag. 286 * \param[in] A the NumericsMatrix, 287 * \param[in] flag a boolean. 288 */ 289 void NM_set_LU_factorized(NumericsMatrix* A, bool flag); 290 void NM_set_Cholesky_factorized(NumericsMatrix* A, bool flag); 291 void NM_set_LDLT_factorized(NumericsMatrix* A, bool flag); 292 293 /** update the size of the matrix based on the matrix data 294 * \param[in,out] A the matrix which size is updated*/ 295 void NM_update_size(NumericsMatrix* A); 296 297 /** Allocate a csc matrix in A 298 * \param A the matrix 299 * \param nzmax number of non-zero elements 300 */ 301 void NM_csc_alloc(NumericsMatrix* A, CS_INT nzmax); 302 303 /** Allocate a csc matrix in A and set the vector of 304 * column pointers to 0 such that the matrix is empty. 305 * \param A the matrix 306 * \param nzmax number of non-zero elements 307 */ 308 void NM_csc_empty_alloc(NumericsMatrix* A, CS_INT nzmax); 309 310 /** Allocate a triplet matrix in A 311 * \param A the matrix 312 * \param nzmax maximum number of non-zero elements 313 */ 314 void NM_triplet_alloc(NumericsMatrix* A, CS_INT nzmax); 315 316 317 /** Free memory for a NumericsMatrix. Warning: call this function only if you are sure that 318 memory has been allocated for the structure in Numerics. This function is assumed that the memory is "owned" by this structure. 319 Note that this function does not free m. 320 \param m the matrix to be deleted. 321 */ 322 323 void NM_clear(NumericsMatrix* m); 324 325 NumericsMatrix * NM_free(NumericsMatrix* m); 326 327 /** Free memory for a NumericsMatrix except the dense matrix that is assumed not to be owned. 328 \param m the matrix to be cleared. 329 */ 330 void NM_clear_not_dense(NumericsMatrix* m); 331 NumericsMatrix * NM_free_not_dense(NumericsMatrix* m); 332 /** Free memory for a NumericsMatrix except the SBM matrix that is assumed not to be owned. 333 Note that this function does not free m. 334 \param m the matrix to be cleared. 335 */ 336 void NM_clear_not_SBM(NumericsMatrix* m); 337 NumericsMatrix * NM_free_not_SBM(NumericsMatrix* m); 338 339 340 341 342 343 /** Free memory for a NumericsMatrix except for a given storage. Warning: call this function only if you are sure that 344 memory has been allocated for the structure in Numerics. This function is assumed that the memory is "owned" by this structure. 345 Note that this function does not free m. 346 \param m the matrix to be deleted. 347 \param storageType to be kept. 348 */ 349 void NM_clear_other_storages(NumericsMatrix* M, NM_types storageType); 350 351 /**************************************************/ 352 /** setters and getters *************/ 353 /**************************************************/ 354 355 /** insert an non zero entry into a NumericsMatrix. 356 * for storageType = NM_SPARSE, a conversion to triplet is done for performing the entry in the 357 * matrix. This method is expensive in terms of memory management. For a lot of entries, use 358 * preferably a triplet matrix. 359 * \param M the NumericsMatrix 360 * \param i row index 361 * \param j column index 362 * \param val the value to be inserted. 363 * \param threshold a threshold to filter the small value in magnitude (useful for dense to sparse conversion) 364 */ 365 void NM_zentry(NumericsMatrix* M, int i, int j, double val, double threshold); 366 367 /** insert an entry into a NumericsMatrix. 368 * for storageType = NM_SPARSE, a conversion to triplet is done for performing the entry in the 369 * matrix. This method is expensive in terms of memory management. For a lot of entries, use 370 * preferably a triplet matrix. 371 * \param M the NumericsMatrix 372 * \param i row index 373 * \param j column index 374 * \param val the value to be inserted. 375 */ 376 void NM_entry(NumericsMatrix* M, int i, int j, double val); 377 378 /** get the value of a NumericsMatrix. 379 * \param M the NumericsMatrix 380 * \param i row index 381 * \param j column index 382 * \return the value to be inserted. 383 */ 384 double NM_get_value(const NumericsMatrix* const M, int i, int j); 385 386 /** compare to NumericsMatrix up to machine accuracy (DBL_EPSILON) 387 * \param A the NumericsMatrix 388 * \param B the NumericsMatrix 389 */ 390 bool NM_equal(NumericsMatrix* A, NumericsMatrix* B); 391 392 /** compare to NumericsMatrix up to a given tolerance 393 * \param A the NumericsMatrix 394 * \param B the NumericsMatrix 395 * \param tol the tolerance 396 */ 397 bool NM_compare(NumericsMatrix* A, NumericsMatrix* B, double tol); 398 399 /** return the number of non-zero element. For a dense matrix, it is the 400 * product of the dimensions (e.g. an upper bound). For a sparse matrix, it is the true number 401 * \param M the matrix 402 * \return the number (or an upper bound) of non-zero elements in the matrix 403 */ 404 size_t NM_nnz(const NumericsMatrix* M); 405 406 /** get the (square) diagonal block of a NumericsMatrix. No allocation is done. 407 * \param[in] M a NumericsMatrix 408 * \param[in] block_row_nb the number of the block Row. Useful only in sparse case 409 * \param[in] start_row the starting row. Useful only in dense case. 410 * \param[in] size of the diag block. Only useful in dense case. 411 * \param[out] Block the target. In the dense and sparse case (*Block) must be allocated by caller. 412 * In case of SBM case **Bout contains the resulting block (from the SBM). 413 */ 414 void NM_extract_diag_block(NumericsMatrix* M, int block_row_nb, size_t start_row, 415 int size, double **Block); 416 417 /** get a 3x3 diagonal block of a NumericsMatrix. No allocation is done. 418 * \param[in] M a NumericsMatrix 419 * \param[in] block_row_nb the number of the block row 420 * \param[out] Block the target. In the dense and sparse case (*Block) must be allocated by caller. 421 * In case of SBM case **Bout contains the resulting block (from the SBM). 422 */ 423 424 void NM_extract_diag_block3(NumericsMatrix* M, int block_row_nb, double **Block); 425 426 /** get a 5x5 diagonal block of a NumericsMatrix. No allocation is done. 427 * \param[in] M a NumericsMatrix 428 * \param[in] block_row_nb the number of the block row 429 * \param[out] Block the target. In the dense and sparse case (*Block) must be allocated by caller. 430 * In case of SBM case **Bout contains the resulting block (from the SBM). 431 */ 432 void NM_extract_diag_block5(NumericsMatrix* M, int block_row_nb, double **Block); 433 434 /** get a 3x3 diagonal block of a NumericsMatrix. No allocation is done. 435 * \param[in] M a NumericsMatrix 436 * \param[in] block_row_nb the number of the block row 437 * \param[out] Block the target. 438 * In all cases (dense, sbm, and sparse) (*Block) must be allocated by caller. 439 * A copy is always performed 440 */ 441 void NM_copy_diag_block3(NumericsMatrix* M, int block_row_nb, double **Block); 442 443 444 /** Set the submatrix B into the matrix A on the position defined in 445 * (start_i, start_j) position. 446 * \param[in] A a pointer to NumerixMatrix 447 * \param[in] B a pointer toNumericsMatrix 448 * \param[in] start_i a start row index 449 * \param[in] start_j a start column index 450 */ 451 void NM_insert(NumericsMatrix* A, const NumericsMatrix* const B, 452 const unsigned int start_i, const unsigned int start_j); 453 454 /**************************************************/ 455 /** Matrix - vector product *************/ 456 /**************************************************/ 457 458 /** Matrix - vector product y = A*x + y 459 \param[in] sizeX dim of the vector x 460 \param[in] sizeY dim of the vector y 461 \param[in] A the matrix to be multiplied 462 \param[in] x the vector to be multiplied 463 \param[in,out] y the resulting vector 464 */ 465 void NM_prod_mv_3x3(int sizeX, int sizeY, NumericsMatrix* A, 466 double* const x, double* y); 467 468 /** Row of a Matrix - vector product y = rowA*x or y += rowA*x, rowA being a submatrix of A (sizeY rows and sizeX columns) 469 \param[in] sizeX dim of the vector x 470 \param[in] sizeY dim of the vector y 471 \param[in] currentRowNumber position of the first row of rowA in A (warning: real row if A is a double*, block-row if A is a SparseBlockStructuredMatrix) 472 \param[in] A the matrix to be multiplied 473 \param[in] x the vector to be multiplied 474 \param[in,out] y the resulting vector 475 \param[in] init = 0 for y += Ax, =1 for y = Ax 476 */ 477 void NM_row_prod(int sizeX, int sizeY, int currentRowNumber, const NumericsMatrix* const A, const double* const x, double* y, int init); 478 479 /** Row of a Matrix - vector product y = rowA*x or y += rowA*x, rowA being a submatrix of A (sizeY rows and sizeX columns) 480 \param[in] sizeX dim of the vector x 481 \param[in] sizeY dim of the vector y 482 \param[in] block_start block number (only used for SBM) 483 \param[in] row_start position of the first row of A (unused if A is SBM) 484 \param[in] A the matrix to be multiplied 485 \param[in] x the vector to be multiplied 486 \param[in,out] y the resulting vector 487 \param[in] xsave storage for saving the part of x set to 0 488 \param[in] init if True y = Ax, else y += Ax 489 */ 490 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); 491 492 /** Row of a Matrix - vector product y = rowA*x or y += rowA*x, rowA being a submatrix of A (3 rows and sizeX columns) 493 \param[in] sizeX dim of the vector x 494 \param[in] block_start block number (only used for SBM) 495 \param[in] row_start position of the first row of A (unused if A is SBM) 496 \param[in] A the matrix to be multiplied 497 \param[in] x the vector to be multiplied 498 \param[in,out] y the resulting vector 499 \param[in] init if True y = Ax, else y += Ax 500 */ 501 void NM_row_prod_no_diag3(size_t sizeX, int block_start, size_t row_start, NumericsMatrix* A, double* x, double* y, bool init); 502 503 504 void NM_row_prod_no_diag1x1(size_t sizeX, int block_start, size_t row_start, NumericsMatrix* A, double* x, double* y, bool init); 505 506 /** Matrix vector multiplication : y = alpha A x + beta y 507 * \param[in] alpha scalar 508 * \param[in] A a NumericsMatrix 509 * \param[in] x pointer on a dense vector of size A->size1 510 * \param[in] beta scalar 511 * \param[in,out] y pointer on a dense vector of size A->size1 512 */ 513 void NM_gemv(const double alpha, NumericsMatrix* A, const double *x, 514 const double beta, 515 double *y); 516 517 /** Matrix matrix multiplication : C = alpha A B + beta C 518 * \param[in] alpha scalar 519 * \param[in] A a NumericsMatrix 520 * \param[in] B a NumericsMatrix 521 * \param[in] beta scalar 522 * \param[in,out] C a NumericsMatrix 523 */ 524 void NM_gemm(const double alpha, NumericsMatrix* A, NumericsMatrix* B, 525 const double beta, NumericsMatrix *C); 526 527 /** Matrix matrix multiplication : C = A B 528 * \param[in] A a NumericsMatrix 529 * \param[in] B a NumericsMatrix 530 * \param[in,out] C a NumericsMatrix 531 */ 532 RawNumericsMatrix * NM_multiply(NumericsMatrix* A, NumericsMatrix* B); 533 534 /** Transposed matrix multiplication : y += alpha transpose(A) x + y 535 * \param[in] alpha scalar 536 * \param[in] A a NumericsMatrix 537 * \param[in] x pointer on a dense vector of size A->size1 538 * \param[in] beta scalar 539 * \param[in,out] y pointer on a dense vector of size A->size1 540 */ 541 void NM_tgemv(const double alpha, NumericsMatrix* A, const double *x, 542 const double beta, 543 double *y); 544 545 /**************************************************/ 546 /** matrix conversion display *********************/ 547 /**************************************************/ 548 549 550 /**************************************************/ 551 /** matrix and vector display *********************/ 552 /**************************************************/ 553 554 void NM_dense_to_sparse(const NumericsMatrix* const A, NumericsMatrix* B, double threshold); 555 556 /** Copy a NumericsMatrix into another with dense storage. 557 \param A source matrix (any kind of storage) 558 \param B targeted matrix, must be dense with the same dimension as A 559 */ 560 int NM_to_dense(const NumericsMatrix* const A, NumericsMatrix* B); 561 562 /** Screen display of the matrix content stored as a double * array in Fortran style 563 \param m the matrix to be displayed 564 \param nRow the number of rows 565 \param nCol the number of columns 566 \param lDim the leading dimension of M 567 */ 568 void NM_dense_display_matlab(double * m, int nRow, int nCol, int lDim); 569 570 /** Screen display of the matrix content stored as a double * array in Fortran style 571 \param m the matrix to be displayed 572 \param nRow the number of rows 573 \param nCol the number of columns 574 \param lDim the leading dimension of M 575 */ 576 void NM_dense_display(double * m, int nRow, int nCol, int lDim); 577 578 /** Screen display of the vector content stored as a double * array 579 \param m the vector to be displayed 580 \param nRow the number of rows 581 */ 582 void NM_vector_display(double * m, int nRow); 583 584 585 /** Screen display of the matrix content 586 \param M the matrix to be displayed 587 */ 588 void NM_display(const NumericsMatrix* const M); 589 590 /** Screen display of the matrix storage 591 \param M the matrix to be displayed 592 */ 593 void NM_display_storageType(const NumericsMatrix* const M); 594 595 596 /** Screen display raw by raw of the matrix content 597 \param m the matrix to be displayed 598 */ 599 void NM_display_row_by_row(const NumericsMatrix* const m); 600 601 /**************************************************/ 602 /** matrix I/O *********************/ 603 /**************************************************/ 604 605 /** PrintInFile of the matrix content 606 \param M the matrix to be printed 607 \param filename the corresponding name of the file 608 */ 609 void NM_write_in_filename(const NumericsMatrix* const M, const char *filename); 610 611 /** Read in file of the matrix content 612 \param M the matrix to be read 613 \param filename the corresponding name of the file 614 */ 615 void NM_read_in_filename(NumericsMatrix* const M, const char *filename); 616 617 /** PrintInFile of the matrix content 618 \param M the matrix to be printed 619 \param file filename the corresponding file 620 */ 621 622 void NM_write_in_file(const NumericsMatrix* const M, FILE* file); 623 624 /** Read in file of the matrix content without performing memory allocation 625 \param M the matrix to be read 626 \param file the corresponding file 627 */ 628 void NM_read_in_file(NumericsMatrix* const M, FILE *file); 629 630 /** Create from file a NumericsMatrix with memory allocation 631 \param file the corresponding file 632 \return 0 if the matrix 633 */ 634 RawNumericsMatrix* NM_new_from_file(FILE *file); 635 RawNumericsMatrix* NM_new_from_filename(const char * filename); 636 637 /** NM_write_in_file_scilab of the matrix content 638 \param M the matrix to be printed 639 \param file the corresponding file 640 */ 641 void NM_write_in_file_scilab(const NumericsMatrix* const M, FILE* file); 642 643 /** NM_write_in_file_python of the matrix content 644 \param M the matrix to be printed 645 \param file the corresponding file 646 */ 647 void NM_write_in_file_python(const NumericsMatrix* const M, FILE* file); 648 649 /** Read in file for scilab of the matrix content 650 \param M the matrix to be read 651 \param file the corresponding file 652 */ 653 void NM_read_in_file_scilab(NumericsMatrix* const M, FILE *file); 654 655 /** Clear dense storage, if it is existent. 656 * \param[in,out] A a Numericsmatrix 657 */ 658 void NM_clearDense(NumericsMatrix* A); 659 660 /** Clear sparse block storage, if it is existent. 661 * \param[in,out] A a Numericsmatrix 662 */ 663 void NM_clearSparseBlock(NumericsMatrix* A); 664 665 /** Clear sparse data, if it is existent. 666 The linear solver parameters are also cleared. 667 * \param[in,out] A a Numericsmatrix 668 */ 669 void NM_clearSparse(NumericsMatrix* A); 670 671 /** Clear triplet storage, if it is existent. 672 * \param[in,out] A a Numericsmatrix 673 */ 674 void NM_clearTriplet(NumericsMatrix* A); 675 676 /** Clear half triplet storage, if it is existent. 677 * \param[in,out] A a Numericsmatrix 678 */ 679 void NM_clearHalfTriplet(NumericsMatrix* A); 680 681 /** Clear compressed column storage, if it is existent. 682 * \param[in,out] A a Numericsmatrix 683 */ 684 void NM_clearCSC(NumericsMatrix* A); 685 686 /** Clear transposed compressed column storage, if it is existent. 687 * \param[in,out] A a Numericsmatrix 688 */ 689 void NM_clearCSCTranspose(NumericsMatrix* A); 690 691 /** Clear compressed row storage, if it is existent. 692 * \param[in,out] A a Numericsmatrix 693 */ 694 void NM_clearCSR(NumericsMatrix* A); 695 696 /** Clear triplet, csc, csc transposed storage, if they are existent. 697 * Linear solver parameters are preserved. 698 * \param[in,out] A a Numericsmatrix 699 */ 700 void NM_clearSparseStorage(NumericsMatrix *A); 701 702 703 704 /** XXXXXX: to be rewritten Direct computation of the solution of a real system of linear 705 * equations: A x = b. The factorized matrix A is kept for future solve. 706 * If A is already factorized, the solve the linear system from it 707 * \warning this is not enable for all the solvers, your mileage may vary 708 * \param[in,out] A a NumericsMatrix. On a dense factorisation 709 * A.iWork is initialized. 710 * \param[in,out] b pointer on a dense vector of size A->size1 711 * \param keep if set to NM_KEEP_FACTORS, keep all the info related to the factorization to 712 * allow for future solves. If A is already factorized, just solve the linear 713 * system. If set to NM_PRESERVE, preserve the original matrix (just used in 714 * the dense case). if NM_NONE, discard everything. 715 * \return 0 if successful, else the error is specific to the backend solver 716 * used 717 */ 718 719 /* LU factorization of the matrix. If the matrix has already been 720 * factorized (i.e if NM_LU_factorized(A) return true), nothing is 721 * done. To force a new factorization one has to set factorization 722 * flag to false : NM_set_LU_factorized(A, false) before the call to 723 * NM_LU_factorize. 724 * If the matrix is preserved, that means that a call to 725 * NM_preserve(A) has been done before the call to NM_LU_factorize, 726 * it is not destroyed, but the factorized part remains accessible for 727 * subsequent calls to NM_LU_solve. 728 * If the matrix is not preserved, then it is replaced by the 729 * factorized part. 730 * \param[in] A the NumericsMatrix 731 * \return an int, 0 means the matrix has been factorized. */ 732 int NM_LU_factorize(NumericsMatrix* A); 733 int NM_Cholesky_factorize(NumericsMatrix* A); 734 int NM_LDLT_factorize(NumericsMatrix* A); 735 736 /* Solve linear system with multiple right hand size. A call to 737 * NM_LU_factorize is done at the beginning. 738 739 * \param[in] A the NumericsMatrix. A is not destroyed if it has 740 * been preserved by a call to NM_preserve(A). 741 742 * \param[in,out] b the right hand size which is a pointer on a 743 * matrix of double. It is replaced by the solutions 744 745 * \param[in] nrhs the number of right hand side. 746 * \return 0 if the solve succeeded. 747 */ 748 int NM_LU_solve(NumericsMatrix* A, double *b, unsigned int nrhs); 749 int NM_LU_solve_matrix_rhs(NumericsMatrix* Ao, NumericsMatrix* B); 750 int NM_Cholesky_solve(NumericsMatrix* A, double *b, unsigned int nrhs); 751 int NM_Cholesky_solve_matrix_rhs(NumericsMatrix* Ao, NumericsMatrix* B); 752 int NM_LDLT_solve(NumericsMatrix* A, double *b, unsigned int nrhs); 753 754 755 int NM_gesv_expert(NumericsMatrix* A, double *b, unsigned keep); 756 int NM_posv_expert(NumericsMatrix* A, double *b, unsigned keep); 757 758 int NM_gesv_expert_multiple_rhs(NumericsMatrix* A, double *b, unsigned int n_rhs, unsigned keep); 759 760 761 762 /** Computation of the inverse of a NumericsMatrix A usinf NM_gesv_expert 763 * \param[in,out] A a NumericsMatrix. 764 * \return the matrix inverse. 765 */ 766 NumericsMatrix* NM_LU_inv(NumericsMatrix* A); 767 768 int NM_inverse_diagonal_block_matrix_in_place(NumericsMatrix* A); 769 770 /** Direct computation of the solution of a real system of linear 771 * equations: A x = b. 772 * \param[in,out] A a NumericsMatrix. On a dense factorisation 773 * A.iWork is initialized. 774 * \param[in,out] b pointer on a dense vector of size A->size1 775 * \param preserve preserve the original matrix data. Only useful in the 776 * dense case, where the LU factorization is done in-place. 777 * \return 0 if successful, else the error is specific to the backend solver 778 * used 779 */ NM_gesv(NumericsMatrix * A,double * b,bool preserve)780 static inline int NM_gesv(NumericsMatrix* A, double *b, bool preserve) 781 { 782 return NM_gesv_expert(A, b, preserve ? NM_PRESERVE : NM_NONE); 783 } 784 785 /** Computation of the inverse of a NumericsMatrix A usinf NM_gesv_expert 786 * \param[in,out] A a NumericsMatrix. 787 * \return the matrix inverse. 788 */ 789 NumericsMatrix* NM_gesv_inv(NumericsMatrix* A); 790 791 /** Set the linear solver 792 * \param A the matrix 793 * \param solver_id id of the solver 794 */ 795 void NM_setSparseSolver(NumericsMatrix* A, NSM_linear_solver solver_id); 796 797 /** Get Matrix internal data with initialization if needed. 798 * \param[in,out] A a NumericsMatrix. 799 * \return a pointer on internal data. 800 */ 801 NumericsMatrixInternalData* NM_internalData(NumericsMatrix* A); 802 803 /** Allocate the internalData structure (but not its content!) 804 * \param M the matrix to modify 805 */ 806 void NM_internalData_new(NumericsMatrix* M); 807 808 /** Copy the internalData structure 809 * \param M the matrix to modify 810 */ 811 void NM_internalData_copy(const NumericsMatrix* const A, NumericsMatrix* B ); 812 813 /** Integer work vector initialization, if needed. 814 * \param[in,out] A pointer on a NumericsMatrix. 815 * \param[in] size number of element to allocate 816 * \param[in] sizeof_elt of an element in bytes (result of sizeof for instance) 817 * \return pointer on A->iWork allocated space of with the right size 818 */ 819 void* NM_iWork(NumericsMatrix *A, size_t size, size_t sizeof_elt); 820 821 /** Double workspace initialization, if needed. 822 * \param[in,out] A pointer on a NumericsMatrix. 823 * \param[in] size the size of needed space. 824 * \return pointer on A->dWork allocated space of with the right size 825 */ 826 double* NM_dWork(NumericsMatrix *A, int size); 827 828 /** Add a constant term to the diagonal elements, when the block of the SBM 829 * are 3x3 830 * \param M the matrix 831 * \param alpha the term to add 832 */ 833 void NM_add_to_diag3(NumericsMatrix* M, double alpha); 834 835 /** Add a constant term to the diagonal elements, when the block of the SBM 836 * are 5x5 837 * \param M the matrix 838 * \param alpha the term to add 839 */ 840 void NM_add_to_diag5(NumericsMatrix* M, double alpha); 841 842 /** Add two matrices with coefficients C = alpha*A + beta*B 843 * \param alpha the first coefficient 844 * \param A the first matrix 845 * \param beta the second coefficient 846 * \param B the second matrix 847 * \return C a new NumericsMatrix 848 */ 849 RawNumericsMatrix * NM_add(double alpha, NumericsMatrix* A, double beta, NumericsMatrix* B); 850 851 /** Multiply a matrix with a double alpha*A --> A 852 * \param alpha the coefficient 853 * \param A the matrix 854 */ 855 void NM_scal(double alpha, NumericsMatrix* A); 856 857 /** assert that a NumericsMatrix has the right structure given its type 858 * \param type expected type 859 * \param M the matrix to check 860 */ NM_assert(NM_types type,NumericsMatrix * M)861 static inline void NM_assert(NM_types type, NumericsMatrix* M) 862 { 863 #ifndef NDEBUG 864 assert(M && "NM_assert :: the matrix is NULL"); 865 assert(M->storageType == type && "NM_assert :: the matrix has the wrong type"); 866 switch(type) 867 { 868 case NM_DENSE: 869 assert(M->matrix0); 870 break; 871 case NM_SPARSE_BLOCK: 872 assert(M->matrix1); 873 break; 874 case NM_SPARSE: 875 assert(M->matrix2); 876 break; 877 default: 878 assert(0 && "NM_assert :: unknown storageType"); 879 } 880 #endif 881 } 882 883 /** Check the matrix (the sparse format for now) 884 * \param A the matrix to check 885 * \return 0 if the matrix storage is fine, 1 if not*/ 886 int NM_check(const NumericsMatrix* const A); 887 888 /** Compute the 1-norm of a sparse matrix = max (sum (abs (A))), largest column sum of a matrix (the sparse format for now) 889 * \param A the matrix 890 * \return the norm*/ 891 double NM_norm_1(NumericsMatrix* const A); 892 893 /** Compute the inf-norm of a sparse matrix = max (sum (abs (A^T))), largest row sum of a matrix (the sparse format for now) 894 * \param A the matrix 895 * \return the norm*/ 896 double NM_norm_inf(NumericsMatrix* const A); 897 898 int NM_is_symmetric(NumericsMatrix* A); 899 double NM_symmetry_discrepancy(NumericsMatrix* A); 900 901 902 /** Pass a NumericsMatrix through swig typemaps. 903 * This is only useful in python. 904 * \param A the matrix 905 * \return a NumericsMatrix 906 */ NM_convert(NumericsMatrix * A)907 static inline NumericsMatrix* NM_convert(NumericsMatrix* A) 908 { 909 return A; 910 } 911 912 /** Compute the maximum eigenvalue with the iterated power method 913 * \param A the matrix 914 * \return the maximum eigenvalue*/ 915 double NM_iterated_power_method(NumericsMatrix* A, double tol, int itermax); 916 917 /* Compute the maximum values by columns 918 * \param A the matrix 919 * \param max the vector of max that must be preallocated 920 * \return info 921 */ 922 int NM_max_by_columns(NumericsMatrix *A, double * max); 923 924 /* Compute the maximum values by rows 925 * \param A the matrix 926 * \param max the vector of max that must be preallocated 927 * \return info 928 */ 929 int NM_max_by_rows(NumericsMatrix *A, double * max); 930 931 /* Compute the maximum absolute values by columns 932 * \param A the matrix 933 * \param max the vector of max that must be preallocated 934 * \return info 935 */ 936 int NM_max_abs_by_columns(NumericsMatrix *A, double * max); 937 938 /* Compute the maximum absolute values by rows 939 * \param A the matrix 940 * \param max the vector of max that must be preallocated 941 * \return info 942 */ 943 int NM_max_abs_by_rows(NumericsMatrix *A, double * max); 944 945 /* Compute the balancing matrices for a given matrix by iteration 946 * \param A the matrix 947 * \param tol tolerance on the balanced matrix 948 * \param itermax max number of iterations 949 * \param alloated structure for the balancing matrices and the balanced matrix 950 * \return 0 if succeed. 951 */ 952 int NM_compute_balancing_matrices(NumericsMatrix* A, double tol, int itermax, BalancingMatrices * B); 953 954 /* Create a Balancing Matrices structure 955 * \param A the matrix to be balanced 956 */ 957 BalancingMatrices * NM_BalancingMatrices_new(NumericsMatrix* A); 958 959 960 /* free a Balancing Matrices structure 961 */ 962 BalancingMatrices * NM_BalancingMatrices_free(BalancingMatrices* A); 963 964 /* Reset the version of a NM_types storage. 965 *\param M the NumericsMatrix, 966 *\param id the NM_types storage 967 */ 968 void NM_reset_version(NumericsMatrix* M, NM_types id); 969 970 /* Reset versions of all storages. 971 *\param M the NumericsMatrix 972 */ 973 void NM_reset_versions(NumericsMatrix* M); 974 975 976 #ifdef WITH_OPENSSL 977 /* Compute sha1 hash of matrix values. Matrices of differents size and same 978 * values have the same hash. 979 * \param[in] A the matrix 980 * \param[in,out] digest an allocated space of size SHA_DIGEST_LENGTH 981 */ 982 void NM_compute_values_sha1(NumericsMatrix* A, unsigned char * digest); 983 984 /* Get stored sha1 hash of a matrix. 985 * \param[in] A the matrix 986 * \return a pointer on the current raw sha1 hash 987 */ 988 unsigned char* NM_values_sha1(NumericsMatrix* A); 989 990 /* Set sha1 hash of a matrix. Matrices of differents size and same 991 * values have the same hash. 992 * \param[in] A the matrix 993 */ 994 void NM_set_values_sha1(NumericsMatrix* A); 995 996 /* Clear sha1 hash of a matrix. 997 * \param[in] A the matrix 998 */ 999 void NM_clear_values_sha1(NumericsMatrix* A); 1000 1001 /* Check if matrix has beend modified after a previous NM_set_values_sha1. 1002 * \param[in] A the NumericsMatrix 1003 * \return true if the matrix is the same 1004 */ 1005 bool NM_check_values_sha1(NumericsMatrix* A); 1006 1007 /* Compare two matrices with sha1. NM_set_values_sha1 must be called 1008 * on the two matrices before. 1009 * \param[in] A a NumericsMatrix 1010 * \param[in] B a NumericsMatrix 1011 * \return true if the matrices have the same values 1012 */ 1013 bool NM_equal_values_sha1(NumericsMatrix* A, NumericsMatrix* B); 1014 1015 #endif 1016 1017 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 1018 } 1019 #endif 1020 1021 #endif 1022