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 / 99
Functions: 0.0% 0 / 0 / 4
Branches: 0.0% 0 / 0 / 38
OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.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 "linearSolverLapack.h"
45
46 extern int dgesv_(int *n, int *nrhs, double *a, int *lda,
47 int *ipiv, double *b, int *ldb, int *info);
48
49 extern int dgetrs_(char* tran, int *n, int *nrhs, double *a, int *lda,
50 int *ipiv, double *b, int *ldb, int *info);
51 /*! \fn allocate memory for linear system solver lapack
52 *
53 */
54 ✗ int allocateLapackData(int size, void** voiddata)
55 {
56 ✗ DATA_LAPACK* data = (DATA_LAPACK*) calloc(1, sizeof(DATA_LAPACK));
57
58 ✗ data->ipiv = (int*) calloc(size, sizeof(int));
59 ✗ assertStreamPrint(NULL, 0 != data->ipiv, "Could not allocate data for linear solver lapack.");
60 ✗ data->nrhs = 1;
61 ✗ data->info = 0;
62 ✗ data->work = _omc_allocateVectorData(size);
63
64 ✗ data->x = _omc_createVector(size, NULL);
65 ✗ data->b = _omc_createVector(size, NULL);
66 ✗ data->A = _omc_createMatrix(size, size, NULL);
67
68 ✗ *voiddata = (void*)data;
69 ✗ return 0;
70 }
71
72 /*! \fn free memory of lapack
73 *
74 */
75 ✗ int freeLapackData(void **voiddata)
76 {
77 ✗ DATA_LAPACK* data = (DATA_LAPACK*) *voiddata;
78
79 ✗ free(data->ipiv);
80 ✗ _omc_deallocateVectorData(data->work);
81
82 ✗ _omc_destroyVector(data->x);
83 ✗ _omc_destroyVector(data->b);
84 ✗ _omc_destroyMatrix(data->A);
85
86 ✗ free(data);
87 ✗ voiddata[0] = NULL;
88
89 ✗ return 0;
90 }
91
92 /*! \fn getAnalyticalJacobian
93 *
94 * function calculates analytical jacobian
95 *
96 * \param [ref] [data]
97 * \param [out] [jac]
98 *
99 * \author wbraun
100 *
101 */
102 ✗ void getAnalyticalJacobianLapack(DATA* data, threadData_t *threadData, LINEAR_SYSTEM_DATA* systemData, double* jac)
103 {
104 int k;
105 ✗ JACOBIAN* jacobian = systemData->jacobian;
106 ✗ JACOBIAN* parentJacobian = systemData->parentJacobian;
107
108 /* call generic dense Jacobian */
109 ✗ evalJacobian(data, threadData, jacobian, parentJacobian, jac, TRUE);
110
111 ✗ for (k = 0; k < (jacobian->sizeRows) * (jacobian->sizeCols); k++)
112 ✗ jac[k] = -jac[k];
113 ✗ }
114
115 /*! \fn wrapper_fvec_lapack for the residual function
116 *
117 */
118 static int wrapper_fvec_lapack(_omc_vector* x, _omc_vector* f, int* iflag, RESIDUAL_USERDATA* resUserData, int sysNumber)
119 {
120 ✗ resUserData->data->simulationInfo->linearSystemData[sysNumber].residualFunc(resUserData, x->data, f->data, iflag);
121 ✗ return 0;
122 }
123
124 /*! \fn solve linear system with lapack method
125 *
126 * \param [in] [data]
127 * [sysNumber] index of the corresponding linear system
128 *
129 * \author wbraun
130 */
131 ✗ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
132 {
133 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
134 ✗ int i, iflag = 1;
135 ✗ LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
136
137 ✗ DATA_LAPACK* solverData = (DATA_LAPACK*) systemData->solverData[0];
138 int success = 1;
139
140 /* We are given the number of the linear system.
141 * We want to look it up among all equations. */
142 ✗ int eqSystemNumber = systemData->equationIndex;
143 ✗ int indexes[2] = {1,eqSystemNumber};
144 _omc_scalar residualNorm = 0;
145 double tmpJacEvalTime;
146 ✗ int reuseMatrixJac = (data->simulationInfo->currentContext == CONTEXT_SYM_JACOBIAN && data->simulationInfo->currentJacobianEval > 0);
147
148 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_LS, omc_dummyFileInfo, 0, indexes,
149 "Start solving Linear System %d (size %d) at time %g with Lapack Solver",
150 ✗ eqSystemNumber, (int) systemData->size, data->localData[0]->timeValue);
151
152 /* set data */
153 ✗ _omc_setVectorData(solverData->x, aux_x);
154 ✗ _omc_setVectorData(solverData->b, systemData->b);
155 ✗ _omc_setMatrixData(solverData->A, systemData->A);
156
157 ✗ rt_ext_tp_tick(&(solverData->timeClock));
158 ✗ if (0 == systemData->method) {
159
160 ✗ if (!reuseMatrixJac) {
161 /* reset matrix A */
162 ✗ memset(systemData->A, 0, (systemData->size)*(systemData->size)*sizeof(double));
163 /* update matrix A */
164 ✗ systemData->setA(data, threadData, systemData);
165 }
166
167 /* update vector b (rhs) */
168 ✗ systemData->setb(data, threadData, systemData);
169 } else {
170 ✗ if (!reuseMatrixJac) {
171 /* calculate jacobian -> matrix A*/
172 ✗ if(systemData->jacobianIndex != -1) {
173 ✗ getAnalyticalJacobianLapack(data, threadData, systemData, solverData->A->data);
174 } else {
175 assertStreamPrint(threadData, 1, "jacobian function pointer is invalid" );
176 }
177 }
178 /* calculate vector b (rhs) */
179 ✗ _omc_copyVector(solverData->work, solverData->x);
180
181 ✗ wrapper_fvec_lapack(solverData->work, solverData->b, &iflag, &resUserData, sysNumber);
182 }
183 ✗ tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
184 ✗ systemData->jacobianTime += tmpJacEvalTime;
185 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
186
187 /* Log A*x=b */
188 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_LS_V)){
189 ✗ _omc_printVector(solverData->x, "Vector old x", OMC_LOG_LS_V);
190 ✗ _omc_printMatrix(solverData->A, "Matrix A", OMC_LOG_LS_V);
191 ✗ _omc_printVector(solverData->b, "Vector b", OMC_LOG_LS_V);
192 }
193
194 ✗ rt_ext_tp_tick(&(solverData->timeClock));
195
196 /* if reuseMatrixJac use also previous factorization */
197 ✗ if (!reuseMatrixJac)
198 {
199 /* Solve system */
200 ✗ dgesv_((int*) &systemData->size,
201 (int*) &solverData->nrhs,
202 ✗ solverData->A->data,
203 (int*) &systemData->size,
204 solverData->ipiv,
205 ✗ solverData->b->data,
206 ✗ (int*) &systemData->size,
207 &solverData->info);
208
209 } /* further Jacobian evaluations */
210 else
211 {
212 ✗ char trans = 'N';
213 /* Solve system */
214 ✗ dgetrs_(&trans,
215 (int*) &systemData->size,
216 (int*) &solverData->nrhs,
217 ✗ solverData->A->data,
218 (int*) &systemData->size,
219 solverData->ipiv,
220 ✗ solverData->b->data,
221 ✗ (int*) &systemData->size,
222 &solverData->info);
223 }
224
225
226 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
227
228 ✗ if(solverData->info < 0)
229 {
230 ✗ warningStreamPrint(OMC_LOG_LS, 0, "Error solving linear system of equations (no. %d) at time %f. Argument %d illegal.", (int)systemData->equationIndex, data->localData[0]->timeValue, (int)solverData->info);
231 success = 0;
232 }
233 ✗ else if(solverData->info > 0)
234 {
235 ✗ warningStreamPrintWithLimit(OMC_LOG_LS, 0, ++(systemData->numberOfFailures) /* Update counter */, data->simulationInfo->maxWarnDisplays,
236 "Failed to solve linear system of equations (no. %d) at time %f, system is singular for U[%d, %d].",
237 ✗ (int)systemData->equationIndex, data->localData[0]->timeValue, (int)solverData->info+1, (int)solverData->info+1);
238
239 success = 0;
240
241 /* debug output */
242 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_LS)){
243 ✗ _omc_printMatrix(solverData->A, "Matrix U", OMC_LOG_LS);
244
245 ✗ _omc_printVector(solverData->b, "Output vector x", OMC_LOG_LS);
246 }
247 }
248
249 if (1 == success){
250
251 ✗ if (1 == systemData->method){
252 /* take the solution */
253 ✗ solverData->x = _omc_addVectorVector(solverData->x, solverData->work, solverData->b); // x = xold(work) + xnew(b)
254
255 /* update inner equations */
256 ✗ wrapper_fvec_lapack(solverData->x, solverData->work, &iflag, &resUserData, sysNumber);
257 ✗ residualNorm = _omc_euclideanVectorNorm(solverData->work);
258
259 ✗ if ((isnan(residualNorm)) || (residualNorm>1e-4)){
260 ✗ warningStreamPrintWithLimit(OMC_LOG_LS, 0, ++(systemData->numberOfFailures) /* Update counter */, data->simulationInfo->maxWarnDisplays,
261 "Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
262 ✗ (int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
263 success = 0;
264 }
265 } else {
266 /* take the solution */
267 ✗ _omc_copyVector(solverData->x, solverData->b);
268 }
269
270 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_LS_V)) {
271 ✗ if (1 == systemData->method) {
272 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
273 } else {
274 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Solution x:");
275 }
276 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
277
278 ✗ for(i = 0; i < systemData->size; ++i) {
279 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "[%d] %s = %.15g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
280 }
281
282 ✗ messageClose(OMC_LOG_LS_V);
283 }
284 }
285
286 ✗ return success;
287 }
288