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 / 522
Functions: 0.0% 0 / 0 / 8
Branches: 0.0% 0 / 0 / 309
OMCompiler/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.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 nonlinearSolverHybrd.c
29 *
30 *
31 */
32
33 #ifdef __cplusplus
34 extern "C" {
35 #endif
36
37 #include <math.h>
38 #include <stdlib.h>
39 #include <string.h> /* memcpy */
40
41 #include "../simulation_info_json.h"
42 #include "../jacobian_util.h"
43 #include "../../util/omc_error.h"
44 #include "../../util/varinfo.h"
45 #include "model_help.h"
46 #include "../../gc/omc_gc.h"
47
48 #include "nonlinearSystem.h"
49 #include "nonlinearSolverHybrd.h"
50
51 extern double enorm_(integer *n, double *x);
52
53 static void wrapper_fvec_hybrj(const integer *n_p, const double* x, double* f, double* fjac, const integer* ldjac, integer* iflag, void* userData);
54
55 /**
56 * @brief Allocate memory for non-linear hybrid solver.
57 *
58 * @param size Size of non-linear system.
59 * @param userData Information about the non-linear system (number, Jacobian, data, threadData, ...)
60 * @return DATA_HYBRD* Pointer to allocated hybrid data.
61 */
62 ✗ DATA_HYBRD* allocateHybrdData(size_t size, NLS_USERDATA* userData)
63 {
64 ✗ DATA_HYBRD* hybrdData = (DATA_HYBRD*) malloc(sizeof(DATA_HYBRD));
65 ✗ assertStreamPrint(NULL, hybrdData != NULL, "allocationHybrdData() failed!");
66
67 /* fjac/fjacobian receive evalJacobian's dense output (strided by the
68 * analytic Jacobian's own sizeCols) and getAnalyticalJacobian's memcpy uses
69 * sizeRows*sizeCols directly -- both can exceed size*(size+1) for a
70 * partial-slice Jacobian with extra addressable-but-not-genuinely-unknown
71 * seed columns (see allocateHomotopyData's identical fix and
72 * NBJacobian.mo's partialSliceSeedCandidates whole-array fallback). Size
73 * those two buffers off the larger of the two; `size`/`n`, r__ (MINPACK's
74 * own internal packed triangular factor, sized purely off the genuine
75 * unknown count), and everything else below stays genuine (the solver
76 * itself must never see phantom unknowns). */
77 ✗ size_t jacCols = size + 1;
78 ✗ if (userData != NULL && userData->analyticJacobian != NULL &&
79 ✗ (size_t)userData->analyticJacobian->sizeCols > jacCols) {
80 jacCols = (size_t)userData->analyticJacobian->sizeCols;
81 }
82
83 ✗ hybrdData->initialized = FALSE;
84 ✗ hybrdData->resScaling = (double*) malloc(size*sizeof(double));
85 ✗ hybrdData->fvecScaled = (double*) malloc(size*sizeof(double));
86 ✗ hybrdData->useXScaling = 1;
87 ✗ hybrdData->xScalefactors = (double*) malloc(size*sizeof(double));
88
89 ✗ hybrdData->n = size;
90 ✗ hybrdData->x = (double*) malloc((size+1)*sizeof(double));
91 ✗ hybrdData->xSave = (double*) malloc((size+1)*sizeof(double));
92 ✗ hybrdData->xScaled = (double*) malloc((size+1)*sizeof(double));
93 ✗ hybrdData->fvec = (double*) calloc(size, sizeof(double));
94 ✗ hybrdData->fvecSave = (double*) calloc(size, sizeof(double));
95 ✗ hybrdData->xtol = 1e-12;
96 ✗ hybrdData->maxfev = size*10000;
97 ✗ hybrdData->ml = size - 1;
98 ✗ hybrdData->mu = size - 1;
99 ✗ hybrdData->epsfcn = 1e-12;
100 ✗ hybrdData->diag = (double*) malloc(size*sizeof(double));
101 ✗ hybrdData->diagres = (double*) malloc(size*sizeof(double));
102 ✗ hybrdData->mode = 1;
103 ✗ hybrdData->factor = 100.0;
104 ✗ hybrdData->nprint = -1;
105 ✗ hybrdData->info = 0;
106 ✗ hybrdData->nfev = 0;
107 ✗ hybrdData->njev = 0;
108 ✗ hybrdData->fjac = (double*) calloc((size*jacCols), sizeof(double));
109 ✗ hybrdData->fjacobian = (double*) calloc((size*jacCols), sizeof(double));
110 ✗ hybrdData->ldfjac = size;
111 ✗ hybrdData->r__ = (double*) malloc(((size*(size+1))/2)*sizeof(double));
112 ✗ hybrdData->lr = (size*(size + 1)) / 2;
113 ✗ hybrdData->qtf = (double*) malloc(size*sizeof(double));
114 ✗ hybrdData->wa1 = (double*) malloc(size*sizeof(double));
115 ✗ hybrdData->wa2 = (double*) malloc(size*sizeof(double));
116 ✗ hybrdData->wa3 = (double*) malloc(size*sizeof(double));
117 ✗ hybrdData->wa4 = (double*) malloc(size*sizeof(double));
118
119 ✗ hybrdData->numberOfIterations = 0;
120 ✗ hybrdData->numberOfFunctionEvaluations = 0;
121
122 ✗ hybrdData->userData = userData;
123
124 ✗ return hybrdData;
125 }
126
127 /**
128 * @brief Free hybrid solver data.
129 *
130 * @param hybrdData Pointer to hybrid data.
131 */
132 ✗ void freeHybrdData(DATA_HYBRD* hybrdData)
133 {
134 ✗ free(hybrdData->resScaling);
135 ✗ free(hybrdData->fvecScaled);
136 ✗ free(hybrdData->xScalefactors);
137 ✗ free(hybrdData->x);
138 ✗ free(hybrdData->xSave);
139 ✗ free(hybrdData->xScaled);
140 ✗ free(hybrdData->fvec);
141 ✗ free(hybrdData->fvecSave);
142 ✗ free(hybrdData->diag);
143 ✗ free(hybrdData->diagres);
144 ✗ free(hybrdData->fjac);
145 ✗ free(hybrdData->fjacobian);
146 ✗ free(hybrdData->r__);
147 ✗ free(hybrdData->qtf);
148 ✗ free(hybrdData->wa1);
149 ✗ free(hybrdData->wa2);
150 ✗ free(hybrdData->wa3);
151 ✗ free(hybrdData->wa4);
152
153 ✗ freeNlsUserData(hybrdData->userData);
154
155 ✗ free(hybrdData);
156 ✗ return;
157 }
158
159 /*! \fn printVector
160 *
161 * \param [in] [vector]
162 * \param [in] [size]
163 * \param [in] [logLevel]
164 * \param [in] [name]
165 *
166 * \author wbraun
167 */
168 ✗ static void printVector(const double *vector, const integer size, const int logLevel, const char *name)
169 {
170 int i;
171 ✗ if (!OMC_ACTIVE_STREAM(logLevel)) return;
172 ✗ infoStreamPrint(logLevel, 1, "%s", name);
173 ✗ for(i=0; i<size; i++)
174 ✗ infoStreamPrint(logLevel, 0, "[%2d] %20.12g", i, vector[i]);
175 ✗ messageClose(logLevel);
176 }
177
178 /*! \fn printStatus
179 *
180 * \param [in] [solverData]
181 * \param [in] [nfunc_evals]
182 * \param [in] [xerror]
183 * \param [in] [xerror_scaled]
184 * \param [in] [logLevel]
185 *
186 * \author wbraun
187 */
188 ✗ static void printStatus(DATA *data, DATA_HYBRD *solverData, int eqSystemNumber, const int *nfunc_evals, const double *xerror, const double *xerror_scaled, const int logLevel)
189 {
190 long i;
191
192 ✗ if (!OMC_ACTIVE_STREAM(logLevel)) return;
193 ✗ infoStreamPrint(logLevel, 1, "nls status");
194
195 ✗ infoStreamPrint(logLevel, 1, "variables");
196 ✗ for(i=0; i<solverData->n; i++)
197 ✗ infoStreamPrint(logLevel, 0, "[%ld] %s = %.20e\n - scaling factor internal = %.16e\n"
198 " - scaling factor external = %.16e", i+1,
199 ✗ modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
200 ✗ solverData->x[i], solverData->diag[i], solverData->xScalefactors[i]);
201 ✗ messageClose(logLevel);
202
203 ✗ infoStreamPrint(logLevel, 1, "functions");
204 ✗ for(i=0; i<solverData->n; i++)
205 ✗ infoStreamPrint(logLevel, 0, "res[%ld] = %.20e [scaling factor = %.16e]", i+1, solverData->fvec[i], solverData->resScaling[i]);
206 ✗ messageClose(logLevel);
207
208 ✗ infoStreamPrint(logLevel, 1, "statistics");
209 ✗ infoStreamPrint(logLevel, 0, "nfunc = %d\nerror = %.20e\nerror_scaled = %.20e", *nfunc_evals, *xerror, *xerror_scaled);
210 ✗ messageClose(logLevel);
211
212 ✗ messageClose(logLevel);
213
214 }
215
216 /**
217 * @brief Calculate numeric Jacobian matrix J(x).
218 *
219 * Using finite differences method.
220 *
221 * @param hybrdUserData Pointer to hybrid solver user data.
222 * @param jac Contains values of Jacobian J(x) on exit.
223 * @param x Vector x.
224 * @param f Residual values f(x).
225 * @return int Return 0 on success.
226 */
227 ✗ static int getNumericalJacobian(NLS_USERDATA* hybrdUserData, double* jac, const double* x, double* f)
228 {
229 ✗ NONLINEAR_SYSTEM_DATA* systemData = hybrdUserData->nlsData;
230 ✗ DATA_HYBRD* solverData = (DATA_HYBRD*) systemData->solverData;
231
232 ✗ double delta_h = sqrt(solverData->epsfcn);
233 double delta_hh, delta_hhh, deltaInv;
234 ✗ integer iflag = 1;
235 int i, j, l;
236
237 ✗ memcpy(solverData->xSave, x, solverData->n*sizeof(double));
238
239 ✗ for(i = 0; i < solverData->n ; ++i)
240 {
241 ✗ delta_hhh = solverData->epsfcn * f[i];
242 ✗ delta_hh = fmax(delta_h * fmax(fabs(x[i]), fabs(delta_hhh)), delta_h);
243 ✗ delta_hh = ((f[i] >= 0) ? delta_hh : -delta_hh);
244 ✗ delta_hh = x[i] + delta_hh - x[i];
245 ✗ deltaInv = 1. / delta_hh;
246 ✗ solverData->xSave[i] = x[i] + delta_hh;
247
248 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC))
249 {
250 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "%d. %s = %f (delta_hh = %f)", i+1, modelInfoGetEquation(&hybrdUserData->data->modelData->modelDataXml, systemData->equationIndex).vars[i], solverData->xSave[i], delta_hh);
251 }
252 ✗ wrapper_fvec_hybrj(&solverData->n, (const double*) solverData->xSave, solverData->fvecSave, solverData->fjacobian, &solverData->ldfjac, &iflag, hybrdUserData);
253
254 ✗ for(j = 0; j < solverData->n; ++j)
255 {
256 ✗ l = i*solverData->n+j;
257 ✗ solverData->fjacobian[l] = jac[l] = (solverData->fvecSave[j] - f[j]) * deltaInv;
258 }
259 ✗ solverData->xSave[i] = x[i];
260 }
261
262 ✗ return 0;
263 }
264
265 /**
266 * @brief Calculate analytic Jacobian J(x).
267 *
268 * Using symbolic Jacobian and sparsity + coloring.
269 * x has to be set before calling this function.
270 *
271 * @param hybrdUserData Pointer to hybrid solver user data.
272 * @param jac Contains values of Jacobian J(x) on exit.
273 * @return int Return 0 on success.
274 */
275 ✗ static int getAnalyticalJacobian(NLS_USERDATA* hybrdUserData, double* jac)
276 {
277 ✗ DATA *data = hybrdUserData->data;
278 ✗ threadData_t *threadData = hybrdUserData->threadData;
279 ✗ NONLINEAR_SYSTEM_DATA* systemData = hybrdUserData->nlsData;
280 ✗ DATA_HYBRD* solverData = (DATA_HYBRD*)(systemData->solverData);
281 ✗ JACOBIAN* jacobian = hybrdUserData->analyticJacobian;
282
283 /* call generic dense Jacobian */
284 ✗ evalJacobian(data, threadData, jacobian, NULL, jac, TRUE);
285
286 ✗ memcpy(solverData->fjacobian, jac, (jacobian->sizeRows) * (jacobian->sizeCols) * sizeof(modelica_real));
287
288 ✗ return 0;
289 }
290
291 /**
292 * @brief Residual and Jacobian function.
293 *
294 * @param n Size of arrays x and f.
295 * @param x Vector x.
296 * @param f Residual vector f(x).
297 * Set to residual vector on exit, if iflag=1.
298 * Needs to be set as input, if iflag=2.
299 * @param fjac Array for Jacobian J(x)
300 * @param ldjac Leading dimension of Jacobian.
301 * @param iflag Flag signaling if residual or Jacobian should be evaluated.
302 * iflag = 1 ==> Residual evaluation
303 * iflag = 2 ==> Jacobian evaluation
304 * @param userDataIn User data. Get's typecasted to NLS_USERDATA
305 */
306 ✗ static void wrapper_fvec_hybrj(const integer *n_p, const double* x, double* f, double* fjac, const integer* ldjac, integer* iflag, void* userDataIn)
307 {
308 int i,j;
309 ✗ int n = *n_p;
310 NLS_USERDATA* userData = (NLS_USERDATA*) userDataIn;
311 ✗ DATA* data = userData->data;
312 ✗ threadData_t* threadData = userData->threadData;
313 ✗ NONLINEAR_SYSTEM_DATA* systemData = userData->nlsData;
314 ✗ DATA_HYBRD* hybrdData = (DATA_HYBRD*)(systemData->solverData);
315 ✗ modelica_boolean continuous = data->simulationInfo->solveContinuous;
316 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=userData->solverData};
317
318 ✗ switch(*iflag)
319 {
320 ✗ case 1:
321 /* re-scaling x vector */
322 ✗ if(hybrdData->useXScaling)
323 ✗ for(i=0; i<n; i++)
324 ✗ hybrdData->xScaled[i] = x[i]*hybrdData->xScalefactors[i];
325
326 /* debug output */
327 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_RES)) {
328 ✗ infoStreamPrint(OMC_LOG_NLS_RES, 0, "-- residual function call %d -- scaling = %d", (int)hybrdData->nfev, hybrdData->useXScaling);
329 ✗ printVector(x, n, OMC_LOG_NLS_RES, "x vector (scaled)");
330 ✗ printVector(hybrdData->xScaled, n, OMC_LOG_NLS_RES, "x vector");
331 }
332
333 /* call residual function */
334 ✗ if(hybrdData->useXScaling){
335 ✗ (systemData->residualFunc)(&resUserData, (const double*) hybrdData->xScaled, f, (const int*)iflag);
336 } else {
337 ✗ (systemData->residualFunc)(&resUserData, x, f, (const int*)iflag);
338 }
339 /* A negative iflag makes MINPACK stop. */
340 ✗ if (OMC_ERROR_RAISED()) {
341 ✗ *iflag = -1;
342 ✗ return;
343 }
344
345 /* debug output */
346 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_RES)) {
347 ✗ printVector(f, n, OMC_LOG_NLS_RES, "residuals");
348 ✗ infoStreamPrint(OMC_LOG_NLS_RES, 0, "-- end of residual function call %d --", (int)hybrdData->nfev);
349 }
350
351 ✗ hybrdData->numberOfFunctionEvaluations++;
352 ✗ break;
353 ✗ case 2:
354 /* set residual function continuous for jacobian calculation */
355 ✗ if(continuous)
356 ✗ data->simulationInfo->solveContinuous = FALSE;
357
358 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_RES))
359 ✗ infoStreamPrint(OMC_LOG_NLS_RES, 0, "-- begin calculating jacobian --");
360
361 /* performance measurement */
362 ✗ rt_ext_tp_tick(&systemData->jacobianTimeClock);
363
364 /* call apropreated jacobian function */
365 ✗ if(systemData->jacobianIndex != -1){
366 ✗ integer iflagtmp = 1;
367 ✗ wrapper_fvec_hybrj(n_p, x, f, fjac, ldjac, &iflagtmp, userData);
368
369 ✗ getAnalyticalJacobian(userData, fjac);
370 }
371 else{
372 ✗ getNumericalJacobian(userData, fjac, x, f);
373 }
374
375 /* debug output */
376 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_RES)) {
377 ✗ infoStreamPrint(OMC_LOG_NLS_RES, 0, "-- end calculating jacobian --");
378
379 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC))
380 {
381 ✗ char *buffer = (char*)malloc(sizeof(char)*n*25);
382 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", n, n);
383 ✗ for(i=0; i<n; i++)
384 {
385 char *p = buffer;
386 ✗ for(j=0; j<n; j++)
387 ✗ p += sprintf(p, "%20.12g ", fjac[i*hybrdData->n+j]);
388 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "%s", buffer);
389 }
390 ✗ messageClose(OMC_LOG_NLS_JAC);
391 ✗ free(buffer);
392 }
393 }
394 /* reset residual function again */
395 ✗ if(continuous)
396 ✗ data->simulationInfo->solveContinuous = TRUE;
397
398 /* performance measurement and statistics */
399 ✗ systemData->jacobianTime += rt_ext_tp_tock(&(systemData->jacobianTimeClock));
400 ✗ systemData->numberOfJEval++;
401
402 ✗ break;
403
404 ✗ default:
405 ✗ throwStreamPrint(NULL, "Well, this is embarrasing. The non-linear solver should never call this case.%d", (int)*iflag);
406 break;
407 }
408 }
409
410 /**
411 * @brief Solve non-linear system with hybrid method.
412 *
413 * @param data Runtime data struct.
414 * @param threadData Thread data for error handling.
415 * @param nlsData Pointer to non-linear system data.
416 * @return NLS_SOLVER_STATUS Return NLS_SOLVED on success and NLS_FAILED otherwise.
417 */
418 ✗ NLS_SOLVER_STATUS solveHybrd(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nlsData)
419 {
420 ✗ DATA_HYBRD* hybrdData = (DATA_HYBRD*)nlsData->solverData;
421 ✗ int eqSystemNumber = nlsData->equationIndex;
422
423 int i, j;
424 ✗ integer iflag = 1;
425 double xerror, xerror_scaled;
426 ✗ NLS_SOLVER_STATUS success = NLS_FAILED;
427 modelica_boolean catchedError;
428 ✗ double local_tol = 1e-12;
429 ✗ double initial_factor = hybrdData->factor;
430 ✗ int nfunc_evals = 0;
431 ✗ modelica_boolean continuous = TRUE;
432 ✗ int nonContinuousCase = 0;
433
434 ✗ int giveUp = 0;
435 ✗ int retries = 0;
436 ✗ int retries2 = 0;
437 ✗ int retries3 = 0;
438 ✗ int assertCalled = 0;
439 ✗ int assertRetries = 0;
440 ✗ int assertMessage = 0;
441
442 modelica_boolean* relationsPreBackup;
443
444 ✗ relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));
445
446 ✗ hybrdData->numberOfFunctionEvaluations = 0;
447
448 // Initialize lambda variable
449 ✗ if (nlsData->homotopySupport) {
450 ✗ hybrdData->x[hybrdData->n] = 1.0;
451 ✗ hybrdData->xSave[hybrdData->n] = 1.0;
452 ✗ hybrdData->xScaled[hybrdData->n] = 1.0;
453 }
454 else {
455 ✗ hybrdData->x[hybrdData->n] = 0.0;
456 ✗ hybrdData->xSave[hybrdData->n] = 0.0;
457 ✗ hybrdData->xScaled[hybrdData->n] = 0.0;
458 }
459
460 /* debug output */
461 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
462 {
463 ✗ int indexes[2] = {1,eqSystemNumber};
464 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_NLS_V, omc_dummyFileInfo, 1, indexes,
465 "Start solving Non-Linear System %d (size %d) at time %g with Hybrd Solver",
466 ✗ eqSystemNumber, (int) nlsData->size, data->localData[0]->timeValue);
467
468 ✗ for(i = 0; i < hybrdData->n; i++) {
469 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "%d. %s = %f", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], nlsData->nlsx[i]);
470 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " nominal = %f\nold = %f\nextrapolated = %f",
471 ✗ nlsData->nominal[i], nlsData->nlsxOld[i], nlsData->nlsxExtrapolation[i]);
472 ✗ messageClose(OMC_LOG_NLS_V);
473 }
474 ✗ messageClose(OMC_LOG_NLS_V);
475 }
476
477 /* set x vector */
478 ✗ if(data->simulationInfo->discreteCall)
479 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
480 else
481 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
482
483 ✗ for(i=0; i<hybrdData->n; i++){
484 ✗ hybrdData->xScalefactors[i] = fmax(fabs(hybrdData->x[i]), nlsData->nominal[i]);
485 }
486
487 /* start solving loop */
488 ✗ while(!giveUp && !success)
489 {
490 /* constrain x */
491 ✗ for(i=0; i<hybrdData->n; i++)
492 ✗ hybrdData->x[i] = fmax(nlsData->min[i], fmin(hybrdData->x[i], nlsData->max[i]));
493
494 ✗ for(i=0; i<hybrdData->n; i++)
495 ✗ hybrdData->xScalefactors[i] = fmax(fabs(hybrdData->x[i]), nlsData->nominal[i]);
496
497 /* debug output */
498 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
499 ✗ printVector(hybrdData->xScalefactors, (hybrdData->n), OMC_LOG_NLS_V, "scaling factors x vector");
500 ✗ printVector(hybrdData->x, (hybrdData->n), OMC_LOG_NLS_V, "Iteration variable values");
501 }
502
503 /* Scaling x vector */
504 ✗ if(hybrdData->useXScaling) {
505 ✗ for(i=0; i<hybrdData->n; i++) {
506 ✗ hybrdData->x[i] = (1.0/hybrdData->xScalefactors[i]) * hybrdData->x[i];
507 }
508 }
509
510 /* debug output */
511 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
512 {
513 ✗ printVector(hybrdData->x, hybrdData->n, OMC_LOG_NLS_V, "Iteration variable values (scaled)");
514 }
515
516 /* set residual function continuous */
517 ✗ data->simulationInfo->solveContinuous = continuous;
518
519 giveUp = 1;
520
521 /* try */
522 {
523 catchedError = TRUE;
524 #ifndef OMC_EMCC
525 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
526 #endif
527 ✗ hybrj_(wrapper_fvec_hybrj, &hybrdData->n, hybrdData->x,
528 ✗ hybrdData->fvec, hybrdData->fjac, &hybrdData->ldfjac, &hybrdData->xtol,
529 ✗ &hybrdData->maxfev, hybrdData->diag, &hybrdData->mode, &hybrdData->factor,
530 ✗ &hybrdData->nprint, &hybrdData->info, &hybrdData->nfev, &hybrdData->njev, hybrdData->r__,
531 &hybrdData->lr, hybrdData->qtf, hybrdData->wa1, hybrdData->wa2,
532 ✗ hybrdData->wa3, hybrdData->wa4, hybrdData->userData);
533
534 /* The residual raised: skip the success tail, so the retry counter
535 below keeps counting. */
536 ✗ if (OMC_ERROR_RAISED()) {
537 ✗ OMC_ERROR_CLEAR();
538 } else {
539 ✗ if(assertCalled)
540 {
541 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "After assertions failed, found a solution for which assertions did not fail.");
542 /* re-scaling x vector */
543 ✗ for(i=0; i<hybrdData->n; i++){
544 ✗ if(hybrdData->useXScaling)
545 ✗ nlsData->nlsxOld[i] = hybrdData->x[i]*hybrdData->xScalefactors[i];
546 else
547 ✗ nlsData->nlsxOld[i] = hybrdData->x[i];
548 }
549 }
550 assertRetries = 0;
551 assertCalled = 0;
552 catchedError = FALSE;
553 }
554 #ifndef OMC_EMCC
555 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
556 #endif
557 /* catch */
558 ✗ if (catchedError)
559 {
560 ✗ if (!assertMessage)
561 {
562 ✗ if (OMC_ACTIVE_WARNING_STREAM(OMC_LOG_STDOUT))
563 {
564 ✗ if(data->simulationInfo->initial)
565 ✗ warningStreamPrint(OMC_LOG_STDOUT, 1, "While solving non-linear system an assertion failed during initialization.");
566 else
567 ✗ warningStreamPrint(OMC_LOG_STDOUT, 1, "While solving non-linear system an assertion failed at time %g.", data->localData[0]->timeValue);
568 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "The non-linear solver tries to solve the problem that could take some time.");
569 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "It could help to provide better start-values for the iteration variables.");
570 ✗ if (!OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
571 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "For more information simulate with -lv LOG_NLS_V");
572 ✗ messageCloseWarning(OMC_LOG_STDOUT);
573 }
574 assertMessage = 1;
575 }
576
577 ✗ hybrdData->info = -1;
578 ✗ xerror_scaled = 1;
579 ✗ xerror = 1;
580 assertCalled = 1;
581 }
582 }
583
584 /* reset residual function continuous */
585 ✗ data->simulationInfo->solveContinuous = !continuous;
586
587 /* re-scaling x vector */
588 ✗ if(hybrdData->useXScaling)
589 ✗ for(i=0; i<hybrdData->n; i++)
590 ✗ hybrdData->x[i] = hybrdData->x[i]*hybrdData->xScalefactors[i];
591
592 /* check for proper inputs */
593 ✗ if(hybrdData->info == 0) {
594 ✗ printErrorEqSyst(IMPROPER_INPUT, modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber),
595 ✗ data->localData[0]->timeValue);
596 }
597
598 ✗ if(hybrdData->info != -1)
599 {
600 /* evaluate with discontinuities */
601 ✗ if(data->simulationInfo->discreteCall){
602 ✗ int scaling = hybrdData->useXScaling;
603 catchedError = TRUE;
604 ✗ if(scaling)
605 ✗ hybrdData->useXScaling = 0;
606
607 ✗ data->simulationInfo->solveContinuous = FALSE;
608
609 /* try */
610 #ifndef OMC_EMCC
611 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
612 #endif
613 ✗ wrapper_fvec_hybrj(&hybrdData->n, hybrdData->x, hybrdData->fvec, hybrdData->fjac, &hybrdData->ldfjac, &iflag, hybrdData->userData);
614 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { catchedError = FALSE; }
615 #ifndef OMC_EMCC
616 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
617 #endif
618 /* catch */
619 ✗ if (catchedError)
620 {
621 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Non-Linear Solver try to handle a problem with a called assert.");
622
623 ✗ hybrdData->info = -1;
624 ✗ xerror_scaled = 1;
625 ✗ xerror = 1;
626 assertCalled = 1;
627 }
628
629 ✗ if(scaling)
630 ✗ hybrdData->useXScaling = 1;
631
632 ✗ updateRelationsPre(data);
633 }
634 }
635
636 ✗ if(hybrdData->info != -1)
637 {
638 /* scaling residual vector */
639 {
640 int l=0;
641 ✗ for(i=0; i<hybrdData->n; i++){
642 ✗ hybrdData->resScaling[i] = 1e-16;
643 ✗ for(j=0; j<hybrdData->n; j++){
644 ✗ hybrdData->resScaling[i] = (fabs(hybrdData->fjacobian[l]) > hybrdData->resScaling[i])
645 ✗ ? fabs(hybrdData->fjacobian[l]) : hybrdData->resScaling[i];
646 ✗ l++;
647 }
648 ✗ hybrdData->fvecScaled[i] = hybrdData->fvec[i] * (1 / hybrdData->resScaling[i]);
649 }
650 /* debug output */
651 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
652 {
653 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "scaling factors for residual vector");
654 ✗ for(i=0; i<hybrdData->n; i++)
655 {
656 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "scaled residual [%d] : %.20e", i, hybrdData->fvecScaled[i]);
657 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "scaling factor [%d] : %.20e", i, hybrdData->resScaling[i]);
658 ✗ messageClose(OMC_LOG_NLS_V);
659 }
660 ✗ messageClose(OMC_LOG_NLS_V);
661 }
662
663 /* debug output */
664 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC))
665 {
666 ✗ char *buffer = (char*)malloc(sizeof(char)*hybrdData->n*15);
667
668 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", (int)hybrdData->n, (int)hybrdData->n);
669 ✗ for(i=0; i<hybrdData->n; i++)
670 {
671 char *p = buffer;
672 ✗ for(j=0; j<hybrdData->n; j++)
673 ✗ p += sprintf(p, "%10g ", hybrdData->fjacobian[i*hybrdData->n+j]);
674 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "%s", buffer);
675 }
676 ✗ messageClose(OMC_LOG_NLS_JAC);
677 ✗ free(buffer);
678 }
679
680 /* check for error */
681 ✗ xerror_scaled = enorm_(&hybrdData->n, hybrdData->fvecScaled);
682 ✗ xerror = enorm_(&hybrdData->n, hybrdData->fvec);
683 }
684 }
685
686 /* reset non-contunuousCase */
687 ✗ if(nonContinuousCase && xerror > local_tol && xerror_scaled > local_tol)
688 {
689 ✗ memcpy(data->simulationInfo->relationsPre, relationsPreBackup, sizeof(modelica_boolean)*data->modelData->nRelations);
690 nonContinuousCase = 0;
691 }
692
693 ✗ if(hybrdData->info < 4 && xerror > local_tol && xerror_scaled > local_tol)
694 ✗ hybrdData->info = 4;
695
696 /* solution found */
697 ✗ if(hybrdData->info == 1 || xerror <= local_tol || xerror_scaled <= local_tol)
698 {
699 int scaling;
700
701 ✗ success = NLS_SOLVED;
702 ✗ nfunc_evals += hybrdData->nfev;
703 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)){
704 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "System solved");
705 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "%d retries\n%d restarts", retries, retries2+retries3);
706 ✗ messageClose(OMC_LOG_NLS_V);
707 }
708 ✗ scaling = hybrdData->useXScaling;
709 ✗ if(scaling)
710 ✗ hybrdData->useXScaling = 0;
711
712 /* take the solution */
713 ✗ memcpy(nlsData->nlsx, hybrdData->x, hybrdData->n*(sizeof(double)));
714
715 /* try */
716 {
717 catchedError = TRUE;
718 #ifndef OMC_EMCC
719 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
720 #endif
721 ✗ wrapper_fvec_hybrj(&hybrdData->n, hybrdData->x, hybrdData->fvec, hybrdData->fjac, &hybrdData->ldfjac, &iflag, hybrdData->userData);
722 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { catchedError = FALSE; }
723 #ifndef OMC_EMCC
724 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
725 #endif
726 /* catch */
727 ✗ if (catchedError) {
728 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Non-Linear Solver try to handle a problem with a called assert.");
729
730 ✗ hybrdData->info = 4;
731 ✗ xerror_scaled = 1;
732 ✗ xerror = 1;
733 assertCalled = 1;
734 success = NLS_FAILED;
735 giveUp = 0;
736 }
737 }
738 ✗ if(scaling)
739 ✗ hybrdData->useXScaling = 1;
740 }
741 ✗ else if((hybrdData->info == 4 || hybrdData->info == 5) && assertRetries < 1+hybrdData->n && assertCalled)
742 {
743 /* case only used, when the Modelica code called an assert
744 * then, we try to modify start values to avoid the assert call.*/
745 int i;
746
747 ✗ memcpy(hybrdData->x, nlsData->nlsxOld, hybrdData->n*(sizeof(double)));
748
749 /* set all zero values to nominal values */
750 ✗ if(assertRetries < 1)
751 {
752 ✗ for(i=0; i<hybrdData->n; i++)
753 {
754 ✗ if(nlsData->nlsx[i] == 0)
755 {
756 ✗ nlsData->nlsx[i] = nlsData->nominal[i];
757 ✗ hybrdData->x[i] = nlsData->nominal[i];
758 }
759 }
760 }
761 /* change initial guess values one by one */
762 ✗ else if(assertRetries < hybrdData->n+1)
763 {
764 ✗ i = assertRetries-1;
765 ✗ hybrdData->x[i] += 0.01*nlsData->nominal[i];
766 }
767
768 ✗ giveUp = 0;
769 ✗ nfunc_evals += hybrdData->nfev;
770 ✗ assertRetries++;
771 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
772 {
773 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - try to handle a problem with a called assert vary initial value a bit. (Retry: %d)",assertRetries);
774 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
775 }
776 }
777 ✗ else if((hybrdData->info == 4 || hybrdData->info == 5) && retries < 3)
778 {
779 /* first try to decrease factor */
780
781 /* set x vector */
782 ✗ if(data->simulationInfo->discreteCall)
783 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
784 else
785 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
786
787 ✗ hybrdData->factor = hybrdData->factor / 10.0;
788
789 ✗ retries++;
790 ✗ giveUp = 0;
791 ✗ nfunc_evals += hybrdData->nfev;
792 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
793 {
794 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t decreasing initial step bound to %f.", hybrdData->factor);
795 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
796 }
797 }
798 ✗ else if((hybrdData->info == 4 || hybrdData->info == 5) && retries < 4)
799 {
800 /* try to vary the initial values */
801
802 ✗ for(i = 0; i < hybrdData->n; i++)
803 ✗ hybrdData->x[i] += nlsData->nominal[i] * 0.1;
804
805 ✗ hybrdData->factor = initial_factor;
806 ✗ retries++;
807 ✗ giveUp = 0;
808 ✗ nfunc_evals += hybrdData->nfev;
809
810 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
811 {
812 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "iteration making no progress:\t vary solution point by 1%%.");
813 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
814 }
815 }
816 ✗ else if((hybrdData->info == 4 || hybrdData->info == 5) && retries < 5)
817 {
818 /* try old values as x-Scaling factors */
819
820 /* set x vector */
821 ✗ if(data->simulationInfo->discreteCall)
822 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
823 else
824 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
825
826
827 ✗ for(i=0; i<hybrdData->n; i++)
828 ✗ hybrdData->xScalefactors[i] = fmax(fabs(nlsData->nlsxOld[i]), nlsData->nominal[i]);
829
830 ✗ retries++;
831 ✗ giveUp = 0;
832 ✗ nfunc_evals += hybrdData->nfev;
833 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
834 {
835 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "iteration making no progress:\t try old values as scaling factors.");
836 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
837 }
838 }
839 ✗ else if((hybrdData->info == 4 || hybrdData->info == 5) && retries < 6)
840 {
841 int scaling = 0;
842 /* try to disable x-Scaling */
843
844 /* set x vector */
845 ✗ if(data->simulationInfo->discreteCall)
846 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
847 else
848 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
849
850 ✗ scaling = hybrdData->useXScaling;
851 ✗ if(scaling)
852 ✗ hybrdData->useXScaling = 0;
853
854 /* reset x-scaling factors */
855 ✗ for(i=0; i<hybrdData->n; i++)
856 ✗ hybrdData->xScalefactors[i] = fmax(fabs(hybrdData->x[i]), nlsData->nominal[i]);
857
858 ✗ retries++;
859 ✗ giveUp = 0;
860 ✗ nfunc_evals += hybrdData->nfev;
861
862 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
863 {
864 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "iteration making no progress:\t try without scaling at all.");
865 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
866 }
867 }
868 ✗ else if((hybrdData->info == 4 || hybrdData->info == 5) && retries < 7 && data->simulationInfo->discreteCall)
869 {
870 /* try to solve non-continuous
871 * work-a-round: since other wise some model does
872 * stuck in event iteration. e.g.: Modelica.Mechanics.Rotational.Examples.HeatLosses
873 */
874
875 ✗ memcpy(hybrdData->x, nlsData->nlsxOld, hybrdData->n*(sizeof(double)));
876 ✗ retries++;
877
878 /* try to solve a discontinuous system */
879 ✗ continuous = FALSE;
880
881 ✗ nonContinuousCase = 1;
882 ✗ memcpy(relationsPreBackup, data->simulationInfo->relationsPre, sizeof(modelica_boolean)*data->modelData->nRelations);
883
884 ✗ giveUp = 0;
885 ✗ nfunc_evals += hybrdData->nfev;
886 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
887 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t try to solve a discontinuous system.");
888 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
889 }
890 /* Then try with old values (instead of extrapolating )*/
891 ✗ } else if((hybrdData->info == 4 || hybrdData->info == 5) && retries2 < 1) {
892 int scaling = 0;
893 /* set x vector */
894 ✗ memcpy(hybrdData->x, nlsData->nlsxOld, hybrdData->n*(sizeof(double)));
895
896 ✗ scaling = hybrdData->useXScaling;
897 ✗ if(!scaling)
898 ✗ hybrdData->useXScaling = 1;
899
900 ✗ continuous = TRUE;
901 ✗ hybrdData->factor = initial_factor;
902
903 ✗ retries = 0;
904 ✗ retries2++;
905 ✗ giveUp = 0;
906 ✗ nfunc_evals += hybrdData->nfev;
907 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
908 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t use old values instead extrapolated.");
909 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
910 }
911 /* try to vary the initial values */
912 ✗ } else if((hybrdData->info == 4 || hybrdData->info == 5) && retries2 < 2) {
913 /* set x vector */
914 ✗ if(data->simulationInfo->discreteCall)
915 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
916 else
917 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
918 ✗ for(i = 0; i < hybrdData->n; i++) {
919 ✗ hybrdData->x[i] *= 1.01;
920 };
921
922 ✗ retries = 0;
923 ✗ retries2++;
924 ✗ giveUp = 0;
925 ✗ nfunc_evals += hybrdData->nfev;
926 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
927 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0,
928 " - iteration making no progress:\t vary initial point by adding 1%%.");
929 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
930 }
931 /* try to vary the initial values */
932 ✗ } else if((hybrdData->info == 4 || hybrdData->info == 5) && retries2 < 3) {
933 /* set x vector */
934 ✗ if(data->simulationInfo->discreteCall)
935 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
936 else
937 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
938 ✗ for(i = 0; i < hybrdData->n; i++) {
939 ✗ hybrdData->x[i] *= 0.99;
940 };
941
942 ✗ retries = 0;
943 ✗ retries2++;
944 ✗ giveUp = 0;
945 ✗ nfunc_evals += hybrdData->nfev;
946 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
947 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t vary initial point by -1%%.");
948 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
949 }
950 /* try to vary the initial values */
951 ✗ } else if((hybrdData->info == 4 || hybrdData->info == 5) && retries2 < 4) {
952 /* set x vector */
953 ✗ memcpy(hybrdData->x, nlsData->nominal, hybrdData->n*(sizeof(double)));
954 ✗ retries = 0;
955 ✗ retries2++;
956 ✗ giveUp = 0;
957 ✗ nfunc_evals += hybrdData->nfev;
958 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
959 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t try scaling factor as initial point.");
960 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
961 }
962 /* try own scaling factors */
963 ✗ } else if((hybrdData->info == 4 || hybrdData->info == 5) && retries2 < 5 && !assertCalled) {
964 /* set x vector */
965 ✗ if(data->simulationInfo->discreteCall)
966 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
967 else
968 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
969
970 ✗ for(i = 0; i < hybrdData->n; i++) {
971 ✗ hybrdData->diag[i] = fabs(hybrdData->resScaling[i]);
972 ✗ if(hybrdData->diag[i] <= 1e-16)
973 ✗ hybrdData->diag[i] = 1e-16;
974 }
975 ✗ retries = 0;
976 ✗ retries2++;
977 ✗ giveUp = 0;
978 ✗ hybrdData->mode = 2;
979 ✗ nfunc_evals += hybrdData->nfev;
980 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
981 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t try with own scaling factors.");
982 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
983 }
984 /* try without internal scaling */
985 ✗ } else if((hybrdData->info == 4 || hybrdData->info == 5) && retries3 < 1) {
986 /* set x vector */
987 ✗ if(data->simulationInfo->discreteCall)
988 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
989 else
990 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
991
992 ✗ for(i = 0; i < hybrdData->n; i++)
993 ✗ hybrdData->diag[i] = 1.0;
994
995 ✗ hybrdData->useXScaling = 1;
996 ✗ retries = 0;
997 ✗ retries2 = 0;
998 ✗ retries3++;
999 ✗ hybrdData->mode = 2;
1000 ✗ giveUp = 0;
1001 ✗ nfunc_evals += hybrdData->nfev;
1002 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
1003 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t disable solver internal scaling.");
1004 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
1005 }
1006 /* try to reduce the tolerance a bit */
1007 ✗ } else if((hybrdData->info == 4 || hybrdData->info == 5) && retries3 < 6) {
1008 /* set x vector */
1009 ✗ if(data->simulationInfo->discreteCall)
1010 ✗ memcpy(hybrdData->x, nlsData->nlsx, hybrdData->n*(sizeof(double)));
1011 else
1012 ✗ memcpy(hybrdData->x, nlsData->nlsxExtrapolation, hybrdData->n*(sizeof(double)));
1013
1014 /* reduce tolarance */
1015 ✗ local_tol = local_tol*10;
1016
1017 ✗ hybrdData->factor = initial_factor;
1018 ✗ hybrdData->mode = 1;
1019
1020 ✗ retries = 0;
1021 ✗ retries2 = 0;
1022 ✗ retries3++;
1023
1024 ✗ giveUp = 0;
1025 ✗ nfunc_evals += hybrdData->nfev;
1026 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
1027 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, " - iteration making no progress:\t reduce the tolerance slightly to %e.", local_tol);
1028 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
1029 }
1030 ✗ } else if(hybrdData->info >= 2 && hybrdData->info <= 5) {
1031
1032 /* while the initialization it's ok to every time a solution */
1033 ✗ if(!data->simulationInfo->initial){
1034 ✗ printErrorEqSyst(ERROR_AT_TIME, modelInfoGetEquation(&data->modelData->modelDataXml, eqSystemNumber), data->localData[0]->timeValue);
1035 }
1036 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
1037 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "### No Solution! ###\n after %d restarts", retries*retries2*retries3);
1038 ✗ printStatus(data, hybrdData, eqSystemNumber, &nfunc_evals, &xerror, &xerror_scaled, OMC_LOG_NLS_V);
1039 }
1040 /* take the best approximation */
1041 ✗ memcpy(nlsData->nlsx, hybrdData->x, hybrdData->n*(sizeof(double)));
1042
1043 giveUp = 1;
1044 success = NLS_FAILED;
1045 ✗ break;
1046 }
1047 }
1048
1049 /* reset some solving data */
1050 ✗ hybrdData->factor = initial_factor;
1051 ✗ hybrdData->mode = 1;
1052
1053 /* write statistics */
1054 ✗ nlsData->numberOfFEval += hybrdData->numberOfFunctionEvaluations;
1055 /* iteration in hybrid are equal to the nfev numbers */
1056 ✗ nlsData->numberOfIterations += nfunc_evals;
1057
1058 ✗ free(relationsPreBackup);
1059
1060 ✗ return success;
1061 }
1062
1063 #ifdef __cplusplus
1064 }
1065 #endif
1066