Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 128
Functions: 0.0% 0 / 0 / 4
Branches: 0.0% 0 / 0 / 76
OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverKlu.c
Line Branch Exec Source
1 /*
2 * This file belongs to the OpenModelica Run-Time System
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings
5 * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights
6 * reserved.
7 *
8 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
9 * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY
10 * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
11 * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL
12 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
13 *
14 * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License
15 * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs:
16 * http://www.openmodelica.org or https://github.com/OpenModelica/ or
17 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU
18 * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW
19 * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY
23 * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF
24 * OSMC-PL.
25 *
26 */
27
28 /*! \file linearSolverKlu.c
29 */
30
31 #include "omc_config.h"
32
33 #ifdef WITH_SUITESPARSE
34 #include <math.h>
35 #include <stdlib.h>
36 #include <string.h>
37
38 #include "simulation_data.h"
39 #include "simulation/simulation_info_json.h"
40 #include "util/omc_error.h"
41 #include "omc_math.h"
42 #include "util/varinfo.h"
43 #include "model_help.h"
44
45 #include "linearSystem.h"
46 #include "linearSolverKlu.h"
47
48 static void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n);
49 static void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n);
50
51 /*! \fn allocate memory for linear system solver Klu
52 *
53 */
54 ✗ int allocateKluData(int n_row, int n_col, int nz, void** voiddata)
55 {
56 ✗ DATA_KLU* data = (DATA_KLU*) malloc(sizeof(DATA_KLU));
57 ✗ assertStreamPrint(NULL, 0 != data, "Could not allocate data for linear solver Klu.");
58
59 ✗ data->symbolic = NULL;
60 ✗ data->numeric = NULL;
61
62 ✗ data->n_col = n_col;
63 ✗ data->n_row = n_row;
64 ✗ data->nnz = nz;
65
66 ✗ data->Ap = (int*) calloc((n_row+1),sizeof(int));
67 ✗ data->Ai = (int*) calloc(nz,sizeof(int));
68 ✗ data->Ax = (double*) calloc(nz,sizeof(double));
69 ✗ data->work = (double*) calloc(n_col,sizeof(double));
70
71 ✗ data->numberSolving = 0;
72 ✗ klu_defaults(&(data->common));
73
74 ✗ *voiddata = (void*)data;
75
76 ✗ return 0;
77 }
78
79
80 /*! \fn free memory for linear system solver Klu
81 *
82 */
83 ✗ int freeKluData(void **voiddata)
84 {
85 ✗ DATA_KLU* data = (DATA_KLU*) *voiddata;
86
87 ✗ free(data->Ap);
88 ✗ free(data->Ai);
89 ✗ free(data->Ax);
90 ✗ free(data->work);
91
92
93 ✗ if(data->symbolic)
94 ✗ klu_free_symbolic(&data->symbolic, &data->common);
95 ✗ if(data->numeric)
96 ✗ klu_free_numeric(&data->numeric, &data->common);
97
98 ✗ return 0;
99 }
100
101 /*! \fn getAnalyticalJacobian
102 *
103 * function calculates analytical jacobian
104 *
105 * \param [ref] [data]
106 * \param [in] [sysNumber]
107 *
108 * \author wbraun
109 *
110 */
111 ✗ static void getAnalyticalJacobian(DATA* data, threadData_t *threadData,
112 LINEAR_SYSTEM_DATA* systemData)
113 {
114 int i,j,l,nth;
115 ✗ JACOBIAN* jacobian = systemData->jacobian;
116 ✗ JACOBIAN* parentJacobian = systemData->parentJacobian;
117 ✗ const SPARSE_PATTERN* sp = jacobian->sparsePattern;
118
119 /* evaluate constant equations of Jacobian */
120 ✗ if (jacobian->constantEqns != NULL) {
121 ✗ jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
122 }
123
124 /* evaluate Jacobian */
125 ✗ for (i = 0; i < sp->maxColors; i++) {
126 /* activate seed variable for the corresponding color */
127 ✗ for (j = 0; j < jacobian->sizeCols; j++)
128 ✗ if (sp->colorCols[j]-1 == i)
129 ✗ jacobian->seedVars[j] = 1.0;
130
131 /* Evaluate Jacobian column */
132 ✗ jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
133
134 ✗ for (j = 0; j < jacobian->sizeCols; j++) {
135 ✗ if (sp->colorCols[j]-1 == i) {
136 ✗ for (nth = sp->leadindex[j]; nth < sp->leadindex[j+1]; nth++) {
137 ✗ l = sp->index[nth];
138 ✗ systemData->setAElement(j, l, -jacobian->resultVars[l], nth, systemData, threadData);
139 }
140 /* de-activate seed variable for the corresponding color */
141 ✗ jacobian->seedVars[j] = 0.0;
142 }
143 }
144 }
145 ✗ }
146
147 /*! \fn residual_wrapper for the residual function
148 *
149 */
150 static int residual_wrapper(double* x, double* f, RESIDUAL_USERDATA* userData, int sysNumber)
151 {
152 ✗ int iflag = 0;
153 ✗ userData->data->simulationInfo->linearSystemData[sysNumber].residualFunc(userData, x, f, &iflag);
154 ✗ return 0;
155 }
156
157 /*! \fn solve linear system with Klu method
158 *
159 * \param [in] [data]
160 * [sysNumber] index of the corresponding linear system
161 *
162 *
163 * author: wbraun
164 */
165 ✗ int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
166 {
167 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
168 ✗ LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
169 ✗ DATA_KLU* solverData = (DATA_KLU*)systemData->solverData[0];
170 _omc_scalar residualNorm = 0;
171
172 ✗ int i, j, status = 0, success = 0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
173 double tmpJacEvalTime;
174 ✗ int reuseMatrixJac = (data->simulationInfo->currentContext == CONTEXT_SYM_JACOBIAN && data->simulationInfo->currentJacobianEval > 0);
175
176 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_LS, omc_dummyFileInfo, 0, indexes,
177 "Start solving Linear System %d (size %d) at time %g with Klu Solver",
178 ✗ eqSystemNumber, (int) systemData->size, data->localData[0]->timeValue);
179
180 ✗ rt_ext_tp_tick(&(solverData->timeClock));
181 ✗ if (0 == systemData->method)
182 {
183 ✗ if (!reuseMatrixJac){
184 /* set A matrix */
185 ✗ solverData->Ap[0] = 0;
186 ✗ systemData->setA(data, threadData, systemData);
187 ✗ solverData->Ap[solverData->n_row] = solverData->nnz;
188 }
189
190 /* set b vector */
191 ✗ systemData->setb(data, threadData, systemData);
192 } else {
193
194 ✗ if (!reuseMatrixJac){
195 ✗ solverData->Ap[0] = 0;
196 /* calculate jacobian -> matrix A*/
197 ✗ if(systemData->jacobianIndex != -1){
198 ✗ getAnalyticalJacobian(data, threadData, systemData);
199 } else {
200 assertStreamPrint(threadData, 1, "jacobian function pointer is invalid" );
201 }
202 ✗ solverData->Ap[solverData->n_row] = solverData->nnz;
203 }
204
205 /* calculate vector b (rhs) */
206 ✗ memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row);
207
208 ✗ residual_wrapper(solverData->work, systemData->b, &resUserData, sysNumber);
209 }
210 ✗ tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
211 ✗ systemData->jacobianTime += tmpJacEvalTime;
212 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
213
214 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_LS_V))
215 {
216 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Old solution x:");
217 ✗ for(i = 0; i < solverData->n_row; ++i)
218 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
219 ✗ messageClose(OMC_LOG_LS_V);
220
221 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Matrix A n_rows = %d", solverData->n_row);
222 ✗ for (i=0; i<solverData->n_row; i++){
223 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "%d. Ap => %d -> %d", i, solverData->Ap[i], solverData->Ap[i+1]);
224 ✗ for (j=solverData->Ap[i]; j<solverData->Ap[i+1]; j++){
225 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "A[%d,%d] = %f", i, solverData->Ai[j], solverData->Ax[j]);
226 }
227 }
228 ✗ messageClose(OMC_LOG_LS_V);
229
230 ✗ for (i=0; i<solverData->n_row; i++)
231 {
232 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "b[%d] = %e", i, systemData->b[i]);
233 }
234 }
235 ✗ rt_ext_tp_tick(&(solverData->timeClock));
236
237 /* symbolic pre-ordering of A to reduce fill-in of L and U */
238 ✗ if (0 == solverData->numberSolving)
239 {
240 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Perform analyze settings:\n - ordering used: %d\n - current status: %d", solverData->common.ordering, solverData->common.status);
241 ✗ solverData->symbolic = klu_analyze(solverData->n_col, solverData->Ap, solverData->Ai, &solverData->common);
242 }
243
244 /* if reuseMatrixJac use also previous factorization */
245 ✗ if (!reuseMatrixJac)
246 {
247 /* compute the LU factorization of A */
248 ✗ if (0 == solverData->common.status){
249 ✗ if(solverData->numeric){
250 /* Just refactor using the same pivots, but check that the refactor is still accurate */
251 ✗ klu_refactor(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, solverData->numeric, &solverData->common);
252 ✗ klu_rgrowth(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, solverData->numeric, &solverData->common);
253 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Klu rgrowth after refactor: %f", solverData->common.rgrowth);
254 /* If rgrowth is small then do a whole factorization with new pivots (What should this tolerance be?) */
255 ✗ if (solverData->common.rgrowth < 1e-3){
256 ✗ klu_free_numeric(&solverData->numeric, &solverData->common);
257 ✗ solverData->numeric = klu_factor(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, &solverData->common);
258 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Klu new factorization performed.");
259 }
260 } else {
261 ✗ solverData->numeric = klu_factor(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, &solverData->common);
262 }
263 }
264 }
265
266 ✗ if (0 == solverData->common.status){
267 ✗ if (1 == systemData->method){
268 ✗ if (klu_solve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->b, &solverData->common)){
269 success = 1;
270 }
271 } else {
272 ✗ if (klu_tsolve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->b, &solverData->common)){
273 success = 1;
274 }
275 }
276 }
277
278 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
279
280 /* print solution */
281 ✗ if (1 == success){
282
283 ✗ if (1 == systemData->method){
284 /* take the solution */
285 ✗ for(i = 0; i < solverData->n_row; ++i)
286 ✗ aux_x[i] += systemData->b[i];
287
288 /* update inner equations */
289 ✗ residual_wrapper(aux_x, solverData->work, &resUserData, sysNumber);
290 ✗ residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
291
292 ✗ if ((isnan(residualNorm)) || (residualNorm>1e-4)) {
293 ✗ warningStreamPrintWithLimit(OMC_LOG_LS, 0, ++(systemData->numberOfFailures) /* Update counter */, data->simulationInfo->maxWarnDisplays,
294 "Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
295 ✗ (int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
296 success = 0;
297 }
298 } else {
299 /* the solution is automatically in x */
300 ✗ memcpy(aux_x, systemData->b, sizeof(double)*systemData->size);
301 }
302
303 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_LS_V))
304 {
305 ✗ if (1 == systemData->method) {
306 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
307 } else {
308 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Solution x:");
309 }
310 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
311
312 ✗ for(i = 0; i < systemData->size; ++i)
313 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
314
315 ✗ messageClose(OMC_LOG_LS_V);
316 }
317 }
318 else
319 {
320 ✗ warningStreamPrintWithLimit(OMC_LOG_STDOUT, 0, ++(systemData->numberOfFailures) /* Update counter */, data->simulationInfo->maxWarnDisplays,
321 "Failed to solve linear system of equations (no. %d) at time %f, system status %d.",
322 ✗ (int)systemData->equationIndex, data->localData[0]->timeValue, status);
323 }
324 ✗ solverData->numberSolving += 1;
325
326 ✗ return success;
327 }
328
329 static
330 void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n)
331 {
332 int i, j, k, l;
333
334 char **buffer = (char**)malloc(sizeof(char*)*n);
335 for (l=0; l<n; l++)
336 {
337 buffer[l] = (char*)malloc(sizeof(char)*n*20);
338 buffer[l][0] = 0;
339 }
340
341 char **p = (char**)malloc(sizeof(char*)*n);
342 for (l=0; l<n; l++)
343 p[l] = buffer[l];
344
345 k = 0;
346 for (i = 0; i < n; i++)
347 {
348 for (j = 0; j < n; j++)
349 {
350 if ((k < Ap[i + 1]) && (Ai[k] == j))
351 {
352 p[j] += sprintf(p[j], " %5g ", Ax[k]);
353 k++;
354 }
355 else
356 {
357 p[j] += sprintf(p[j], " %5g ", 0.0);
358 }
359 }
360 }
361 for (l = 0; l < n; l++)
362 {
363 infoStreamPrint(OMC_LOG_LS_V, 0, "%s", buffer[l]);
364 free(buffer[l]);
365 }
366 free(p);
367 free(buffer);
368 }
369
370 static
371 void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n)
372 {
373 int i, j, k;
374 char *buffer = (char*)malloc(sizeof(char)*n*15);
375 char *q;
376 k = 0;
377 for (i = 0; i < n; i++)
378 {
379 q = buffer;
380 for (j = 0; j < n; j++)
381 {
382 if ((k < Ap[i + 1]) && (Ai[k] == j))
383 {
384 q += sprintf(q, " %5.2g ", Ax[k]);
385 k++;
386 }
387 else
388 {
389 q += sprintf(q, " %5.2g ", 0.0);
390 }
391 }
392 infoStreamPrint(OMC_LOG_LS_V, 0, "%s", buffer);
393 }
394 free(buffer);
395 }
396
397 #endif
398