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 / 314
Functions: 0.0% 0 / 0 / 12
Branches: 0.0% 0 / 0 / 144
OMCompiler/SimulationRuntime/c/simulation/solver/newtonIteration.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 newtonIteration.c
29 */
30
31 #ifdef __cplusplus
32 extern "C" {
33 #endif
34
35 #include <math.h>
36 #include <stdlib.h>
37 #include <string.h> /* memcpy */
38
39 #include "simulation/simulation_info_json.h"
40 #include "model_help.h"
41 #include "omc_math.h"
42 #include "util/omc_error.h"
43 #include "util/varinfo.h"
44
45 #include "nonlinearSystem.h"
46 #include "newtonIteration.h"
47
48 #include "external_input.h"
49
50 /* Private function prototypes */
51
52 int solveLinearSystem(int n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData);
53 void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f,
54 double* scaledError_f, int n, double* x, double* fvec);
55 void scaling_residual_vector(DATA_NEWTON* solverData);
56 void damping_heuristic(double* x, genericResidualFunc f,
57 double current_fvec_enorm, int n, double* fvec, double* lambda, int* k,
58 DATA_NEWTON* solverData, NLS_USERDATA* userData);
59 void damping_heuristic2(double damping_parameter, double* x, genericResidualFunc f,
60 double current_fvec_enorm, int n, double* fvec, int* k,
61 DATA_NEWTON* solverData, NLS_USERDATA* userdata);
62 void LineSearch(double* x, genericResidualFunc f,
63 double current_fvec_enorm, int n, double* fvec, int* k,
64 DATA_NEWTON* solverData, NLS_USERDATA* userdata);
65 void Backtracking(double* x, genericResidualFunc f, double current_fvec_enorm,
66 int n, double* fvec, DATA_NEWTON* solverData,
67 NLS_USERDATA* userdata);
68 void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps);
69
70 /* Extern function prototypes */
71
72 extern double enorm_(int *n, double *x);
73 extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
74 extern void dgetrf_(int *m, int *n, doublereal *fjac, int *lda, int* iwork, int *info);
75 extern void dgetrs_(char *trans, int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
76
77 /**
78 * @brief Allocate NLS Newton data.
79 *
80 * @param size Size of non-linear system.
81 * @param userData Pointer to set NLS user data.
82 * @return DATA_NEWTON* Allocated memory.
83 */
84 ✗ DATA_NEWTON* allocateNewtonData(int size, NLS_USERDATA* userData)
85 {
86 ✗ DATA_NEWTON* newtonData = (DATA_NEWTON*) malloc(sizeof(DATA_NEWTON));
87 ✗ assertStreamPrint(NULL, NULL != newtonData, "allocationNewtonData() failed. Out of memory.");
88
89 ✗ newtonData->resScaling = (double*) malloc(size*sizeof(double));
90 ✗ newtonData->fvecScaled = (double*) malloc(size*sizeof(double));
91
92 ✗ newtonData->n = size;
93 ✗ newtonData->x = (double*) malloc((size+1)*sizeof(double));
94 ✗ newtonData->fvec = (double*) calloc(size,sizeof(double));
95 ✗ newtonData->xtol = 1e-6;
96 ✗ newtonData->ftol = 1e-6;
97 ✗ newtonData->maxfev = size*100;
98 ✗ newtonData->epsfcn = DBL_EPSILON;
99 ✗ newtonData->fjac = (double*) malloc((size*(size+1))*sizeof(double));
100
101 ✗ newtonData->rwork = (double*) malloc((size)*sizeof(double));
102 ✗ newtonData->iwork = (int*) malloc(size*sizeof(int));
103
104 /* damped newton */
105 ✗ newtonData->x_new = (double*) malloc((size+1)*sizeof(double));
106 ✗ newtonData->x_increment = (double*) malloc(size*sizeof(double));
107 ✗ newtonData->f_old = (double*) calloc(size,sizeof(double));
108 ✗ newtonData->fvec_minimum = (double*) calloc(size,sizeof(double));
109 ✗ newtonData->delta_f = (double*) calloc(size,sizeof(double));
110 ✗ newtonData->delta_x_vec = (double*) calloc(size,sizeof(double));
111
112 ✗ newtonData->factorization = 0;
113 ✗ newtonData->calculate_jacobian = 1;
114 ✗ newtonData->numberOfIterations = 0;
115 ✗ newtonData->numberOfFunctionEvaluations = 0;
116
117 ✗ newtonData->userData = userData;
118
119 ✗ return newtonData;
120 }
121
122 /**
123 * @brief Free NLS Newton data.
124 *
125 * @param newtonData Pointer to Newton data.
126 */
127 ✗ void freeNewtonData(DATA_NEWTON* newtonData)
128 {
129 ✗ free(newtonData->resScaling);
130 ✗ free(newtonData->fvecScaled);
131 ✗ free(newtonData->x);
132 ✗ free(newtonData->fvec);
133 ✗ free(newtonData->fjac);
134 ✗ free(newtonData->rwork);
135 ✗ free(newtonData->iwork);
136
137 /* damped newton */
138 ✗ free(newtonData->x_new);
139 ✗ free(newtonData->x_increment);
140 ✗ free(newtonData->f_old);
141 ✗ free(newtonData->fvec_minimum);
142 ✗ free(newtonData->delta_f);
143 ✗ free(newtonData->delta_x_vec);
144
145 ✗ freeNlsUserData(newtonData->userData);
146 ✗ free(newtonData);
147 ✗ }
148
149 /**
150 * @brief Solve system with Newton-Raphson.
151 *
152 * @param f Residual function.
153 * @param solverData Solver data for containing information for Newton solver.
154 * @param userData Void pointer containing user data for supplied function f and damping heuristics.
155 * @return int Returns 0.
156 */
157
158 ✗ int _omc_newton(genericResidualFunc f, DATA_NEWTON* solverData, void* userData)
159 {
160 ✗ int i, j, k = 0, l = 0, nrsh = 1;
161 ✗ int n = solverData->n; /* size of equation */
162 ✗ double *x = solverData->x;
163 ✗ double *fvec = solverData->fvec;
164 ✗ double *eps = &(solverData->ftol); /* tolerance for x */
165 double *fdeps = &(solverData->epsfcn);
166 int * maxfev = &(solverData->maxfev);
167 ✗ double *fjac = solverData->fjac;
168 double *work = solverData->rwork;
169 ✗ int *iwork = solverData->iwork;
170 int *info = &(solverData->info);
171 int calc_jac = 1;
172
173 ✗ double error_f = 1.0 + *eps, scaledError_f = 1.0 + *eps, delta_x = 1.0 + *eps, delta_f = 1.0 + *eps, delta_x_scaled = 1.0 + *eps, lambda = 1.0;
174 double current_fvec_enorm, enorm_new;
175
176 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
177 {
178 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "######### Start Newton maxfev: %d #########", (int)*maxfev);
179
180 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "x vector");
181 ✗ for(i=0; i<n; i++)
182 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "x[%d]: %e ", i, x[i]);
183 ✗ messageClose(OMC_LOG_NLS_V);
184
185 ✗ messageClose(OMC_LOG_NLS_V);
186 }
187
188 ✗ *info = 1;
189
190 /* calculate the function values */
191 ✗ (*f)(n, x, fvec, userData, 1);
192
193 ✗ solverData->nfev++;
194
195 /* save current fvec in f_old*/
196 ✗ memcpy(solverData->f_old, fvec, n*sizeof(double));
197
198 ✗ error_f = current_fvec_enorm = enorm_(&n, fvec);
199
200 ✗ memcpy(solverData->fvecScaled, solverData->fvec, n*sizeof(double));
201
202 ✗ while(error_f > *eps && scaledError_f > *eps && delta_x > *eps && delta_f > *eps && delta_x_scaled > *eps)
203 {
204 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
205 {
206 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "\n**** start Iteration: %d *****", (int) l);
207
208 /* Debug output */
209 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "function values");
210 ✗ for(i=0; i<n; i++)
211 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "fvec[%d]: %e ", i, fvec[i]);
212 ✗ messageClose(OMC_LOG_NLS_V);
213 }
214
215 /* calculate jacobian if no matrix is given */
216 ✗ if (calc_jac == 1 && solverData->calculate_jacobian >= 0)
217 {
218 ✗ (*f)(n, x, fvec, userData, 0);
219 ✗ solverData->factorization = 0;
220 ✗ calc_jac = solverData->calculate_jacobian;
221 }
222 else
223 {
224 ✗ solverData->factorization = 1;
225 ✗ calc_jac--;
226 }
227
228
229 /* debug output */
230 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_JAC))
231 {
232 ✗ char *buffer = (char*)malloc(sizeof(char)*solverData->n*15);
233
234 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", n, n);
235 ✗ for(i=0; i<solverData->n;i++)
236 {
237 char *p = buffer;
238 ✗ for(j=0; j<solverData->n; j++)
239 ✗ p += sprintf(p, "%10g ", fjac[i*n+j]);
240 ✗ infoStreamPrint(OMC_LOG_NLS_JAC, 0, "%s", buffer);
241 }
242 ✗ messageClose(OMC_LOG_NLS_JAC);
243 ✗ free(buffer);
244 }
245
246 ✗ if (solveLinearSystem(n, iwork, fvec, fjac, solverData) != 0)
247 {
248 ✗ *info=-1;
249 ✗ break;
250 }
251 else
252 {
253 ✗ for (i = 0; i < n; i++)
254 ✗ solverData->x_new[i] = x[i]-solverData->x_increment[i];
255
256 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V)) {
257 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "x_increment");
258 ✗ for(i = 0; i < n; i++) {
259 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "x_increment[%d] = %e ", i, solverData->x_increment[i]);
260 }
261 ✗ messageClose(OMC_LOG_NLS_V);
262 }
263
264 ✗ if (solverData->newtonStrategy == NEWTON_DAMPED)
265 {
266 ✗ damping_heuristic(x, f, current_fvec_enorm, n, fvec, &lambda, &k, solverData, userData);
267 }
268 ✗ else if (solverData->newtonStrategy == NEWTON_DAMPED2)
269 {
270 ✗ damping_heuristic2(0.75, x, f, current_fvec_enorm, n, fvec, &k, solverData, userData);
271 }
272 ✗ else if (solverData->newtonStrategy == NEWTON_DAMPED_LS)
273 {
274 ✗ LineSearch(x, f, current_fvec_enorm, n, fvec, &k, solverData, userData);
275 }
276 ✗ else if (solverData->newtonStrategy == NEWTON_DAMPED_BT)
277 {
278 ✗ Backtracking(x, f, current_fvec_enorm, n, fvec, solverData, userData);
279 }
280 else
281 {
282 /* calculate the function values */
283 ✗ (*f)(n, solverData->x_new, fvec, userData, 1);
284 ✗ solverData->nfev++;
285 }
286
287 ✗ calculatingErrors(solverData, &delta_x, &delta_x_scaled, &delta_f, &error_f, &scaledError_f, n, x, fvec);
288
289 /* updating x */
290 ✗ memcpy(x, solverData->x_new, n*sizeof(double));
291
292 /* updating f_old */
293 ✗ memcpy(solverData->f_old, fvec, n*sizeof(double));
294
295 ✗ current_fvec_enorm = error_f;
296
297 /* check if maximum iteration is reached */
298 ✗ if (++l > *maxfev)
299 {
300 ✗ *info = -1;
301 ✗ if (solverData->initial) {
302 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Newton iteration: Maximal number of iteration reached at initialization, but no root found.");
303 } else {
304 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Newton iteration: Maximal number of iteration reached at time %f, but no root found.", solverData->time);
305 }
306 break;
307 }
308 /* check if maximum iteration is reached */
309 ✗ if (k > 5)
310 {
311 ✗ *info = -1;
312 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Newton iteration: Maximal number of iterations reached.");
313 ✗ break;
314 }
315 }
316
317 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_NLS_V))
318 {
319 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "x vector");
320 ✗ for(i = 0; i < n; i++)
321 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "x[%d] = %e ", i, x[i]);
322 ✗ messageClose(OMC_LOG_NLS_V);
323 ✗ printErrors(delta_x, delta_x_scaled, delta_f, error_f, scaledError_f, eps);
324 }
325 }
326
327 ✗ solverData->numberOfIterations += l;
328 ✗ solverData->numberOfFunctionEvaluations += solverData->nfev;
329
330 ✗ return 0;
331 }
332
333 /**
334 * @brief Print errors.
335 *
336 * Print if tolerance is reached.
337 * Errors computed by calculatingErrors.
338 *
339 * @param delta_x delta_x := ||x_new - x_old||
340 * @param delta_x_scaled delta_x_scaled := delta_x / scaling_factor
341 * @param delta_f delta_f := || f_old - f_new ||
342 * @param error_f enorm_(n,fvec)
343 * @param scaledError_f
344 * @param eps
345 */
346 ✗ void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps)
347 {
348 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "errors ");
349 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "delta_x = %e \ndelta_x_scaled = %e \ndelta_f = %e \nerror_f = %e \nscaledError_f = %e", delta_x, delta_x_scaled, delta_f, error_f, scaledError_f);
350
351 ✗ if (delta_x < *eps)
352 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "delta_x reached eps");
353 ✗ if (delta_x_scaled < *eps)
354 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "delta_x_scaled reached eps");
355 ✗ if (delta_f < *eps)
356 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "delta_f reached eps");
357 ✗ if (error_f < *eps)
358 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "error_f reached eps");
359 ✗ if (scaledError_f < *eps)
360 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "scaledError_f reached eps");
361
362 ✗ messageClose(OMC_LOG_NLS_V);
363 ✗ }
364
365 /*! \fn solveLinearSystem
366 *
367 * function solves linear system J*(x_{n+1} - x_n) = f using lapack
368 */
369 ✗ int solveLinearSystem(int n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData)
370 {
371 ✗ int i, nrsh=1, lapackinfo;
372 ✗ char trans = 'N';
373
374 /* if no factorization is given, calculate it */
375 ✗ if (solverData->factorization == 0)
376 {
377 /* solve J*(x_{n+1} - x_n)=f */
378 ✗ dgetrf_(&n, &n, fjac, &n, iwork, &lapackinfo);
379 ✗ solverData->factorization = 1;
380 ✗ dgetrs_(&trans, &n, &nrsh, fjac, &n, iwork, fvec, &n, &lapackinfo);
381 }
382 else
383 {
384 ✗ dgetrs_(&trans, &n, &nrsh, fjac, &n, iwork, fvec, &n, &lapackinfo);
385 }
386
387 ✗ if(lapackinfo > 0)
388 {
389 ✗ warningStreamPrint(OMC_LOG_NLS, 0, "Newton iteration linear solver: Jacobian matrix singular.");
390 ✗ return -1;
391 }
392 ✗ else if(lapackinfo < 0)
393 {
394 ✗ warningStreamPrint(OMC_LOG_NLS, 0, "illegal input in argument %d", (int)lapackinfo);
395 ✗ return -1;
396 }
397 else
398 {
399 /* save solution of J*(x_{n+1} - x_n)=f */
400 ✗ memcpy(solverData->x_increment, fvec, n*sizeof(double));
401 }
402
403 ✗ return 0;
404 }
405
406 /**
407 * @brief Calculate delta and error.
408 *
409 * Current value of x from input `x`, old value from `solverData->x_new`.
410 * Current value of f(x) from input `fvec`, old value from `solverData->fvecScaled`
411 *
412 * @param solverData Newton solver data.
413 * @param delta_x delta_x := ||x_new - x_old||
414 * @param delta_x_scaled delta_x_scaled := delta_x / scaling_factor, where
415 * scaling_factor := ||x||
416 * @param delta_f delta_f := ||f_old - f_new||
417 * @param error_f error_f := ||fvec||
418 * @param scaledError_f scaledError_f := || fvec ./ resScaling||, where
419 * resScaling is from solverData.
420 * @param n Length of arrays x and fvec.
421 * @param x New vector x.
422 * @param fvec New vector f(x).
423 */
424 ✗ void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f,
425 double* scaledError_f, int n, double* x, double* fvec)
426 {
427 int i=0;
428 double scaling_factor;
429
430 /* delta_x = || x_new-x_old || */
431 ✗ for (i=0; i<n; i++)
432 ✗ solverData->delta_x_vec[i] = x[i]-solverData->x_new[i];
433
434 ✗ *delta_x = enorm_(&n,solverData->delta_x_vec);
435
436 ✗ scaling_factor = enorm_(&n,x);
437 ✗ if (scaling_factor > 1) {
438 ✗ *delta_x_scaled = *delta_x * 1./ scaling_factor;
439 } else {
440 ✗ *delta_x_scaled = *delta_x;
441 }
442
443 /* delta_f = || f_old - f_new || */
444 ✗ for (i=0; i<n; i++)
445 ✗ solverData->delta_f[i] = solverData->f_old[i]-fvec[i];
446
447 ✗ *delta_f=enorm_(&n, solverData->delta_f);
448
449 ✗ *error_f = enorm_(&n,fvec);
450
451 /* scaling residual vector */
452 ✗ scaling_residual_vector(solverData);
453
454 ✗ for (i=0; i<n; i++) {
455 ✗ solverData->fvecScaled[i]=fvec[i]/solverData->resScaling[i];
456 }
457 ✗ *scaledError_f = enorm_(&n,solverData->fvecScaled);
458 ✗ }
459
460 /**
461 * @brief Compute residual scaling vector.
462 *
463 * scalingVector[i] = 1 / ||Jac(i,:)||
464 * Warn if Jacobian row is all zeros i.e. the Jacobian is singular.
465 *
466 * @param solverData Newton solver data.
467 * @param scalingVector Residual scaling vector.
468 */
469 ✗ void compute_scaling_vector(DATA_NEWTON* solverData, double* scalingVector) {
470 int i;
471 int jac_row_start;
472
473 ✗ for(i=0; i<solverData->n; i++)
474 {
475 ✗ jac_row_start = i*solverData->n;
476 ✗ scalingVector[i] = _omc_gen_maximumVectorNorm(&(solverData->fjac[jac_row_start]), solverData->n);
477 ✗ if(scalingVector[i] <= 0.0) {
478 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Jacobian matrix is singular. Scaling of residual entry is set to 1e-16.");
479 ✗ scalingVector[i] = 1e-16;
480 }
481 ✗ else if (!isfinite(scalingVector[i]))
482 {
483 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Jacobian entry is inf or nan. Scaling of residual entry will be set to 1.0.");
484 ✗ scalingVector[i] = 1.0;
485 }
486 }
487 ✗ }
488
489 /**
490 * @brief Scale residual vector.
491 *
492 * Save result in solverData->fvecScaled.
493 *
494 * @param solverData Newton solver data.
495 */
496 ✗ void scaling_residual_vector(DATA_NEWTON* solverData)
497 {
498 int i;
499
500 ✗ compute_scaling_vector(solverData, solverData->resScaling);
501 ✗ for(i=0; i<solverData->n; i++)
502 {
503 ✗ solverData->fvecScaled[i] = solverData->fvec[i] / solverData->resScaling[i];
504 }
505 ✗ }
506
507 /*! \fn damping_heuristic
508 *
509 * first damping heuristic:
510 * x_increment will be halved until the Euclidean norm of the residual function
511 * is smaller than the Euclidean norm of the current point
512 *
513 * treshold for damping = 0.01
514 * compiler flag: -newton = damped
515 */
516 ✗ void damping_heuristic(double* x, genericResidualFunc f,
517 double current_fvec_enorm, int n, double* fvec, double* lambda, int* k,
518 DATA_NEWTON* solverData, NLS_USERDATA* userData)
519 {
520 int i;
521 double enorm_new, treshold = 1e-2;
522 modelica_boolean startDamping = FALSE; /* remember to close log message */
523
524 /* calculate new function values */
525 ✗ (*f)(n, solverData->x_new, fvec, userData, 1);
526 ✗ solverData->nfev++;
527
528 ✗ enorm_new=enorm_(&n,fvec);
529
530 ✗ if (enorm_new >= current_fvec_enorm) {
531 startDamping = TRUE;
532 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "Start Damping: enorm_new : %e; current_fvec_enorm: %e ", enorm_new, current_fvec_enorm);
533 }
534
535 ✗ while (enorm_new >= current_fvec_enorm)
536 {
537 ✗ *lambda*=0.5;
538
539
540 ✗ for (i=0; i<n; i++)
541 ✗ solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i];
542
543
544 /* calculate new function values */
545 ✗ (*f)(n, solverData->x_new, fvec, userData, 1);
546 ✗ solverData->nfev++;
547
548 ✗ enorm_new=enorm_(&n,fvec);
549
550 ✗ if (*lambda <= treshold)
551 {
552 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Warning: lambda reached a threshold.");
553
554 /* if damping is without success, trying full newton step;
555 after 5 full newton steps try a very little step */
556 ✗ if (*k >= 5)
557 ✗ for (i=0; i<n; i++)
558 ✗ solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i];
559 else
560 ✗ for (i=0; i<n; i++)
561 ✗ solverData->x_new[i]=x[i]-solverData->x_increment[i];
562
563 /* calculate new function values */
564 ✗ (*f)(n, solverData->x_new, fvec, userData, 1);
565 ✗ solverData->nfev++;
566
567 ✗ (*k)++;
568
569 ✗ break;
570 }
571 }
572
573 ✗ *lambda = 1;
574
575 ✗ if (startDamping)
576 ✗ messageClose(OMC_LOG_NLS_V);
577 ✗ }
578
579 /*! \fn damping_heuristic2
580 *
581 * second (default) damping heuristic:
582 * x_increment will be multiplied by 3/4 until the Euclidean norm of the
583 * residual function is smaller than the Euclidean norm of the current point
584 *
585 * treshold for damping = 0.0001
586 * compiler flag: -newton = damped2
587 */
588 ✗ void damping_heuristic2(double damping_parameter, double* x, genericResidualFunc f,
589 double current_fvec_enorm, int n, double* fvec, int* k,
590 DATA_NEWTON* solverData, NLS_USERDATA* userdata)
591 {
592 int i;
593 double enorm_new, treshold = 1e-4, lambda=1;
594 modelica_boolean startDamping = FALSE; /* remember to close log message */
595
596 /* calculate new function values */
597 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
598 ✗ solverData->nfev++;
599
600 ✗ enorm_new=enorm_(&n,fvec);
601
602 ✗ if (enorm_new >= current_fvec_enorm) {
603 startDamping = TRUE;
604 ✗ infoStreamPrint(OMC_LOG_NLS_V, 1, "StartDamping:");
605 }
606
607 ✗ while (enorm_new >= current_fvec_enorm)
608 {
609 ✗ lambda*=damping_parameter;
610
611 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "lambda = %e, k = %d", lambda, *k);
612
613 ✗ for (i=0; i<n; i++)
614 ✗ solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i];
615
616
617 /* calculate new function values */
618 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
619 ✗ solverData->nfev++;
620
621 ✗ enorm_new=enorm_(&n,fvec);
622
623 ✗ if (lambda <= treshold)
624 {
625 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Warning: lambda reached a threshold.");
626
627 /* if damping is without success, trying full newton step;
628 after 5 full newton steps try a very little step */
629 ✗ if (*k >= 5)
630 ✗ for (i=0; i<n; i++)
631 ✗ solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i];
632 else
633 ✗ for (i=0; i<n; i++)
634 ✗ solverData->x_new[i]=x[i]-solverData->x_increment[i];
635
636 /* calculate new function values */
637 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
638 ✗ solverData->nfev++;
639
640 ✗ (*k)++;
641
642 ✗ break;
643 }
644 }
645
646 ✗ if (startDamping)
647 ✗ messageClose(OMC_LOG_NLS_V);
648 ✗ }
649
650 /*! \fn LineSearch
651 *
652 * third damping heuristic:
653 * Along the tangent 5 five points are selected. For every point the Euclidean
654 * norm of the residual function will be calculated and the minimum is chosen
655 * for the further iteration.
656 *
657 * compiler flag: -newton = damped_ls
658 */
659 ✗ void LineSearch(double* x, genericResidualFunc f,
660 double current_fvec_enorm, int n, double* fvec, int* k,
661 DATA_NEWTON* solverData, NLS_USERDATA* userdata)
662 {
663 int i,j;
664 double enorm_new, enorm_minimum=current_fvec_enorm, lambda_minimum=0;
665 ✗ double lambda[5]={1.25,1,0.75,0.5,0.25};
666
667
668 ✗ for (j=0; j<5; j++)
669 {
670 ✗ for (i=0; i<n; i++)
671 ✗ solverData->x_new[i]=x[i]-lambda[j]*solverData->x_increment[i];
672
673 /* calculate new function values */
674 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
675 ✗ solverData->nfev++;
676
677 ✗ enorm_new=enorm_(&n,fvec);
678
679 /* searching minimal enorm */
680 ✗ if (enorm_new < enorm_minimum)
681 {
682 enorm_minimum = enorm_new;
683 ✗ lambda_minimum = lambda[j];
684 ✗ memcpy(solverData->fvec_minimum, fvec,n*sizeof(double));
685 }
686 }
687
688 ✗ infoStreamPrint(OMC_LOG_NLS_V,0,"lambda_minimum = %e", lambda_minimum);
689
690 ✗ if (lambda_minimum == 0)
691 {
692 ✗ warningStreamPrint(OMC_LOG_NLS_V, 0, "Warning: lambda_minimum = 0 ");
693
694 /* if damping is without success, trying full newton step;
695 after 5 full newton steps try a very little step */
696 ✗ if (*k >= 5)
697 {
698 lambda_minimum = 0.125;
699
700 /* calculate new function values */
701 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
702 ✗ solverData->nfev++;
703 }
704 else
705 {
706 lambda_minimum = 1;
707
708 /* calculate new function values */
709 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
710 ✗ solverData->nfev++;
711 }
712
713 ✗ (*k)++;
714 }
715 else
716 {
717 /* save new function values */
718 ✗ memcpy(fvec, solverData->fvec_minimum, n*sizeof(double));
719 }
720
721 ✗ for (i=0; i<n; i++)
722 ✗ solverData->x_new[i]=x[i]-lambda_minimum*solverData->x_increment[i];
723 ✗ }
724
725 /*! \fn Backtracking
726 *
727 * forth damping heuristic:
728 * Calculate new function h:R^n->R ; h(x) = 1/2 * ||f(x)|| ^2
729 * g(lambda) = h(x_old + lambda * x_increment)
730 * find minimum of g with golden ratio method
731 * tau = golden ratio
732 *
733 * compiler flag: -newton = damped_bt
734 */
735 ✗ void Backtracking(double* x,
736 genericResidualFunc f,
737 double current_fvec_enorm,
738 int n,
739 double* fvec,
740 DATA_NEWTON* solverData,
741 NLS_USERDATA* userdata)
742 {
743 int i,j;
744 double enorm_new, enorm_f, lambda, a1, b1, a, b, tau, g1, g2;
745 double tolerance = 1e-3;
746
747 /* saving current function values in f_old */
748 ✗ memcpy(solverData->f_old, fvec, n*sizeof(double));
749
750 ✗ for (i=0; i<n; i++)
751 ✗ solverData->x_new[i]=x[i]-solverData->x_increment[i];
752
753 /* calculate new function values */
754 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
755 ✗ solverData->nfev++;
756
757
758 /* calculate new enorm */
759 ✗ enorm_new = enorm_(&n,fvec);
760
761 /* Backtracking only if full newton step is useless */
762 ✗ if (enorm_new >= current_fvec_enorm)
763 {
764 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "Start Backtracking\n enorm_new= %f \t current_fvec_enorm=%f", enorm_new, current_fvec_enorm);
765
766 /* h(x) = 1/2 * ||f(x)|| ^2
767 * g(lambda) = h(x_old + lambda * x_increment)
768 * find minimum of g with golden ratio method
769 * tau = golden ratio
770 * */
771
772 a = 0;
773 b = 1;
774 tau = 0.618033988749894848;
775
776 a1 = a + (1-tau)*(b-a);
777 /* g1 = g(a1) = h(x_old - a1 * x_increment) = 1/2 * ||f(x_old- a1 * x_increment)||^2 */
778 ✗ solverData->x_new[i] = x[i]- a1 * solverData->x_increment[i];
779 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
780 ✗ solverData->nfev++;
781 ✗ enorm_f= enorm_(&n,fvec);
782 ✗ g1 = 0.5 * enorm_f * enorm_f;
783
784
785 b1 = a + tau * (b-a);
786 /* g2 = g(b1) = h(x_old - b1 * x_increment) = 1/2 * ||f(x_old- b1 * x_increment)||^2 */
787 ✗ solverData->x_new[i] = x[i]- b1 * solverData->x_increment[i];
788 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
789 ✗ solverData->nfev++;
790 ✗ enorm_f= enorm_(&n,fvec);
791 ✗ g2 = 0.5 * enorm_f * enorm_f;
792
793 ✗ while ( (b - a) > tolerance)
794 {
795 ✗ if (g1<g2)
796 {
797 b = b1;
798 b1 = a1;
799 ✗ a1 = a + (1-tau)*(b-a);
800 g2 = g1;
801
802 /* g1 = g(a1) = h(x_old - a1 * x_increment) = 1/2 * ||f(x_old- a1 * x_increment)||^2 */
803 ✗ solverData->x_new[i] = x[i]- a1 * solverData->x_increment[i];
804 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
805 ✗ solverData->nfev++;
806 ✗ enorm_f= enorm_(&n,fvec);
807 ✗ g1 = 0.5 * enorm_f * enorm_f;
808 }
809 else
810 {
811 a = a1;
812 a1 = b1;
813 ✗ b1 = a + tau * (b-a);
814 g1 = g2;
815
816 /* g2 = g(b1) = h(x_old - b1 * x_increment) = 1/2 * ||f(x_old- b1 * x_increment)||^2 */
817 ✗ solverData->x_new[i] = x[i]- b1 * solverData->x_increment[i];
818 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
819 ✗ solverData->nfev++;
820 ✗ enorm_f= enorm_(&n,fvec);
821 ✗ g2 = 0.5 * enorm_f * enorm_f;
822 }
823 }
824
825 ✗ lambda = (a+b)/2;
826
827 /* print lambda */
828 ✗ infoStreamPrint(OMC_LOG_NLS_V, 0, "Backtracking - lambda = %e", lambda);
829
830 ✗ for (i=0; i<n; i++)
831 ✗ solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i];
832
833 /* calculate new function values */
834 ✗ (*f)(n, solverData->x_new, fvec, userdata, 1);
835 ✗ solverData->nfev++;
836 }
837 ✗ }
838
839 #ifdef __cplusplus
840 }
841 #endif
842