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