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 SparseBlockMatrix_H 20 #define SparseBlockMatrix_H 21 22 #include <stdio.h> // for size_t, FILE 23 #include "CSparseMatrix.h" // for CSparseMatrix 24 #include "NumericsFwd.h" // for SparseBlockStructuredMatrix, SparseBlockC... 25 26 #include "SiconosConfig.h" // for BUILD_AS_CPP // IWYU pragma: keep 27 #include "NumericsDataVersion.h" // versioning 28 29 /*!\file SparseBlockMatrix.h 30 \brief Structure definition and functions related to 31 SparseBlockStructuredMatrix 32 33 */ 34 35 /** Structure to store sparse block matrices with square diagonal 36 blocks. 37 38 Note: the sparse format is the same as the one used by Boost C++ 39 library to store compressed sparse row matrices. The same member 40 names have been adopted in order to simplify usage from Siconos 41 Kernel : filled1, filled2, index1_data, index2_data. 42 Reference : 43 http://ublas.sourceforge.net/refdoc/classboost_1_1numeric_1_1ublas_1_1compressed__matrix.html 44 45 \param nbblocks the total number of non null blocks 46 \param **block : *block contains the double values of one block in 47 Fortran storage (column by column) **block is 48 the list of non null blocks 49 \param blocknumber0 the first dimension of the block matrix 50 (number of block rows) 51 \param blocknumber1 the second dimension of the block matrix 52 (number of block columns) 53 \param *blocksize0 the list of sums of the number of rows of the 54 first column of blocks of M: blocksize0[i] = blocksize0[i-1] + 55 ni, ni being the number of rows of the block at row i 56 \param *blocksize1 the list of sums of the number of columns of the 57 first row of blocks of M: blocksize1[i] = blocksize1[i-1] + ni, 58 ni being the number of columns of the block at column i 59 \param filled1 index of the last non empty line + 1 60 \param filled2 number of non null blocks 61 \param index1_data index1_data is of size equal to number of non 62 empty lines + 1. A block with number blockNumber inside a row 63 numbered rowNumber verify index1_data[rowNumber]<= blockNumber 64 <index1_data[rowNumber+1]` 65 66 \param index2_data index2_data is of size filled2 67 index2_data[blockNumber] -> columnNumber. 68 69 70 Related functions: SBM_gemv(), SBM_row_prod(), SBM_clear(), 71 SBM_print, SBM_diagonal_block_index() 72 * If we consider the matrix M and the right-hand-side q defined as 73 * 74 * \f$ 75 * M=\left[\begin{array}{cccc|cc|cc} 76 * 1 & 2 & 0 & 4 & 3 &-1 & 0 & 0\\ 77 * 2 & 1 & 0 & 0 & 4 & 1 & 0 & 0\\ 78 * 0 & 0 & 1 &-1 & 0 & 0 & 0 & 0\\ 79 * 5 & 0 &-1 & 6 & 0 & 6 & 0 & 0\\ 80 * \hline 81 * 0 & 0 & 0 & 0 & 1 & 0 & 0 & 5\\ 82 * 0 & 0 & 0 & 0 & 0 & 2 & 0 & 2\\ 83 * \hline 84 * 0 & 0 & 2 & 1 & 0 & 0 & 2 & 2\\ 85 * 0 & 0 & 2 & 2 & 0 & 0 & -1 & 2\\ 86 * \end{array}\right] \quad, q=\left[\begin{array}{c}-1\\-1\\0\\-1\\\hline 1\\0\\\hline -1\\2\end{array}\right]. 87 * \f$ 88 * 89 * then 90 * - the number of non null blocks is 6 (nbblocks=6) 91 * - the number of rows of blocks is 3 (blocknumber0 =3) and the 92 number of columns of blocks is 3 (blocknumber1 =3) 93 * - the vector blocksize0 is equal to {4,6,8} and the vector 94 blocksize1 is equal to {4,6,8} 95 * - the integer filled1 is equal to 4 96 * - the integer filled2 is equal to 6 97 * - the vector index1_data is equal to {0,2,4,6} 98 * - the vector index2_data is equal to {0,1,1,2,0,2} 99 * - the block contains all non null block matrices stored in Fortran 100 order (column by column) as 101 * block[0] = {1,2,0,5,2,1,0,0,0,0,1,-1,4,0,-1,6} 102 * block[1] = {3,4,0,0,-1,1,0,6} 103 * ... 104 * block[5] = {2,-1,2,2} 105 */ 106 107 108 struct SparseBlockStructuredMatrix 109 { 110 /* the number of non null blocks */ 111 unsigned int nbblocks; 112 double **block; 113 /* the number of rows of blocks */ 114 unsigned int blocknumber0; 115 /* the number of columns of blocks */ 116 unsigned int blocknumber1; 117 /* the vector of cumulated row sizes of blocks */ 118 unsigned int *blocksize0; 119 /* the vector of cumulated column sizes of blocks */ 120 unsigned int *blocksize1; 121 122 /* the index of the last non empty line + 1 */ 123 size_t filled1; 124 /* the size of index2_data that corresponds to the number of non null blocks*/ 125 size_t filled2; 126 127 size_t *index1_data; 128 size_t *index2_data; 129 130 /* the indices of the diagonal blocks */ 131 unsigned int * diagonal_blocks; 132 133 NumericsDataVersion version; /**< version of storage */ 134 }; 135 136 struct SparseBlockCoordinateMatrix 137 { 138 /** number of blocks */ 139 unsigned int nbblocks; 140 141 /** number of rows */ 142 unsigned int blocknumber0; 143 144 /** number of columns */ 145 unsigned int blocknumber1; 146 147 /** block pointers */ 148 double **block; 149 150 /** cumulative number of rows in blocks */ 151 unsigned int *blocksize0; 152 153 /** cumulative number of columns in blocks */ 154 155 unsigned int *blocksize1; 156 157 /** row indices */ 158 unsigned int *row; 159 160 /** column indices */ 161 unsigned int *column; 162 }; 163 164 struct SparseBlockStructuredMatrixPred 165 { 166 int nbbldiag; 167 int **indic; 168 int **indicop; 169 double **submatlcp; 170 double **submatlcpop; 171 int **ipiv; 172 int *sizesublcp; 173 int *sizesublcpop; 174 double **subq; 175 double **bufz; 176 double **newz; 177 double **workspace; 178 }; 179 180 struct SBM_index_by_column 181 { 182 /* the index of the last non empty line + 1 */ 183 size_t filled3; 184 /* the size of index2_data that corresponds of the number of non null blocks*/ 185 size_t filled4; 186 187 size_t * index3_data; 188 size_t * index4_data; 189 size_t * blockMap; 190 }; 191 192 193 194 #define NUMERICS_SBM_FREE_BLOCK 4 195 #define NUMERICS_SBM_FREE_SBM 8 196 197 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 198 extern "C" 199 { 200 #endif 201 202 /** Creation of an empty Sparse Block Matrix. 203 * \return a pointer on allocated and initialized space 204 */ 205 SparseBlockStructuredMatrix* SBM_new(void); 206 207 /** set Sparse Block Matrix. fields to NULL 208 * \param sbm a matrix 209 */ 210 void SBM_null(SparseBlockStructuredMatrix* sbm); 211 212 /** SparseMatrix - vector product y = alpha*A*x + beta*y 213 \param[in] sizeX dim of the vectors x 214 \param[in] sizeY dim of the vectors y 215 \param[in] alpha coefficient 216 \param[in] A the matrix to be multiplied 217 \param[in] x the vector to be multiplied 218 \param[in] beta coefficient 219 \param[in,out] y the resulting vector 220 */ 221 void SBM_gemv(unsigned int sizeX, unsigned int sizeY, 222 double alpha, const SparseBlockStructuredMatrix* const A, 223 const double* x, double beta, double* y); 224 225 /** SparseMatrix - vector product y = A*x + y for block of size 3x3 226 \param[in] sizeX dim of the vectors x 227 \param[in] sizeY dim of the vectors y 228 \param[in] A the matrix to be multiplied 229 \param[in] x the vector to be multiplied 230 \param[in,out] y the resulting vector 231 */ 232 void SBM_gemv_3x3(unsigned int sizeX, unsigned int sizeY, 233 const SparseBlockStructuredMatrix* const A, 234 double* const x, double* y); 235 236 /** SparseBlockStructuredMatrix - SparseBlockStructuredMatrix product C = alpha*A*B + beta*C 237 * The routine has to be used with precaution. The allocation of C is not done 238 * since we want to add beta*C. We assume that the structure and the allocation 239 * of the matrix C are right. Especially: 240 * - the blocks C(i,j) must exists 241 * - the sizes of blocks must be consistent 242 * - no extra block must be present in C 243 * \param[in] alpha coefficient 244 * \param[in] A the matrix to be multiplied 245 * \param[in] B the matrix to be multiplied 246 * \param[in] beta coefficient 247 * \param[in,out] C the resulting matrix 248 */ 249 void SBM_gemm_without_allocation(double alpha, const SparseBlockStructuredMatrix* const A, 250 const SparseBlockStructuredMatrix* const B, 251 double beta, SparseBlockStructuredMatrix* C); 252 253 /** SparseBlockStructuredMatrix - SparseBlockStructuredMatrix multiplication C = A *B 254 Correct allocation is performed 255 \param[in] A the matrix to be multiplied 256 \param[in] B the matrix to be multiplied 257 \return C the resulting matrix 258 */ 259 SparseBlockStructuredMatrix* SBM_multiply(const SparseBlockStructuredMatrix* const A, const SparseBlockStructuredMatrix* const B); 260 261 /** Perform the allocation of a zero matrix that is compatible qith multiplication 262 */ 263 SparseBlockStructuredMatrix* SBM_zero_matrix_for_multiply(const SparseBlockStructuredMatrix* const A, const SparseBlockStructuredMatrix* const B); 264 265 /** SparseBlockStructuredMatrix - SparseBlockStructuredMatrix addition C = alpha*A + beta*B 266 \param[in] A the matrix to be added 267 \param[in] B the matrix to be added 268 \param[in] alpha coefficient 269 \param[in] beta coefficient 270 \return C the resulting matrix 271 */ 272 SparseBlockStructuredMatrix * SBM_add(SparseBlockStructuredMatrix * A, SparseBlockStructuredMatrix * B, double alpha, double beta); 273 274 /** SparseBlockStructuredMatrix - SparseBlockStructuredMatrix addition C = alpha*A + beta*B + gamma*C without allocation 275 We assume that C has the correct structure 276 \param[in] A the matrix to be added 277 \param[in] B the matrix to be added 278 \param[in] alpha coefficient 279 \param[in] beta coefficient 280 \param[in] gamma coefficient 281 \param[in,out] C the resulting matrix 282 */ 283 284 void SBM_add_without_allocation(SparseBlockStructuredMatrix * A, SparseBlockStructuredMatrix * B, 285 double alpha, double beta, 286 SparseBlockStructuredMatrix * C, 287 double gamma); 288 289 290 /** Multiply a matrix with a double alpha*A --> A 291 * \param alpha the coefficient 292 * \param A the matrix 293 */ 294 void SBM_scal(double alpha, SparseBlockStructuredMatrix * A); 295 296 297 /** Row of a SparseMatrix - vector product y = rowA*x or y += rowA*x, rowA being a row of blocks of A 298 \param[in] sizeX dim of the vector x 299 \param[in] sizeY dim of the vector y 300 \param[in] currentRowNumber number of the required row of blocks 301 \param[in] A the matrix to be multiplied 302 \param[in] x the vector to be multiplied 303 \param[in,out] y the resulting vector 304 \param[in] init = 0 for y += Ax, =1 for y = Ax 305 */ 306 void SBM_row_prod(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix* const A, const double* const x, double* y, int init); 307 308 /** Row of a SparseMatrix - vector product y = rowA*x or y += rowA*x, rowA being a row of blocks of A 309 \param[in] sizeX dim of the vector x 310 \param[in] sizeY dim of the vector y 311 \param[in] currentRowNumber number of the required row of blocks 312 \param[in] A the matrix to be multiplied 313 \param[in] x the vector to be multiplied 314 \param[in,out] y the resulting vector 315 \param[in] init = 0 for y += Ax, =1 for y = Ax 316 */ 317 void SBM_row_prod_no_diag(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix* const A, const double* const x, double* y, int init); 318 319 /** Row of a SparseMatrix - vector product y = rowA*x or y += rowA*x, rowA being a row of blocks of A of size 3x3 320 \param[in] sizeX dim of the vector x 321 \param[in] sizeY dim of the vector y 322 \param[in] currentRowNumber number of the required row of blocks 323 \param[in] A the matrix to be multiplied 324 \param[in] x the vector to be multiplied 325 \param[in,out] y the resulting vector 326 */ 327 void SBM_row_prod_no_diag_3x3(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix* const A, double* const x, double* y); 328 void SBM_row_prod_no_diag_1x1(unsigned int sizeX, unsigned int sizeY, unsigned int currentRowNumber, const SparseBlockStructuredMatrix* const A, double* const x, double* y); 329 330 331 void SBM_extract_component_3x3(const SparseBlockStructuredMatrix* const A, 332 SparseBlockStructuredMatrix* B, 333 unsigned int *row_components, unsigned int row_components_size, 334 unsigned int *col_components, unsigned int col_components_size); 335 336 /** Destructor for SparseBlockStructuredMatrix objects 337 \param blmat SparseBlockStructuredMatrix the matrix to be destroyed. 338 */ 339 void SBM_clear(SparseBlockStructuredMatrix * blmat); 340 341 /** To free a SBM matrix (for example allocated by NM_new_from_file). 342 * \param[in] A the SparseBlockStructuredMatrix that mus be de-allocated. 343 * \param[in] level use NUMERICS_SBM_FREE_BLOCK | NUMERICS_SBM_FREE_SBM 344 */ 345 void SBMfree(SparseBlockStructuredMatrix* A, unsigned int level); 346 347 /** Screen display of the matrix content 348 \param m the matrix to be displayed 349 */ 350 void SBM_print(const SparseBlockStructuredMatrix* const m); 351 352 /** print in file of the matrix content 353 \param m the matrix to be displayed 354 \param file the corresponding file 355 */ 356 void SBM_write_in_file(const SparseBlockStructuredMatrix* const m, FILE* file); 357 358 /** read in file of the matrix content without performing memory allocation 359 \param M the matrix to be displayed 360 \param file the corresponding name of the file 361 */ 362 void SBM_read_in_file(SparseBlockStructuredMatrix* const M, FILE *file); 363 364 /** Create from file a SparseBlockStructuredMatrix with memory allocation 365 \param file the corresponding name of the file 366 \return the matrix to be displayed 367 */ 368 SparseBlockStructuredMatrix* SBM_new_from_file(FILE *file); 369 370 /** print in file of the matrix content in Scilab format for each block 371 \param M the matrix to be displayed 372 \param file the corresponding file 373 */ 374 void SBM_write_in_fileForScilab(const SparseBlockStructuredMatrix* const M, FILE* file); 375 376 377 /** print in file of the matrix content 378 \param M the matrix to be displayed 379 \param filename the corresponding file 380 */ 381 void SBM_write_in_filename(const SparseBlockStructuredMatrix* const M, const char *filename); 382 383 /** read in file of the matrix content 384 \param M the matrix to be displayed 385 \param filename the corresponding name of the file 386 */ 387 void SBM_read_in_filename(SparseBlockStructuredMatrix* const M, const char *filename); 388 389 /** Destructor for SparseBlockStructuredMatrixPred objects 390 * \param blmatpred SparseBlockStructuredMatrix, the matrix to be destroyed. 391 */ 392 void SBM_clear_pred(SparseBlockStructuredMatrixPred *blmatpred); 393 394 /** Compute the indices of blocks of the diagonal block 395 \param M the SparseBlockStructuredMatrix matrix 396 \return the indices for all the rows 397 */ 398 unsigned int * SBM_diagonal_block_indices(SparseBlockStructuredMatrix* const M); 399 400 /** Find index of the diagonal block in a row 401 \param M the SparseBlockStructuredMatrix matrix 402 \param row the row of the required block 403 \return pos the position of the block 404 */ 405 unsigned int SBM_diagonal_block_index(SparseBlockStructuredMatrix* const M, unsigned int row); 406 407 /** insert an entry into a SparseBlockStructuredMatrix. 408 * This method is expensive in terms of memory management. For a lot of entries, use 409 * an alternative 410 * \param M the SparseBlockStructuredMatrix 411 * \param i row index 412 * \param j column index 413 * \param val the value to be inserted. 414 */ 415 int SBM_entry(SparseBlockStructuredMatrix* M, unsigned int row, unsigned int col, double val); 416 417 /** get the element of row i and column j of the matrix M 418 \param M the SparseBlockStructuredMatrix matrix 419 \param row the row index 420 \param col the column index 421 \return the value 422 */ 423 double SBM_get_value(const SparseBlockStructuredMatrix* const M, unsigned int row, unsigned int col); 424 425 /** Copy of a SBM A into B 426 \param[in] A the SparseBlockStructuredMatrix matrix to be copied 427 \param[out] B the SparseBlockStructuredMatrix matrix copy of A 428 \param[in] copyBlock if copyBlock then the content of block are copied, else only the pointers are copied. 429 \return 0 if ok 430 */ 431 int SBM_copy(const SparseBlockStructuredMatrix* const A, SparseBlockStructuredMatrix* B, unsigned int copyBlock); 432 433 434 /** Transpose by copy of a SBM A into B 435 \param[in] A the SparseBlockStructuredMatrix matrix to be copied 436 \param[out] B the SparseBlockStructuredMatrix matrix copy of transpose A 437 \return 0 if ok 438 */ 439 int SBM_transpose(const SparseBlockStructuredMatrix* const A, SparseBlockStructuredMatrix* B); 440 441 /** Inverse (in place) a square diagonal block matrix 442 \param[in,out] M the SparseBlockStructuredMatrix matrix to be inversed 443 \param ipiv worksapce for storign ipiv 444 \return 0 ik ok 445 */ 446 int SBM_inverse_diagonal_block_matrix_in_place(const SparseBlockStructuredMatrix* M, int* ipiv); 447 448 449 /** Copy a SBM into a Dense Matrix 450 \param[in] A the SparseBlockStructuredMatrix matrix 451 \param[in] denseMat pointer on the filled dense Matrix 452 */ 453 void SBM_to_dense(const SparseBlockStructuredMatrix* const A, double *denseMat); 454 455 /** Copy a SBM into a Sparse (CSR) Matrix 456 \param[in] A the SparseBlockStructuredMatrix matrix 457 \param[in] outSparseMat pointer on the filled sparse Matrix 458 \return 0 if ok 459 */ 460 int SBM_to_sparse(const SparseBlockStructuredMatrix* const A, CSparseMatrix *outSparseMat); 461 462 /** initMemory of a Sparse (CSR) Matrix form a SBM matrix 463 \param[in] A the SparseBlockStructuredMatrix matrix 464 \param[in] sparseMat pointer on the initialized sparse Matrix 465 \return 0 if ok 466 */ 467 int SBM_to_sparse_init_memory(const SparseBlockStructuredMatrix* const A, CSparseMatrix *sparseMat); 468 469 /**Copy a block row of the SBM into a Dense Matrix 470 \param[in] A the SparseBlockStructuredMatrix matrix to be inversed. 471 \param[in] row the block row index copied. 472 \param[in] denseMat pointer on the filled dense Matrix. 473 \param[in] rowPos line pos in the dense matrix. 474 \param[in] rowNb total number of line of the dense matrix. 475 The number of line copied is contained in M. 476 \ 477 */ 478 void SBM_row_to_dense(const SparseBlockStructuredMatrix* const A, int row, double *denseMat, int rowPos, int rowNb); 479 480 /* 481 * \param [in] rowIndex: permutation: the row numC of C is the row rowIndex[numC] of A. 482 * \param [in] A The source SBM. 483 * \param [out] C The target SBM. It assumes the structure SBM has been allocated. 484 * The memory allocation for its menber is done inside. 485 * NB : The blocks are not copied. 486 */ 487 void SBM_row_permutation(unsigned int *rowIndex, SparseBlockStructuredMatrix* A, SparseBlockStructuredMatrix* C); 488 489 /* 490 * \param [in] colIndex: permutation: the col numC of C is the col colIndex[numC] of A. 491 * \param [in] A The source SBM. 492 * \param [out] C The target SBM. It assumes the structure SBM has been allocated. 493 * The memory allocation for its menber is done inside. 494 * NB : The blocks are not copied. 495 */ 496 void SBM_column_permutation(unsigned int *colIndex, SparseBlockStructuredMatrix* A, SparseBlockStructuredMatrix* C); 497 498 void SBCM_null(SparseBlockCoordinateMatrix* MC); 499 500 SparseBlockCoordinateMatrix* SBCM_new(void); 501 502 503 504 505 /** allocate a SparseBlockCoordinateMatrix from a list of 3x3 506 * blocks 507 * \param[in] m the number of rows 508 * \param[in] n the number of colums 509 * \param[in] nbblocks the number of blocks 510 * \param[in] row a pointer to row of each block 511 * \param[in] column a pointer to column of each block 512 * \param[in] block a pointer to each block 513 * \return a pointer to a SparseBlockCoordinateMatrix structure 514 */ 515 SparseBlockCoordinateMatrix* SBCM_new_3x3(unsigned int m, unsigned int n, 516 unsigned int nbblocks, 517 unsigned int *row, unsigned int *column, double *block); 518 519 520 /** free allocated memory in newSparseBlockCoordinateMatrix functions 521 * \param[in] MC matrix pointer */ 522 void SBCM_free_3x3(SparseBlockCoordinateMatrix *MC); 523 524 /** copy a SparseBlockCoordinateMatrix to a SparseBlockStructuredMatrix 525 * \param[in] MC the SparseBlockCoordinateMatrix matrix 526 * \return a pointer to a SparseBlockCoordinateMatrix structure 527 */ 528 SparseBlockStructuredMatrix* SBCM_to_SBM(SparseBlockCoordinateMatrix* MC); 529 530 531 /** free a SparseBlockStructuredMatrix created with SBCM_to_SBM 532 * \param[in,out] M a SparseBlockStructuredMatrix to free*/ 533 void SBM_free_from_SBCM(SparseBlockStructuredMatrix* M); 534 535 536 /** Copy a Sparse Matrix into a SBM, with fixed blocksize 537 \param[in] blocksize the blocksize 538 \param[in] sparseMat pointer on the Sparse Matrix 539 \param[in,out] outSBM pointer on an empty SparseBlockStructuredMatrix 540 \return 0 in ok 541 */ 542 int SBM_from_csparse(int blocksize, const CSparseMatrix* const sparseMat, SparseBlockStructuredMatrix* outSBM); 543 544 545 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 546 } 547 #endif 548 549 #endif /* NSSPACK_H */ 550