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 / 184
Functions: 0.0% 0 / 0 / 12
Branches: 0.0% 0 / 0 / 244
OMCompiler/SimulationRuntime/cpp/Solver/LinearSolver/LinearSolver.cpp
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 #include <Core/ModelicaDefine.h>
29 #include <Core/Modelica.h>
30 /** @addtogroup solverLinearSolver
31 *
32 * @{
33 */
34
35 #include <Core/Math/ILapack.h>
36 #include <Solver/LinearSolver/FactoryExport.h>
37 #include <Core/Utils/extension/logger.hpp>
38 #include <Solver/LinearSolver/LinearSolver.h>
39
40
41 ✗ LinearSolver::LinearSolver(ILinSolverSettings* settings,shared_ptr<ILinearAlgLoop> algLoop)
42 :AlgLoopSolverDefaultImplementation()
43 , _algLoop (algLoop)
44
45
46 ✗ , _yNames (NULL)
47 ✗ , _yNominal (NULL)
48 ✗ , _y (NULL)
49 ✗ , _y0 (NULL)
50 ✗ , _y_old (NULL)
51 ✗ , _y_new (NULL)
52 ✗ , _b (NULL)
53 ✗ , _A (NULL)
54 ✗ , _ihelpArray (NULL)
55 ✗ , _jhelpArray (NULL)
56 ✗ , _zeroVec (NULL)
57
58 #if defined(klu)
59 , _kluSymbolic (NULL)
60 , _kluNumeric (NULL)
61 , _kluCommon (NULL)
62 , _Ai (NULL)
63 , _Ap (NULL)
64 , _Ax (NULL)
65 #endif
66
67 ✗ , _iterationStatus (CONTINUE)
68 ✗ , _firstCall (true)
69 ✗ , _hasDgesvFactors (false)
70 ✗ , _hasDgetc2Factors (false)
71 ✗ , _scale (NULL)
72 ✗ , _generateoutput (false)
73 ✗ , _fNominal (NULL)
74
75 {
76 ✗ _max_dimSys = 100;
77 ✗ _max_dimZeroFunc=50;
78 ✗ if (_algLoop)
79 {
80 ✗ _single_instance = false;
81 ✗ AlgLoopSolverDefaultImplementation::initialize(_algLoop->getDimZeroFunc(),_algLoop->getDimReal());
82 }
83 else
84 {
85 ✗ _single_instance = true;
86 ✗ AlgLoopSolverDefaultImplementation::initialize(_max_dimZeroFunc,_max_dimSys);
87 }
88
89 ✗ }
90
91 ✗ LinearSolver::~LinearSolver()
92 {
93 ✗ if (_yNames) delete [] _yNames;
94 ✗ if (_yNominal) delete [] _yNominal;
95 ✗ if (_y) delete [] _y;
96 ✗ if (_y0) delete [] _y0;
97 ✗ if (_y_old) delete [] _y_old;
98 ✗ if (_y_new) delete [] _y_new;
99 ✗ if (_b) delete [] _b;
100 ✗ if (_A) delete [] _A;
101 ✗ if (_ihelpArray) delete [] _ihelpArray;
102 ✗ if (_jhelpArray) delete [] _jhelpArray;
103 ✗ if (_zeroVec) delete [] _zeroVec;
104 ✗ if (_scale) delete [] _scale;
105 ✗ if (_fNominal) delete [] _fNominal;
106
107 #if defined(klu)
108 if (_sparse == true) {
109 if (_kluCommon) {
110 if (_kluSymbolic)
111 klu_free_symbolic(&_kluSymbolic, _kluCommon);
112 if (_kluNumeric)
113 klu_free_numeric(&_kluNumeric, _kluCommon);
114 delete _kluCommon;
115 }
116 if (_Ap)
117 delete [] _Ap;
118 if (_Ai)
119 delete [] _Ai;
120 }
121 #endif
122 ✗ }
123
124 ✗ void LinearSolver::initialize()
125 {
126
127 ✗ if(_firstCall)
128 ✗ _algLoop->initialize();
129
130 ✗ _firstCall = false;
131 //(Re-) Initialization of algebraic loop
132 ✗ if(!_algLoop)
133 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
134
135 ✗ _sparse = _algLoop->getUseSparseFormat();
136 ✗ _dimSys =_algLoop->getDimReal();
137 ✗ if (_dimSys>0) {
138 // Initialization of vector of unknowns
139 ✗ if (_yNames) delete [] _yNames;
140 ✗ if (_yNominal) delete [] _yNominal;
141 ✗ if (_y) delete [] _y;
142 ✗ if (_y0) delete [] _y0;
143 ✗ if (_y_old) delete [] _y_old;
144 ✗ if (_y_new) delete [] _y_new;
145 ✗ if (_b) delete [] _b;
146 ✗ if (_A) delete [] _A;
147 ✗ if (_ihelpArray) delete [] _ihelpArray;
148 ✗ if (_jhelpArray) delete [] _jhelpArray;
149 ✗ if (_zeroVec) delete [] _zeroVec;
150 ✗ if (_scale) delete [] _scale;
151 ✗ if (_fNominal) delete [] _fNominal;
152
153 ✗ _yNames = new const char* [_dimSys];
154 ✗ _yNominal = new double[_dimSys];
155 ✗ _y = new double[_dimSys];
156 ✗ _y0 = new double[_dimSys];
157 ✗ _y_old = new double[_dimSys];
158 ✗ _y_new = new double[_dimSys];
159 ✗ _b = new double[_dimSys];
160 ✗ _A = new double[_dimSys*_dimSys];
161 ✗ _ihelpArray = new long int[_dimSys];
162 ✗ _jhelpArray = new long int[_dimSys];
163 ✗ _zeroVec = new double[_dimSys];
164 ✗ _scale = new double[_dimSys];
165 ✗ _fNominal = new double[_dimSys];
166
167 ✗ _algLoop->getNamesReal(_yNames);
168 ✗ _algLoop->getNominalReal(_yNominal);
169 ✗ _algLoop->getReal(_y);
170 ✗ _algLoop->getReal(_y0);
171 ✗ _algLoop->getReal(_y_new);
172 ✗ _algLoop->getReal(_y_old);
173 ✗ memset(_b, 0, _dimSys*sizeof(double));
174 ✗ memset(_ihelpArray, 0, _dimSys*sizeof(long int));
175 ✗ memset(_jhelpArray, 0, _dimSys*sizeof(long int));
176 ✗ memset(_A, 0, _dimSys*_dimSys*sizeof(double));
177 ✗ memset(_zeroVec, 0, _dimSys*sizeof(double));
178 ✗ memset(_scale, 0, _dimSys*sizeof(double));
179
180 #if defined(klu)
181 if (_sparse) {
182 _kluCommon = new klu_common;
183 ok = klu_defaults(_kluCommon);
184 if (ok != 1)
185 throw ModelicaSimulationError(ALGLOOP_SOLVER,"error initializing Sparse Solver KLU");
186
187 sparsematrix_t& A = _algLoop->getSparseAMatrix();
188
189 _nonzeros = A.nnz();
190
191 _Ap = new int[(_dimSys + 1)];
192 _Ai = new int[_nonzeros];
193
194 int const* Ti= A.index1_data().begin();
195 int const* Tj= A.index2_data().begin();
196
197 _Ax= A.value_data().begin();
198
199 memcpy(_Ap,Ti, sizeof(int)*(_dimSys + 1));
200 memcpy(_Ai,Tj, sizeof(int)*(_nonzeros));
201
202 _kluSymbolic = klu_analyze(_dimSys, _Ap, _Ai, _kluCommon);
203 _kluNumeric = klu_factor(_Ap, _Ai, _Ax, _kluSymbolic, _kluCommon);
204 if (_kluNumeric == NULL)
205 throw ModelicaSimulationError(ALGLOOP_SOLVER, "error during numerical factorization with Sparse Solver KLU");
206 }
207 #endif
208
209 }
210
211 ✗ LOGGER_WRITE_BEGIN("LinearSolver: eq" + to_string(_algLoop->getEquationIndex()) +
212 " initialized", LC_LS, LL_DEBUG);
213 ✗ LOGGER_WRITE_VECTOR("yNames", _yNames, _dimSys, LC_LS, LL_DEBUG);
214 ✗ LOGGER_WRITE_VECTOR("yNominal", _yNominal, _dimSys, LC_LS, LL_DEBUG);
215 ✗ LOGGER_WRITE_END(LC_LS, LL_DEBUG);
216 ✗ }
217
218
219 ✗ void LinearSolver::solve(shared_ptr<ILinearAlgLoop> algLoop, bool first_solve)
220 {
221 ✗ if (first_solve)
222 {
223 _algLoop = algLoop;
224 ✗ _firstCall = true;
225 }
226 ✗ if (_algLoop != algLoop)
227 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
228 ✗ solve();
229 ✗ }
230
231 ✗ void LinearSolver::solve()
232 {
233 ✗ if (_firstCall)
234 {
235 ✗ initialize();
236 }
237 ✗ if(!_algLoop)
238 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
239 ✗ _iterationStatus = CONTINUE;
240
241 ✗ LOGGER_WRITE_BEGIN("LinearSolver: eq" + to_string(_algLoop->getEquationIndex()) +
242 " at time " + to_string(_algLoop->getSimTime()) + ":",
243 LC_LS, LL_DEBUG);
244
245 ✗ if (_algLoop->isLinearTearing())
246 ✗ _algLoop->setReal(_zeroVec); //if the system is linear tearing it means that the system is of the form Ax-b=0, so plugging in x=0 yields -b for the left hand side
247
248 ✗ _algLoop->evaluate();
249 ✗ _algLoop->getb(_b);
250
251 //if !_sparse, we use LAPACK routines, otherwise we use KLU to solve the linear system
252 ✗ if (!_sparse) {
253 //use lapack
254 ✗ long int dimRHS = 1; // Dimension of right hand side of linear system (=_b)
255 ✗ long int info = 0; // Return-flag of Fortran code
256
257 ✗ if (!_algLoop->getFreeVariablesLock()) {
258 ✗ const matrix_t& A = _algLoop->getAMatrix();
259 const double* Atemp = A.data().begin();
260
261 ✗ memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double));
262 ✗ _hasDgesvFactors = false;
263 ✗ _hasDgetc2Factors = false;
264
265 // scale Jacobian
266 ✗ std::fill(_fNominal, _fNominal + _dimSys, 1e-6);
267 ✗ for (int j = 0, idx = 0; j < _dimSys; j++) {
268 ✗ for (int i = 0; i < _dimSys; i++, idx++) {
269 ✗ _fNominal[i] = std::max(std::abs(Atemp[idx]), _fNominal[i]);
270 }
271 }
272
273 ✗ LOGGER_WRITE_VECTOR("fNominal", _fNominal, _dimSys, LC_LS, LL_DEBUG);
274
275 ✗ for (int j = 0, idx = 0; j < _dimSys; j++)
276 ✗ for (int i = 0; i < _dimSys; i++, idx++)
277 ✗ _A[idx] /= _fNominal[i];
278 }
279
280 ✗ for (int i = 0; i < _dimSys; i++)
281 ✗ _b[i] /= _fNominal[i];
282
283 ✗ if (_generateoutput) {
284 std::cout << std::endl;
285 std::cout << "We solve a linear system with coefficient matrix" << std::endl;
286 ✗ for (int i=0; i<_dimSys; i++) {
287 ✗ for (int j=0; j<_dimSys; j++) {
288 ✗ std::cout << _A[i+j*_dimSys] << " ";
289 }
290 std::cout << std::endl;
291 }
292 std::cout << "and right hand side" << std::endl;
293 ✗ for (int i=0; i<_dimSys; i++) {
294 ✗ std::cout << _b[i] << " ";
295 }
296 std::cout << std::endl;
297 }
298
299 ✗ if (!_hasDgesvFactors && !_hasDgetc2Factors) {
300 ✗ dgesv_(&_dimSys, &dimRHS, _A, &_dimSys, _ihelpArray, _b, &_dimSys, &info);
301 ✗ _hasDgesvFactors = true;
302 }
303 ✗ else if (_hasDgesvFactors) {
304 // solve using previously obtained dgesv factors
305 ✗ char trans = 'N';
306 ✗ dgetrs_(&trans, &_dimSys, &dimRHS, _A, &_dimSys, _ihelpArray, _b, &_dimSys, &info);
307 }
308 else {
309 // solve using previously obtained dgetc2 factors
310 ✗ dgesc2_(&_dimSys, _A, &_dimSys, _b, _ihelpArray, _jhelpArray, _scale);
311 ✗ info = 0;
312 }
313
314 ✗ if (info != 0) {
315 ✗ dgetc2_(&_dimSys, _A, &_dimSys, _ihelpArray, _jhelpArray, &info);
316 ✗ dgesc2_(&_dimSys, _A, &_dimSys, _b, _ihelpArray, _jhelpArray, _scale);
317 ✗ _hasDgetc2Factors = true;
318 ✗ LOGGER_WRITE("LinearSolver: Linear system singular, using perturbed system matrix.", LC_LS, LL_DEBUG);
319 ✗ _iterationStatus = DONE;
320 }
321 else
322 ✗ _iterationStatus = DONE;
323 }
324 else {
325 #if defined(klu)
326 //writing entries of A
327 sparsematrix_t& A = _algLoop->getSparseAMatrix();
328 _Ax = A.value_data().begin();
329
330 if (_generateoutput) {
331
332 std::cout << std::endl;
333
334 std::cout << "_Ap=(";
335 for (int i=0; i<_dimSys+1; i++) {
336 std::cout << " " << _Ap[i];
337 }
338 std::cout << ")" << std::endl;
339
340 std::cout << "_Ai=(";
341 for (int i=0; i<_nonzeros; i++) {
342 std::cout << " " << _Ai[i];
343 }
344 std::cout << ")" << std::endl;
345
346 std::cout << "_Ax=(";
347 for (int i=0; i<_nonzeros; i++) {
348 std::cout << " " << _Ax[i];
349 }
350 std::cout << ")" << std::endl;
351
352
353 double* a = new double[_dimSys*_dimSys];
354 memset(a, 0, _dimSys*_dimSys*sizeof(double));
355
356 for (int i=0; i<_dimSys; i++) {
357 for (int j=0; j<_dimSys; j++) {
358 for (int k=_Ap[j]; k<_Ap[j+1]; k++)
359 if (i == _Ai[k])
360 a[i+j*_dimSys] = _Ax[k];
361 }
362 }
363
364 std::cout << std::endl;
365 std::cout << "We solve a linear system with coefficient matrix" << std::endl;
366 for (int i=0; i<_dimSys; i++) {
367 for (int j=0; j<_dimSys; j++) {
368 std::cout << a[i+j*_dimSys] << " ";
369 }
370 std::cout << std::endl;
371 }
372
373 delete [] a;
374
375
376 std::cout << "and right hand side" << std::endl;
377 for (int i=0; i<_dimSys; i++) {
378 std::cout << _b[i] << " ";
379 }
380 std::cout << std::endl;
381 }
382
383 int ok = klu_refactor(_Ap, _Ai, _Ax, _kluSymbolic, _kluNumeric, _kluCommon) ;
384
385 //checking for accuracy of refactorization
386 ok = klu_rgrowth(_Ap, _Ai, _Ax, _kluSymbolic, _kluNumeric, _kluCommon);
387 if (ok != 1)
388 throw ModelicaSimulationError(ALGLOOP_SOLVER,"Sparse Solver KLU: error checking accuracy of refactorization by computing reciprocal pivot growth");
389 if (_kluCommon->rgrowth < 1e-3) {
390 klu_free_numeric(&_kluNumeric, _kluCommon);
391 _kluNumeric = klu_factor(_Ap, _Ai, _Ax, _kluSymbolic, _kluCommon);
392 if (_kluNumeric == NULL)
393 throw ModelicaSimulationError(ALGLOOP_SOLVER,"error during numerical factorization with Sparse Solver KLU");
394 }
395
396 ok = klu_solve(_kluSymbolic, _kluNumeric, _dimSys, 1, _b, _kluCommon) ;
397 if (ok != 1)
398 throw ModelicaSimulationError(ALGLOOP_SOLVER,"error solving Sparse Solver KLU");
399 _iterationStatus = DONE;
400
401 #else
402 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER,"error solving linear system with klu not implemented");
403 #endif
404 }
405
406 //we need to revert the sign of y, because the sign of b was changed before.
407 ✗ if (_algLoop->isLinearTearing()) {
408 ✗ for (int i=0; i<_dimSys; i++)
409 ✗ _y[i] = -_b[i];
410 }
411 else {
412 ✗ memcpy(_y, _b, _dimSys*sizeof(double));
413 }
414
415 ✗ if (_generateoutput) {
416 std::cout << "The solution of the linear system is given by" << std::endl;
417 ✗ for (int i=0; i<_dimSys; i++) {
418 ✗ std::cout << _y[i] << " ";
419 }
420 std::cout << std::endl;
421 }
422
423 ✗ _algLoop->setReal(_y);
424 ✗ if (_algLoop->isLinearTearing())
425 ✗ _algLoop->evaluate();//resets the right hand side to zero in the case of linear tearing. Otherwise, the b vector on the right hand side needs no update.
426
427 ✗ LOGGER_WRITE_VECTOR("y*", _y, _dimSys, LC_LS, LL_DEBUG);
428 ✗ LOGGER_WRITE_END(LC_LS, LL_DEBUG);
429 ✗ }
430
431 ✗ ILinearAlgLoopSolver::ITERATIONSTATUS LinearSolver::getIterationStatus()
432 {
433 ✗ return _iterationStatus;
434 }
435
436 ✗ bool* LinearSolver::getConditionsWorkArray()
437 {
438 ✗ return AlgLoopSolverDefaultImplementation::getConditionsWorkArray();
439
440 }
441 ✗ bool* LinearSolver::getConditions2WorkArray()
442 {
443
444 ✗ return AlgLoopSolverDefaultImplementation::getConditions2WorkArray();
445 }
446
447
448 ✗ double* LinearSolver::getVariableWorkArray()
449 {
450
451 ✗ return AlgLoopSolverDefaultImplementation::getVariableWorkArray();
452
453 }
454
455
456
457 ✗ void LinearSolver::stepCompleted(double time)
458 {
459 ✗ memcpy(_y0, _y, _dimSys*sizeof(double));
460 ✗ memcpy(_y_old, _y_new, _dimSys*sizeof(double));
461 ✗ memcpy(_y_new, _y, _dimSys*sizeof(double));
462 ✗ }
463
464 /**
465 * \brief Restores all algloop variables for a output step
466 * \return Return_Description
467 * \details Details
468 */
469 ✗ void LinearSolver::restoreOldValues()
470 {
471 ✗ memcpy(_y, _y_old, _dimSys*sizeof(double));
472 ✗ }
473
474
475 /**
476 * \brief Restores all algloop variables for last output step
477 * \return Return_Description
478 * \details Details
479 */
480 ✗ void LinearSolver::restoreNewValues()
481 {
482 ✗ memcpy(_y, _y_new, _dimSys*sizeof(double));
483 ✗ }
484
485
486 /** @} */ // end of solverLinearSolver
487