Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 75.6% 130 / 0 / 172
Functions: 58.3% 7 / 0 / 12
Branches: 26.9% 80 / 0 / 297
OMCompiler/SimulationRuntime/cpp/Solver/Dgesv/DgesvSolver.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 solverDgesvSolver
31 *
32 * @{
33 */
34
35 #include <Core/Math/ILapack.h>
36 #include <Solver/Dgesv/FactoryExport.h>
37 #include <Core/Utils/extension/logger.hpp>
38 #include <Solver/Dgesv/DgesvSolver.h>
39
40 #include <iostream>
41
42
43 69 DgesvSolver::DgesvSolver(ILinSolverSettings* settings,shared_ptr<ILinearAlgLoop> algLoop)
44 :AlgLoopSolverDefaultImplementation()
45 , _algLoop (algLoop)
46
47 69 , _yNames (NULL)
48 69 , _yNominal (NULL)
49 69 , _y (NULL)
50 69 , _y0 (NULL)
51 69 , _y_old (NULL)
52 69 , _y_new (NULL)
53 69 , _b (NULL)
54 69 , _A (NULL)
55 69 , _iHelp (NULL)
56 69 , _jHelp (NULL)
57 69 , _zeroVec (NULL)
58 69 , _iterationStatus (CONTINUE)
59 69 , _firstCall (true)
60 69 , _hasDgesvFactors (false)
61 69 , _hasDgetc2Factors (false)
62
1/2
✓ Branch 2 taken 69 times.
✗ Branch 3 not taken.
 69 , _fNominal (NULL)
63 {
64
1/2
✓ Branch 0 taken 69 times.
✗ Branch 1 not taken.
 69 if (_algLoop)
65 {
66
3/3
✓ Branch 1 taken 69 times.
✓ Branch 4 taken 69 times.
✓ Branch 7 taken 69 times.
 69 AlgLoopSolverDefaultImplementation::initialize(_algLoop->getDimZeroFunc(),_algLoop->getDimReal());
67 }
68 else
69 {
70 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "solve for single instance is not supported");
71 }
72 69 }
73
74 138 DgesvSolver::~DgesvSolver()
75 {
76
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_yNames) delete [] _yNames;
77
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_yNominal) delete [] _yNominal;
78
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_y) delete [] _y;
79
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_y0) delete [] _y0;
80
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_y_old) delete [] _y_old;
81
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_y_new) delete [] _y_new;
82
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_b) delete [] _b;
83
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_A) delete [] _A;
84
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_iHelp) delete [] _iHelp;
85
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_jHelp) delete [] _jHelp;
86
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_zeroVec) delete [] _zeroVec;
87
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 1 time.
 69 if (_fNominal) delete [] _fNominal;
88 138 }
89
90 68 void DgesvSolver::initialize()
91 {
92
1/2
✓ Branch 0 taken 68 times.
✗ Branch 1 not taken.
 68 _firstCall = false;
93 //(Re-) Initialization of algebraic loop
94
1/2
✓ Branch 0 taken 68 times.
✗ Branch 1 not taken.
 68 if(_algLoop)
95 68 _algLoop->initialize();
96 else
97 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
98
99 68 int _dimSys = _algLoop->getDimReal();
100
1/2
✓ Branch 0 taken 68 times.
✗ Branch 1 not taken.
 68 if (_dimSys > 0) {
101 // Initialization of vector of unknowns
102
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_yNames) delete [] _yNames;
103
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_yNominal) delete [] _yNominal;
104
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_y) delete [] _y;
105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_y0) delete [] _y0;
106
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_y_old) delete [] _y_old;
107
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_y_new) delete [] _y_new;
108
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_b) delete [] _b;
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_A) delete [] _A;
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_iHelp) delete [] _iHelp;
111
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_jHelp) delete [] _jHelp;
112
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_zeroVec) delete [] _zeroVec;
113
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 68 times.
 68 if (_fNominal) delete [] _fNominal;
