1 /* Copyright 2021 INRIA. 2 * Siconos is a program dedicated to modeling, simulation and control 3 * of non smooth dynamical systems. 4 * Siconos is a free software; you can redistribute it and/or modify 5 * it under the terms of the GNU Lesser General Public License as published by 6 * the Free Software Foundation; either version 2 of the License, or 7 * (at your option) any later version. 8 * Siconos is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU Lesser General Public License for more details. 12 * 13 * You should have received a copy of the GNU Lesser General Public License 14 * along with Siconos; if not, write to the Free Software 15 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 16 * 17 * Contact: Vincent ACARY, siconos-team@lists.gforge.inria.fr 18 */ 19 20 #ifndef SparseMatrix_H 21 #define SparseMatrix_H 22 23 #include <stdio.h> 24 25 /*!\file CSparseMatrix.h 26 \brief Structure definition and functions related to sparse matrix storage in Numerics 27 */ 28 29 #include "SiconosConfig.h" 30 31 /* Compile-time assertion: users of SparseMatrix.h must have CS_LONG 32 * set if and only if SICONOS_INT64 is also set. If it is unset, it is 33 * set here. */ 34 #ifdef SICONOS_INT64 35 #ifndef CS_LONG 36 #define CS_LONG 37 #endif 38 #else 39 #ifdef CS_LONG 40 #error "CS_LONG (set) does not correspond with SICONOS_INT64 (unset)" 41 #endif 42 #endif 43 44 #ifndef CS_INT 45 46 /* From cs.h */ 47 #ifdef CS_LONG 48 #define CS_INT long 49 #else 50 #define CS_INT int 51 #endif 52 53 /* Treat CXSparse structs as opaque types. Users may #include "cs.h" 54 * to use them outside Siconos. */ 55 struct cs_dl_sparse; 56 struct cs_di_sparse; 57 struct cs_dl_symbolic; 58 struct cs_di_symbolic; 59 struct cs_dl_numeric; 60 struct cs_di_numeric; 61 typedef struct cs_dl_symbolic cs_dls; 62 typedef struct cs_di_symbolic cs_dis; 63 typedef struct cs_dl_numeric cs_dln; 64 typedef struct cs_di_numeric cs_din; 65 66 #ifdef SICONOS_INT64 // SWIG gives syntax error for CS_NAME(_sparse) 67 #ifndef css 68 #define css cs_dls 69 #endif 70 #ifndef csn 71 #define csn cs_dln 72 #endif 73 #else 74 #ifndef css 75 #define css cs_dis 76 #endif 77 #ifndef csn 78 #define csn cs_din 79 #endif 80 #endif 81 82 #endif 83 84 #ifdef SICONOS_INT64 // SWIG gives syntax error for CS_NAME(_sparse) 85 typedef struct cs_dl_sparse CSparseMatrix; 86 #else 87 typedef struct cs_di_sparse CSparseMatrix; 88 #endif 89 90 91 /* we use csparse from Timothy Davis 92 93 Timothy Davis, 94 Direct Methods for Sparse Linear Systems, 95 SIAM, 2006, 96 ISBN: 0898716136, 97 LC: QA188.D386. 98 99 matrix in compressed row/column or triplet form : 100 { 101 csi nzmax ; : maximum number of entries 102 csi m ; : number of rows 103 csi n ; : number of columns 104 csi *p ; : compressed: row (size m+1) or column (size n+1) pointers; triplet: row indices (size nz) 105 csi *i ; : compressed: column or row indices, size nzmax; triplet: column indices (size nz) 106 double *x ; : numerical values, size nzmax 107 csi nz ; : # of entries in triplet matrix; 108 -1 for compressed columns; 109 -2 for compressed rows 110 111 } 112 113 csi is either int64_t or int32_t and this is controlled at compile time*/ 114 115 116 #define NSM_UNKNOWN_ERR(func, orig) \ 117 fprintf(stderr, #func ": unknown origin %d for sparse matrix\n", orig); 118 119 120 #define NSM_NROW_CSR(mat) mat->n 121 #define NSM_NCOL_CSR(mat) mat->m 122 123 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 124 extern "C" 125 { 126 #endif 127 128 /** \struct CSparseMatrix_factors 129 * Information used and produced by CSparse for an LU factorization*/ 130 typedef struct { 131 CS_INT n; /**< size of linear system */ 132 css* S; /**< symbolic analysis */ 133 csn* N; /**< numerics factorization */ 134 } CSparseMatrix_factors; 135 136 /** compute a LU factorization of A and store it in a workspace 137 * \param order control if ordering is used 138 * \param A the sparse matrix 139 * \param tol the tolerance 140 * \param cs_lu_A the parameter structure that eventually holds the factors 141 * \return 1 if the factorization was successful, 1 otherwise 142 */ 143 int CSparseMatrix_lu_factorization(CS_INT order, const CSparseMatrix *A, double tol, CSparseMatrix_factors * cs_lu_A); 144 145 /** compute a Cholesky factorization of A and store it in a workspace 146 * \param order control if ordering is used 147 * \param A the sparse matrix 148 * \param cs_chol_A the parameter structure that eventually holds the factors 149 * \return 1 if the factorization was successful, 1 otherwise 150 */ 151 int CSparseMatrix_chol_factorization(CS_INT order, const CSparseMatrix *A, CSparseMatrix_factors * cs_chol_A); 152 153 /** compute a LDLT factorization of A and store it in a workspace 154 * \param order control if ordering is used 155 * \param A the sparse matrix 156 * \param cs_ldlt_A the parameter structure that eventually holds the factors 157 * \return 1 if the factorization was successful, 1 otherwise 158 */ 159 int CSparseMatrix_ldlt_factorization(CS_INT order, const CSparseMatrix *A, CSparseMatrix_factors * cs_ldlt_A); 160 161 /** reuse a LU factorization (stored in the cs_lu_A) to solve a linear system Ax = b 162 * \param cs_lu_A contains the LU factors of A, permutation information 163 * \param x workspace 164 * \param[in,out] b on input RHS of the linear system; on output the solution 165 * \return 0 if failed, 1 otherwise*/ 166 CS_INT CSparseMatrix_solve(CSparseMatrix_factors* cs_lu_A, double* x, double *b); 167 168 /** reuse a LU factorization (stored in the cs_lu_A) to solve a linear system Ax = B 169 * with a sparse r.h.s 170 * \param cs_lu_A contains the LU factors of A, permutation information 171 * \param X a csc sparse matrix workspace 172 * \param[in,out] B on input sparse RHS of the linear system; on output the solution 173 * \return 0 if failed, 1 otherwise*/ 174 CS_INT CSparseMatrix_spsolve(CSparseMatrix_factors* cs_lu_A, CSparseMatrix* X, CSparseMatrix* B); 175 176 /** reuse a Cholesky factorization (stored in the cs_chol_A) to solve a linear system Ax = b 177 * \param cs_chol_A contains the Cholesky factors of A, permutation information 178 * \param x workspace 179 * \param[in,out] b on input RHS of the linear system; on output the solution 180 * \return 0 if failed, 1 otherwise*/ 181 CS_INT CSparseMatrix_chol_solve(CSparseMatrix_factors* cs_chol_A, double* x, double *b); 182 183 /** reuse a Cholesky factorization (stored in the cs_chol_A) to solve a linear system Ax = B 184 * with a sparse r.h.s 185 * \param cs_chol_A contains the Cholesky factors of A, permutation information 186 * \param X a csc sparse matrix workspace 187 * \param[in,out] b on input sparse RHS of the linear system; on output the solution 188 * \return 0 if failed, 1 otherwise*/ 189 CS_INT CSparseMatrix_chol_spsolve(CSparseMatrix_factors* cs_chol_A, CSparseMatrix* X, CSparseMatrix* B); 190 191 /** reuse a LDLT factorization (stored in the cs_ldlt_A) to solve a linear system Ax = b 192 * \param cs_ldlt_A contains the LDLT factors of A, permutation information 193 * \param x workspace 194 * \param[in,out] b on input RHS of the linear system; on output the solution 195 * \return 0 if failed, 1 otherwise*/ 196 CS_INT CSparseMatrix_ldlt_solve(CSparseMatrix_factors* cs_ldlt_A, double* x, double *b); 197 /** Free a workspace related to a LU factorization 198 * \param cs_lu_A the structure to free 199 */ 200 void CSparseMatrix_free_lu_factors(CSparseMatrix_factors* cs_lu_A); 201 202 /** Matrix vector multiplication : y = alpha*A*x+beta*y 203 * \param[in] alpha matrix coefficient 204 * \param[in] A the sparse matrix 205 * \param[in] x pointer on a dense vector of size A->n 206 * \param[in] beta vector coefficient 207 * \param[in, out] y pointer on a dense vector of size A->n 208 * \return 0 if A x or y is NULL else 1 209 */ 210 int CSparseMatrix_aaxpby(const double alpha, const CSparseMatrix *A, const double *x, 211 const double beta, double *y); 212 213 /** Allocate a CSparse matrix for future copy (as in NSM_copy) 214 * \param m the matrix used as model 215 * \return an newly allocated matrix 216 */ 217 CSparseMatrix* CSparseMatrix_alloc_for_copy(const CSparseMatrix* const m); 218 219 CS_INT CSparseMatrix_to_dense(const CSparseMatrix* const A, double * B); 220 221 /** print a matrix to std output 222 * \param A matrix to print 223 * \param brief if positive, print only a portion of the matrix 224 */ 225 int CSparseMatrix_print(const CSparseMatrix *A, int brief); 226 227 /** print a matrix to a text file 228 * \param A matrix to print 229 * \param brief if positive, print only a portion of the matrix 230 * \param file file descriptor*/ 231 int CSparseMatrix_print_in_file(const CSparseMatrix *A, int brief, FILE* file); 232 233 CSparseMatrix * CSparseMatrix_new_from_file(FILE* file); 234 235 /** Add an entry to a triplet matrix only if the absolute value is 236 * greater than threshold 237 * \param T the sparse matrix 238 * \param i row index 239 * \param j column index 240 * \param x the value 241 * \return integer value : 1 if the absolute value is less than 242 * DBL_EPSILON, otherwise the return value of cs_entry. 243 */ 244 CS_INT CSparseMatrix_zentry(CSparseMatrix *T, CS_INT i, CS_INT j, double x, double threshold); 245 246 /** Add an entry to a symmetric triplet matrix only if the absolute value is 247 * greater than threshold 248 * \param T the sparse matrix 249 * \param i row index 250 * \param j column index 251 * \param x the value 252 * \return integer value : 1 if the absolute value is less than 253 * DBL_EPSILON, otherwise the return value of cs_entry. 254 */ 255 CS_INT CSparseMatrix_symmetric_zentry(CSparseMatrix *T, CS_INT i, CS_INT j, double x, double threshold); 256 257 /** Add an entry to a triplet matrix 258 * \param T the sparse matrix 259 * \param i row index 260 * \param j column index 261 * \param x the value 262 */ 263 CS_INT CSparseMatrix_entry(CSparseMatrix *T, CS_INT i, CS_INT j, double x); 264 265 /** Add an entry to a symmetric triplet matrix 266 * \param T the sparse matrix 267 * \param i row index 268 * \param j column index 269 * \param x the value 270 */ 271 CS_INT CSparseMatrix_symmetric_entry(CSparseMatrix *T, CS_INT i, CS_INT j, double x); 272 273 /** Check if the given triplet matrix is properly constructed (col and row indices are correct) 274 * \param T the sparse matrix to check 275 * \return 0 if the matrix is fine, 1 otherwise 276 * */ 277 int CSparseMatrix_check_triplet(CSparseMatrix *T); 278 279 /** Check if the given triplet matrix is properly constructed (col and row indices are correct) 280 * \param T the sparse matrix to check 281 * \return 0 if the matrix is fine, 1 otherwise 282 * */ 283 int CSparseMatrix_check_csc(CSparseMatrix *T); 284 285 /** Free space allocated for a SparseMatrix. note : cs_spfree also 286 * free the cs_struct this fails when the struct is allocated on 287 * the stack. 288 * \param A the sparse matrix 289 * \return NULL on success 290 */ 291 CSparseMatrix* CSparseMatrix_spfree_on_stack(CSparseMatrix* A); 292 293 /** Copy a CSparseMatrix inside another CSparseMatrix. 294 * Reallocations are performed if B cannot hold a copy of A 295 * \param[in] A a CSparseMatrix 296 * \param[in,out] B a CSparseMatrix 297 */ 298 void CSparseMatrix_copy(const CSparseMatrix* const A, CSparseMatrix* B); 299 300 /** Multiply a matrix with a double alpha*A --> A 301 * \param alpha the coefficient 302 * \param A the matrix 303 */ 304 int CSparseMatrix_scal(const double alpha, const CSparseMatrix *A); 305 306 307 /** Return the element A(i,j) 308 * \param A the sparse matrix 309 * \param i the row index 310 * \param j the column index 311 */ 312 double CSparseMatrix_get_value(const CSparseMatrix *A, CS_INT i, CS_INT j); 313 314 /** print a matrix to a text file in pyhton format 315 * \param m matrix to print 316 * \param file file descriptor*/ 317 void CSparseMatrix_write_in_file_python(const CSparseMatrix* const m, FILE* file); 318 319 /** Compute the max by columns of a sparse matrix 320 * \param A the sparse matrix 321 * \param max, the vector of maximum by columns 322 * \param j the column index 323 */ 324 int CSparseMatrix_max_by_columns(const CSparseMatrix *A, double * max); 325 326 /** Compute the max in absolute value by columns of a sparse matrix 327 * \param A the sparse matrix 328 * \param max, the vector of maximum by columns 329 * \param j the column index 330 */ 331 332 int CSparseMatrix_max_abs_by_columns(const CSparseMatrix *A, double * max); 333 334 #if defined(__cplusplus) && !defined(BUILD_AS_CPP) 335 } 336 #endif 337 338 #endif 339