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 #include "SiconosConfig.h" // for WITH_SUPERLU, SUPERLU_MAJOR_VERSION ... // IWYU pragma: keep
20 #include "siconos_debug.h"
21 #include "numerics_verbose.h"
22 
23 #ifdef WITH_SUPERLU
24 
25 #include <slu_ddefs.h>
26 #include "CSparseMatrix_internal.h"
27 #include "NumericsMatrix_internal.h"
28 #include "NumericsSparseMatrix.h"
29 
30 
31 /** \struct NM_SuperLU_WS NumericsMatrix_internal.h
32  * Structure for holding the data SuperLU needs
33  */
34 struct NM_SuperLU_WS
35 {
36   SuperMatrix* L;
37   SuperMatrix* U;
38   int_t* perm_r;
39   int_t* perm_c;
40   superlu_options_t* options;
41 #ifdef SUPERLU_MAJOR_VERSION
42   GlobalLU_t* Glu;
43 #endif
44 };
45 
NM_SuperLU_factorize(NumericsMatrix * A)46 NM_SuperLU_WS* NM_SuperLU_factorize(NumericsMatrix* A)
47 {
48   SuperMatrix SA, SAC;
49   SuperLUStat_t stat;
50   int_t *etree;
51 
52   int status;
53 
54   NSM_linear_solver_params* params = NSM_linearSolverParams(A);
55 
56   if(params->linear_solver_data)
57   {
58     return (NM_SuperLU_WS*) params->linear_solver_data;
59   }
60 
61   params->linear_solver_data = calloc(1, sizeof(NM_SuperLU_WS));
62   NM_SuperLU_WS* superlu_ws = (NM_SuperLU_WS*) params->linear_solver_data;
63 
64   if(!superlu_ws->options) superlu_ws->options = (superlu_options_t*)malloc(sizeof(superlu_options_t));
65 
66 #ifdef SUPERLU_MAJOR_VERSION
67   if(!superlu_ws->Glu) superlu_ws->Glu = (GlobalLU_t*)malloc(sizeof(GlobalLU_t));
68 #endif
69 
70   set_default_options(superlu_ws->options);
71 
72   if(verbose > 1)
73     superlu_ws->options->PrintStat = YES;
74   /* TODO SuperLU_PIVOT_TOLERANCE, SuperLU_ORDERING, SuperLU_SCALE
75    * SuperLU_DROPTOL, SuperLU_STRATEGY, SuperLU_IRSTEP*/
76 
77   StatInit(&stat);
78 
79   CSparseMatrix* C = NM_csc(A);
80 
81   superlu_ws->L = (SuperMatrix *) SUPERLU_MALLOC(sizeof(SuperMatrix));
82   superlu_ws->U = (SuperMatrix *) SUPERLU_MALLOC(sizeof(SuperMatrix));
83   if(!(superlu_ws->perm_r = intMalloc(C->m))) ABORT("Malloc fails for perm_r[].");
84   if(!(superlu_ws->perm_c = intMalloc(C->n))) ABORT("Malloc fails for perm_c[].");
85   if(!(etree = intMalloc(C->n))) ABORT("Malloc fails for etree[].");
86 
87   /* Symbolic part  */
88   int_t* indices;
89   int_t* pointers;
90   size_t nnz = NSM_nnz(C);
91 
92   if(sizeof(*C->i) != sizeof(*indices))
93   {
94     int_t* iWork = (int_t*)NM_iWork(A, (size_t)(nnz + C->n) + 1, sizeof(int_t));
95     indices = iWork;
96     pointers = &iWork[nnz];
97 
98     for(size_t i = 0; i < nnz  ; ++i) indices[i] = (int_t)C->i[i];
99     for(size_t i = 0; i <= C->n; ++i) pointers[i] = (int_t)C->p[i];
100   }
101   else
102   {
103     indices = (int_t*)C->i;
104     pointers = (int_t*)C->p;
105   }
106 
107   dCreate_CompCol_Matrix(&SA, C->m, C->n, nnz, C->x, indices, pointers, SLU_NC, SLU_D, SLU_GE);
108 
109   int permc_spec = 3;
110   get_perm_c(permc_spec, &SA, superlu_ws->perm_c);
111 
112   sp_preorder(superlu_ws->options, &SA, superlu_ws->perm_c, etree, &SAC);
113 
114   /* Numerical part */
115   int panel_size = sp_ienv(1);
116   int relax = sp_ienv(2);
117   double drop_tol = 0.0;
118 
119   /* Versions 4.x and earlier do not include a #define'd version numbers  */
120 #ifndef SUPERLU_MAJOR_VERSION
121   dgstrf(superlu_ws->options, &SAC, drop_tol, relax, &panel_size, etree, NULL, 0, superlu_ws->perm_c, superlu_ws->perm_r, superlu_ws->L, superlu_ws->U, &stat, &status);
122 #else
123   dgstrf(superlu_ws->options, &SAC, relax, panel_size, etree, NULL, 0, superlu_ws->perm_c, superlu_ws->perm_r, superlu_ws->L, superlu_ws->U, superlu_ws->Glu, &stat, &status);
124 #endif
125 
126   if(status)
127   {
128     fprintf(stderr, "dgstrf() error returns INFO= %d\n", status);
129     return NULL;
130   }
131 
132   SUPERLU_FREE(etree);
133   Destroy_SuperMatrix_Store(&SA);
134   Destroy_CompCol_Permuted(&SAC);
135   StatFree(&stat);
136 
137   return superlu_ws;
138 }
139 
NM_SuperLU_solve(NumericsMatrix * A,double * b,NM_SuperLU_WS * superlu_ws)140 int NM_SuperLU_solve(NumericsMatrix* A, double* b, NM_SuperLU_WS* superlu_ws)
141 {
142   SuperMatrix B;
143   SuperLUStat_t stat;
144 
145   int status;
146 
147   CSparseMatrix* C = NM_csc(A);
148 
149   dCreate_Dense_Matrix(&B, C->n, 1, b, C->n, SLU_DN, SLU_D, SLU_GE);
150 
151   StatInit(&stat);
152 
153   dgstrs(NOTRANS, superlu_ws->L, superlu_ws->U, superlu_ws->perm_c, superlu_ws->perm_r, &B, &stat, &status);
154 
155   Destroy_SuperMatrix_Store(&B);
156   StatFree(&stat);
157 
158   return status;
159 }
160 
NM_SuperLU_free(void * p)161 void NM_SuperLU_free(void* p)
162 {
163   assert(p);
164   NSM_linear_solver_params* params = (NSM_linear_solver_params*) p;
165   assert(params);
166   NM_SuperLU_WS* superlu_ws = (NM_SuperLU_WS*) params->linear_solver_data;
167   assert(superlu_ws);
168 
169   SUPERLU_FREE(superlu_ws->perm_r);
170   SUPERLU_FREE(superlu_ws->perm_c);
171   Destroy_SuperNode_Matrix(superlu_ws->L);
172   Destroy_CompCol_Matrix(superlu_ws->U);
173   SUPERLU_FREE(superlu_ws->L);
174   SUPERLU_FREE(superlu_ws->U);
175 
176   superlu_ws->perm_r = NULL;
177   superlu_ws->perm_c = NULL;
178   superlu_ws->L = NULL;
179   superlu_ws->U = NULL;
180 
181   if(superlu_ws->options)
182   {
183     free(superlu_ws->options);
184     superlu_ws->options = NULL;
185   }
186 
187 #ifdef SUPERLU_MAJOR_VERSION
188   if(superlu_ws->Glu)
189   {
190     free(superlu_ws->Glu);
191     superlu_ws->Glu = NULL;
192   }
193 #endif
194   /* Here we free superlu_ws ...  */
195   free(superlu_ws);
196   params->linear_solver_data = NULL;
197 
198 }
199 
NM_SuperLU_extra_display(NM_SuperLU_WS * superlu_ws)200 void NM_SuperLU_extra_display(NM_SuperLU_WS* superlu_ws)
201 {
202 #if 0
203   if(verbose > 2)
204   {
205     SuperLU_FN(report_info)(superlu_ws->control, superlu_ws->info);
206 
207     if(verbose > 3)
208     {
209       SuperLU_FN(report_control)(superlu_ws->control);
210     }
211   }
212   else if(verbose > 1)
213   {
214     if(superlu_ws->control[SuperLU_IRSTEP] > 0)
215     {
216       printf("SuperLU : backward error estimate omega1 %g\n", superlu_ws->info[SuperLU_OMEGA1]);
217       printf("SuperLU : backward error estimate omega2 %g\n", superlu_ws->info[SuperLU_OMEGA2]);
218     }
219     printf("SuperLU : solve FLOPS %g\n", superlu_ws->info[SuperLU_SOLVE_FLOPS]);
220     printf("SuperLU : solve time %g\n", superlu_ws->info[SuperLU_SOLVE_TIME]);
221     printf("SuperLU : wall time %g\n", superlu_ws->info[SuperLU_SOLVE_WALLTIME]);
222 
223   }
224 #endif
225 }
226 
227 #endif /* WITH_SuperLU */
228