114
115 68 _yNames = new const char* [_dimSys];
116 68 _yNominal = new double[_dimSys];
117 68 _y = new double[_dimSys];
118 68 _y0 = new double[_dimSys];
119 68 _y_old = new double[_dimSys];
120 68 _y_new = new double[_dimSys];
121 68 _b = new double[_dimSys];
122 68 _A = new double[_dimSys*_dimSys];
123 68 _iHelp = new long int[_dimSys];
124 68 _jHelp = new long int[_dimSys];
125 68 _zeroVec = new double[_dimSys];
126 68 _fNominal = new double[_dimSys];
127
128 68 _algLoop->getNamesReal(_yNames);
129 68 _algLoop->getNominalReal(_yNominal);
130 68 _algLoop->getReal(_y);
131 68 _algLoop->getReal(_y0);
132 68 _algLoop->getReal(_y_new);
133 68 _algLoop->getReal(_y_old);
134 68 memset(_b, 0, _dimSys*sizeof(double));
135 68 memset(_iHelp, 0, _dimSys*sizeof(long int));
136 68 memset(_jHelp, 0, _dimSys*sizeof(long int));
137 68 memset(_A, 0, _dimSys*_dimSys*sizeof(double));
138 68 memset(_zeroVec, 0, _dimSys*sizeof(double));
139 }
140 else {
141 ✗ _iterationStatus = SOLVERERROR;
142 }
143
144
145 ✗ LOGGER_WRITE_BEGIN("DgesvSolver: eq" + to_string(_algLoop->getEquationIndex()) +
146 " initialized", LC_LS, LL_DEBUG);
147 ✗ LOGGER_WRITE_VECTOR("yNames", _yNames, _dimSys, LC_LS, LL_DEBUG);
148 ✗ LOGGER_WRITE_VECTOR("yNominal", _yNominal, _dimSys, LC_LS, LL_DEBUG);
149 ✗ LOGGER_WRITE_END(LC_LS, LL_DEBUG);
150 68 }
151
152 ✗ void DgesvSolver::solve(shared_ptr<ILinearAlgLoop> algLoop,bool first_solve)
153 {
154 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "solve for single instance is not supported");
155 }
156 24014046 bool* DgesvSolver::getConditionsWorkArray()
157 {
158 24014046 return AlgLoopSolverDefaultImplementation::getConditionsWorkArray();
159
160 }
161 24014046 bool* DgesvSolver::getConditions2WorkArray()
162 {
163
164 24014046 return AlgLoopSolverDefaultImplementation::getConditions2WorkArray();
165 }
166
167
168 24014046 double* DgesvSolver::getVariableWorkArray()
169 {
170
171 24014046 return AlgLoopSolverDefaultImplementation::getVariableWorkArray();
172
173 }
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24015868 times.
 24015868 void DgesvSolver::solve()
175 {
176
177
178
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24015868 times.
 24015868 if(!_algLoop)
179 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "algloop system is not initialized");
180
181
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 24015800 times.
 24015868 if (_firstCall)
182 {
183 68 initialize();
184 }
185
1/2
✓ Branch 0 taken 24015868 times.
✗ Branch 1 not taken.
 24015868 _iterationStatus = CONTINUE;
186
187
188 ✗ LOGGER_WRITE_BEGIN("DgesvSolver: eq" + to_string(_algLoop->getEquationIndex()) +
189 " at time " + to_string(_algLoop->getSimTime()) + ":",
190 LC_LS, LL_DEBUG);
191
192 //use lapack
193 24015868 long int dimRHS = 1; // Dimension of right hand side of linear system (=_b)
194 24015868 long int info = 0; // Return flag of dgesv
195 24015868 double scale = 0.0; // Scale factor of dgesc2
196
197
2/2
✓ Branch 1 taken 24014046 times.
✓ Branch 2 taken 1822 times.
 24015868 if (_algLoop->isLinearTearing())
198 24014046 _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
199
200 24015868 _algLoop->evaluate();
201 24015868 _algLoop->getb(_b);
202
203
1/2
✓ Branch 1 taken 24015868 times.
✗ Branch 2 not taken.
 24015868 if (!_algLoop->getFreeVariablesLock()) {
204 24015868 const matrix_t& A = _algLoop->getAMatrix();
205 const double* Atemp = A.data().begin();
206
207 24015868 memcpy(_A, Atemp, _dimSys*_dimSys*sizeof(double));
208 24015868 _hasDgesvFactors = false;
209 24015868 _hasDgetc2Factors = false;
210
211 24015868 std::fill(_fNominal, _fNominal + _dimSys, 1e-6);
212
2/2
✓ Branch 0 taken 89702262 times.
✓ Branch 1 taken 24015868 times.
 113718130 for (int j = 0, idx = 0; j < _dimSys; j++)
213
2/2
✓ Branch 0 taken 352455108 times.
✓ Branch 1 taken 89702262 times.
 442157370 for (int i = 0; i < _dimSys; i++, idx++)
214
2/2
✓ Branch 0 taken 194109714 times.
✓ Branch 1 taken 158345394 times.
 546564822 _fNominal[i] = std::max(std::abs(Atemp[idx]), _fNominal[i]);
215
216 ✗ LOGGER_WRITE_VECTOR("fNominal", _fNominal, _dimSys, LC_LS, LL_DEBUG);
217
218
2/2
✓ Branch 0 taken 89702262 times.
✓ Branch 1 taken 24015868 times.
 113718130 for (int j = 0, idx = 0; j < _dimSys; j++)
219
2/2
✓ Branch 0 taken 352455108 times.
✓ Branch 1 taken 89702262 times.
 442157370 for (int i = 0; i < _dimSys; i++, idx++)
220 352455108 _A[idx] /= _fNominal[i];
221 }
222
223
2/2
✓ Branch 0 taken 89702262 times.
✓ Branch 1 taken 24015868 times.
 113718130 for (int i = 0; i < _dimSys; i++)
