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 / 187
Functions: 0.0% 0 / 0 / 16
Branches: 0.0% 0 / 0 / 74
OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.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 nonlinear_solver.c
29 */
30
31 #include <math.h>
32 #include <stdlib.h>
33 #include <string.h> /* memcpy */
34
35 #include "../../simulation_data.h"
36 #include "../simulation_info_json.h"
37 #include "../jacobian_util.h"
38 #include "../../util/omc_error.h"
39 #include "omc_math.h"
40 #include "../../util/varinfo.h"
41 #include "model_help.h"
42
43 #include "linearSystem.h"
44 #include "linearSolverTotalPivot.h"
45
46 ✗ void debugMatrixDoubleLS(int logName, char* matrixName, double* matrix, int n, int m)
47 {
48 ✗ if(OMC_ACTIVE_STREAM(logName))
49 {
50 int i, j;
51 int sparsity = 0;
52 ✗ char *buffer = (char*)malloc(sizeof(char)*m*18);
53
54 ✗ infoStreamPrint(logName, 1, "%s [%dx%d-dim]", matrixName, n, m);
55 ✗ for(i=0; i<n;i++)
56 {
57 char *p = buffer;
58 ✗ for(j=0; j<m; j++)
59 {
60 if (sparsity)
61 {
62 if (fabs(matrix[i + j*(m-1)])<1e-12)
63 p += sprintf(p, " 0");
64 else
65 p += sprintf(p, " *");
66 }
67 else
68 {
69 ✗ p += sprintf(p, " %12.4g", matrix[i + j*(m-1)]);
70 }
71 }
72 ✗ infoStreamPrint(logName, 0, "%s", buffer);
73 }
74 ✗ messageClose(logName);
75 ✗ free(buffer);
76 }
77 ✗ }
78
79 ✗ void debugVectorDoubleLS(int logName, char* vectorName, double* vector, int n)
80 {
81 ✗ if(OMC_ACTIVE_STREAM(logName))
82 {
83 int i;
84 ✗ char *buffer = (char*)malloc(sizeof(char)*n*22);
85
86 ✗ infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
87 {
88 char *p = buffer;
89 ✗ if (vector[0]<-1e+300)
90 ✗ p += sprintf(p, " -INF");
91 ✗ else if (vector[0]>1e+300)
92 ✗ p += sprintf(p, " +INF");
93 else
94 ✗ p += sprintf(p, " %16.8g", vector[0]);
95 ✗ for(i=1; i<n;i++)
96 {
97 ✗ if (vector[i]<-1e+300)
98 ✗ p += sprintf(p, " -INF");
99 ✗ else if (vector[i]>1e+300)
100 ✗ p += sprintf(p, " +INF");
101 else
102 ✗ p += sprintf(p, " %16.8g", vector[i]);
103 }
104 }
105 ✗ infoStreamPrint(logName, 0, "%s", buffer);
106 ✗ free(buffer);
107 ✗ messageClose(logName);
108 }
109 ✗ }
110
111 ✗ void debugStringLS(int logName, char* message)
112 {
113 ✗ infoStreamPrint(logName, 0, "%s", message);
114 ✗ }
115
116 ✗ void debugIntLS(int logName, char* message, int value)
117 {
118 ✗ infoStreamPrint(logName, 1, "%s %d", message, value);
119 ✗ }
120
121 ✗ void vecMultScalingLS(int n, double *a, double *b, double *c)
122 {
123 int i;
124 ✗ for (i=0;i<n;i++)
125 ✗ c[i] = a[i]*fabs(b[i]);
126 ✗ }
127
128 ✗ void vecAddScalLS(int n, double *a, double *b, double s, double *c)
129 {
130 int i;
131 ✗ for (i=0;i<n;i++)
132 ✗ c[i] = a[i] + s*b[i];
133 ✗ }
134
135 ✗ void vecAddLS(int n, double *a, double *b, double *c)
136 {
137 int i;
138 ✗ for (i=0;i<n;i++)
139 ✗ c[i] = a[i] + b[i];
140 ✗ }
141
142 ✗ void vecCopyLS(int n, double *a, double *b)
143 {
144 ✗ memcpy(b, a, n*(sizeof(double)));
145 ✗ }
146
147 ✗ void vecConstLS(int n, double value, double *a)
148 {
149 int i;
150 ✗ for (i=0;i<n;i++)
151 ✗ a[i] = value;
152 ✗ }
153
154 ✗ void vecScalarMultLS(int n, double *a, double s, double *b)
155 {
156 int i;
157 ✗ for (i=0;i<n;i++)
158 ✗ b[i] = s*a[i];
159 ✗ }
160
161 ✗ void getIndicesOfPivotElementLS(int *n, int *m, int *l, double* A, int *indRow, int *indCol, int *pRow, int *pCol, double *absMax)
162 {
163 int i, j;
164
165 ✗ *absMax = fabs(A[indRow[*l] + indCol[*l]* *n]);
166 ✗ *pCol = *l;
167 ✗ *pRow = *l;
168 ✗ for (i = *l; i < *n; i++) {
169 ✗ for (j = *l; j < *m; j++) {
170 ✗ if (fabs(A[indRow[i] + indCol[j]* *n]) > *absMax) {
171 ✗ *absMax = fabs(A[indRow[i] + indCol[j]* *n]);
172 ✗ *pCol = j;
173 ✗ *pRow = i;
174 }
175 }
176 }
177 ✗ }
178
179 /**
180 * @brief Linear solver for A*x = b based on a total pivot search.
181 *
182 * \author bbachmann
183 *
184 * @param data Simulation data.
185 * @param n Size of matrix a
186 * @param x On return: Solution dim n+1, last column is 1 for solvable systems.
187 * @param Ab Matrix A|b: first n columns are matrix A, last column is -b
188 * @param indRow Work array for row indices, used for coloring.
189 * @param indCol Work array for column indices, used for coloring.
190 * @param rank On return: Rank of matrix A|b.
191 * @return int Return 0 on success, -1 if system is under-determined.
192 */
193 ✗ int solveSystemWithTotalPivotSearchLS(DATA* data, int n, double* x, double* Ab, int* indRow, int* indCol, int *rank)
194 {
195 ✗ int i, k, j, l, m=n+1, nrsh=1, singular=0;
196 int pCol, pRow;
197 double hValue;
198 double hInt;
199 double absMax;
200 int r,s;
201 int permutation = 1;
202
203 /* assume full rank of matrix A|b [n x (n+1)] */
204 ✗ *rank = n;
205
206 ✗ for (i=0; i<n; i++) {
207 ✗ indRow[i] = i;
208 }
209 ✗ for (i=0; i<m; i++) {
210 ✗ indCol[i] = i;
211 }
212
213 ✗ for (i = 0; i < n; i++) {
214 ✗ getIndicesOfPivotElementLS(&n, &n, &i, Ab, indRow, indCol, &pRow, &pCol, &absMax);
215 /* this criteria should be evaluated and may be improved in future */
216 ✗ if (absMax < DBL_EPSILON) {
217 ✗ *rank = i;
218 ✗ if (data->simulationInfo->initial) {
219 ✗ warningStreamPrint(OMC_LOG_LS, 1, "Total Pivot: Matrix (nearly) singular at initialization.");
220 } else {
221 ✗ warningStreamPrint(OMC_LOG_LS, 1, "Total Pivot: Matrix (nearly) singular at time %f.", data->localData[0]->timeValue);
222 }
223 ✗ warningStreamPrint(OMC_LOG_LS, 0, "Continuing anyway. For more information please use -lv %s.", OMC_LOG_STREAM_NAME[OMC_LOG_LS]);
224 ✗ messageCloseWarning(OMC_LOG_LS);
225 ✗ infoStreamPrint(OMC_LOG_LS, 0, "rank = %u", *rank);
226 ✗ break;
227 }
228 /* swap row indices */
229 ✗ if (pRow!=i) {
230 ✗ hInt = indRow[i];
231 ✗ indRow[i] = indRow[pRow];
232 ✗ indRow[pRow] = hInt;
233 }
234 /* swap column indices */
235 ✗ if (pCol!=i) {
236 ✗ hInt = indCol[i];
237 ✗ indCol[i] = indCol[pCol];
238 ✗ indCol[pCol] = hInt;
239 }
240
241 /* Gauss elimination of row indRow[i] */
242 ✗ for (k=i+1; k<n; k++) {
243 ✗ hValue = -Ab[indRow[k] + indCol[i]*n]/Ab[indRow[i] + indCol[i]*n];
244 ✗ for (j=i+1; j<m; j++) {
245 ✗ Ab[indRow[k] + indCol[j]*n] = Ab[indRow[k] + indCol[j]*n] + hValue*Ab[indRow[i] + indCol[j]*n];
246 }
247 ✗ Ab[indRow[k] + indCol[i]*n] = 0;
248 }
249 }
250
251 ✗ debugMatrixDoubleLS(OMC_LOG_LS_V,"LGS: matrix Ab manipulated",Ab, n, n+1);
252 /* solve even singular matrix */
253 ✗ for (i=n-1;i>=0; i--) {
254 ✗ if (i>=*rank) {
255 /* this criteria should be evaluated and may be improved in future */
256 ✗ if (fabs(Ab[indRow[i] + n*n])>1e-12) {
257 ✗ warningStreamPrint(OMC_LOG_LS, 0, "under-determined linear system not solvable!");
258 ✗ return -1;
259 } else {
260 ✗ x[indCol[i]] = 0.0;
261 }
262 } else {
263 ✗ x[indCol[i]] = -Ab[indRow[i] + n*n];
264 ✗ for (j=n-1; j>i; j--) {
265 ✗ x[indCol[i]] = x[indCol[i]] - Ab[indRow[i] + indCol[j]*n]*x[indCol[j]];
266 }
267 ✗ x[indCol[i]]=x[indCol[i]]/Ab[indRow[i] + indCol[i]*n];
268 }
269 }
270 ✗ x[n]=1.0;
271 ✗ debugVectorDoubleLS(OMC_LOG_LS_V,"LGS: solution vector x",x, n+1);
272
273 ✗ return 0;
274 }
275
276 /*! \fn allocate memory for linear system solver totalpivot
277 *
278 * \author bbachmann
279 */
280 ✗ int allocateTotalPivotData(int size, void** voiddata)
281 {
282 ✗ DATA_TOTALPIVOT* data = (DATA_TOTALPIVOT*) malloc(sizeof(DATA_TOTALPIVOT));
283
284 /* memory for linear system */
285 ✗ data->Ab = (double*) calloc((size*(size+1)),sizeof(double));
286 ✗ data->b = (double*) malloc(size*sizeof(double));
287 ✗ data->x = (double*) calloc(size+1,sizeof(double));
288
289 /* used for pivot strategy */
290 ✗ data->indRow =(int*) calloc(size,sizeof(int));
291 ✗ data->indCol =(int*) calloc(size+1,sizeof(int));
292
293 ✗ voiddata[1] = (void*)data;
294 ✗ return 0;
295 }
296
297 /*! \fn free memory for nonlinear solver totalpivot
298 *
299 * \author bbachmann
300 */
301 ✗ int freeTotalPivotData(void** voiddata)
302 {
303 ✗ DATA_TOTALPIVOT* data = (DATA_TOTALPIVOT*) voiddata[1];
304
305 /* memory for linear system */
306 ✗ free(data->Ab);
307 ✗ free(data->b);
308 ✗ free(data->x);
309
310 /* used for pivot strategy */
311 ✗ free(data->indRow);
312 ✗ free(data->indCol);
313
314 ✗ free(voiddata[1]);
315 ✗ voiddata[1] = NULL;
316
317 ✗ return 0;
318 }
319
320 /*! \fn getAnalyticalJacobian
321 *
322 * function calculates analytical jacobian
323 *
324 * \param [ref] [data]
325 * \param [out] [jac]
326 *
327 * \author wbraun
328 *
329 */
330 ✗ void getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, LINEAR_SYSTEM_DATA* systemData, modelica_real* jac)
331 {
332 ✗ JACOBIAN* jacobian = systemData->jacobian;
333 ✗ JACOBIAN* parentJacobian = systemData->parentJacobian;
334
335 /* call generic dense Jacobian */
336 ✗ evalJacobian(data, threadData, jacobian, parentJacobian, jac, TRUE);
337 ✗ }
338
339 /*! \fn wrapper_fvec_hybrd for the residual Function
340 * calls for the subroutine fcn(n, x, fvec, iflag, data)
341 *
342 *
343 */
344 static int wrapper_fvec_totalpivot(double* x, double* f, RESIDUAL_USERDATA* resUserData, int sysNumber)
345 {
346 int currentSys = sysNumber;
347 ✗ int iflag = 0;
348
349 ✗ resUserData->data->simulationInfo->linearSystemData[currentSys].residualFunc(resUserData, x, f, &iflag);
350 ✗ return 0;
351 }
352
353 /**
354 * @brief Solve linear system with total pivot method.
355 *
356 * \author bbachmann
357 *
358 * @param data Runtime data struct.
359 * @param threadData Thread data for error handling.
360 * @param sysNumber Index of the corresponding non-linear system.
361 * @param aux_x Work array with old values of x. Will be overwritten with solution.
362 * @return int Return 1 on success and 0 on failure.
363 */
364 ✗ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
365 {
366 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
367 int i, j;
368 ✗ LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
369 ✗ DATA_TOTALPIVOT* solverData = (DATA_TOTALPIVOT*) systemData->solverData[1];
370
371 ✗ int n = systemData->size, status;
372 double fdeps = 1e-8;
373 double xTol = 1e-8;
374 ✗ int eqSystemNumber = systemData->equationIndex;
375 ✗ int indexes[2] = {1,eqSystemNumber};
376 int rank;
377 _omc_scalar residualNorm = 0;
378
379 /* We are given the number of the linear system.
380 * We want to look it up among all equations. */
381 /* int eqSystemNumber = systemData->equationIndex; */
382 int success = 1;
383 double tmpJacEvalTime;
384
385 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_LS, omc_dummyFileInfo, 0, indexes,
386 "Start solving Linear System %d (size %d) at time %g with Total Pivot Solver",
387 ✗ eqSystemNumber, (int) systemData->size, data->localData[0]->timeValue);
388
389 ✗ debugVectorDoubleLS(OMC_LOG_LS_V,"SCALING",systemData->nominal,n);
390 ✗ debugVectorDoubleLS(OMC_LOG_LS_V,"Old VALUES",aux_x,n);
391
392 ✗ rt_ext_tp_tick(&(solverData->timeClock));
393 ✗ if (0 == systemData->method) {
394
395 /* reset matrix A */
396 ✗ vecConstLS(n*n, 0.0, systemData->A);
397 /* update matrix A -> first n columns of matrix Ab*/
398 ✗ systemData->setA(data, threadData, systemData);
399 ✗ vecCopyLS(n*n, systemData->A, solverData->Ab);
400
401 /* update vector b (rhs) -> -b is last column of matrix Ab*/
402 ✗ rt_ext_tp_tick(&(solverData->timeClock));
403 ✗ systemData->setb(data, threadData, systemData);
404 ✗ vecScalarMultLS(n, systemData->b, -1.0, solverData->Ab + n*n);
405 } else {
406
407 /* calculate jacobian -> first n columns of matrix Ab*/
408 ✗ if(systemData->jacobianIndex != -1){
409 ✗ getAnalyticalJacobianTotalPivot(data, threadData, systemData, solverData->Ab);
410 } else {
411 assertStreamPrint(threadData, 1, "jacobian function pointer is invalid" );
412 }
413 /* calculate vector b (rhs) -> -b is last column of matrix Ab */
414 ✗ wrapper_fvec_totalpivot(aux_x, solverData->Ab + n*n, &resUserData, sysNumber);
415 }
416 ✗ tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
417 ✗ systemData->jacobianTime += tmpJacEvalTime;
418 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
419 ✗ debugMatrixDoubleLS(OMC_LOG_LS_V,"LGS: matrix Ab",solverData->Ab, n, n+1);
420
421 ✗ rt_ext_tp_tick(&(solverData->timeClock));
422 ✗ status = solveSystemWithTotalPivotSearchLS(data, n, solverData->x, solverData->Ab, solverData->indRow, solverData->indCol, &rank);
423 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
424
425 ✗ if (status != 0) {
426 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Error solving linear system of equations (no. %d) at time %f.", (int)systemData->equationIndex, data->localData[0]->timeValue);
427 success = 0;
428 } else {
429 ✗ debugVectorDoubleLS(OMC_LOG_LS_V, "SOLUTION:", solverData->x, n+1);
430 ✗ if (1 == systemData->method) {
431 /* add the solution to old solution vector*/
432 ✗ vecAddLS(n, aux_x, solverData->x, aux_x);
433 ✗ wrapper_fvec_totalpivot(aux_x, solverData->b, &resUserData, sysNumber);
434 } else {
435 /* take the solution */
436 ✗ vecCopyLS(n, solverData->x, aux_x);
437 }
438
439 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_LS_V)) {
440 ✗ if (1 == systemData->method) {
441 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
442 } else {
443 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Solution x:");
444 }
445 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
446 ✗ for(i=0; i<systemData->size; ++i)
447 {
448 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
449 }
450 ✗ messageClose(OMC_LOG_LS_V);
451 }
452 }
453 ✗ return success;
454 }
455