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