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 /*!\file lumod_wrapper.h
20  * \brief Wrapper around lumod */
21 
22 #ifndef LUMOD_WRAPPER_H
23 #define LUMOD_WRAPPER_H
24 
25 #include "NumericsFwd.h"    // for NumericsMatrix
26 #include "SiconosLapack.h"  // for lapack_int
27 #include "assert.h"         // for assert
28 
29 #define SN_LUMOD_NEED_REFACTORIZATION 1
30 
31 /**\struct SN_lumod_dense_data lumod_wrapper.h
32  * Data structure for the LUMOD (successive rank-one update of a matrix) */
33 typedef struct {
34   unsigned n; /**< size of the matrix H*/
35   unsigned maxmod; /**< maximum number of changes */
36   unsigned k; /**< number of rows (or columns) of C */
37   double* LU_H; /**< LU factors of the initial matrix H */
38   lapack_int* ipiv_LU_H; /**< pivot for the LU factorization of H*/
39   unsigned* factorized_basis; /**< basis when H was factorized and storing for the info where the columns of the non basic variables are in U and when a basic variable exited */
40   int* row_col_indx; /**< Store the information to which column or row the variable correspond */
41   double* Uk; /**< matrix which keeps track of modified columns in Hk*/
42   double* Yk; /**< matrix of changed columns in Hk (row major)*/
43   double* L_C; /**< L matrix of the LU factorization for Ck*/
44   double* U_C; /**< U matrix of the LU factorization  for Ck*/
45   double* y; /**< y vector for LUMOD */
46   double* z; /**< z vector for LUMOD */
47   double* w; /**< work vector for LUMOD */
48 } SN_lumod_dense_data;
49 
50 
51 
SN_lumod_need_refactorization(int info)52 static inline unsigned SN_lumod_need_refactorization(int info) { return (info == SN_LUMOD_NEED_REFACTORIZATION ? 1 : 0); }
53 
SN_lumod_find_arg_var(int * array,int p,unsigned n)54 static inline unsigned SN_lumod_find_arg_var(int* array, int p, unsigned n)
55 {
56   for (unsigned i = 0; i < 2*n + 1; ++i)
57   {
58     if (array[i] == p) return i;
59   }
60   assert(0 && "We should not be here");
61   return 0;
62 }
63 
64 SN_lumod_dense_data* SN_lumod_dense_allocate(unsigned n, unsigned maxmod);
65 void SM_lumod_dense_free(SN_lumod_dense_data* lumod_data);
66 int SN_lumod_dense_solve(SN_lumod_dense_data* lumod_data, double* x, double* col_tilde);
67 void SN_lumod_add_row_col(SN_lumod_dense_data* lumod_data, unsigned leaving_indx, double* col);
68 void SN_lumod_replace_col(SN_lumod_dense_data* lumod_data, unsigned index_col, double* col);
69 void SN_lumod_replace_row(SN_lumod_dense_data* lumod_data, unsigned index_row, unsigned leaving_indx);
70 void SN_lumod_delete_row_col(SN_lumod_dense_data* lumod_data, unsigned index_row, unsigned index_col);
71 int SN_lumod_factorize(SN_lumod_dense_data* lumod_data, unsigned* basis, NumericsMatrix* M, double* covering_vector);
72 
73 #endif
74