224 89702262 _b[i] /= _fNominal[i];
225
226 double det;
227
3/4
✓ Branch 0 taken 2121007 times.
✓ Branch 1 taken 21894861 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2121007 times.
 24015868 if (_dimSys == 1 && (det = _A[0]) != 0.0) {
228 2121007 _b[0] /= det;
229 2121007 info = 0;
230 }
231
3/4
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 21894857 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
 21894861 else if (_dimSys == 2 && (det = _A[0]*_A[3] - _A[1]*_A[2]) != 0.0) {
232 4 double b0 = (_b[0]*_A[3] - _b[1]*_A[2]) / det;
233 4 _b[1] = (_A[0]*_b[1] - _A[1]*_b[0]) / det;
234 4 _b[0] = b0;
235 4 info = 0;
236 }
237
2/4
✓ Branch 0 taken 21894857 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 21894857 times.
✗ Branch 3 not taken.
 21894857 else if (!_hasDgesvFactors && !_hasDgetc2Factors) {
238 21894857 dgesv_(&_dimSys, &dimRHS, _A, &_dimSys, _iHelp, _b, &_dimSys, &info);
239 21894857 _hasDgesvFactors = true;
240 }
241 ✗ else if (_hasDgesvFactors) {
242 // solve using previously obtained dgesv factors
243 ✗ char trans = 'N';
244 ✗ dgetrs_(&trans, &_dimSys, &dimRHS, _A, &_dimSys, _iHelp, _b, &_dimSys, &info);
245 }
246 else {
247 // solve using previously obtained dgetc2 factors
248 ✗ dgesc2_(&_dimSys, _A, &_dimSys, _b, _iHelp, _jHelp, &scale);
249 ✗ info = 0;
250 }
251
252
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24015868 times.
 24015868 if (info > 0) {
253 ✗ long int info2 = 0;
254 ✗ dgetc2_(&_dimSys, _A, &_dimSys, _iHelp, _jHelp, &info2);
255 ✗ dgesc2_(&_dimSys, _A, &_dimSys, _b, _iHelp, _jHelp, &scale);
256 ✗ _hasDgetc2Factors = true;
257 ✗ LOGGER_WRITE("total pivoting: dgesv/dgetc2 infos: " + to_string(info) + "/" + to_string(info2) +
258 ", dgesc2 scale: " + to_string(scale) + ")", LC_LS, LL_DEBUG);
259 }
260
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24015868 times.
 24015868 else if (info < 0) {
261 ✗ _iterationStatus = SOLVERERROR;
262 ✗ LOGGER_WRITE_END(LC_LS, LL_DEBUG);
263 ✗ if (_algLoop->isLinearTearing())
264 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "error solving linear tearing system (dgesv info: " + to_string(info) + ")");
265 else
266 ✗ throw ModelicaSimulationError(ALGLOOP_SOLVER, "error solving linear system (dgesv info: " + to_string(info) + ")");
267 }
268 24015868 _iterationStatus = DONE;
269
270 //we need to revert the sign of y, because the sign of b was changed before.
271
2/2
✓ Branch 1 taken 24014046 times.
✓ Branch 2 taken 1822 times.
 24015868 if (_algLoop->isLinearTearing()) {
272
2/2
✓ Branch 0 taken 89694978 times.
✓ Branch 1 taken 24014046 times.
 113709024 for (int i = 0; i < _dimSys; i++)
273 89694978 _y[i] = -_b[i];
274 }
275 else {
276 1822 memcpy(_y, _b, _dimSys*sizeof(double));
277 }
278
279 24015868 _algLoop->setReal(_y);
280
2/2
✓ Branch 1 taken 24014046 times.
✓ Branch 2 taken 1822 times.
 24015868 if (_algLoop->isLinearTearing())
281 24014046 _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.
282
283 ✗ LOGGER_WRITE_VECTOR("y*", _y, _dimSys, LC_LS, LL_DEBUG);
284 ✗ LOGGER_WRITE_END(LC_LS, LL_DEBUG);
285 24015868 }
286
287 ✗ ILinearAlgLoopSolver::ITERATIONSTATUS DgesvSolver::getIterationStatus()
288 {
289 ✗ return _iterationStatus;
290 }
291
292 ✗ void DgesvSolver::stepCompleted(double time)
293 {
294 ✗ memcpy(_y0, _y, _dimSys*sizeof(double));
295 ✗ memcpy(_y_old, _y_new, _dimSys*sizeof(double));
296 ✗ memcpy(_y_new, _y, _dimSys*sizeof(double));
297 ✗ }
298
299 /**
300 * \brief Restores all algloop variables for a output step
301 * \return Return_Description
302 * \details Details
303 */
304 ✗ void DgesvSolver::restoreOldValues()
305 {
306 ✗ memcpy(_y, _y_old, _dimSys*sizeof(double));
307 ✗ }
308
309
310 /**
311 * \brief Restores all algloop variables for last output step
312 * \return Return_Description
313 * \details Details
314 */
315 ✗ void DgesvSolver::restoreNewValues()
316 {
317 ✗ memcpy(_y, _y_new, _dimSys*sizeof(double));
318 ✗ }
319
320
321 /** @} */ // end of solverDgesvSolver
322