Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 7.4% 23 / 0 / 309
Functions: 30.8% 4 / 0 / 13
Branches: 3.8% 5 / 0 / 133
OMCompiler/SimulationRuntime/c/simulation/solver/linearSystem.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 linearSystem.c
29 */
30
31 #include <math.h>
32 #include <string.h>
33
34 #include "model_help.h"
35 #include "../../util/omc_error.h"
36 #include "../jacobian_util.h"
37 #include "../../util/rtclock.h"
38 #include "nonlinearSystem.h"
39 #include "linearSystem.h"
40 #include "linearSolverLapack.h"
41 #if !defined(OMC_MINIMAL_RUNTIME)
42 #include "linearSolverKlu.h"
43 #include "linearSolverLis.h"
44 #include "linearSolverUmfpack.h"
45 #endif
46 #include "linearSolverTotalPivot.h"
47 #include "../options.h"
48 #include "../simulation_info_json.h"
49
50 static void setAElement(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t* threadData);
51 static void setAElementLis(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t* threadData);
52 static void setAElementUmfpack(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t* threadData);
53 static void setAElementKlu(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t* threadData);
54 static void setBElement(int row, double value, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t* threadData);
55 static void setBElementLis(int row, double value, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t* threadData);
56
57 int check_linear_solution(DATA *data, int printFailingSystems, int sysNumber);
58
59 /*! \fn int initializeLinearSystems(DATA *data)
60 *
61 * This function allocates memory for all linear systems.
62 *
63 * \param [ref] [data]
64 */
65 1 int initializeLinearSystems(DATA *data, threadData_t *threadData)
66 {
67 int i, nnz;
68 int size;
69 1 LINEAR_SYSTEM_DATA *linsys = data->simulationInfo->linearSystemData;
70
71 1 infoStreamPrint(OMC_LOG_LS, 1, "initialize linear system solvers");
72 1 infoStreamPrint(OMC_LOG_LS, 0, "%ld linear systems", data->modelData->nLinearSystems);
73
74
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (LSS_DEFAULT == data->simulationInfo->lssMethod) {
75 #ifdef WITH_SUITESPARSE
76 1 data->simulationInfo->lssMethod = LSS_KLU;
77 #elif !defined(OMC_MINIMAL_RUNTIME)
78 data->simulationInfo->lssMethod = LSS_LIS;
79 #endif
80 }
81
82
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(i=0; i<data->modelData->nLinearSystems; ++i)
83 {
84 ✗ size = linsys[i].size;
85 ✗ nnz = linsys[i].nnz;
86 ✗ linsys[i].totalTime = 0;
87 ✗ linsys[i].failed = 0;
88
89 /* allocate system data */
90 ✗ linsys[i].b = (double*) malloc(size*sizeof(double));
91
92 /* check if analytical jacobian is created */
93 ✗ if (1 == linsys[i].method)
94 {
95 ✗ if(linsys[i].jacobianIndex != -1)
96 {
97 ✗ assertStreamPrint(threadData, 0 != linsys[i].analyticalJacobianColumn, "jacobian function pointer is invalid" );
98 }
99 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[linsys[i].jacobianIndex]);
100 ✗ if(linsys[i].initialAnalyticalJacobian(data, threadData, jacobian))
101 {
102 ✗ linsys[i].jacobianIndex = -1;
103 ✗ throwStreamPrint(threadData, "Failed to initialize the jacobian for torn linear system %d.", (int)linsys[i].equationIndex);
104 }
105 /* evalJacobian() fills sizeRows*sizeCols entries, the solvers hand it a
106 * size x (size+1) buffer. A wider Jacobian would run over it. */
107 ✗ if(jacobian->sizeRows != size || jacobian->sizeCols != size)
108 {
109 ✗ linsys[i].jacobianIndex = -1;
110 ✗ throwStreamPrint(threadData, "Jacobian of torn linear system %d is %ux%u, but the system has size %d.", (int)linsys[i].equationIndex, jacobian->sizeRows, jacobian->sizeCols, size);
111 }
112 ✗ nnz = jacobian->sparsePattern->nnz;
113 ✗ linsys[i].nnz = nnz;
114
115 ✗ linsys[i].jacobian = jacobian;
116 }
117
118 /* -ls reaches klu and umfpack too, so naming a solver can override the format. */
119 ✗ if (omc_flag[FLAG_LSS]) {
120 ✗ linsys[i].useSparseSolver = 1;
121 ✗ } else if (omc_flag[FLAG_LS]) {
122 ✗ linsys[i].useSparseSolver = 0;
123 ✗ } else if (linsys[i].matrixFormat == OMC_MATRIX_SPARSE) {
124 ✗ linsys[i].useSparseSolver = 1;
125 }
126
127 /* Allocate nominal, min and max */
128 ✗ linsys[i].nominal = (double*) malloc(size*sizeof(double));
129 ✗ linsys[i].min = (double*) malloc(size*sizeof(double));
130 ✗ linsys[i].max = (double*) malloc(size*sizeof(double));
131
132 /* Init sparsitiy pattern */
133 ✗ linsys[i].initializeStaticLSData(data, threadData, &linsys[i], 1 /* true */);
134
135 /* allocate solver data */
136 /* the implementation of matrix A is solver-specific */
137 ✗ if(linsys[i].useSparseSolver == 1)
138 {
139 ✗ switch(data->simulationInfo->lssMethod)
140 {
141 #ifdef WITH_SUITESPARSE
142 ✗ case LSS_UMFPACK:
143 ✗ linsys[i].setAElement = setAElementUmfpack;
144 ✗ linsys[i].setBElement = setBElement;
145 ✗ allocateUmfPackData(size, size, nnz, linsys[i].solverData);
146 ✗ break;
147 ✗ case LSS_KLU:
148 ✗ linsys[i].setAElement = setAElementKlu;
149 ✗ linsys[i].setBElement = setBElement;
150 ✗ allocateKluData(size, size, nnz, linsys[i].solverData);
151 ✗ break;
152 #else
153 ✗ case LSS_KLU:
154 case LSS_UMFPACK:
155 ✗ throwStreamPrint(threadData, "OMC is compiled without UMFPACK, if you want use klu or umfpack please compile OMC with UMFPACK.");
156 break;
157 #endif
158 #if !defined(OMC_MINIMAL_RUNTIME)
159 ✗ case LSS_LIS:
160 ✗ linsys[i].setAElement = setAElementLis;
161 ✗ linsys[i].setBElement = setBElementLis;
162 ✗ allocateLisData(size, size, nnz, linsys[i].solverData);
163 ✗ break;
164 #else
165 ✗ case LSS_LIS_NOT_AVAILABLE:
166 ✗ throwStreamPrint(threadData, "OMC is compiled without sparse linear solver Lis.");
167 break;
168 #endif
169 #if defined(OMC_MINIMAL_RUNTIME) && !defined(WITH_SUITESPARSE)
170 ✗ case LSS_DEFAULT:
171 {
172 ✗ int indexes[2] = {1, linsys[i].equationIndex};
173 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_STDOUT, omc_dummyFileInfo, 0, indexes, "The simulation runtime does not have access to sparse solvers. Defaulting to a dense linear system solver instead.");
174 ✗ linsys[i].useSparseSolver = 0;
175 break;
176 }
177 #endif
178 ✗ default:
179 ✗ throwStreamPrint(threadData, "unrecognized sparse linear solver (%d)", data->simulationInfo->lssMethod);
180 }
181 }
182 ✗ if(linsys[i].useSparseSolver == 0) /* Not an else-statement because there might not be a sparse linear solver available */
183 {
184 ✗ switch(data->simulationInfo->lsMethod)
185 {
186 ✗ case LS_LAPACK:
187 ✗ linsys[i].setAElement = setAElement;
188 ✗ linsys[i].setBElement = setBElement;
189 ✗ linsys[i].A = (double*) malloc(size*size*sizeof(double));
190 ✗ allocateLapackData(size, linsys[i].solverData);
191 ✗ break;
192
193 #if !defined(OMC_MINIMAL_RUNTIME)
194 ✗ case LS_LIS:
195 ✗ linsys[i].setAElement = setAElementLis;
196 ✗ linsys[i].setBElement = setBElementLis;
197 ✗ allocateLisData(size, size, nnz, linsys[i].solverData);
198 ✗ break;
199 #endif
200 #ifdef WITH_SUITESPARSE
201 ✗ case LS_UMFPACK:
202 ✗ linsys[i].setAElement = setAElementUmfpack;
203 ✗ linsys[i].setBElement = setBElement;
204 ✗ allocateUmfPackData(size, size, nnz, linsys[i].solverData);
205 ✗ break;
206 ✗ case LS_KLU:
207 ✗ linsys[i].setAElement = setAElementKlu;
208 ✗ linsys[i].setBElement = setBElement;
209 ✗ allocateKluData(size, size, nnz, linsys[i].solverData);
210 ✗ break;
211 #else
212 ✗ case LS_UMFPACK:
213 ✗ throwStreamPrint(threadData, "OMC is compiled without UMFPACK, if you want use umfpack please compile OMC with UMFPACK.");
214 break;
215 #endif
216
217 ✗ case LS_TOTALPIVOT:
218 ✗ linsys[i].setAElement = setAElement;
219 ✗ linsys[i].setBElement = setBElement;
220 ✗ linsys[i].A = (double*) malloc(size*size*sizeof(double));
221 ✗ allocateTotalPivotData(size, linsys[i].solverData);
222 ✗ break;
223
224 ✗ case LS_DEFAULT:
225 ✗ linsys[i].setAElement = setAElement;
226 ✗ linsys[i].setBElement = setBElement;
227 ✗ linsys[i].A = (double*) malloc(size*size*sizeof(double));
228 ✗ allocateLapackData(size, linsys[i].solverData);
229 ✗ allocateTotalPivotData(size, linsys[i].solverData);
230 ✗ break;
231
232 ✗ default:
233 ✗ throwStreamPrint(threadData, "unrecognized dense linear solver (%d)", data->simulationInfo->lsMethod);
234 }
235 }
236 }
237
238 1 messageClose(OMC_LOG_LS);
239
240 1 return 0;
241 }
242
243 /**
244 * @brief Set min, max, nominal for linear systems.
245 *
246 * This function allocates memory for sparsity pattern and
247 * initialized nominal, min, max and spsarsity pattern.
248 *
249 * @param data Pointer to data.
250 * @param threadData Thread data for error handling.
251 * @return int Return 0.
252 */
253 1 int updateStaticDataOfLinearSystems(DATA *data, threadData_t *threadData)
254 {
255 int i, nnz;
256 int size;
257 1 LINEAR_SYSTEM_DATA *linsys = data->simulationInfo->linearSystemData;
258
259 1 infoStreamPrint(OMC_LOG_LS_V, 1, "update static data of linear system solvers");
260
261
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(i=0; i<data->modelData->nLinearSystems; ++i)
262 {
263 // check if LS is initialized
264 ✗ if (linsys[i].nominal == NULL || linsys[i].min == NULL || linsys[i].max==NULL)
265 {
266 ✗ throwStreamPrint(threadData, "Static data of Linear system not initialized for linear system %i",i);
267 }
268 ✗ linsys[i].initializeStaticLSData(data, threadData, &linsys[i], 0 /* false */);
269 }
270
271 1 messageClose(OMC_LOG_LS_V);
272
273 1 return 0;
274 }
275
276 /*! \fn int printLinearSystemStatistics(DATA *data)
277 *
278 * This function print memory for all linear systems.
279 *
280 * \param [ref] [data]
281 * [in] [sysNumber] index of corresponding linear System
282 */
283 ✗ void printLinearSystemSolvingStatistics(DATA *data, int sysNumber, int logLevel)
284 {
285 ✗ LINEAR_SYSTEM_DATA* linsys = data->simulationInfo->linearSystemData;
286 ✗ infoStreamPrint(logLevel, 1, "Linear system %d with (size = %d, nonZeroElements = %d, density = %.2f %%) solver statistics:",
287 ✗ (int)linsys[sysNumber].equationIndex, (int)linsys[sysNumber].size, (int)linsys[sysNumber].nnz,
288 ✗ (((double) linsys[sysNumber].nnz) / ((double)(linsys[sysNumber].size*linsys[sysNumber].size)))*100 );
289 ✗ infoStreamPrint(logLevel, 0, " number of calls : %ld", linsys[sysNumber].numberOfCall);
290 ✗ infoStreamPrint(logLevel, 0, " average time per call : %g", linsys[sysNumber].totalTime/linsys[sysNumber].numberOfCall);
291 ✗ infoStreamPrint(logLevel, 0, " time of jacobian evaluations : %g", linsys[sysNumber].jacobianTime);
292 ✗ infoStreamPrint(logLevel, 0, " total time : %g", linsys[sysNumber].totalTime);
293 ✗ messageClose(logLevel);
294 ✗ }
295
296 /*! \fn freeLinearSystems
297 *
298 * This function frees memory of linear systems.
299 *
300 * \param [ref] [data]
301 */
302 1 int freeLinearSystems(DATA *data, threadData_t *threadData)
303 {
304 int i;
305 1 LINEAR_SYSTEM_DATA* linsys = data->simulationInfo->linearSystemData;
306
307 1 infoStreamPrint(OMC_LOG_LS_V, 1, "free linear system solvers");
308
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for(i=0; i<data->modelData->nLinearSystems; ++i)
309 {
310 /* free system and solver data */
311 ✗ free(linsys[i].nominal); linsys[i].nominal = NULL;
312 ✗ free(linsys[i].min); linsys[i].min = NULL;
313 ✗ free(linsys[i].max); linsys[i].max = NULL;
314
315 ✗ free(linsys[i].b); linsys[i].b = NULL;
316
317 /* ToDo Implement unique function to free a JACOBIAN */
318 ✗ if (1 == linsys[i].method) {
319 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[linsys[i].jacobianIndex]);
320 ✗ freeJacobian(jacobian);
321 /* Note: The Jacobian of data->simulationInfo itself will be free later. */
322 /* jacobian only aliases the Jacobian freed above. */
323 ✗ linsys[i].jacobian = NULL;
324 }
325
326 ✗ if(linsys[i].useSparseSolver == 1)
327 {
328 ✗ switch(data->simulationInfo->lssMethod)
329 {
330 #if !defined(OMC_MINIMAL_RUNTIME)
331 ✗ case LSS_LIS:
332 ✗ freeLisData(linsys[i].solverData);
333 ✗ break;
334 #endif
335
336 #ifdef WITH_SUITESPARSE
337 ✗ case LSS_UMFPACK:
338 ✗ freeUmfPackData(linsys[i].solverData);
339 ✗ break;
340 ✗ case LSS_KLU:
341 ✗ freeKluData(linsys[i].solverData);
342 ✗ break;
343 #else
344 ✗ case LSS_UMFPACK:
345 ✗ throwStreamPrint(threadData, "OMC is compiled without UMFPACK, if you want use umfpack please compile OMC with UMFPACK.");
346 break;
347 #endif
348
349 ✗ default:
350 ✗ throwStreamPrint(threadData, "unrecognized sparse linear solver (%d)", data->simulationInfo->lssMethod);
351 }
352 }
353 else { // useSparseSolver == 0
354 ✗ switch(data->simulationInfo->lsMethod) {
355 ✗ case LS_LAPACK:
356 ✗ free(linsys[i].A);
357 ✗ freeLapackData(linsys[i].solverData);
358 ✗ break;
359
360 #if !defined(OMC_MINIMAL_RUNTIME)
361 ✗ case LS_LIS:
362 ✗ freeLisData(linsys[i].solverData);
363 ✗ break;
364 #endif
365
366 #ifdef WITH_SUITESPARSE
367 ✗ case LS_UMFPACK:
368 ✗ freeUmfPackData(linsys[i].solverData);
369 ✗ break;
370 ✗ case LS_KLU:
371 ✗ freeKluData(linsys[i].solverData);
372 ✗ break;
373 #else
374 ✗ case LS_UMFPACK:
375 ✗ throwStreamPrint(threadData, "OMC is compiled without UMFPACK, if you want use umfpack please compile OMC with UMFPACK.");
376 break;
377 #endif
378
379 ✗ case LS_TOTALPIVOT:
380 ✗ free(linsys[i].A);
381 ✗ freeTotalPivotData(linsys[i].solverData);
382 ✗ break;
383
384 ✗ case LS_DEFAULT:
385 ✗ free(linsys[i].A);
386 ✗ freeLapackData(linsys[i].solverData);
387 ✗ freeTotalPivotData(linsys[i].solverData);
388 ✗ break;
389
390 ✗ default:
391 ✗ throwStreamPrint(threadData, "unrecognized dense linear solver (%d)", data->simulationInfo->lsMethod);
392 }
393 }
394 }
395
396 1 messageClose(OMC_LOG_LS_V);
397
398 1 return 0;
399 }
400
401 /*! \fn solve linear system
402 *
403 * \param [ref] data
404 * [ref] threadData
405 * [in] sysNumber index of corresponding linear System
406 * [in/out] aux_x Auxiliary vector with values for x. Contains solution on output.
407 */
408 ✗ int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
409 {
410 int retVal;
411 int success;
412 int logLevel;
413 ✗ LINEAR_SYSTEM_DATA* linsys = &(data->simulationInfo->linearSystemData[sysNumber]);
414
415 ✗ if (!linsys->logActive) {
416 ✗ deactivateLogging();
417 }
418
419 ✗ rt_ext_tp_tick(&(linsys->totalTimeClock));
420
421 /* enable to avoid division by zero */
422 ✗ data->simulationInfo->noThrowDivZero = 1;
423
424 ✗ if(linsys->useSparseSolver == 1)
425 {
426 ✗ switch(data->simulationInfo->lssMethod)
427 {
428 #if !defined(OMC_MINIMAL_RUNTIME)
429 ✗ case LSS_LIS:
430 ✗ success = solveLis(data, threadData, sysNumber, aux_x);
431 ✗ break;
432 #else
433 ✗ case LSS_LIS_NOT_AVAILABLE:
434 ✗ throwStreamPrint(threadData, "OMC is compiled without UMFPACK, if you want use umfpack please compile OMC with UMFPACK.");
435 break;
436 #endif
437 #ifdef WITH_SUITESPARSE
438 ✗ case LSS_KLU:
439 ✗ success = solveKlu(data, threadData, sysNumber, aux_x);
440 ✗ break;
441 ✗ case LSS_UMFPACK:
442 ✗ success = solveUmfPack(data, threadData, sysNumber, aux_x);
443 ✗ if (!success && linsys->strictTearingFunctionCall != NULL){
444 debugString(OMC_LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
445 ✗ success = linsys->strictTearingFunctionCall(data, threadData);
446 ✗ if (success) success=2;
447 }
448 break;
449 #else
450 ✗ case LSS_KLU:
451 case LSS_UMFPACK:
452 ✗ throwStreamPrint(threadData, "OMC is compiled without UMFPACK, if you want use umfpack please compile OMC with UMFPACK.");
453 break;
454 #endif
455 ✗ default:
456 ✗ throwStreamPrint(threadData, "unrecognized sparse linear solver (%d)", data->simulationInfo->lssMethod);
457 }
458 }
459
460 else{
461 ✗ switch(data->simulationInfo->lsMethod)
462 {
463 ✗ case LS_LAPACK:
464 ✗ success = solveLapack(data, threadData, sysNumber, aux_x);
465 ✗ break;
466
467 #if !defined(OMC_MINIMAL_RUNTIME)
468 ✗ case LS_LIS:
469 ✗ success = solveLis(data, threadData, sysNumber, aux_x);
470 ✗ break;
471 #endif
472 #ifdef WITH_SUITESPARSE
473 ✗ case LS_KLU:
474 ✗ success = solveKlu(data, threadData, sysNumber, aux_x);
475 ✗ break;
476 ✗ case LS_UMFPACK:
477 ✗ success = solveUmfPack(data, threadData, sysNumber, aux_x);
478 ✗ if (!success && linsys->strictTearingFunctionCall != NULL){
479 debugString(OMC_LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
480 ✗ success = linsys->strictTearingFunctionCall(data, threadData);
481 ✗ if (success) success=2;
482 }
483 break;
484 #else
485 ✗ case LS_UMFPACK:
486 ✗ throwStreamPrint(threadData, "OMC is compiled without UMFPACK, if you want use umfpack please compile OMC with UMFPACK.");
487 break;
488 #endif
489
490 ✗ case LS_TOTALPIVOT:
491 ✗ success = solveTotalPivot(data, threadData, sysNumber, aux_x);
492 ✗ break;
493
494 ✗ case LS_DEFAULT:
495 ✗ success = solveLapack(data, threadData, sysNumber, aux_x);
496
497 /* check if solution process was successful, if not use alternative tearing set if available (dynamic tearing)*/
498 ✗ if (!success && linsys->strictTearingFunctionCall != NULL){
499 debugString(OMC_LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used.");
500 ✗ success = linsys->strictTearingFunctionCall(data, threadData);
501 ✗ if (success){
502 success=2;
503 ✗ linsys->failed = 0;
504 }
505 else {
506 ✗ linsys->failed = 1;
507 }
508 }
509 else{
510 /* if there is no alternative tearing set, use fallback solver */
511 ✗ if (!success){
512 ✗ if (linsys->failed){
513 logLevel = OMC_LOG_LS;
514 } else {
515 logLevel = OMC_LOG_STDOUT;
516 }
517 ✗ warningStreamPrintWithLimit(logLevel, 0, linsys->numberOfFailures, data->simulationInfo->maxWarnDisplays,
518 ✗ "The default linear solver fails, the fallback solver with total pivoting is started at time %f. That might raise performance issues, for more information use -lv LOG_LS.", data->localData[0]->timeValue);
519 ✗ success = solveTotalPivot(data, threadData, sysNumber, aux_x);
520 ✗ linsys->failed = 1;
521 } else {
522 ✗ linsys->failed = 0;
523 }
524 }
525 break;
526
527 ✗ default:
528 ✗ throwStreamPrint(threadData, "unrecognized dense linear solver (%d)", data->simulationInfo->lsMethod);
529 }
530 }
531 ✗ linsys->solved = success;
532
533 ✗ linsys->totalTime += rt_ext_tp_tock(&(linsys->totalTimeClock));
534 ✗ linsys->numberOfCall++;
535
536 ✗ retVal = check_linear_solution(data, 1, sysNumber);
537
538 ✗ if (!linsys->logActive) {
539 ✗ reactivateLogging();
540 }
541
542 ✗ return retVal;
543 }
544
545 /*! \fn check_linear_solutions
546 *
547 * This function check whether some of linear systems
548 * are failed to solve. If one is failed it returns 1 otherwise 0.
549 *
550 * \param [in] [data]
551 * \param [in] [printFailingSystems]
552 * \param [out] [returnValue] It returns >0 if fail otherwise 0.
553 *
554 * \author wbraun
555 */
556 2 int check_linear_solutions(DATA *data, int printFailingSystems)
557 {
558 long i;
559
560
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 for(i=0; i<data->modelData->nLinearSystems; ++i)
561 {
562 ✗ if(check_linear_solution(data, printFailingSystems, i))
563 {
564 return 1;
565 }
566 }
567
568 return 0;
569 }
570
571 /*! \fn check_linear_solution
572 * This function check whether some of linear systems
573 * are failed to solve. If one is failed it returns 1 otherwise 0.
574 *
575 * \param [in] [data]
576 * \param [in] [printFailingSystems]
577 * \param [in] [sysNumber] index of corresponding linear System
578 * \param [out] [returnValue] It returns 1 if fail otherwise 0.
579 *
580 * \author wbraun
581 */
582 ✗ int check_linear_solution(DATA *data, int printFailingSystems, int sysNumber)
583 {
584 ✗ LINEAR_SYSTEM_DATA* linsys = data->simulationInfo->linearSystemData;
585 long j, i = sysNumber;
586
587 const size_t buff_size = 2048;
588 char *start_buffer;
589 char *nominal_buffer;
590
591 ✗ if(linsys[i].solved == 0)
592 {
593 ✗ int index = linsys[i].equationIndex, indexes[2] = {1,index};
594 ✗ if (!printFailingSystems)
595 {
596 return 1;
597 }
598 ✗ warningStreamPrintWithEquationIndexes(OMC_LOG_STDOUT, omc_dummyFileInfo, 1, indexes, "Solving linear system %d fails at time %g. For more information use -lv LOG_LS.", index, data->localData[0]->timeValue);
599
600 ✗ start_buffer = (char*) malloc(buff_size * sizeof(char));
601 ✗ assertStreamPrint(NULL, start_buffer != NULL, "Out of memory.");
602 ✗ nominal_buffer = (char*) malloc(buff_size * sizeof(char));
603 ✗ assertStreamPrint(NULL, nominal_buffer != NULL, "Out of memory.");
604
605 ✗ for(j=0; j<modelInfoGetEquation(&data->modelData->modelDataXml, (linsys[i]).equationIndex).numVar; ++j)
606 {
607 int done=0;
608 long k;
609 ✗ const MODEL_DATA *mData = data->modelData;
610 ✗ for(k=0; k<mData->nVariablesRealArray && !done; ++k)
611 {
612 ✗ if (!strcmp(mData->realVarsData[k].info.name, modelInfoGetEquation(&data->modelData->modelDataXml, (linsys[i]).equationIndex).vars[j]))
613 {
614 done = 1;
615 ✗ real_vector_to_string(&mData->realVarsData[k].attribute.start, mData->realVarsData[k].dimension.numberOfDimensions == 0, start_buffer, buff_size);
616 ✗ real_vector_to_string(&mData->realVarsData[k].attribute.nominal, mData->realVarsData[k].dimension.numberOfDimensions == 0, nominal_buffer, buff_size);
617 ✗ warningStreamPrint(OMC_LOG_LS, 0, "[%ld] Real %s(start=%s, nominal=%s)",
618 j+1,
619 ✗ mData->realVarsData[k].info.name,
620 start_buffer,
621 nominal_buffer);
622 }
623 }
624 ✗ if (!done)
625 {
626 ✗ warningStreamPrint(OMC_LOG_LS, 0, "[%ld] Real %s(start=?, nominal=?)", j+1, modelInfoGetEquation(&data->modelData->modelDataXml, (linsys[i]).equationIndex).vars[j]);
627 }
628 }
629 ✗ messageCloseWarning(OMC_LOG_STDOUT);
630 ✗ free(start_buffer);
631 ✗ free(nominal_buffer);
632
633 ✗ return 1;
634 }
635
636 ✗ if(linsys[i].solved == 2)
637 {
638 ✗ linsys[i].solved = 1;
639 ✗ return 2;
640 }
641
642 return 0;
643 }
644
645 /*! \fn setAElement
646 * This function sets the (col, row)-value of linsys->A.
647 *
648 * \param [in] [row]
649 * \param [in] [col]
650 * \param [in] [value]
651 * \param [in] [nth] number element in matrix,
652 * is ignored here, used only for sparse
653 * \param [ref] [data]
654 *
655 */
656 ✗ static void setAElement(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t* threadData)
657 {
658 ✗ linearSystemData->A[row + col * linearSystemData->size] = value;
659 ✗ }
660
661 /*! \fn setBElement
662 * This function sets the row-th value of linsys->b[row] = value.
663 *
664 * \param [in] [row]
665 * \param [in] [value]
666 * \param [ref] [data]
667 */
668 ✗ static void setBElement(int row, double value, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t *threadData)
669 {
670 ✗ linearSystemData->b[row] = value;
671 ✗ }
672
673 #if !defined(OMC_MINIMAL_RUNTIME)
674 ✗ static void setAElementLis(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t *threadData)
675 {
676 ✗ DATA_LIS* sData = (DATA_LIS*) linearSystemData->solverData[0];
677 ✗ lis_matrix_set_value(LIS_INS_VALUE, row, col, value, sData->A);
678 ✗ }
679
680 ✗ static void setBElementLis(int row, double value, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t *threadData)
681 {
682 ✗ DATA_LIS* sData = (DATA_LIS*) linearSystemData->solverData[0];
683 ✗ lis_vector_set_value(LIS_INS_VALUE, row, value, sData->b);
684 ✗ }
685 #endif
686
687 #ifdef WITH_SUITESPARSE
688 ✗ static void setAElementUmfpack(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t *threadData)
689 {
690 ✗ DATA_UMFPACK* sData = (DATA_UMFPACK*) linearSystemData->solverData[0];
691
692 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, " set %d. -> (%d,%d) = %f", nth, row, col, value);
693 ✗ if (row > 0) {
694 ✗ if (sData->Ap[row] == 0) {
695 ✗ sData->Ap[row] = nth;
696 }
697 }
698
699 ✗ sData->Ai[nth] = col;
700 ✗ sData->Ax[nth] = value;
701 ✗ }
702
703 ✗ static void setAElementKlu(int row, int col, double value, int nth, LINEAR_SYSTEM_DATA* linearSystemData, threadData_t *threadData)
704 {
705 ✗ DATA_KLU* sData = (DATA_KLU*) linearSystemData->solverData[0];
706
707 ✗ if (row > 0) {
708 ✗ if (sData->Ap[row] == 0) {
709 ✗ sData->Ap[row] = nth;
710 }
711 }
712
713 ✗ sData->Ai[nth] = col;
714 ✗ sData->Ax[nth] = value;
715 ✗ }
716
717 #endif
718