Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/newtonIteration.c
Warning:line 170, column 9
Value stored to 'data' during its initialization is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2014, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
11 * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
14 * ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the OSMC (Open Source Modelica Consortium)
17 * Public License (OSMC-PL) are obtained from OSMC, either from the above
18 * address, from the URLs: http://www.openmodelica.org or
19 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica
20 * distribution. GNU version 3 is obtained from:
21 * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from:
22 * http://www.opensource.org/licenses/BSD-3-Clause.
23 *
24 * This program is distributed WITHOUT ANY WARRANTY; without even the implied
25 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS
26 * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE
27 * CONDITIONS OF OSMC-PL.
28 *
29 */
30
31/*! \file newtonIteration.c
32 */
33
34#ifdef __cplusplus
35extern "C" {
36#endif
37
38#include <math.h>
39#include <stdlib.h>
40#include <string.h> /* memcpy */
41
42#include "simulation/simulation_info_json.h"
43#include "util/omc_error.h"
44#include "util/varinfo.h"
45#include "model_help.h"
46
47#include "nonlinearSystem.h"
48#include "newtonIteration.h"
49
50#include "external_input.h"
51
52
53extern double enorm_(int *n, double *x);
54int solveLinearSystem(int* n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData);
55void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f,
56 double* scaledError_f, int* n, double* x, double* fvec);
57void scaling_residual_vector(DATA_NEWTON* solverData);
58void damping_heuristic(double* x, int(*f)(int*, double*, double*, void*, int),
59 double current_fvec_enorm, int* n, double* fvec, double* lambda, int* k, DATA_NEWTON* solverData, void* userdata);
60void damping_heuristic2(double damping_parameter, double* x, int(*f)(int*, double*, double*, void*, int),
61 double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata);
62void LineSearch(double* x, int(*f)(int*, double*, double*, void*, int),
63 double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata);
64void Backtracking(double* x, int(*f)(int*, double*, double*, void*, int),
65 double current_fvec_enorm, int* n, double* fvec, DATA_NEWTON* solverData, void* userdata);
66void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps);
67
68
69#ifdef __cplusplus
70extern "C" {
71#endif
72
73extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
74extern void dgetrf_(int *m, int *n, doublereal *fjac, int *lda, int* iwork, int *info);
75extern void dgetrs_(char *trans, int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
76
77#ifdef __cplusplus
78}
79#endif
80
81
82/*! \fn allocateNewtonData
83 * allocate memory for nonlinear system solver
84 */
85int allocateNewtonData(int size, void** voiddata)
86{
87 DATA_NEWTON* data = (DATA_NEWTON*) malloc(sizeof(DATA_NEWTON));
88
89 *voiddata = (void*)data;
90 assertStreamPrint(NULL, NULL != data, "allocationNewtonData() failed!")if (!(((void*)0) != data)) {throwStreamPrint((((void*)0)), "allocationNewtonData() failed!"
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "./simulation/solver/newtonIteration.c",
90, __extension__ __PRETTY_FUNCTION__); }));}
;
91
92 data->resScaling = (double*) malloc(size*sizeof(double));
93 data->fvecScaled = (double*) malloc(size*sizeof(double));
94
95 data->n = size;
96 data->x = (double*) malloc((size+1)*sizeof(double));
97 data->fvec = (double*) calloc(size,sizeof(double));
98 data->xtol = 1e-6;
99 data->ftol = 1e-6;
100 data->maxfev = size*100;
101 data->epsfcn = DBL_EPSILON2.2204460492503131e-16;
102 data->fjac = (double*) malloc((size*(size+1))*sizeof(double));
103
104 data->rwork = (double*) malloc((size)*sizeof(double));
105 data->iwork = (int*) malloc(size*sizeof(int));
106
107 /* damped newton */
108 data->x_new = (double*) malloc((size+1)*sizeof(double));
109 data->x_increment = (double*) malloc(size*sizeof(double));
110 data->f_old = (double*) calloc(size,sizeof(double));
111 data->fvec_minimum = (double*) calloc(size,sizeof(double));
112 data->delta_f = (double*) calloc(size,sizeof(double));
113 data->delta_x_vec = (double*) calloc(size,sizeof(double));
114
115 data->factorization = 0;
116 data->calculate_jacobian = 1;
117 data->numberOfIterations = 0;
118 data->numberOfFunctionEvaluations = 0;
119
120 return 0;
121}
122
123/*! \fn freeNewtonData
124 *
125 * free memory for nonlinear solver newton
126 *
127 */
128int freeNewtonData(void **voiddata)
129{
130 DATA_NEWTON* data = (DATA_NEWTON*) *voiddata;
131
132 free(data->resScaling);
133 free(data->fvecScaled);
134 free(data->x);
135 free(data->fvec);
136 free(data->fjac);
137 free(data->rwork);
138 free(data->iwork);
139
140 /* damped newton */
141 free(data->x_new);
142 free(data->x_increment);
143 free(data->f_old);
144 free(data->fvec_minimum);
145 free(data->delta_f);
146 free(data->delta_x_vec);
147
148 return 0;
149}
150
151/*! \fn solve system with Newton-Raphson
152 *
153 * \param [in] [n] size of equation
154 * [eps] tolerance for x
155 * [h] tolerance for f'
156 * [k] maximum number of iterations
157 * [work] work array size of (n*X)
158 * [f] user provided function
159 * [data] userdata
160 * [info]
161 * [calculate_jacobian] flag which decides whether Jacobian is calculated
162 * (0) once for the first calculation
163 * (i) every i steps (=1 means original newton method)
164 * (-1) never, factorization has to be given in A
165 *
166 */
167int _omc_newton(int(*f)(int*, double*, double*, void*, int), DATA_NEWTON* solverData, void* userdata)
168{
169 DATA_USER* uData = (DATA_USER*) userdata;
170 DATA* data = (DATA*) uData->data;
Value stored to 'data' during its initialization is never read
171 int i, j, k = 0, l = 0, nrsh = 1;
172 int *n = &(solverData->n);
173 double *x = solverData->x;
174 double *fvec = solverData->fvec;
175 double *eps = &(solverData->ftol);
176 double *fdeps = &(solverData->epsfcn);
177 int * maxfev = &(solverData->maxfev);
178 double *fjac = solverData->fjac;
179 double *work = solverData->rwork;
180 int *iwork = solverData->iwork;
181 int *info = &(solverData->info);
182 int calc_jac = 1;
183
184 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;
185 double current_fvec_enorm, enorm_new;
186
187 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
188 {
189 infoStreamPrint(LOG_NLS_V, 1, "######### Start Newton maxfev: %d #########", (int)*maxfev);
190
191 infoStreamPrint(LOG_NLS_V, 1, "x vector");
192 for(i=0; i<*n; i++)
193 infoStreamPrint(LOG_NLS_V, 0, "x[%d]: %e ", i, x[i]);
194 messageClose(LOG_NLS_V);
195
196 messageClose(LOG_NLS_V);
197 }
198
199 *info = 1;
200
201 /* calculate the function values */
202 (*f)(n, x, fvec, userdata, 1);
203
204 solverData->nfev++;
205
206 /* save current fvec in f_old*/
207 memcpy(solverData->f_old, fvec, *n*sizeof(double));
208
209 error_f = current_fvec_enorm = enorm_(n, fvec);
210
211 memcpy(solverData->fvecScaled, solverData->fvec, *n*sizeof(double));
212
213 while(error_f > *eps && scaledError_f > *eps && delta_x > *eps && delta_f > *eps && delta_x_scaled > *eps)
214 {
215 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
216 {
217 infoStreamPrint(LOG_NLS_V, 0, "\n**** start Iteration: %d *****", (int) l);
218
219 /* Debug output */
220 infoStreamPrint(LOG_NLS_V, 1, "function values");
221 for(i=0; i<*n; i++)
222 infoStreamPrint(LOG_NLS_V, 0, "fvec[%d]: %e ", i, fvec[i]);
223 messageClose(LOG_NLS_V);
224 }
225
226 /* calculate jacobian if no matrix is given */
227 if (calc_jac == 1 && solverData->calculate_jacobian >= 0)
228 {
229 (*f)(n, x, fvec, userdata, 0);
230 solverData->factorization = 0;
231 calc_jac = solverData->calculate_jacobian;
232 }
233 else
234 {
235 solverData->factorization = 1;
236 calc_jac--;
237 }
238
239
240 /* debug output */
241 if(ACTIVE_STREAM(LOG_NLS_JAC)(useStream[LOG_NLS_JAC]))
242 {
243 char *buffer = (char*)malloc(sizeof(char)*solverData->n*15);
244
245 infoStreamPrint(LOG_NLS_JAC, 1, "jacobian matrix [%dx%d]", (int)*n, (int)*n);
246 for(i=0; i<solverData->n;i++)
247 {
248 buffer[0] = 0;
249 for(j=0; j<solverData->n; j++)
250 sprintf(buffer, "%s%10g ", buffer, fjac[i*(*n)+j]);
251 infoStreamPrint(LOG_NLS_JAC, 0, "%s", buffer);
252 }
253 messageClose(LOG_NLS_JAC);
254 free(buffer);
255 }
256
257 if (solveLinearSystem(n, iwork, fvec, fjac, solverData) != 0)
258 {
259 *info=-1;
260 break;
261 }
262 else
263 {
264 for (i =0; i<*n; i++)
265 solverData->x_new[i]=x[i]-solverData->x_increment[i];
266
267 infoStreamPrint(LOG_NLS_V,1,"x_increment");
268 for(i=0; i<*n; i++)
269 infoStreamPrint(LOG_NLS_V, 0, "x_increment[%d] = %e ", i, solverData->x_increment[i]);
270 messageClose(LOG_NLS_V);
271
272 if (solverData->newtonStrategy == NEWTON_DAMPED)
273 {
274 damping_heuristic(x, f, current_fvec_enorm, n, fvec, &lambda, &k, solverData, userdata);
275 }
276 else if (solverData->newtonStrategy == NEWTON_DAMPED2)
277 {
278 damping_heuristic2(0.75, x, f, current_fvec_enorm, n, fvec, &k, solverData, userdata);
279 }
280 else if (solverData->newtonStrategy == NEWTON_DAMPED_LS)
281 {
282 LineSearch(x, f, current_fvec_enorm, n, fvec, &k, solverData, userdata);
283 }
284 else if (solverData->newtonStrategy == NEWTON_DAMPED_BT)
285 {
286 Backtracking(x, f, current_fvec_enorm, n, fvec, solverData, userdata);
287 }
288 else
289 {
290 /* calculate the function values */
291 (*f)(n, solverData->x_new, fvec, userdata, 1);
292 solverData->nfev++;
293 }
294
295 calculatingErrors(solverData, &delta_x, &delta_x_scaled, &delta_f, &error_f, &scaledError_f, n, x, fvec);
296
297 /* updating x */
298 memcpy(x, solverData->x_new, *n*sizeof(double));
299
300 /* updating f_old */
301 memcpy(solverData->f_old, fvec, *n*sizeof(double));
302
303 current_fvec_enorm = error_f;
304
305 /* check if maximum iteration is reached */
306 if (++l > *maxfev)
307 {
308 *info = -1;
309 warningStreamPrint(LOG_NLS_V, 0, "Warning: maximal number of iteration reached but no root found");
310 break;
311 }
312 /* check if maximum iteration is reached */
313 if (k > 5)
314 {
315 *info = -1;
316 warningStreamPrint(LOG_NLS_V, 0, "Warning: maximal number threshold reached");
317 break;
318 }
319 }
320
321 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
322 {
323 infoStreamPrint(LOG_NLS_V,1,"x vector");
324 for(i=0; i<*n; i++)
325 infoStreamPrint(LOG_NLS_V, 0, "x[%d] = %e ", i, x[i]);
326 messageClose(LOG_NLS_V);
327 printErrors(delta_x, delta_x_scaled, delta_f, error_f, scaledError_f, eps);
328 }
329 }
330
331 solverData->numberOfIterations += l;
332 solverData->numberOfFunctionEvaluations += solverData->nfev;
333
334 return 0;
335}
336
337
338/*! \fn printErrors
339 *
340 * function prints errors, that reached tolerance
341 */
342
343void printErrors(double delta_x, double delta_x_scaled, double delta_f, double error_f, double scaledError_f, double* eps)
344{
345 infoStreamPrint(LOG_NLS_V, 1, "errors ");
346 infoStreamPrint(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);
347
348 if (delta_x < *eps)
349 infoStreamPrint(LOG_NLS_V, 0, "delta_x reached eps");
350 if (delta_x_scaled < *eps)
351 infoStreamPrint(LOG_NLS_V, 0, "delta_x_scaled reached eps");
352 if (delta_f < *eps)
353 infoStreamPrint(LOG_NLS_V, 0, "delta_f reached eps");
354 if (error_f < *eps)
355 infoStreamPrint(LOG_NLS_V, 0, "error_f reached eps");
356 if (scaledError_f < *eps)
357 infoStreamPrint(LOG_NLS_V, 0, "scaledError_f reached eps");
358
359 messageClose(LOG_NLS_V);
360}
361
362/*! \fn solveLinearSystem
363 *
364 * function solves linear system J*(x_{n+1} - x_n) = f using lapack
365 */
366int solveLinearSystem(int* n, int* iwork, double* fvec, double *fjac, DATA_NEWTON* solverData)
367{
368 int i, nrsh=1, lapackinfo;
369 char trans = 'N';
370
371 /* if no factorization is given, calculate it */
372 if (solverData->factorization == 0)
373 {
374 /* solve J*(x_{n+1} - x_n)=f */
375 dgetrf_(n, n, fjac, n, iwork, &lapackinfo);
376 solverData->factorization = 1;
377 dgetrs_(&trans, n, &nrsh, fjac, n, iwork, fvec, n, &lapackinfo);
378 }
379 else
380 {
381 dgetrs_(&trans, n, &nrsh, fjac, n, iwork, fvec, n, &lapackinfo);
382 }
383
384 if(lapackinfo > 0)
385 {
386 warningStreamPrint(LOG_NLS, 0, "Jacobian Matrix singular!");
387 return -1;
388 }
389 else if(lapackinfo < 0)
390 {
391 warningStreamPrint(LOG_NLS, 0, "illegal input in argument %d", (int)lapackinfo);
392 return -1;
393 }
394 else
395 {
396 /* save solution of J*(x_{n+1} - x_n)=f */
397 memcpy(solverData->x_increment, fvec, *n*sizeof(double));
398 }
399
400 return 0;
401}
402
403/*! \fn calculatingErrors
404 *
405 * function calculates the errors
406 */
407void calculatingErrors(DATA_NEWTON* solverData, double* delta_x, double* delta_x_scaled, double* delta_f, double* error_f,
408 double* scaledError_f, int* n, double* x, double* fvec)
409{
410 int i=0;
411 double scaling_factor;
412
413 /* delta_x = || x_new-x_old || */
414 for (i=0; i<*n; i++)
415 solverData->delta_x_vec[i] = x[i]-solverData->x_new[i];
416
417 *delta_x = enorm_(n,solverData->delta_x_vec);
418
419 scaling_factor = enorm_(n,x);
420 if (scaling_factor > 1)
421 *delta_x_scaled = *delta_x * 1./ scaling_factor;
422 else
423 *delta_x_scaled = *delta_x;
424
425 /* delta_f = || f_old - f_new || */
426 for (i=0; i<*n; i++)
427 solverData->delta_f[i] = solverData->f_old[i]-fvec[i];
428
429 *delta_f=enorm_(n, solverData->delta_f);
430
431 *error_f = enorm_(n,fvec);
432
433 /* scaling residual vector */
434 scaling_residual_vector(solverData);
435
436 for (i=0; i<*n; i++)
437 solverData->fvecScaled[i]=fvec[i]/solverData->resScaling[i];
438 *scaledError_f = enorm_(n,solverData->fvecScaled);
439}
440
441/*! \fn calculatingErrors
442 *
443 * function scales the residual vector using the jacobian (heuristic)
444 */
445void scaling_residual_vector(DATA_NEWTON* solverData)
446{
447 int i,j,k;
448 for(i=0, k=0; i<solverData->n; i++)
449 {
450 solverData->resScaling[i] = 0.0;
451 for(j=0; j<solverData->n; j++, ++k)
452 {
453 solverData->resScaling[i] = fmax(fabs(solverData->fjac[k]), solverData->resScaling[i]);
454 }
455 if(solverData->resScaling[i] <= 0.0){
456 warningStreamPrint(LOG_NLS_V, 1, "Jacobian matrix is singular.");
457 solverData->resScaling[i] = 1e-16;
458 }
459 solverData->fvecScaled[i] = solverData->fvec[i] / solverData->resScaling[i];
460 }
461}
462
463/*! \fn damping_heuristic
464 *
465 * first damping heuristic:
466 * x_increment will be halved until the Euclidean norm of the residual function
467 * is smaller than the Euclidean norm of the current point
468 *
469 * treshold for damping = 0.01
470 * compiler flag: -newton = damped
471 */
472void damping_heuristic(double* x, int(*f)(int*, double*, double*, void*, int),
473 double current_fvec_enorm, int* n, double* fvec, double* lambda, int* k, DATA_NEWTON* solverData, void* userdata)
474{
475 int i,j=0;
476 double enorm_new, treshold = 1e-2;
477
478 /* calculate new function values */
479 (*f)(n, solverData->x_new, fvec, userdata, 1);
480 solverData->nfev++;
481
482 enorm_new=enorm_(n,fvec);
483
484 if (enorm_new >= current_fvec_enorm)
485 infoStreamPrint(LOG_NLS_V, 1, "Start Damping: enorm_new : %e; current_fvec_enorm: %e ", enorm_new, current_fvec_enorm);
486
487 while (enorm_new >= current_fvec_enorm)
488 {
489 j++;
490
491 *lambda*=0.5;
492
493
494 for (i=0; i<*n; i++)
495 solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i];
496
497
498 /* calculate new function values */
499 (*f)(n, solverData->x_new, fvec, userdata, 1);
500 solverData->nfev++;
501
502 enorm_new=enorm_(n,fvec);
503
504 if (*lambda <= treshold)
505 {
506 warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda reached a threshold.");
507
508 /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */
509 if (*k >= 5)
510 for (i=0; i<*n; i++)
511 solverData->x_new[i]=x[i]-*lambda*solverData->x_increment[i];
512 else
513 for (i=0; i<*n; i++)
514 solverData->x_new[i]=x[i]-solverData->x_increment[i];
515
516 /* calculate new function values */
517 (*f)(n, solverData->x_new, fvec, userdata, 1);
518 solverData->nfev++;
519
520 (*k)++;
521
522 break;
523 }
524 }
525
526 *lambda = 1;
527
528 messageClose(LOG_NLS_V);
529}
530
531/*! \fn damping_heuristic2
532 *
533 * second (default) damping heuristic:
534 * x_increment will be multiplied by 3/4 until the Euclidean norm of the residual function
535 * is smaller than the Euclidean norm of the current point
536 *
537 * treshold for damping = 0.0001
538 * compiler flag: -newton = damped2
539 */
540void damping_heuristic2(double damping_parameter, double* x, int(*f)(int*, double*, double*, void*, int),
541 double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata)
542{
543 int i,j=0;
544 double enorm_new, treshold = 1e-4, lambda=1;
545
546 /* calculate new function values */
547 (*f)(n, solverData->x_new, fvec, userdata, 1);
548 solverData->nfev++;
549
550 enorm_new=enorm_(n,fvec);
551
552 if (enorm_new >= current_fvec_enorm)
553 infoStreamPrint(LOG_NLS_V, 1, "StartDamping: ");
554
555 while (enorm_new >= current_fvec_enorm)
556 {
557 j++;
558
559 lambda*=damping_parameter;
560
561 infoStreamPrint(LOG_NLS_V, 0, "lambda = %e, k = %d", lambda, *k);
562
563 for (i=0; i<*n; i++)
564 solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i];
565
566
567 /* calculate new function values */
568 (*f)(n, solverData->x_new, fvec, userdata, 1);
569 solverData->nfev++;
570
571 enorm_new=enorm_(n,fvec);
572
573 if (lambda <= treshold)
574 {
575 warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda reached a threshold.");
576
577 /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */
578 if (*k >= 5)
579 for (i=0; i<*n; i++)
580 solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i];
581 else
582 for (i=0; i<*n; i++)
583 solverData->x_new[i]=x[i]-solverData->x_increment[i];
584
585 /* calculate new function values */
586 (*f)(n, solverData->x_new, fvec, userdata, 1);
587 solverData->nfev++;
588
589 (*k)++;
590
591 break;
592 }
593 }
594
595 messageClose(LOG_NLS_V);
596}
597
598/*! \fn LineSearch
599 *
600 * third damping heuristic:
601 * Along the tangent 5 five points are selected. For every point the Euclidean norm of
602 * the residual function will be calculated and the minimum is chosen for the further iteration.
603 *
604 * compiler flag: -newton = damped_ls
605 */
606void LineSearch(double* x, int(*f)(int*, double*, double*, void*, int),
607 double current_fvec_enorm, int* n, double* fvec, int* k, DATA_NEWTON* solverData, void* userdata)
608{
609 int i,j;
610 double enorm_new, enorm_minimum=current_fvec_enorm, lambda_minimum=0;
611 double lambda[5]={1.25,1,0.75,0.5,0.25};
612
613
614 for (j=0; j<5; j++)
615 {
616 for (i=0; i<*n; i++)
617 solverData->x_new[i]=x[i]-lambda[j]*solverData->x_increment[i];
618
619 /* calculate new function values */
620 (*f)(n, solverData->x_new, fvec, userdata, 1);
621 solverData->nfev++;
622
623 enorm_new=enorm_(n,fvec);
624
625 /* searching minimal enorm */
626 if (enorm_new < enorm_minimum)
627 {
628 enorm_minimum = enorm_new;
629 lambda_minimum = lambda[j];
630 memcpy(solverData->fvec_minimum, fvec,*n*sizeof(double));
631 }
632 }
633
634 infoStreamPrint(LOG_NLS_V,0,"lambda_minimum = %e", lambda_minimum);
635
636 if (lambda_minimum == 0)
637 {
638 warningStreamPrint(LOG_NLS_V, 0, "Warning: lambda_minimum = 0 ");
639
640 /* if damping is without success, trying full newton step; after 5 full newton steps try a very little step */
641 if (*k >= 5)
642 {
643 lambda_minimum = 0.125;
644
645 /* calculate new function values */
646 (*f)(n, solverData->x_new, fvec, userdata, 1);
647 solverData->nfev++;
648 }
649 else
650 {
651 lambda_minimum = 1;
652
653 /* calculate new function values */
654 (*f)(n, solverData->x_new, fvec, userdata, 1);
655 solverData->nfev++;
656 }
657
658 (*k)++;
659 }
660 else
661 {
662 /* save new function values */
663 memcpy(fvec, solverData->fvec_minimum, *n*sizeof(double));
664 }
665
666 for (i=0; i<*n; i++)
667 solverData->x_new[i]=x[i]-lambda_minimum*solverData->x_increment[i];
668}
669
670/*! \fn Backtracking
671 *
672 * forth damping heuristic:
673 * Calculate new function h:R^n->R ; h(x) = 1/2 * ||f(x)|| ^2
674 * g(lambda) = h(x_old + lambda * x_increment)
675 * find minimum of g with golden ratio method
676 * tau = golden ratio
677 *
678 * compiler flag: -newton = damped_bt
679 */
680void Backtracking(double* x, int(*f)(int*, double*, double*, void*, int),
681 double current_fvec_enorm, int* n, double* fvec, DATA_NEWTON* solverData, void* userdata)
682{
683 int i,j;
684 double enorm_new, enorm_f, lambda, a1, b1, a, b, tau, g1, g2;
685 double tolerance = 1e-3;
686
687 /* saving current function values in f_old */
688 memcpy(solverData->f_old, fvec, *n*sizeof(double));
689
690 for (i=0; i<*n; i++)
691 solverData->x_new[i]=x[i]-solverData->x_increment[i];
692
693 /* calculate new function values */
694 (*f)(n, solverData->x_new, fvec, userdata, 1);
695 solverData->nfev++;
696
697
698 /* calculate new enorm */
699 enorm_new = enorm_(n,fvec);
700
701 /* Backtracking only if full newton step is useless */
702 if (enorm_new >= current_fvec_enorm)
703 {
704 infoStreamPrint(LOG_NLS_V, 0, "Start Backtracking\n enorm_new= %f \t current_fvec_enorm=%f",enorm_new, current_fvec_enorm);
705
706 /* h(x) = 1/2 * ||f(x)|| ^2
707 * g(lambda) = h(x_old + lambda * x_increment)
708 * find minimum of g with golden ratio method
709 * tau = golden ratio
710 * */
711
712 a = 0;
713 b = 1;
714 tau = 0.61803398875;
715
716 a1 = a + (1-tau)*(b-a);
717 /* g1 = g(a1) = h(x_old - a1 * x_increment) = 1/2 * ||f(x_old- a1 * x_increment)||^2 */
718 solverData->x_new[i] = x[i]- a1 * solverData->x_increment[i];
719 (*f)(n, solverData->x_new, fvec, userdata, 1);
720 solverData->nfev++;
721 enorm_f= enorm_(n,fvec);
722 g1 = 0.5 * enorm_f * enorm_f;
723
724
725 b1 = a + tau * (b-a);
726 /* g2 = g(b1) = h(x_old - b1 * x_increment) = 1/2 * ||f(x_old- b1 * x_increment)||^2 */
727 solverData->x_new[i] = x[i]- b1 * solverData->x_increment[i];
728 (*f)(n, solverData->x_new, fvec, userdata, 1);
729 solverData->nfev++;
730 enorm_f= enorm_(n,fvec);
731 g2 = 0.5 * enorm_f * enorm_f;
732
733 while ( (b - a) > tolerance)
734 {
735 if (g1<g2)
736 {
737 b = b1;
738 b1 = a1;
739 a1 = a + (1-tau)*(b-a);
740 g2 = g1;
741
742 /* g1 = g(a1) = h(x_old - a1 * x_increment) = 1/2 * ||f(x_old- a1 * x_increment)||^2 */
743 solverData->x_new[i] = x[i]- a1 * solverData->x_increment[i];
744 (*f)(n, solverData->x_new, fvec, userdata, 1);
745 solverData->nfev++;
746 enorm_f= enorm_(n,fvec);
747 g1 = 0.5 * enorm_f * enorm_f;
748 }
749 else
750 {
751 a = a1;
752 a1 = b1;
753 b1 = a + tau * (b-a);
754 g1 = g2;
755
756 /* g2 = g(b1) = h(x_old - b1 * x_increment) = 1/2 * ||f(x_old- b1 * x_increment)||^2 */
757 solverData->x_new[i] = x[i]- b1 * solverData->x_increment[i];
758 (*f)(n, solverData->x_new, fvec, userdata, 1);
759 solverData->nfev++;
760 enorm_f= enorm_(n,fvec);
761 g2 = 0.5 * enorm_f * enorm_f;
762 }
763 }
764
765
766 lambda = (a+b)/2;
767
768
769 /* print lambda */
770 infoStreamPrint(LOG_NLS_V, 0, "Backtracking - lambda = %e", lambda);
771
772 for (i=0; i<*n; i++)
773 solverData->x_new[i]=x[i]-lambda*solverData->x_increment[i];
774
775 /* calculate new function values */
776 (*f)(n, solverData->x_new, fvec, userdata, 1);
777 solverData->nfev++;
778 }
779}
780
781
782#ifdef __cplusplus
783}
784#endif