Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/nonlinearSolverHomotopy.c
Warning:line 1355, column 7
Value stored to 'assert' 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#if !defined(OMC_NUM_LINEAR_SYSYTEMS) || OMC_NUM_NONLINEAR_SYSTEMS>0
32
33#if defined(__AVR__)
34#warning "AVR CPUs are not suitable for non-linear solvers"
35#endif
36
37/*! \file nonlinearSolverHomotopy.c
38* \author bbachmann
39*/
40
41#include <math.h>
42#include <stdlib.h>
43#include <string.h> /* memcpy */
44
45#include "../options.h"
46#include "../simulation_info_json.h"
47#include "../../util/omc_error.h"
48#include "../../util/omc_file.h"
49#include "../../util/varinfo.h"
50#include "model_help.h"
51#include "../../meta/meta_modelica.h"
52#if !defined(OMC_MINIMAL_RUNTIME)
53#include "../../util/write_csv.h"
54#endif
55
56#include "nonlinearSystem.h"
57#include "nonlinearSolverHomotopy.h"
58#include "nonlinearSolverHybrd.h"
59
60#ifdef __cplusplus
61extern "C" {
62#endif
63
64extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
65
66#ifdef __cplusplus
67}
68#endif
69
70/*! \typedef DATA_HOMOTOPY
71 * define memory structure for nonlinear system solver
72 * \author bbachmann
73 */
74typedef struct DATA_HOMOTOPY
75{
76 int initialized; /* 1 = initialized, else = 0*/
77
78 int n; /* dimension; n == size */
79 int m; /* dimension: m == size+1 */
80
81 double xtol_sqrd; /* tolerance for updating solution vector */
82 double ftol_sqrd; /* tolerance for accepting accuracy */
83
84 double error_f_sqrd;
85
86 double* resScaling; /* residual scaling */
87 double* fvecScaled; /* function values scaled */
88 double* hvecScaled; /* function values scaled */
89 double* dxScaled; /* scaled solution vector */
90
91 double* minValue; /* min-attribute of variable, only pointer */
92 double* maxValue; /* max-attribute of variable, only pointer */
93 double* xScaling; /* nominal-attrbute [x.nominal,lambda.nominal] with lambda.nominal=1.0 */
94
95 /* used in wrapper_*/
96 double* f1;
97 double* f2;
98 /* used for steepest descent method */
99 double* gradFx;
100
101 /* return value, if success info == 1 */
102 int info;
103 int numberOfIterations; /* over the whole simulation time */
104 int numberOfFunctionEvaluations; /* over the whole simulation time */
105 int maxNumberOfIterations; /* number of Newton steps */
106
107 /* strict tearing set or casual tearing set */
108 int casualTearingSet;
109
110 /* newton algorithm*/
111 double* x;
112 double* x0;
113 double* xStart;
114 double* x1;
115 double* finit;
116 double* fx0;
117 double* fJac;
118 double* fJacx0;
119
120 /* debug arrays */
121 double* debug_fJac;
122 double* debug_dx;
123
124 /* homotopy parameters */
125 int initHomotopy; /* homotopy method used for the initialization with lambda from the homotopy()-operator */
126 double startDirection;
127 double tau;
128 double* y0;
129 double* y1;
130 double* y2;
131 double* yt;
132 double* dy0;
133 double* dy1;
134 double* dy2;
135 double* hvec;
136 double* hJac;
137 double* hJac2;
138 double* hJacInit;
139 double* ones;
140
141 /* linear system */
142 int* indRow;
143 int* indCol;
144
145 int (*f) (struct DATA_HOMOTOPY*, double*, double*);
146 int (*f_con) (struct DATA_HOMOTOPY*, double*, double*);
147 int (*fJac_f) (struct DATA_HOMOTOPY*, double*, double*);
148 int (*h_function)(struct DATA_HOMOTOPY*, double*, double*);
149 int (*hJac_dh) (struct DATA_HOMOTOPY*, double*, double*);
150
151 DATA* data;
152 threadData_t *threadData;
153 int sysNumber;
154 int eqSystemNumber;
155 double timeValue;
156 int mixedSystem;
157
158 void* dataHybrid;
159
160} DATA_HOMOTOPY;
161
162/*! \fn allocateHomotopyData
163 * allocate memory for nonlinear system solver
164 * \author bbachmann
165 */
166int allocateHomotopyData(int size, void** voiddata)
167{
168 DATA_HOMOTOPY* data = (DATA_HOMOTOPY*) malloc(sizeof(DATA_HOMOTOPY));
169
170 *voiddata = (void*)data;
171 assertStreamPrint(NULL, 0 != data, "allocationHomotopyData() failed!")if (!(0 != data)) {throwStreamPrint((((void*)0)), "allocationHomotopyData() failed!"
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "./simulation/solver/nonlinearSolverHomotopy.c"
, 171, __extension__ __PRETTY_FUNCTION__); }));}
;
172
173 data->initialized = 0;
174 data->n = size;
175 data->m = size + 1;
176 data->xtol_sqrd = newtonXTol*newtonXTol;
177 data->ftol_sqrd = newtonFTol*newtonFTol;
178
179 data->error_f_sqrd = 0;
180
181 data->maxNumberOfIterations = size*100;
182 data->numberOfIterations = 0;
183 data->numberOfFunctionEvaluations = 0;
184
185 data->resScaling = (double*) calloc(size,sizeof(double));
186 data->fvecScaled = (double*) calloc(size,sizeof(double));
187 data->hvecScaled = (double*) calloc(size,sizeof(double));
188 data->dxScaled = (double*) calloc(size,sizeof(double));
189
190 data->xScaling = (double*) calloc((size+1),sizeof(double));
191
192 data->f1 = (double*) calloc(size,sizeof(double));
193 data->f2 = (double*) calloc(size,sizeof(double));
194 data->gradFx = (double*) calloc(size,sizeof(double));
195
196 /* damped newton */
197 data->x = (double*) calloc((size+1),sizeof(double));
198 data->x0 = (double*) calloc((size+1),sizeof(double));
199 data->xStart = (double*) calloc(size,sizeof(double));
200 data->x1 = (double*) calloc((size+1),sizeof(double));
201 data->finit = (double*) calloc(size,sizeof(double));
202 data->fx0 = (double*) calloc(size,sizeof(double));
203 data->fJac = (double*) calloc((size*(size+1)),sizeof(double));
204 data->fJacx0 = (double*) calloc((size*(size+1)),sizeof(double));
205
206 /* debug arrays */
207 data->debug_dx = (double*) calloc(size,sizeof(double));
208 data->debug_fJac = (double*) calloc((size*(size+1)),sizeof(double));
209
210 /* homotopy */
211 data->y0 = (double*) calloc((size+1),sizeof(double));
212 data->y1 = (double*) calloc((size+1),sizeof(double));
213 data->y2 = (double*) calloc((size+1),sizeof(double));
214 data->yt = (double*) calloc((size+1),sizeof(double));
215 data->dy0 = (double*) calloc((size+1),sizeof(double));
216 data->dy1 = (double*) calloc((size+homBacktraceStrategy),sizeof(double));
217 data->dy2 = (double*) calloc((size+1),sizeof(double));
218 data->hvec = (double*) calloc(size,sizeof(double));
219 data->hJac = (double*) calloc(size*(size+1),sizeof(double));
220 data->hJac2 = (double*) calloc((size+1)*(size+2),sizeof(double));
221 data->hJacInit = (double*) calloc(size*(size+1),sizeof(double));
222 data->ones = (double*) calloc(size+1,sizeof(double));
223
224 /* linear system */
225 data->indRow =(int*) calloc(size+homBacktraceStrategy-1,sizeof(int));
226 data->indCol =(int*) calloc(size+homBacktraceStrategy,sizeof(int));
227
228 allocateHybrdData(size, &data->dataHybrid);
229
230 assertStreamPrint(NULL, 0 != *voiddata, "allocationHomotopyData() voiddata failed!")if (!(0 != *voiddata)) {throwStreamPrint((((void*)0)), "allocationHomotopyData() voiddata failed!"
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "./simulation/solver/nonlinearSolverHomotopy.c"
, 230, __extension__ __PRETTY_FUNCTION__); }));}
;
231 return 0;
232}
233
234/*! \fn freeHomotopyData
235 *
236 * free memory for nonlinear system solver
237 * \author bbachmann
238 */
239int freeHomotopyData(void **voiddata)
240{
241 DATA_HOMOTOPY* data = (DATA_HOMOTOPY*) *voiddata;
242
243 free(data->resScaling);
244 free(data->fvecScaled);
245 free(data->hvecScaled);
246 free(data->x);
247 free(data->debug_dx);
248 free(data->finit);
249 free(data->f1);
250 free(data->f2);
251 free(data->gradFx);
252 free(data->fJac);
253 free(data->fJacx0);
254 free(data->debug_fJac);
255
256 /* damped newton */
257 free(data->x0);
258 free(data->xStart);
259 free(data->x1);
260 free(data->dxScaled);
261
262 /* homotopy */
263 free(data->fx0);
264 free(data->hvec);
265 free(data->hJac);
266 free(data->hJac2);
267 free(data->hJacInit);
268 free(data->y0);
269 free(data->y1);
270 free(data->y2);
271 free(data->yt);
272 free(data->dy0);
273 free(data->dy1);
274 free(data->dy2);
275 free(data->xScaling);
276 free(data->ones);
277
278 /* linear system */
279 free(data->indRow);
280 free(data->indCol);
281
282 freeHybrdData(&data->dataHybrid);
283
284 return 0;
285}
286
287/* Prototypes for debug functions
288 * \author bbachmann
289 */
290
291void printUnknowns(int logName, DATA_HOMOTOPY *solverData)
292{
293 long i;
294 int eqSystemNumber = solverData->eqSystemNumber;
295 DATA *data = solverData->data;
296
297 if (!ACTIVE_STREAM(logName)(useStream[logName])) return;
298 infoStreamPrint(logName, 1, "nls status");
299 infoStreamPrint(logName, 1, "variables");
300 messageClose(logName);
301
302 for(i=0; i<solverData->n; i++)
303 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g\t\t min = %16.8g\t\t max = %16.8g", i+1,
304 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
305 solverData->x[i], solverData->xScaling[i], solverData->minValue[i], solverData->maxValue[i]);
306 messageClose(logName);
307}
308
309void printNewtonStep(int logName, DATA_HOMOTOPY *solverData)
310{
311 long i;
312 int eqSystemNumber = solverData->eqSystemNumber;
313 DATA *data = solverData->data;
314
315 if (!ACTIVE_STREAM(logName)(useStream[logName])) return;
316 infoStreamPrint(logName, 1, "newton step");
317 infoStreamPrint(logName, 1, "variables");
318 messageClose(logName);
319
320 for(i=0; i<solverData->n; i++)
321 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t step = %16.8g\t\t old = %16.8g", i+1,
322 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
323 solverData->x1[i], solverData->dy0[i], solverData->x[i]);
324 messageClose(logName);
325}
326
327void printHomotopyUnknowns(int logName, DATA_HOMOTOPY *solverData)
328{
329 long i;
330 int eqSystemNumber = solverData->eqSystemNumber;
331 DATA *data = solverData->data;
332
333 if (!ACTIVE_STREAM(logName)(useStream[logName])) return;
334 infoStreamPrint(logName, 1, "homotopy status");
335 infoStreamPrint(logName, 1, "variables");
336 messageClose(logName);
337
338 for(i=0; i<solverData->n; i++)
339 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g\t\t min = %16.8g\t\t max = %16.8g", i+1,
340 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
341 solverData->y0[i], solverData->xScaling[i], solverData->minValue[i], solverData->maxValue[i]);
342 if (solverData->initHomotopy) {
343 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g\t\t min = %16.8g\t\t max = %16.8g", i+1,
344 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
345 solverData->y0[i], solverData->xScaling[i], solverData->minValue[i], solverData->maxValue[i]);
346 }
347 else {
348 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t nom = %16.8g", i+1,
349 "LAMBDA",
350 solverData->y0[solverData->n], solverData->xScaling[solverData->n]);
351 }
352 messageClose(logName);
353}
354
355void printHomotopyPredictorStep(int logName, DATA_HOMOTOPY *solverData)
356{
357 long i;
358 int eqSystemNumber = solverData->eqSystemNumber;
359 DATA *data = solverData->data;
360
361 if (!ACTIVE_STREAM(logName)(useStream[logName])) return;
362 infoStreamPrint(logName, 1, "predictor status");
363 infoStreamPrint(logName, 1, "variables");
364 messageClose(logName);
365
366 for(i=0; i<solverData->n; i++)
367 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
368 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
369 solverData->yt[i], solverData->dy0[i], solverData->y0[i], solverData->tau);
370 if (solverData->initHomotopy) {
371 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
372 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
373 solverData->yt[i], solverData->dy0[i], solverData->y0[i], solverData->tau);
374 } else {
375 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
376 "LAMBDA",
377 solverData->yt[solverData->n], solverData->dy0[i], solverData->y0[i], solverData->tau);
378 }
379 messageClose(logName);
380}
381
382void printHomotopyCorrectorStep(int logName, DATA_HOMOTOPY *solverData)
383{
384 long i;
385 int eqSystemNumber = solverData->eqSystemNumber;
386 DATA *data = solverData->data;
387
388 if (!ACTIVE_STREAM(logName)(useStream[logName])) return;
389 infoStreamPrint(logName, 1, "corrector status");
390 infoStreamPrint(logName, 1, "variables");
391 messageClose(logName);
392
393 for(i=0; i<solverData->n; i++)
394 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
395 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
396 solverData->y1[i], solverData->dy1[i], solverData->yt[i], solverData->tau);
397 if (solverData->initHomotopy) {
398 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
399 modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i],
400 solverData->y1[i], solverData->dy1[i], solverData->yt[i], solverData->tau);
401 } else {
402 infoStreamPrint(logName, 0, "[%2ld] %30s = %16.8g\t\t dy = %16.8g\t\t old = %16.8g\t\t tau = %16.8g", i+1,
403 "LAMBDA",
404 solverData->y1[solverData->n], solverData->dy1[i], solverData->yt[i], solverData->tau);
405 }
406 messageClose(logName);
407}
408
409void debugMatrixPermutedDouble(int logName, char* matrixName, double* matrix, int n, int m, int* indRow, int* indCol)
410{
411 if(ACTIVE_STREAM(logName)(useStream[logName]))
412 {
413 int i, j;
414 int sparsity = 0;
415 char *buffer = (char*)malloc(sizeof(char)*m*20);
416
417 infoStreamPrint(logName, 1, "%s [%dx%d-dim]", matrixName, n, m);
418 for(i=0; i<n;i++)
419 {
420 buffer[0] = 0;
421 for(j=0; j<m; j++)
422 {
423 if (sparsity)
424 {
425 if (fabs(matrix[indRow[i] + indCol[j]*(m-1)])<1e-12)
426 sprintf(buffer, "%s 0", buffer);
427 else
428 sprintf(buffer, "%s *", buffer);
429 }
430 else
431 {
432 sprintf(buffer, "%s%16.8g ", buffer, matrix[indRow[i] + indCol[j]*(m-1)]);
433 }
434 }
435 infoStreamPrint(logName, 0, "%s", buffer);
436 }
437 messageClose(logName);
438 free(buffer);
439 }
440}
441
442void debugMatrixDouble(int logName, char* matrixName, double* matrix, int n, int m)
443{
444 if(ACTIVE_STREAM(logName)(useStream[logName]))
445 {
446 int i, j;
447 int sparsity = 0;
448 char *buffer = (char*)malloc(sizeof(char)*m*20);
449
450 infoStreamPrint(logName, 1, "%s [%dx%d-dim]", matrixName, n, m);
451 for(i=0; i<n;i++)
452 {
453 buffer[0] = 0;
454 for(j=0; j<m; j++)
455 {
456 if (sparsity)
457 {
458 if (fabs(matrix[i + j*(m-1)])<1e-12)
459 sprintf(buffer, "%s 0", buffer);
460 else
461 sprintf(buffer, "%s *", buffer);
462 }
463 else
464 {
465 sprintf(buffer, "%s%16.8g ", buffer, matrix[i + j*(m-1)]);
466 }
467 }
468 infoStreamPrint(logName, 0, "%s", buffer);
469 }
470 messageClose(logName);
471 free(buffer);
472 }
473}
474
475void debugVectorDouble(int logName, char* vectorName, double* vector, int n)
476{
477 if(ACTIVE_STREAM(logName)(useStream[logName]))
478 {
479 int i;
480 char *buffer = (char*)malloc(sizeof(char)*n*20);
481
482 infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
483 buffer[0] = 0;
484 for(i=0; i<n;i++)
485 {
486 if (vector[i]<-1e+300)
487 sprintf(buffer, "%s -INF ", buffer);
488 else if (vector[i]>1e+300)
489 sprintf(buffer, "%s +INF ", buffer);
490 else
491 sprintf(buffer, "%s%16.8g ", buffer, vector[i]);
492 }
493 infoStreamPrint(logName, 0, "%s", buffer);
494 messageClose(logName);
495 free(buffer);
496 }
497}
498
499void debugVectorBool(int logName, char* vectorName, modelica_boolean* vector, int n)
500{
501 if(ACTIVE_STREAM(logName)(useStream[logName]))
502 {
503 int i;
504 char *buffer = (char*)malloc(sizeof(char)*n*20);
505
506 infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
507 buffer[0] = 0;
508 for(i=0; i<n;i++)
509 {
510 if (vector[i]<-1e+300)
511 sprintf(buffer, "%s -INF ", buffer);
512 else if (vector[i]>1e+300)
513 sprintf(buffer, "%s +INF ", buffer);
514 else
515 sprintf(buffer, "%s %d", buffer, vector[i]);
516 }
517 infoStreamPrint(logName, 0, "%s", buffer);
518 messageClose(logName);
519 free(buffer);
520 }
521}
522
523void debugVectorInt(int logName, char* vectorName, int* vector, int n)
524{
525 if(ACTIVE_STREAM(logName)(useStream[logName]))
526 {
527 int i;
528 char *buffer = (char*)malloc(sizeof(char)*n*20);
529
530 infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
531 buffer[0] = 0;
532 for(i=0; i<n;i++)
533 {
534 if (vector[i]<-1e+300)
535 sprintf(buffer, "%s -INF ", buffer);
536 else if (vector[i]>1e+300)
537 sprintf(buffer, "%s +INF ", buffer);
538 else
539 sprintf(buffer, "%s%d ", buffer, vector[i]);
540 }
541 infoStreamPrint(logName, 0, "%s", buffer);
542 messageClose(logName);
543 free(buffer);
544 }
545}
546
547
548/* Prototypes for linear algebra functions
549 * \author bbachmann
550 */
551
552double vec2Norm(int n, double *x)
553{
554 int i;
555 double norm=0.0;
556 for (i=0;i<n;i++)
557 norm+=x[i]*x[i];
558 return sqrt(norm);
559}
560
561double vec2NormSqrd(int n, double *x)
562{
563 int i;
564 double norm=0.0;
565 for (i=0;i<n;i++)
566 norm+=x[i]*x[i];
567 return norm;
568}
569
570double vecMaxNorm(int n, double *x)
571{
572 int i;
573 double norm=fabs(x[0]);
574 for (i=1;i<n;i++)
575 if (fabs(x[i])>norm)
576 norm=fabs(x[i]);
577 return norm;
578}
579
580void vecAdd(int n, double *a, double *b, double *c)
581{
582 int i;
583 for (i=0;i<n;i++)
584 c[i] = a[i] + b[i];
585}
586
587void vecAddScal(int n, double *a, double *b, double s, double *c)
588{
589 int i;
590 for (i=0;i<n;i++)
591 c[i] = a[i] + s*b[i];
592}
593
594void vecScalarMult(int n, double *a, double s, double *b)
595{
596 int i;
597 for (i=0;i<n;i++)
598 b[i] = s*a[i];
599}
600
601void vecLinearComb(int n, double *a, double r, double *b, double s, double *c)
602{
603 int i;
604 for (i=0;i<n;i++)
605 c[i] = r*a[i] + s*b[i];
606}
607
608void vecCopy(int n, double *a, double *b)
609{
610 memcpy(b, a, n*(sizeof(double)));
611}
612
613void vecCopyBool(int n, modelica_boolean *a, modelica_boolean *b)
614{
615 memcpy(b, a, n*(sizeof(modelica_boolean)));
616}
617
618void vecAddInv(int n, double *a, double *b)
619{
620 int i;
621 for (i=0;i<n;i++)
622 b[i] = -a[i];
623}
624
625void vecDiff(int n, double *a, double *b, double *c)
626{
627 int i;
628 for (i=0;i<n;i++)
629 c[i] = a[i] - b[i];
630}
631
632int isNotEqualVectorInt(int n, modelica_boolean *a, modelica_boolean *b)
633{
634 int i, isNotEqual = 0;
635 for (i=0;i<n;i++)
636 isNotEqual += abs(a[i] - b[i]);
637 return isNotEqual;
638}
639
640void vecMultScaling(int n, double *a, double *b, double *c)
641{
642 int i;
643 for (i=0;i<n;i++)
644 c[i] = (fabs(b[i])>0 ? a[i]*fabs(b[i]):a[i]);
645}
646
647void vecDivScaling(int n, double *a, double *b, double *c)
648{
649 int i;
650 for (i=0;i<n;i++)
651 c[i] = (fabs(b[i])>0 ? a[i]/fabs(b[i]):a[i]);
652}
653
654void vecNormalize(int n, double *a, double *b)
655{
656 int i;
657 double norm = vec2Norm(n,a);
658 for (i=0;i<n;i++)
659 b[i] = (norm>0 ? a[i]/norm:a[i]);
660}
661
662void vecConst(int n, double value, double *a)
663{
664 int i;
665 for (i=0;i<n;i++)
666 a[i] = value;
667}
668
669double vecScalarProd(int n, double *a, double *b)
670{
671 int i;
672 double prod;
673
674 for (i=0,prod=0;i<n;i++)
675 prod = prod + a[i]*b[i];
676
677 return prod;
678}
679
680/* Matrix has dimension [n x m], vector [m] */
681void matVecMult(int n, int m, double *A, double *b, double *c)
682{
683 int i, j;
684 for (i=0;i<n;i++) {
685 c[i] = 0.0;
686 for (j=0;j<m;j++)
687 c[i] += A[i+j*(m-1)]*b[j];
688 }
689}
690
691/* Matrix has dimension [n x m], vector [m] */
692void matVecMultAbs(int n, int m, double *A, double *b, double *c)
693{
694 int i, j;
695 for (i=0;i<n;i++) {
696 c[i] = 0.0;
697 for (j=0;j<m;j++)
698 c[i] += fabs(A[i+j*(m-1)]*b[j]);
699 }
700}
701
702/* Matrix has dimension [n x (n+1)] */
703void matVecMultBB(int n, double *A, double *b, double *c)
704{
705 int i, j;
706 for (i=0;i<n;i++) {
707 c[i] = 0.0;
708 for (j=0;j<n;j++)
709 c[i] += A[i+j*n]*b[j];
710 }
711}
712
713/* Matrix has dimension [n x (n+1)] */
714void matVecMultAbsBB(int n, double *A, double *b, double *c)
715{
716 int i, j;
717 for (i=0;i<n;i++) {
718 c[i] = 0.0;
719 for (j=0;j<n;j++)
720 c[i] += fabs(A[i+j*n]*b[j]);
721 }
722}
723
724/* Matrix has dimension [n x (n+1)] */
725void matAddBB(int n, double* A, double* B, double* C)
726{
727 int i, j;
728
729 for (i=0;i<n;i++) {
730 for (j=0;j<n+1;j++)
731 C[i + j*n] = A[i + j*n] + B[i + j*n];
732 }
733}
734
735/* Matrix has dimension [n x (n+1)] */
736void matDiffBB(int n, double* A, double* B, double* C)
737{
738 int i, j;
739
740 for (i=0;i<n;i++) {
741 for (j=0;j<n;j++)
742 C[i + j*n] = A[i + j*n] - B[i + j*n];
743 }
744}
745
746/* Matrix has dimension [n x m] */
747void scaleMatrixRows(int n, int m, double *A)
748{
749 const double delta = sqrt(DBL_EPSILON2.2204460492503131e-16);
750 int i, j;
751 double rowMax;
752 for (i=0;i<n;i++) {
753 rowMax = 0; /* This might be changed to delta */
754 for (j=0;j<n;j++) {
755 if (fabs(A[i+j*(m-1)]) > rowMax) {
756 rowMax = fabs(A[i+j*(m-1)]);
757 }
758 }
759 if (rowMax>0) {
760 for (j=0;j<m;j++)
761 A[i+j*(m-1)] /= rowMax;
762 }
763 }
764}
765
766/* Build the newton matrix for the corrector step with orthogonal backtrace strategy */
767void orthogonalBacktraceMatrix(DATA_HOMOTOPY* solverData, double* hJac, double* hvec, double* v, double* hJac2, int n, int m)
768{
769 int i, j;
770 for (i=0; i<n; i++)
771 {
772 for (j=0; j<m; j++)
773 {
774 hJac2[i + j*m] = hJac[i + j*(m-1)];
775 }
776 hJac2[i + m*m] = hvec[i];
777 }
778 for (j=0; j<m; j++)
779 {
780 hJac2[n + j*m] = v[j];
781 }
782 hJac2[n + m*m] = 0;
783}
784
785void swapPointer(double* *p1, double* *p2)
786{
787 double* help;
788 help = *p1;
789 *p1 = *p2;
790 *p2 = help;
791}
792
793/*! \fn getAnalyticalJacobian
794 *
795 * function calculates analytical jacobian
796 *
797 * \param [ref] [data]
798 * \param [out] [jac]
799 *
800 * \author wbraun
801 * bbachmann: introduce scaling factor
802 *
803 */
804int getAnalyticalJacobianHomotopy(DATA_HOMOTOPY* solverData, double* jac)
805{
806 DATA* data = solverData->data;
807 threadData_t *threadData = solverData->threadData;
808 int i,j,k,l,ii;
809 NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[solverData->sysNumber]);
810 const int index = systemData->jacobianIndex;
811 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]);
812
813 memset(jac, 0, (solverData->n)*(solverData->n)*sizeof(double));
814
815 if (jacobian->constantEqns != NULL((void*)0)) {
816 jacobian->constantEqns(data, threadData, jacobian, NULL((void*)0));
817 }
818
819 for(i=0; i < jacobian->sparsePattern->maxColors; i++)
820 {
821 /* activate seed variable for the corresponding color */
822 for(ii=0; ii < jacobian->sizeCols; ii++)
823 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
824 jacobian->seedVars[ii] = 1;
825
826 ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, NULL((void*)0));
827
828 for(j = 0; j < jacobian->sizeCols; j++)
829 {
830 if(jacobian->seedVars[j] == 1)
831 {
832 ii = jacobian->sparsePattern->leadindex[j];
833 while(ii < jacobian->sparsePattern->leadindex[j+1])
834 {
835 l = jacobian->sparsePattern->index[ii];
836 k = j*jacobian->sizeRows + l;
837 /* Calculate scaled difference quotient */
838 jac[k] = jacobian->resultVars[l] * solverData->xScaling[j];
839 ii++;
840 };
841 }
842 /* de-activate seed variable for the corresponding color */
843 if(jacobian->sparsePattern->colorCols[j]-1 == i)
844 jacobian->seedVars[j] = 0;
845 }
846 }
847
848 return 0;
849}
850
851/*! \fn getNumericalJacobianHomotopy
852 *
853 * function calculates a jacobian matrix by
854 * numerical method finite differences
855 * \author bbachmann
856 *
857*/
858static int getNumericalJacobianHomotopy(DATA_HOMOTOPY* solverData, double *x, double *fJac)
859{
860 const double delta_h = sqrt(DBL_EPSILON2.2204460492503131e-16*2e1);
861 double delta_hh;
862 double xsave;
863 int i,j,l;
864
865 if (solverData->initHomotopy) {
866 /* Use the homotopy function values hvec and also calculate the lambda column */
867 for(i = 0; i < solverData->n+1; i++) {
868 xsave = x[i];
869 delta_hh = delta_h * (fabs(xsave) + 1.0);
870 if ((xsave + delta_hh >= solverData->maxValue[i]))
871 delta_hh *= -1;
872 x[i] += delta_hh;
873 /* Calculate scaled difference quotient */
874 delta_hh = 1. / delta_hh * solverData->xScaling[i];
875 solverData->h_function(solverData, x, solverData->f2);
876
877 for(j = 0; j < solverData->n; j++) {
878 l = i * solverData->n + j;
879 fJac[l] = (solverData->f2[j] - solverData->hvec[j]) * delta_hh; /* solverData->hvec must be set outside this function based on x */
880 }
881 x[i] = xsave;
882 }
883 } else {
884 /* Use the normal function values f2 and calculate jacobian without the lambda column */
885 for(i = 0; i < solverData->n; i++) {
886 xsave = x[i];
887 delta_hh = delta_h * (fabs(xsave) + 1.0);
888 if ((xsave + delta_hh >= solverData->maxValue[i]))
889 delta_hh *= -1;
890 x[i] += delta_hh;
891 /* Calculate scaled difference quotient */
892 delta_hh = 1. / delta_hh * solverData->xScaling[i];
893 if (solverData->casualTearingSet)
894 solverData->f_con(solverData, x, solverData->f2);
895 else
896 solverData->f(solverData, x, solverData->f2);
897
898 for(j = 0; j < solverData->n; j++) {
899 l = i * solverData->n + j;
900 fJac[l] = (solverData->f2[j] - solverData->f1[j]) * delta_hh; /* solverData->f1 must be set outside this function based on x */
901 }
902 x[i] = xsave;
903 }
904 }
905 return 0;
906}
907
908/*! \fn wrapper_fvec for the residual Function
909 * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data)
910 *
911 * \author bbachmann
912 *
913 */
914static int wrapper_fvec(DATA_HOMOTOPY* solverData, double* x, double* f)
915{
916 void *dataAndThreadData[2] = {solverData->data, solverData->threadData};
917 int iflag = 0;
918
919 /*TODO: change input to residualFunc from data to systemData */
920 (solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].residualFunc(dataAndThreadData, x, f, &iflag);
921 solverData->numberOfFunctionEvaluations++;
922
923 return 0;
924}
925
926/*! \fn wrapper_fvec_constraints for the residual Function
927 * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data)
928 *
929 * \author ptaeuber
930 *
931 */
932int wrapper_fvec_constraints(DATA_HOMOTOPY* solverData, double* x, double* f)
933{
934 void *dataAndThreadData[2] = {solverData->data, solverData->threadData};
935 int iflag = 0;
936 int retVal;
937
938 /*TODO: change input to residualFunc from data to systemData */
939 retVal = (solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].residualFuncConstraints(dataAndThreadData, x, f, &iflag);
940 solverData->numberOfFunctionEvaluations++;
941
942 return retVal;
943}
944
945/*! \fn wrapper_fvec_der for the residual Function
946 * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data)
947 *
948 * \author bbachmann
949 *
950 */
951static int wrapper_fvec_der(DATA_HOMOTOPY* solverData, double* x, double* fJac)
952{
953 int i;
954 int jacobianIndex = (&(solverData->data->simulationInfo->nonlinearSystemData[solverData->sysNumber]))->jacobianIndex;
955 NONLINEAR_SYSTEM_DATA* nonlinsys = &(solverData->data->simulationInfo->nonlinearSystemData[solverData->sysNumber]);
956
957 /* performance measurement */
958 rt_ext_tp_tick(&nonlinsys->jacobianTimeClock);
959
960 /* calculate jacobian */
961 if(jacobianIndex != -1)
962 {
963 /* !!!!!!!!!!! Be sure that actual x is used !!!!!!!!!!! */
964 getAnalyticalJacobianHomotopy(solverData, fJac);
965 }
966 else
967 {
968 getNumericalJacobianHomotopy(solverData, x, fJac);
969 }
970
971 if(ACTIVE_STREAM(LOG_NLS_JAC_TEST)(useStream[LOG_NLS_JAC_TEST]))
972 {
973 int n = solverData->n;
974 /* debugMatrixDouble(LOG_NLS_JAC_TEST,"analytical jacobian:",fJac, n, n+1); */
975 getNumericalJacobianHomotopy(solverData, x, solverData->debug_fJac);
976 /* debugMatrixDouble(LOG_NLS_JAC_TEST,"numerical jacobian:",solverData->debug_fJac, n, n+1); */
977 matDiffBB(n, fJac, solverData->debug_fJac, solverData->debug_fJac);
978 /* debugMatrixDouble(LOG_NLS_JAC_TEST,"Difference of jacobians:",solverData->debug_fJac, n, n+1); */
979 debugDouble(LOG_NLS_JAC_TEST,"error between analytical and numerical jacobian = ", vecMaxNorm(n*n, solverData->debug_fJac));
980 vecDivScaling(n*(n+1), solverData->debug_fJac , fJac, solverData->debug_fJac);
981 debugDouble(LOG_NLS_JAC_TEST,"relative error between analytical and numerical jacobian = ", vecMaxNorm(n*n, solverData->debug_fJac));
982 messageClose(LOG_NLS_JAC_TEST);
983 }
984 /* performance measurement and statistics */
985 nonlinsys->jacobianTime += rt_ext_tp_tock(&(nonlinsys->jacobianTimeClock));
986 nonlinsys->numberOfJEval++;
987
988 return 0;
989}
990
991/*! \fn wrapper_fvec_homotopy_newton for the residual Function
992 *
993 * \author bbachmann
994 *
995 */
996static int wrapper_fvec_homotopy_newton(DATA_HOMOTOPY* solverData, double* x, double* h)
997{
998 int i;
999 int n = solverData->n;
1000
1001 /* Newton homotopy */
1002 wrapper_fvec(solverData, x, solverData->f1);
1003 vecAddScal(solverData->n, solverData->f1, solverData->fx0, - (1-x[n]), h);
1004
1005 return 0;
1006}
1007
1008/*! \fn wrapper_fvec_homotopy_newton_der for the residual Function
1009 *
1010 * \author bbachmann
1011 *
1012 */
1013static int wrapper_fvec_homotopy_newton_der(DATA_HOMOTOPY* solverData, double* x, double* hJac)
1014{
1015 int i, j;
1016 int n = solverData->n;
1017
1018 /* Newton homotopy */
1019 wrapper_fvec_der(solverData, x, hJac);
1020
1021 /* add f(x0) as the last column of the Jacobian*/
1022 vecCopy(n, solverData->fx0, hJac + n*n);
1023
1024 return 0;
1025}
1026
1027/*! \fn wrapper_fvec_homotopy_fixpoint for the residual Function
1028 *
1029 * \author bbachmann
1030 *
1031 */
1032static int wrapper_fvec_homotopy_fixpoint(DATA_HOMOTOPY* solverData, double* x, double* h)
1033{
1034 int i;
1035 int n = solverData->n;
1036
1037 /* Fixpoint homotopy */
1038 wrapper_fvec(solverData, x, solverData->f1);
1039 for (i=0; i<n; i++){
1040 h[i] = x[n]*solverData->f1[i] + (1-x[n]) * (x[i]-solverData->x0[i]);
1041 }
1042
1043 return 0;
1044}
1045
1046/*! \fn wrapper_fvec_homotopy_fixpoint_der for the residual Function
1047 *
1048 * \author bbachmann
1049 *
1050 */
1051static int wrapper_fvec_homotopy_fixpoint_der(DATA_HOMOTOPY* solverData, double* x, double* hJac)
1052{
1053 int i, j;
1054 int n = solverData->n;
1055
1056 /* Fixpoint homotopy */
1057 wrapper_fvec_der(solverData, x, hJac);
1058 for (i=0; i<n; i++){
1059 for (j=0; j<n; j++) {
1060 hJac[i+ j * n] = x[n]*hJac[i+ j * n];
1061 }
1062 hJac[i+ i * n] = hJac[i+ i * n] + (1-x[n]);
1063 hJac[i+ n * n] = solverData->f1[i]-(x[i] - solverData->x0[i]);
1064 }
1065 return 0;
1066}
1067
1068/*! \fn getIndicesOfPivotElement for calculating pivot element
1069 *
1070 * \author bbachmann
1071 *
1072 */
1073 void getIndicesOfPivotElement(int *n, int *m, int *l, double* A, int *indRow, int *indCol, int *pRow, int *pCol, double *absMax)
1074{
1075 int i, j;
1076
1077 *absMax = fabs(A[indRow[*l] + indCol[*l]* *n]);
1078 *pCol = *l;
1079 *pRow = *l;
1080 for (i = *l; i < *n; i++) {
1081 for (j = *l; j < *m; j++) {
1082 if (fabs(A[indRow[i] + indCol[j]* *n]) > *absMax) {
1083 *absMax = fabs(A[indRow[i] + indCol[j]* *n]);
1084 *pCol = j;
1085 *pRow = i;
1086 }
1087 }
1088 }
1089}
1090
1091
1092/*! \fn solveSystemWithTotalPivotSearch for solution of overdetermined linear system
1093 * used for the homotopy solver, for calculating the direction
1094 * used for the newton solver, for calculating the Newton step
1095 *
1096 * \author bbachmann
1097 *
1098 */
1099int solveSystemWithTotalPivotSearch(int n, double* x, double* A, int* indRow, int* indCol, int *pos, int *rank, int casualTearingSet)
1100{
1101 int i, k, j, m=n+1, nPivot=n;
1102 int pCol, pRow;
1103 double hValue;
1104 double hInt;
1105 double absMax, detJac;
1106 int returnValue = 0;
1107
1108 debugMatrixDouble(LOG_NLS_JAC,"Linear System Matrix [Jac res]:",A, n, m);
1109 debugVectorDouble(LOG_NLS_JAC,"vector b:", A+n*n, n);
1110
1111 /* assume full rank of matrix [n x (n+1)] */
1112 *rank = n;
1113
1114 for (i=0; i<n; i++) {
1115 indRow[i] = i;
1116 }
1117 for (i=0; i<m; i++) {
1118 indCol[i] = i;
1119 }
1120 if (*pos>=0) {
1121 indCol[n] = *pos;
1122 indCol[*pos] = n;
1123 } else {
1124 nPivot = n+1;
1125 }
1126
1127 for (i = 0; i < n; i++) {
1128 getIndicesOfPivotElement(&n, &nPivot, &i, A, indRow, indCol, &pRow, &pCol, &absMax);
1129 if (absMax<DBL_EPSILON2.2204460492503131e-16) {
1130 *rank = i;
1131 warningStreamPrint(LOG_NLS_V, 0, "Matrix singular!");
1132 debugInt(LOG_NLS_V,"rank = ", *rank);
1133 debugInt(LOG_NLS_V,"position = ", *pos);
1134 break;
1135 }
1136 /* swap row indices */
1137 if (pRow!=i) {
1138 hInt = indRow[i];
1139 indRow[i] = indRow[pRow];
1140 indRow[pRow] = hInt;
1141 }
1142 /* swap column indices */
1143 if (pCol!=i) {
1144 hInt = indCol[i];
1145 indCol[i] = indCol[pCol];
1146 indCol[pCol] = hInt;
1147 }
1148
1149 /* Gauss elimination of row indRow[i] */
1150 for (k=i+1; k<n; k++) {
1151 hValue = -A[indRow[k] + indCol[i]*n]/A[indRow[i] + indCol[i]*n];
1152 for (j=i+1; j<m; j++) {
1153 A[indRow[k] + indCol[j]*n] = A[indRow[k] + indCol[j]*n] + hValue*A[indRow[i] + indCol[j]*n];
1154 }
1155 A[indRow[k] + indCol[i]*n] = 0;
1156 }
1157 }
1158
1159 for (detJac=1.0,k=0; k<n; k++) detJac *= A[indRow[k] + indCol[k]*n];
1160
1161 debugMatrixPermutedDouble(LOG_NLS_JAC,"Linear System Matrix [Jac res] after decomposition",A, n, m, indRow, indCol);
1162 debugDouble(LOG_NLS_JAC,"Determinant = ", detJac);
1163 if (isnan(detJac)(sizeof ((detJac)) == sizeof (float) ? __isnanf (detJac) : sizeof
((detJac)) == sizeof (double) ? __isnan (detJac) : __isnanl (
detJac))
){
1164 warningStreamPrint(LOG_NLS_V, 0, "Jacobian determinant is NaN.");
1165 return -1;
1166 }
1167 else if (fabs(detJac) < 1e-9 && casualTearingSet)
1168 {
1169 debugString(LOG_DT, "The determinant of the casual tearing set is vanishing, let's fail if this is not the solution...");
1170 returnValue = 1;
1171 }
1172
1173 /* Solve even singular matrices !!! */
1174 for (i=n-1;i>=0; i--) {
1175 if (i>=*rank) {
1176 /* this criteria should be evaluated and may be improved in future */
1177 if (fabs(A[indRow[i] + indCol[n]*n])>1e-6) {
1178 warningStreamPrint(LOG_NLS_V, 0, "under-determined linear system not solvable!");
1179 return -1;
1180 } else {
1181 x[indCol[i]] = 0.0;
1182 }
1183 } else {
1184 x[indCol[i]] = -A[indRow[i] + indCol[n]*n];
1185 for (j=n-1; j>i; j--) {
1186 x[indCol[i]] = x[indCol[i]] - A[indRow[i] + indCol[j]*n]*x[indCol[j]];
1187 }
1188 x[indCol[i]]=x[indCol[i]]/A[indRow[i] + indCol[i]*n];
1189 }
1190 }
1191 x[indCol[n]]=1.0;
1192 debugVectorInt(LOG_NLS_V,"indRow:", indRow, n);
1193 debugVectorInt(LOG_NLS_V,"indCol:", indCol, n+1);
1194 debugVectorDouble(LOG_NLS_V,"vector x (solution):", x, n+1);
1195
1196 /* Return position of largest value (1.0) */
1197 if (*pos<0) {
1198 *pos=indCol[n];
1199 debugInt(LOG_NLS_V,"position of largest value = ", *pos);
1200 }
1201
1202 return returnValue;
1203}
1204
1205
1206/*! \fn linearSolverWrapper
1207 */
1208int linearSolverWrapper(int n, double* x, double* A, int* indRow, int* indCol, int *pos, int *rank, int method, int casualTearingSet)
1209{
1210 /* First try to use lapack and if it fails then
1211 * use solveSystemWithTotalPivotSearch */
1212 int returnValue = -1;
1213 int solverinfo;
1214 int nrhs = 1;
1215 int lda = n;
1216 int k;
1217 double detJac;
1218
1219 debugMatrixDouble(LOG_NLS_JAC,"Linear System Matrix [Jac res]:", A, n, n+1);
1220 debugVectorDouble(LOG_NLS_JAC,"vector b:", x, n);
1221
1222 switch(method){
1223 case NLS_LS_TOTALPIVOT:
1224
1225 solverinfo = solveSystemWithTotalPivotSearch(n, x, A, indRow, indCol, pos, rank, casualTearingSet);
1226 /* in case of failing */
1227 if (solverinfo == -1)
1228 {
1229 /* debug information */
1230 debugString(LOG_NLS_V, "Linear total pivot solver failed!!!");
1231 debugString(LOG_NLS_V, "******************************************************");
1232 }
1233 else if (solverinfo == 1)
1234 {
1235 returnValue = 1;
1236 }
1237 else
1238 {
1239 returnValue = 0;
1240 }
1241 break;
1242 case NLS_LS_LAPACK:
1243 /* Solve system with lapack */
1244 dgesv_((int*) &n,
1245 (int*) &nrhs,
1246 A,
1247 (int*) &lda,
1248 indRow,
1249 x,
1250 (int*) &n,
1251 &solverinfo);
1252
1253 for (detJac=1.0, k=0; k<n; k++) detJac *= A[k + k*n];
1254
1255 debugMatrixDouble(LOG_NLS_JAC,"Linear system matrix [Jac res] after decomposition:", A, n, n+1);
1256 debugDouble(LOG_NLS_JAC,"Determinant = ", detJac);
1257
1258 /* in case of failing */
1259 if (solverinfo != 0)
1260 {
1261 /* debug information */
1262 debugString(LOG_NLS_V, "Linear lapack solver failed!!!");
1263 debugString(LOG_NLS_V, "******************************************************");
1264 }
1265 else if (fabs(detJac) < 1e-9 && casualTearingSet)
1266 {
1267 debugString(LOG_DT, "The determinant of the casual tearing set is vanishing, let's fail if this is not the solution...");
1268 returnValue = 1;
1269 }
1270 else
1271 {
1272 vecScalarMult(n, x, -1, x);
1273 returnValue = 0;
1274 }
1275 break;
1276 default:
1277 throwStreamPrint(0, "Non-Linear solver try to run with a unknown linear solver (%d).", method);
1278 }
1279
1280 /* Debugging error of linear system */
1281 if(ACTIVE_STREAM(LOG_NLS_JAC)(useStream[LOG_NLS_JAC]))
1282 {
1283 double* res = (double*) calloc(n,sizeof(double));
1284 debugVectorDouble(LOG_NLS_JAC,"solution:", x, n);
1285 matVecMult(n, n, A, x, res);
1286 debugVectorDouble(LOG_NLS_JAC,"test solution:", res, n);
1287 debugDouble(LOG_NLS_JAC,"error of linear system = ", vec2Norm(n, res));
1288 free(res);
1289 messageClose(LOG_NLS_JAC);
1290 }
1291
1292 return returnValue;
1293}
1294
1295
1296/*! \fn solve system with damped Newton-Raphson
1297 *
1298 * \author bbachmann
1299 *
1300 */
1301static int newtonAlgorithm(DATA_HOMOTOPY* solverData, double* x)
1302{
1303 int numberOfIterations = 0 ,i, j, n=solverData->n, m=solverData->m;
1304 int pos = solverData->n, rank;
1305 double error_f_sqrd, error_f1_sqrd, error_f2_sqrd, error_f_sqrd_scaled, error_f1_sqrd_scaled;
1306 double delta_x_sqrd, delta_x_sqrd_scaled, grad_f, grad_f_scaled;
1307 int numberOfSmallSteps = 0;
1308 double error_f_old = 1e100;
1309 int countNegativeSteps = 0;
1310 double lambda;
1311 double lambda1, lambda2;
1312 double lambdaMin = 1e-4;
1313 double a2, a3, rhs1, rhs2, D;
1314 double alpha = 1e-1;
1315 int firstrun;
1316 int constraintViolated;
1317 int solverinfo = 0;
1318
1319 int assert = 1;
1320 threadData_t *threadData = solverData->threadData;
1321 NONLINEAR_SYSTEM_DATA* nonlinsys = &(solverData->data->simulationInfo->nonlinearSystemData[solverData->data->simulationInfo->currentNonlinearSystemIndex]);
1322 int linearSolverMethod = solverData->data->simulationInfo->nlsLinearSolver;
1323
1324 /* debug information */
1325 debugString(LOG_NLS_V, "******************************************************");
1326 debugInt(LOG_NLS_V, "NEWTON SOLVER STARTED! equation number: ",solverData->eqSystemNumber);
1327 debugInt(LOG_NLS_V, "maximum number of function evaluation: ", solverData->maxNumberOfIterations);
1328 printUnknowns(LOG_NLS_V, solverData);
1329
1330 /* set default solver message */
1331 solverData->info = 0;
1332
1333 /* calculated error of function values */
1334 error_f_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1335 vecDivScaling(solverData->n, solverData->f1, solverData->resScaling, solverData->fvecScaled);
1336 error_f_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->fvecScaled);
1337
1338 while(1)
1339 {
1340 numberOfIterations++;
1341 /* debug information */
1342 debugInt(LOG_NLS_V, "Iteration:", numberOfIterations);
1343
1344 /* solve jacobian and function value (both stored in hJac, last column is fvec), side effects: jacobian matrix is changed */
1345 if (numberOfIterations>1)
1346 solverinfo = linearSolverWrapper(solverData->n, solverData->dy0, solverData->fJac, solverData->indRow, solverData->indCol, &pos, &rank, linearSolverMethod, solverData->casualTearingSet);
1347
1348 if (solverinfo == -1)
1349 {
1350 /* report solver abortion */
1351 solverData->info=-1;
1352 /* debug information */
1353 debugString(LOG_NLS_V, "NEWTON SOLVER DID ---NOT--- CONVERGE TO A SOLUTION!!!");
1354 debugString(LOG_NLS_V, "******************************************************");
1355 assert = 0;
Value stored to 'assert' is never read
1356 break;
1357 }
1358 else
1359 {
1360 /* Scaling back to original variables */
1361 vecMultScaling(solverData->m, solverData->dy0, solverData->xScaling, solverData->dy0);
1362 /* try full Newton step */
1363 vecAdd(solverData->n, x, solverData->dy0, solverData->x1);
1364 printNewtonStep(LOG_NLS_V, solverData);
1365
1366 /* Damping strategy, performance is very sensitive on the value of lambda */
1367 lambda1 = 1.0;
1368 assert = 1;
1369 firstrun = 1;
1370 while (assert && (lambda1 > lambdaMin))
1371 {
1372 if (!firstrun){
1373 lambda1 *= 0.655;
1374 vecAddScal(solverData->n, x, solverData->dy0, lambda1, solverData->x1);
1375 assert = 1;
1376 }
1377#ifndef OMC_EMCC
1378 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1379#endif
1380 if (solverData->casualTearingSet){
1381 constraintViolated = solverData->f_con(solverData, solverData->x1, solverData->f1);
1382 if (constraintViolated){
1383 lambda1 = lambdaMin-1;
1384 break;
1385 }
1386 }
1387 else
1388 solverData->f(solverData, solverData->x1, solverData->f1);
1389
1390 assert = 0;
1391#ifndef OMC_EMCC
1392 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1393#endif
1394 firstrun = 0;
1395 if (assert){
1396 debugDouble(LOG_NLS_V,"Assert of Newton step: lambda1 =", lambda1);
1397 }
1398 }
1399
1400 if (lambda1 < lambdaMin)
1401 {
1402 debugDouble(LOG_NLS_V,"UPS! MUST HANDLE A PROBLEM (Newton method), time : ", solverData->timeValue);
1403 solverData->info = -1;
1404 break;
1405 }
1406
1407 /* Damping (see Numerical Recipes) */
1408 /* calculate gradient of quadratic function for damping strategy */
1409 grad_f = -2.0*error_f_sqrd;
1410 grad_f_scaled = -2.0*error_f_sqrd_scaled;
1411 error_f1_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1412 vecDivScaling(solverData->n, solverData->f1, solverData->resScaling, solverData->f2);
1413 error_f1_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->f2);
1414 debugDouble(LOG_NLS_V,"Need to damp, grad_f = ", grad_f);
1415 debugDouble(LOG_NLS_V,"Need to damp, error_f = ", sqrt(error_f_sqrd));
1416 debugDouble(LOG_NLS_V,"Need to damp this!! lambda1 = ", lambda1);
1417 debugDouble(LOG_NLS_V,"Need to damp, error_f1 = ", sqrt(error_f1_sqrd));
1418 debugDouble(LOG_NLS_V,"Need to damp, forced error = ", error_f_sqrd + alpha*lambda1*grad_f);
1419 if ((error_f1_sqrd > error_f_sqrd + alpha*lambda1*grad_f)
1420 && (error_f1_sqrd_scaled > error_f_sqrd_scaled + alpha*lambda1*grad_f_scaled)
1421 && (error_f_sqrd > 1e-12) && (error_f_sqrd_scaled > 1e-12))
1422 {
1423 lambda2 = fmax(-lambda1*lambda1*grad_f/(2*(error_f1_sqrd-error_f_sqrd-lambda1*grad_f)),lambdaMin);
1424 debugDouble(LOG_NLS_V,"Need to damp this!! lambda2 = ", lambda2);
1425 vecAddScal(solverData->n, x, solverData->dy0, lambda2, solverData->x1);
1426 assert= 1;
1427#ifndef OMC_EMCC
1428 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1429#endif
1430 if (solverData->casualTearingSet){
1431 constraintViolated = solverData->f_con(solverData, solverData->x1, solverData->f1);
1432 if (constraintViolated){
1433 solverData->info = -1;
1434 break;
1435 }
1436 }
1437 else
1438 solverData->f(solverData, solverData->x1, solverData->f1);
1439
1440 error_f2_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1441 debugDouble(LOG_NLS_V,"Need to damp, error_f2 = ", sqrt(error_f2_sqrd));
1442 assert = 0;
1443#ifndef OMC_EMCC
1444 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1445#endif
1446 if (assert)
1447 {
1448 debugDouble(LOG_NLS_V,"UPS! MUST HANDLE A PROBLEM (Newton method), time : ", solverData->timeValue);
1449 solverData->info = -1;
1450 break;
1451 }
1452 if ((error_f1_sqrd > error_f_sqrd + alpha*lambda2*grad_f) && (error_f_sqrd > 1e-12) && (error_f_sqrd_scaled > 1e-12))
1453 {
1454 rhs1 = error_f1_sqrd - grad_f*lambda1 - error_f_sqrd;
1455 rhs2 = error_f2_sqrd - grad_f*lambda2 - error_f_sqrd;
1456 a3 = (rhs1/(lambda1*lambda1) - rhs2/(lambda2*lambda2))/(lambda1 - lambda2);
1457 a2 = (-lambda2*rhs1/(lambda1*lambda1) + lambda1*rhs2/(lambda2*lambda2))/(lambda1 - lambda2);
1458 if (a3==0.0)
1459 lambda = -grad_f/(2.0*a2);
1460 else
1461 {
1462 D = a2*a2 - 3.0*a3*grad_f;
1463 if (D <= 0.0)
1464 lambda = 0.5*lambda1;
1465 else
1466 if (a2 <= 0.0)
1467 lambda = (-a2+sqrt(D))/(3.0*a3);
1468 else
1469 lambda = -grad_f/(a2+sqrt(D));
1470 }
1471 lambda = fmax(lambda, lambdaMin);
1472 debugDouble(LOG_NLS_V,"Need to damp this!! lambda = ", lambda);
1473 vecAddScal(solverData->n, x, solverData->dy0, lambda, solverData->x1);
1474 assert= 1;
1475#ifndef OMC_EMCC
1476 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1477#endif
1478 if (solverData->casualTearingSet){
1479 constraintViolated = solverData->f_con(solverData, solverData->x1, solverData->f1);
1480 if (constraintViolated){
1481 solverData->info = -1;
1482 break;
1483 }
1484 }
1485 else
1486 solverData->f(solverData, solverData->x1, solverData->f1);
1487
1488 error_f1_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1489 debugDouble(LOG_NLS_V,"Need to damp, error_f1 = ", sqrt(error_f1_sqrd));
1490 assert = 0;
1491#ifndef OMC_EMCC
1492 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1493#endif
1494 if (assert)
1495 {
1496 debugDouble(LOG_NLS_V,"UPS! MUST HANDLE A PROBLEM (Newton method), time : ", solverData->timeValue);
1497 solverData->info = -1;
1498 break;
1499 }
1500 }
1501 }else{
1502 lambda = lambda1;
1503 }
1504 }
1505 /* updating x, fvec, error_f_sqrd */
1506 /* event. swapPointer(&x, &(solverData->x1)); */
1507 vecCopy(solverData->n, solverData->x1, x);
1508
1509 /* Calculate different error measurements */
1510 vecDivScaling(solverData->n, solverData->f1, solverData->resScaling, solverData->fvecScaled);
1511 debugVectorDouble(LOG_NLS_V,"function values:",solverData->f1, n);
1512 debugVectorDouble(LOG_NLS_V,"scaled function values:",solverData->fvecScaled, n);
1513
1514 vecDivScaling(solverData->n, solverData->dy0, solverData->xScaling, solverData->dxScaled);
1515 delta_x_sqrd = vec2NormSqrd(solverData->n, solverData->dy0);
1516 delta_x_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->dxScaled);
1517 error_f_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
1518 error_f_sqrd_scaled = vec2NormSqrd(solverData->n, solverData->fvecScaled);
1519
1520
1521 /* debug information */
1522 debugString(LOG_NLS_V, "error measurements:");
1523 debugDouble(LOG_NLS_V, "delta_x =", sqrt(delta_x_sqrd));
1524 debugDouble(LOG_NLS_V, "delta_x_scaled =", sqrt(delta_x_sqrd_scaled));
1525 debugDouble(LOG_NLS_V, "newtonXTol =", sqrt(solverData->xtol_sqrd));
1526 debugDouble(LOG_NLS_V, "error_f =", sqrt(error_f_sqrd));
1527 debugDouble(LOG_NLS_V, "error_f_scaled =", sqrt(error_f_sqrd_scaled));
1528 debugDouble(LOG_NLS_V, "newtonFTol =", sqrt(solverData->ftol_sqrd));
1529
1530 countNegativeSteps += (error_f_sqrd > 10*error_f_old);
1531 error_f_old = error_f_sqrd;
1532
1533#if !defined(OMC_MINIMAL_RUNTIME)
1534 if (solverData->data->simulationInfo->nlsCsvInfomation){
1535 print_csvLineIterStats(((struct csvStats*) nonlinsys->csvData)->iterStats,
1536 nonlinsys->size,
1537 nonlinsys->numberOfCall+1,
1538 numberOfIterations,
1539 solverData->x,
1540 solverData->f1,
1541 delta_x_sqrd,
1542 delta_x_sqrd_scaled,
1543 error_f_sqrd,
1544 error_f_sqrd_scaled,
1545 lambda
1546 );
1547 }
1548#endif
1549 if (countNegativeSteps > 20)
1550 {
1551 debugInt(LOG_NLS_V,"UPS! Something happened, NegativeSteps = ", countNegativeSteps);
1552 solverData->info = -1;
1553 break;
1554 }
1555
1556 /* solution found */
1557 if (((error_f_sqrd < solverData->ftol_sqrd) || (error_f_sqrd_scaled < solverData->ftol_sqrd)) && ((delta_x_sqrd_scaled < solverData->xtol_sqrd) || (delta_x_sqrd < solverData->xtol_sqrd)))
1558 {
1559 solverData->info = 1;
1560
1561 /* debug information */
1562 debugString(LOG_NLS_V, "NEWTON SOLVER DID CONVERGE TO A SOLUTION!!!");
1563 printUnknowns(LOG_NLS_V, solverData);
1564 debugString(LOG_NLS_V, "******************************************************");
1565
1566 /* update statistics */
1567 solverData->numberOfIterations += numberOfIterations;
1568 solverData->error_f_sqrd = error_f_sqrd;
1569
1570 break;
1571 }
1572 else if (solverinfo == 1){
1573 solverData->info = -1;
1574 debugString(LOG_DT, "It is not the solution.");
1575 break;
1576 }
1577
1578 /* check if maximum iteration is reached */
1579 if (numberOfIterations > solverData->maxNumberOfIterations)
1580 {
1581 solverData->info = -1;
1582 warningStreamPrint(LOG_NLS_V, 0, "Warning: maximal number of iteration reached but no root found");
1583 /* debug information */
1584 debugString(LOG_NLS_V, "NEWTON SOLVER DID ---NOT--- CONVERGE TO A SOLUTION!!!");
1585 debugString(LOG_NLS_V, "******************************************************");
1586
1587 /* update statistics */
1588 solverData->numberOfIterations += numberOfIterations;
1589 break;
1590 }
1591
1592 numberOfSmallSteps += (delta_x_sqrd < solverData->xtol_sqrd*1e4) || (delta_x_sqrd_scaled < solverData->xtol_sqrd*1e4);
1593 /* check changes in unknown vector */
1594 if ((delta_x_sqrd < solverData->xtol_sqrd) || (delta_x_sqrd_scaled < solverData->xtol_sqrd) || (numberOfSmallSteps > 20))
1595 {
1596 if ((error_f_sqrd < solverData->ftol_sqrd*1e6) || (error_f_sqrd_scaled < solverData->ftol_sqrd*1e6))
1597 {
1598 solverData->info = 1;
1599
1600 /* debug information */
1601 debugString(LOG_NLS_V, "NEWTON SOLVER DID CONVERGE TO A SOLUTION WITH LESS ACCURACY!!!");
1602 printUnknowns(LOG_NLS_V, solverData);
1603 debugString(LOG_NLS_V, "******************************************************");
1604 solverData->error_f_sqrd = 0;
1605
1606 } else
1607 {
1608 solverData->info = -1;
1609 debugString(LOG_NLS_V, "Warning: newton solver gets stuck!!!");
1610 /* debug information */
1611 debugString(LOG_NLS_V, "NEWTON SOLVER DID ---NOT--- CONVERGE TO A SOLUTION!!!");
1612 debugString(LOG_NLS_V, "******************************************************");
1613 }
1614 /* update statistics */
1615 solverData->numberOfIterations += numberOfIterations;
1616 break;
1617 }
1618 assert = 1;
1619#ifndef OMC_EMCC
1620 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1621#endif
1622 /* calculate jacobian and function values (both stored in fJac, last column is fvec) */
1623 solverData->fJac_f(solverData, x, solverData->fJac);
1624 assert = 0;
1625#ifndef OMC_EMCC
1626 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1627#endif
1628 if (assert)
1629 {
1630 /* report solver abortion */
1631 solverData->info=-1;
1632 debugString(LOG_NLS_V,"UPS! assert when calculating Jacobian!!!");
1633 break;
1634 }
1635 vecCopy(n, solverData->f1, solverData->fJac + n*n);
1636 /* calculate scaling factor of residuals */
1637 matVecMultAbsBB(solverData->n, solverData->fJac, solverData->ones, solverData->resScaling);
1638 debugVectorDouble(LOG_NLS_JAC, "residuum scaling:", solverData->resScaling, solverData->n);
1639 scaleMatrixRows(solverData->n, solverData->m, solverData->fJac);
1640 vecCopy(n, solverData->fJac + n*n, solverData->dy0);
1641 }
1642 return 0;
1643}
1644
1645/*! \fn solve system with homotopy method
1646 *
1647 * \author bbachmann
1648 */
1649static int homotopyAlgorithm(DATA_HOMOTOPY* solverData, double *x)
1650{
1651 int i, j;
1652 double error_h, error_h_scaled, delta_x;
1653 double vecScalarProduct;
1654
1655 int pos, rank;
1656 int iter = 0;
1657 int maxiter = homMaxNewtonSteps;
1658 int maxTries = homMaxTries;
1659 int numSteps = 0;
1660 int stepAccept = 0;
1661 int correctorStrategy = homBacktraceStrategy; /* 1: go back to the path by fixing one coordinate, 2: go back to the path in an orthogonal direction to the tangent vector */
1662 double bend = 0;
1663 double tau = homTauStart, tauMax = homTauMax, tauMin = homTauMin, hEps = homHEps, adaptBend = homAdaptBend;
1664 double tauDecreasingFactor = homTauDecreasingFactor, tauDecreasingFactorPredictor = homTauDecreasingFactorPredictor;
1665 double tauIncreasingFactor = homTauIncreasingFactor, tauIncreasingThreshold = homTauIncreasingThreshold;
1666 double preTau;
1667 int m = solverData->m;
1668 int n = solverData->n;
1669 int initialStep = 1;
1670 int maxLambdaSteps = homMaxLambdaSteps ? homMaxLambdaSteps : solverData->maxNumberOfIterations;
1671
1672 int assert = 1;
1673 threadData_t *threadData = solverData->threadData;
1674 FILE *pFile = NULL((void*)0);
1675 char buffer[4096];
1676
1677#if !defined(OMC_NO_FILESYSTEM)
1678 const char sep[] = ",";
1679 if(solverData->initHomotopy && ACTIVE_STREAM(LOG_INIT_HOMOTOPY)(useStream[LOG_INIT_HOMOTOPY]))
1680 {
1681 sprintf(buffer, "%s_nonlinsys%d_adaptive_%s_homotopy_%s.csv", solverData->data->modelData->modelFilePrefix, solverData->sysNumber, solverData->data->callback->useHomotopy == 2 ? "global" : "local", solverData->startDirection > 0 ? "pos" : "neg");
1682 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "The homotopy path will be exported to %s.", buffer);
1683 pFile = omc_fopen(buffer, "wt");
1684 fprintf(pFile, "\"sep=%s\"\n%s", sep, "\"lambda\"");
1685 for(i=0; i<n; ++i)
1686 fprintf(pFile, "%s\"%s\"", sep, modelInfoGetEquation(&solverData->data->modelData->modelDataXml,solverData->eqSystemNumber).vars[i]);
1687 fprintf(pFile, "\n");
1688 fprintf(pFile, "0.0");
1689 for(i=0; i<n; ++i)
1690 fprintf(pFile, "%s%.16g", sep, x[i]);
1691 fprintf(pFile, "\n");
1692 }
1693#endif
1694
1695 /* Initialize vector dy2 using chosen startDirection */
1696 /* set start vector, lambda = 0.0 */
1697 vecCopy(solverData->n, x, solverData->y0);
1698 solverData->y0[solverData->n] = 0.0;
1699
1700 vecConst(solverData->n, 0.0, solverData->dy2);
1701 solverData->dy2[solverData->n]= solverData->startDirection;
1702 printHomotopyUnknowns(LOG_NLS_V, solverData);
1703 assert = 1;
1704#ifndef OMC_EMCC
1705 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1706#endif
1707 solverData->h_function(solverData, solverData->y0, solverData->hvec);
1708 assert = 0;
1709#ifndef OMC_EMCC
1710 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1711#endif
1712 /* start iteration; stop, if lambda = solverData->y0[solverData->n] == 1 */
1713 while (solverData->y0[solverData->n]<1)
1714 {
1715 if (solverData->initHomotopy)
1716 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
1717 else
1718 infoStreamPrint(LOG_NLS_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
1719 /* Break loop, iff algorithm gets stuck or lambda accelerates to the wrong direction */
1720 if (iter>=maxTries)
1721 {
1722 if (solverData->initHomotopy) {
1723 if (preTau == tau)
1724 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nNo solution for current step size tau found and tau cannot be decreased any further.\nYou can set the minimum step size tau with:\n\t-homTauMin=<value>\nYou can also try to allow more newton steps in the corrector step with:\n\t-homMaxNewtonSteps=<value>\nor change the tolerance for the solution with:\n\t-homHEps=<value>\nYou can also try to use another backtrace stategy in the corrector step with:\n\t-homBacktraceStrategy=<fix|orthogonal>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
1725 else
1726 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe maximum number of tries for one lambda is reached (%d).\nYou can change the number of tries with:\n\t-homMaxTries=<value>\nYou can also try to allow more newton steps in the corrector step with:\n\t-homMaxNewtonSteps=<value>\nor change the tolerance for the solution with:\n\t-homHEps=<value>\nYou can also try to use another backtrace stategy in the corrector step with:\n\t-homBacktraceStrategy=<fix|orthogonal>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", iter);
1727 }
1728 else
1729 debugInt(LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge: iter = ", iter);
1730 debugString(LOG_NLS_HOMOTOPY, "======================================================");
1731 return -1;
1732 }
1733 if (solverData->y0[solverData->n]<(-1))
1734 {
1735 if (solverData->initHomotopy)
1736 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nlambda is smaller than -1: lambda=%g\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", solverData->y0[solverData->n]);
1737 else
1738 debugDouble(LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge: lambda = ", solverData->y0[solverData->n]);
1739 debugString(LOG_NLS_HOMOTOPY, "======================================================");
1740 return -1;
1741 }
1742 if (numSteps >= maxLambdaSteps)
1743 {
1744 if (solverData->initHomotopy)
1745 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe maximum number of lambda steps is reached (%d).\nYou can change the maximum number of lambda steps with:\n\t-homMaxLambdaSteps=<value>\nYou can also try to influence the step size tau with the following flags:\n\t-homTauDecFac=<value>\n\t-homTauDecFacPredictor=<value>\n\t-homTauIncFac=<value>\n\t-homTauIncThreshold=<value>\n\t-homTauMax=<value>\n\t-homTauMin=<value>\n\t-homTauStart=<value>\nor you can also set the threshold for accepting the current bending with:\n\t-homAdaptBend=<value>\nYou can also try to use another backtrace stategy in the corrector step with:\n\t-homBacktraceStrategy=<fix|orthogonal>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", maxLambdaSteps);
1746 else
1747 debugInt(LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge: numSteps = ", numSteps);
1748 debugString(LOG_NLS_HOMOTOPY, "======================================================");
1749 return -1;
1750 }
1751
1752 stepAccept = 0;
1753
1754 /****************************************************************************
1755 * Predictor step: Calculation of tangent vector! *
1756 ****************************************************************************/
1757 /* If a step succeeded, calculate the homotopy function and corresponding jacobian */
1758 if (iter==0)
1759 {
1760 /* Handle asserts of function calls, mainly necessary for fluid stuff */
1761 assert = 1;
1762#ifndef OMC_EMCC
1763 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1764#endif
1765 solverData->hJac_dh(solverData, solverData->y0, solverData->hJac);
1766 debugMatrixDouble(LOG_NLS_JAC,"Jacobian hJac:",solverData->hJac, solverData->n, solverData->n+1);
1767 scaleMatrixRows(solverData->n, solverData->m, solverData->hJac);
1768 debugMatrixDouble(LOG_NLS_JAC,"Jacobian hJac after scaling:",solverData->hJac, solverData->n, solverData->n+1);
1769 assert = 0;
1770 pos = -1; /* stable solution algorithm for solving a generalized over-determined linear system */
1771#ifndef OMC_EMCC
1772 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1773#endif
1774
1775 if (assert || (solveSystemWithTotalPivotSearch(solverData->n, solverData->dy0, solverData->hJac, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1))
1776 {
1777 /* report solver abortion */
1778 solverData->info=-1;
1779 /* debug information */
1780 if (assert) {
1781 if (solverData->initHomotopy)
1782 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nIt was not possible to calculate the jacobian.\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
1783 else {
1784 debugString(LOG_NLS_HOMOTOPY, "Assert, when calculating Jacobian!");
1785 debugString(LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge");
1786 }
1787 } else {
1788 if (solverData->initHomotopy)
1789 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe system is singular and not solvable.\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
1790 else {
1791 debugString(LOG_NLS_HOMOTOPY, "System singular and not solvable!");
1792 debugString(LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge");
1793 }
1794 }
1795 debugString(LOG_NLS_HOMOTOPY, "======================================================");
1796 /* update statistics */
1797 return -1;
1798 }
1799 /* Scaling back to original variables */
1800 vecMultScaling(solverData->m, solverData->dy0, solverData->xScaling, solverData->dy0);
1801 debugVectorDouble(LOG_NLS_HOMOTOPY, "tangent vector with original scaling: ", solverData->dy0, solverData->m);
1802 debugDouble(LOG_NLS_HOMOTOPY,"length of tangent vector with original scaling: ", vec2Norm(solverData->m, solverData->dy0));
1803 // vecNormalize(solverData->m, solverData->dy0, solverData->dy0);
1804 // debugVectorDouble(LOG_NLS_HOMOTOPY, "normalized tangent vector: ", solverData->dy0, solverData->m);
1805 // debugDouble(LOG_NLS_HOMOTOPY,"length of normalized tangent vector: ", vec2Norm(solverData->m, solverData->dy0));
1806
1807 /* Correct search direction, depending on the last direction (angle < 90 degree) */
1808 vecScalarProduct = vecScalarProd(solverData->m,solverData->dy0,solverData->dy2);
1809 debugDouble(LOG_NLS_HOMOTOPY,"scalar product ", vecScalarProduct);
1810 if (vecScalarProduct<0 || ((fabs(vecScalarProduct)<DBL_EPSILON2.2204460492503131e-16) && (solverData->startDirection == -1) && initialStep))
1811 {
1812 debugInt(LOG_NLS_HOMOTOPY,"initialStep = ", initialStep);
1813 debugInt(LOG_NLS_HOMOTOPY,"solverData->startDirection = ", solverData->startDirection);
1814 debugVectorDouble(LOG_NLS_HOMOTOPY,"step:",solverData->dy0, m);
1815 vecAddInv(solverData->m, solverData->dy0, solverData->dy0);
1816 debugVectorDouble(LOG_NLS_HOMOTOPY,"corrected step:",solverData->dy0, m);
1817 }
1818 /* adapt tau, if lambda + tau*delta_lambda > 1 */
1819 if (fabs(solverData->dy0[solverData->n])>1e-8)
1820 {
1821 tau = fmin(tau,(1-solverData->y0[solverData->n])/fabs(solverData->dy0[solverData->n]));
1822 }
1823 }
1824
1825 assert = 1;
1826 do {
1827 /* do update and store approximated vector in yt */
1828 vecAddScal(solverData->m, solverData->y0, solverData->dy0, tau, solverData->y1);
1829
1830 /* update function value */
1831#ifndef OMC_EMCC
1832 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1833#endif
1834 debugVectorDouble(LOG_NLS_HOMOTOPY,"y1 (predictor step):",solverData->y1, m);
1835 solverData->h_function(solverData, solverData->y1, solverData->hvec);
1836 debugVectorDouble(LOG_NLS_HOMOTOPY,"hvec (predictor step):",solverData->hvec, n);
1837 assert = 0;
1838#ifndef OMC_EMCC
1839 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1840#endif
1841 if (assert){
1842 debugString(LOG_NLS_HOMOTOPY, "Assert, when calculating function value!");
1843 debugString(LOG_NLS_HOMOTOPY, "--- decreasing step size tau in predictor step!");
1844 debugDouble(LOG_NLS_HOMOTOPY, "old tau =", tau);
1845 tau = tau/tauDecreasingFactorPredictor;
1846 debugDouble(LOG_NLS_HOMOTOPY, "new tau =", tau);
1847 }
1848 } while (assert && (tau > tauMin));
1849
1850 if (assert)
1851 {
1852 /* report solver abortion */
1853 solverData->info=-1;
1854 /* debug information */
1855 if (solverData->initHomotopy)
1856 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe step size tau cannot be decreased anymore and current tau=%g already failed.\nYou can influence the calculation of tau with the following flags:\n\t-homTauDecFac=<value>\n\t-homTauDecFacPredictor=<value>\n\t-homTauIncFac=<value>\n\t-homTauIncThreshold=<value>\n\t-homTauMax=<value>\n\t-homTauMin=<value>\n\t-homTauStart=<value>\nYou can also set the threshold for accepting the current bending with:\n\t-homAdaptBend=<value>\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.", tau);
1857 else {
1858 debugString(LOG_NLS_HOMOTOPY, "Assert, because tau cannot be decreased anymore and current tau already failed!");
1859 debugString(LOG_NLS_HOMOTOPY, "Homotopy algorithm did not converge");
1860 }
1861 debugString(LOG_NLS_HOMOTOPY, "======================================================");
1862 /* update statistics */
1863 return -1;
1864 }
1865 vecCopy(solverData->m, solverData->y1, solverData->y2);
1866 vecCopy(solverData->m, solverData->y1, solverData->yt);
1867 vecCopy(solverData->n, solverData->hvec, solverData->hvecScaled);
1868
1869 solverData->tau = tau;
1870 printHomotopyPredictorStep(LOG_NLS_HOMOTOPY, solverData);
1871
1872 /****************************************************************************
1873 * Corrector step: Newton iteration! *
1874 ****************************************************************************/
1875 debugString(LOG_NLS_HOMOTOPY, "Newton iteration for corrector step begins!");
1876
1877 /* If this is the last step, use backtrace strategy with one fixed coordinate and fix lambda */
1878 if (solverData->yt[solverData->n] == 1)
1879 {
1880 debugString(LOG_NLS_HOMOTOPY, "Force '-homBacktraceStrategy=fix' and fix lambda, because this is the last step!");
1881 debugDouble(LOG_NLS_HOMOTOPY, "Set tolerance homHEps to newtonFTol =", newtonFTol);
1882 correctorStrategy = 1;
1883 pos = solverData->n;
1884 hEps = newtonFTol;
1885 }
1886
1887 if (correctorStrategy==1)
1888 debugString(LOG_NLS_HOMOTOPY, "Using backtrace strategy with one fixed coordinate! To change this use: '-homBacktraceStrategy=orthogonal'");
1889 else
1890 debugString(LOG_NLS_HOMOTOPY, "Using backtrace strategy orthogonal to the tangent vector! To change this use: '-homBacktraceStrategy=fix'");
1891
1892
1893 for(j=0;j<maxiter;j++)
1894 {
1895 debugInt(LOG_NLS_HOMOTOPY, "Iteration: ", j+1);
1896 if (vec2Norm(solverData->n, solverData->hvec)<hEps || vec2Norm(solverData->n, solverData->hvecScaled)<hEps)
1897 {
1898 debugString(LOG_NLS_HOMOTOPY, "step accepted!");
1899 stepAccept = 1;
1900 break;
1901 }
1902 assert = 1;
1903#ifndef OMC_EMCC
1904 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1905#endif
1906 /* calculate homotopy jacobian */
1907 solverData->hJac_dh(solverData, solverData->y1, solverData->hJac);
1908 debugMatrixDouble(LOG_NLS_JAC,"Jacobian hJac:",solverData->hJac, solverData->n, solverData->n+1);
1909
1910 if (correctorStrategy==2)
1911 {
1912 /* calculate the newton matrix hJac2 for the orthogonal backtrace strategy */
1913 orthogonalBacktraceMatrix(solverData, solverData->hJac, solverData->hvec, solverData->dy0, solverData->hJac2, solverData->n, solverData->m);
1914 debugMatrixDouble(LOG_NLS_JAC,"Enhanced Jacobian hJac2 (orthogonal backtrace strategy):",solverData->hJac2, solverData->n+1, solverData->m+1);
1915 }
1916
1917 assert = 0;
1918#ifndef OMC_EMCC
1919 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1920#endif
1921 if (assert)
1922 {
1923 debugString(LOG_NLS_HOMOTOPY, "step NOT accepted, because hJac_dh could not be calculated!");
1924 stepAccept = 0;
1925 break;
1926 }
1927 matVecMultAbs(solverData->n, solverData->m, solverData->hJac, solverData->ones, solverData->resScaling);
1928 debugVectorDouble(LOG_NLS_HOMOTOPY, "residuum scaling of function h:", solverData->resScaling, solverData->n);
1929
1930 if (correctorStrategy==1) // fix one coordinate
1931 {
1932 /* copy vector h to column "pos" of the jacobian */
1933 debugVectorDouble(LOG_NLS_HOMOTOPY, "copy vector hvec to column 'pos' of the jacobian: ", solverData->hvec, solverData->n);
1934 vecCopy(solverData->n, solverData->hvec, solverData->hJac + pos*solverData->n);
1935 scaleMatrixRows(solverData->n, solverData->m, solverData->hJac);
1936 if (solveSystemWithTotalPivotSearch(solverData->n, solverData->dy1, solverData->hJac, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1)
1937 {
1938 debugString(LOG_NLS_HOMOTOPY, "step NOT accepted, because solveSystemWithTotalPivotSearch failed!");
1939 stepAccept = 0;
1940 break;
1941 }
1942 solverData->dy1[pos] = 0.0;
1943 }
1944 else // go back in orthogonal direction to tangent vector
1945 {
1946 scaleMatrixRows(solverData->n+1, solverData->m+1, solverData->hJac2);
1947 pos = solverData->n+1;
1948 if (solveSystemWithTotalPivotSearch(solverData->n+1, solverData->dy1, solverData->hJac2, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1)
1949 {
1950 debugString(LOG_NLS_HOMOTOPY, "step NOT accepted, because solveSystemWithTotalPivotSearch failed!");
1951 stepAccept = 0;
1952 break;
1953 }
1954 }
1955
1956 /* Scaling back to original variables */
1957 vecMultScaling(solverData->m, solverData->dy1, solverData->xScaling, solverData->dy1);
1958 debugVectorDouble(LOG_NLS_HOMOTOPY, "solution (original scaling): ", solverData->dy1, solverData->m);
1959
1960 vecAdd(solverData->m, solverData->y1, solverData->dy1, solverData->y2);
1961 vecCopy(solverData->m, solverData->y2, solverData->y1);
1962 debugVectorDouble(LOG_NLS_HOMOTOPY, "new y in newton: ", solverData->y1, solverData->m);
1963 assert = 1;
1964#ifndef OMC_EMCC
1965 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
1966#endif
1967 /* calculate homotopy function */
1968 solverData->h_function(solverData, solverData->y1, solverData->hvec);
1969 assert = 0;
1970#ifndef OMC_EMCC
1971 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
1972#endif
1973 if (assert)
1974 {
1975 debugString(LOG_NLS_HOMOTOPY, "step NOT accepted, because h_function could not be calculated!");
1976 stepAccept = 0;
1977 break;
1978 }
1979 /* Calculate different error measurements */
1980 vecDivScaling(solverData->n, solverData->hvec, solverData->resScaling, solverData->hvecScaled);
1981
1982 delta_x = vec2Norm(solverData->m, solverData->dy1);
1983 error_h = vec2Norm(solverData->n, solverData->hvec);
1984 error_h_scaled = vec2Norm(solverData->n, solverData->hvecScaled);
1985
1986
1987 /* debug information */
1988 debugVectorDouble(LOG_NLS_HOMOTOPY,"function values:",solverData->hvec, n);
1989 debugVectorDouble(LOG_NLS_HOMOTOPY,"scaled function values:",solverData->hvecScaled, n);
1990
1991 debugString(LOG_NLS_HOMOTOPY, "error measurements:");
1992 debugDouble(LOG_NLS_HOMOTOPY, "delta_x =", delta_x);
1993 debugDouble(LOG_NLS_HOMOTOPY, "error_h =", error_h);
1994 debugDouble(LOG_NLS_HOMOTOPY, "error_h_scaled =", error_h_scaled);
1995 debugDouble(LOG_NLS_HOMOTOPY, "hEps =", hEps);
1996
1997 }
1998 debugString(LOG_NLS_HOMOTOPY, "Newton iteration for corrector step finished!");
1999
2000 if (!assert)
2001 {
2002 vecDiff(solverData->m, solverData->y1, solverData->yt, solverData->dy1);
2003 vecDiff(solverData->m, solverData->yt, solverData->y0, solverData->dy2);
2004 printHomotopyCorrectorStep(LOG_NLS_HOMOTOPY, solverData);
2005 bend = vec2Norm(solverData->m,solverData->dy1)/vec2Norm(solverData->m,solverData->dy2);
2006
2007 debugDouble(LOG_NLS_HOMOTOPY, "vector length of predictor step =", vec2Norm(solverData->m,solverData->dy2));
2008 debugDouble(LOG_NLS_HOMOTOPY, "vector length of corrector step =", vec2Norm(solverData->m,solverData->dy1));
2009 debugDouble(LOG_NLS_HOMOTOPY, "bend =", bend);
2010 debugDouble(LOG_NLS_HOMOTOPY, "adaptBend =", adaptBend);
2011 }
2012 if ((bend > adaptBend) || !stepAccept)
2013 {
2014 if (bend<DBL_EPSILON2.2204460492503131e-16)
2015 {
2016 /* debug information */
2017 if (solverData->initHomotopy)
2018 warningStreamPrint(LOG_ASSERT, 0, "Homotopy algorithm did not converge.\nThe value specifying the bending of the homotopy curve is smaller than DBL_EPSILON (increment zero).\nYou can use -lv=LOG_INIT_HOMOTOPY,LOG_NLS_HOMOTOPY to get more information.");
2019 else
2020 debugString(LOG_NLS_HOMOTOPY, "\nINCREMENT ZERO: Homotopy algorithm did not converge\n");
2021 debugString(LOG_NLS_HOMOTOPY, "======================================================");
2022 /* update statistics */
2023 return -1;
2024 }
2025 debugString(LOG_NLS_HOMOTOPY, "The relation between the vector length of corrector step and predictor step is too big:");
2026 debugDouble(LOG_NLS_HOMOTOPY, "bend/adaptBend =", bend/adaptBend);
2027 debugString(LOG_NLS_HOMOTOPY, "--- decreasing step size tau in corrector step!");
2028 preTau = tau;
2029 debugDouble(LOG_NLS_HOMOTOPY, "old tau =", preTau);
2030 tau = fmax(tauMin,tau/tauDecreasingFactor);
2031 debugDouble(LOG_NLS_HOMOTOPY, "new tau =", tau);
2032 if (tau==preTau)
2033 iter = maxTries;
2034 else
2035 iter++;
2036 } else
2037 {
2038 initialStep = 0;
2039 iter = 0;
2040 numSteps++;
2041 if (bend < adaptBend/tauIncreasingThreshold)
2042 {
2043 debugString(LOG_NLS_HOMOTOPY, "--- increasing step size tau in corrector step!");
2044 debugDouble(LOG_NLS_HOMOTOPY, "old tau =", tau);
2045 tau = fmin(tauMax, tau*tauIncreasingFactor);
2046 debugDouble(LOG_NLS_HOMOTOPY, "new tau =", tau);
2047 }
2048 vecCopy(solverData->m, solverData->y1, solverData->y0);
2049 vecCopy(solverData->m, solverData->dy0, solverData->dy2);
2050 debugString(LOG_NLS_HOMOTOPY, "Successfull homotopy step!\n======================================================");
2051 printHomotopyUnknowns(LOG_NLS_HOMOTOPY, solverData);
2052
2053#if !defined(OMC_NO_FILESYSTEM)
2054 if(solverData->initHomotopy && ACTIVE_STREAM(LOG_INIT_HOMOTOPY)(useStream[LOG_INIT_HOMOTOPY]))
2055 {
2056 fprintf(pFile, "%.16g", solverData->y0[n]);
2057 for(i=0; i<n; ++i)
2058 fprintf(pFile, "%s%.16g", sep, solverData->y0[i]);
2059 fprintf(pFile, "\n");
2060 }
2061#endif
2062 }
2063 }
2064 if (solverData->initHomotopy)
2065 infoStreamPrint(LOG_INIT_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
2066 else
2067 infoStreamPrint(LOG_NLS_HOMOTOPY, 0, "homotopy parameter lambda = %g", solverData->y0[solverData->n]);
2068 /* copy solution back to vector x */
2069 vecCopy(solverData->n, solverData->y1, x);
2070
2071 debugString(LOG_NLS_HOMOTOPY, "HOMOTOPY ALGORITHM SUCCEEDED");
2072 if (solverData->initHomotopy) {
2073 solverData->data->simulationInfo->homotopySteps += numSteps+1;
2074 debugInt(LOG_INIT_HOMOTOPY, "Total number of lambda steps for this homotopy loop:", numSteps+1);
2075 }
2076 debugString(LOG_NLS_HOMOTOPY, "======================================================");
2077 solverData->info = 1;
2078
2079#if !defined(OMC_NO_FILESYSTEM)
2080 if(solverData->initHomotopy && ACTIVE_STREAM(LOG_INIT_HOMOTOPY)(useStream[LOG_INIT_HOMOTOPY]))
2081 fclose(pFile);
2082#endif
2083
2084 return 0;
2085}
2086
2087/*! \fn solve non-linear system with a damped Newton method combined with a homotopy approach
2088
2089 *
2090 * \param [in] [data]
2091* [sysNumber] index of the corresponding non-linear system
2092 *
2093 * \author bbachmann
2094 */
2095int solveHomotopy(DATA *data, threadData_t *threadData, int sysNumber)
2096{
2097 NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
2098 DATA_HOMOTOPY* solverData = (DATA_HOMOTOPY*)(systemData->solverData);
2099 DATA_HYBRD* solverDataHybrid;
2100
2101 /*
2102 * Get non-linear equation system
2103 */
2104 int eqSystemNumber = systemData->equationIndex;
2105 int mixedSystem = systemData->mixedSystem;
2106
2107 int i, j;
2108 int success = 0;
2109 double error_f_sqrd, error_f1_sqrd;
2110
2111 int assert = 1;
2112 int giveUp = 0;
2113 int alreadyTested = 0;
2114 int pos;
2115 int rank;
2116 int tries = 0;
2117 int runHomotopy = 0;
2118 int skipNewton = 0;
2119 int numberOfFunctionEvaluationsOld = solverData->numberOfFunctionEvaluations;
2120 solverData->casualTearingSet = systemData->strictTearingFunctionCall != NULL((void*)0);
2121 int constraintViolated;
2122 solverData->initHomotopy = systemData->initHomotopy;
2123
2124 modelica_boolean* relationsPreBackup;
2125 relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));
2126
2127 solverData->f = wrapper_fvec;
2128 solverData->f_con = wrapper_fvec_constraints;
2129 solverData->fJac_f = wrapper_fvec_der;
2130
2131 solverData->data = data;
2132 solverData->threadData = threadData;
2133 solverData->sysNumber = sysNumber;
2134 solverData->eqSystemNumber = systemData->equationIndex;
2135 solverData->mixedSystem = mixedSystem;
2136 solverData->timeValue = data->localData[0]->timeValue;
2137 solverData->minValue = systemData->min;
2138 solverData->maxValue = systemData->max;
2139 solverData->info = 0;
2140
2141 vecConst(solverData->m,1.0,solverData->ones);
2142
2143 debugString(LOG_NLS_V, "------------------------------------------------------");
2144 if (!solverData->initHomotopy)
2145 debugString(LOG_NLS_V, "SOLVING NON-LINEAR SYSTEM USING MIXED SOLVER (Newton/Homotopy solver)");
2146 else
2147 debugString(LOG_NLS_V, "SOLVING HOMOTOPY INITIALIZATION PROBLEM WITH THE HOMOTOPY SOLVER");
2148 debugInt(LOG_NLS_V, "EQUATION NUMBER:", eqSystemNumber);
2149 debugDouble(LOG_NLS_V, "TIME:", solverData->timeValue);
2150 debugInt(LOG_NLS_V, "number of function calls (so far!): ",numberOfFunctionEvaluationsOld);
2151
2152 /* set x vector */
2153 if(data->simulationInfo->discreteCall)
2154 {
2155 vecCopy(solverData->n, systemData->nlsx, solverData->xStart);
2156 debugVectorDouble(LOG_NLS_V,"System values", solverData->xStart, solverData->n);
2157 } else
2158 {
2159 vecCopy(solverData->n, systemData->nlsxExtrapolation, solverData->xStart);
2160 debugVectorDouble(LOG_NLS_V,"System extrapolation", solverData->xStart, solverData->n);
2161 }
2162 vecCopy(solverData->n, solverData->xStart, solverData->x0);
2163 // Initialize lambda variable
2164 if ((solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].homotopySupport && !solverData->initHomotopy && (solverData->data)->simulationInfo->nonlinearSystemData[solverData->sysNumber].size > solverData->n) {
2165 solverData->x0[solverData->n] = 1.0;
2166 solverData->x[solverData->n] = 1.0;
2167 solverData->x1[solverData->n] = 1.0;
2168 } else {
2169 solverData->x0[solverData->n] = 0.0;
2170 solverData->x[solverData->n] = 0.0;
2171 solverData->x1[solverData->n] = 0.0;
2172 }
2173 /* Use actual working point for scaling */
2174 for (i=0;i<solverData->n;i++){
2175 solverData->xScaling[i] = fmax(systemData->nominal[i],fabs(solverData->x0[i]));
2176 }
2177 solverData->xScaling[solverData->n] = 1.0;
2178
2179 debugVectorDouble(LOG_NLS_V,"Nominal values", systemData->nominal, solverData->n);
2180 debugVectorDouble(LOG_NLS_V,"Scaling values", solverData->xScaling, solverData->m);
2181
2182
2183 if (!solverData->initHomotopy) {
2184 /* Handle asserts of function calls, mainly necessary for fluid stuff */
2185 assert = 1;
2186 giveUp = 1;
2187 while (tries<=2)
2188 {
2189 debugVectorDouble(LOG_NLS_V,"x0", solverData->x0, solverData->n);
2190 /* evaluate with discontinuities */
2191 if(data->simulationInfo->discreteCall)
2192 {
2193 ((DATA*)data)->simulationInfo->solveContinuous = 0;
2194 }
2195 /* evaluate with discontinuities */
2196#ifndef OMC_EMCC
2197 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
2198#endif
2199 if (mixedSystem)
2200 memcpy(relationsPreBackup, data->simulationInfo->relations, sizeof(modelica_boolean)*data->modelData->nRelations);
2201
2202 if (solverData->casualTearingSet){
2203 constraintViolated = solverData->f_con(solverData, solverData->x0, solverData->f1);
2204 if (constraintViolated){
2205 giveUp = 1;
2206 break;
2207 }
2208 }
2209 else
2210 solverData->f(solverData, solverData->x0, solverData->f1);
2211
2212 /* Try to get out of here!!! */
2213 error_f_sqrd = vec2NormSqrd(solverData->n, solverData->f1);
2214 if ((error_f_sqrd - solverData->error_f_sqrd)<=0)
2215 {
2216 //infoStreamPrint(LOG_STDOUT, 0, "No Iteration at time %g needed new f = %g and old f1 = %g", solverData->timeValue, error_f_sqrd, solverData->error_f_sqrd);
2217 if (mixedSystem && data->simulationInfo->discreteCall && isNotEqualVectorInt(((DATA*)data)->modelData->nRelations, ((DATA*)data)->simulationInfo->relations, relationsPreBackup)){}
2218 else
2219 {
2220 success = 1;
2221
2222 debugString(LOG_NLS_V, "NO ITERATION NECESSARY!!!");
2223 debugString(LOG_NLS_V, "******************************************************");
2224 debugString(LOG_NLS_V,"SYSTEM SOLVED");
2225 debugInt(LOG_NLS_V, "number of function calls: ",solverData->numberOfFunctionEvaluations-numberOfFunctionEvaluationsOld);
2226 debugString(LOG_NLS_V, "------------------------------------------------------");
2227
2228 vecCopy(solverData->n, solverData->x0, systemData->nlsx);
2229 debugVectorDouble(LOG_NLS_V,"Solution", solverData->x0, solverData->n);
2230
2231 ((DATA*)data)->simulationInfo->solveContinuous = 0;
2232
2233 free(relationsPreBackup);
2234
2235 systemData->numberOfFEval = solverData->numberOfFunctionEvaluations;
2236
2237 return success;
2238 }
2239 }
2240 solverData->fJac_f(solverData, solverData->x0, solverData->fJac);
2241 vecCopy(solverData->n, solverData->f1, solverData->fJac + solverData->n*solverData->n);
2242 vecCopy(solverData->n*solverData->m, solverData->fJac, solverData->fJacx0);
2243 if (mixedSystem)
2244 memcpy(relationsPreBackup, data->simulationInfo->relations, sizeof(modelica_boolean)*data->modelData->nRelations);
2245 /* calculate scaling factor of residuals */
2246 matVecMultAbsBB(solverData->n, solverData->fJac, solverData->ones, solverData->resScaling);
2247 debugVectorDouble(LOG_NLS_JAC, "residuum scaling:", solverData->resScaling, solverData->n);
2248 scaleMatrixRows(solverData->n, solverData->m, solverData->fJac);
2249
2250 pos = solverData->n;
2251 assert = (solveSystemWithTotalPivotSearch(solverData->n, solverData->dy0, solverData->fJac, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1);
2252 if (!assert)
2253 debugString(LOG_NLS_V, "regular initial point!!!");
2254 giveUp = 0;
2255#ifndef OMC_EMCC
2256 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
2257#endif
2258 if (assert && solverData->casualTearingSet)
2259 {
2260 giveUp = 1;
2261 break;
2262 }
2263 if (assert)
2264 {
2265 tries += 1;
2266 }
2267 else
2268 break;
2269 /* break symmetry, when varying start values */
2270 /* try to find regular initial point, if necessary */
2271 if (tries == 1)
2272 {
2273 debugString(LOG_NLS_V, "assert handling:\t vary initial guess by +1%.");
2274 for(i = 0; i < solverData->n; i++)
2275 solverData->x0[i] = solverData->xStart[i] + solverData->xScaling[i]*i/solverData->n*0.01;
2276 }
2277 if (tries == 2)
2278 {
2279 debugString(LOG_NLS_V,"assert handling:\t vary initial guess by +10%.");
2280 for(i = 0; i < solverData->n; i++)
2281 solverData->x0[i] = solverData->xStart[i] + solverData->xScaling[i]*i/solverData->n*0.1;
2282 }
2283 }
2284 ((DATA*)data)->simulationInfo->solveContinuous = 1;
2285 vecCopy(solverData->n, solverData->x0, solverData->x);
2286 vecCopy(solverData->n, solverData->f1, solverData->fx0);
2287 }
2288
2289 /* start solving loop */
2290 while(!giveUp && !success)
2291 {
2292 giveUp = 1;
2293
2294 if (!skipNewton && !solverData->initHomotopy){
2295
2296 /* set x vector */
2297 if(data->simulationInfo->discreteCall){
2298 memcpy(systemData->nlsx, solverData->x, solverData->n*(sizeof(double)));
2299 }
2300 else{
2301 memcpy(systemData->nlsxExtrapolation, solverData->x, solverData->n*(sizeof(double)));
2302 }
2303
2304 newtonAlgorithm(solverData, solverData->x);
2305
2306 // If this is the casual tearing set (only exists for dynamic tearing), break after first try
2307 if (solverData->info == -1 && solverData->casualTearingSet){
2308 infoStreamPrint(LOG_NLS_V, 0, "### No Solution for the casual tearing set at the first try! ###");
2309 break;
2310 }
2311
2312 if (solverData->info == -1){
2313 solverDataHybrid = (DATA_HYBRD*)(solverData->dataHybrid);
2314 systemData->solverData = solverDataHybrid;
2315
2316 solverData->info = solveHybrd(data, threadData, sysNumber);
2317
2318 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
2319 systemData->solverData = solverData;
2320 }
2321 }
2322
2323 /* solution found */
2324 if(solverData->info == 1)
2325 {
2326 success = 1;
2327 /* This case may be switched off, because of event chattering!!!*/
2328 if(mixedSystem && data->simulationInfo->discreteCall && (alreadyTested<1))
2329 {
2330 debugVectorBool(LOG_NLS_V,"Relations Pre vector ", ((DATA*)data)->simulationInfo->relationsPre, ((DATA*)data)->modelData->nRelations);
2331 debugVectorBool(LOG_NLS_V,"Relations Backup vector ", relationsPreBackup, ((DATA*)data)->modelData->nRelations);
2332 ((DATA*)data)->simulationInfo->solveContinuous = 0;
2333
2334 if (solverData->casualTearingSet){
2335 constraintViolated = solverData->f_con(solverData, solverData->x, solverData->f1);
2336 if (constraintViolated){
2337 success = 0;
2338 break;
2339 }
2340 }
2341 else
2342 solverData->f(solverData, solverData->x, solverData->f1);
2343
2344 debugVectorBool(LOG_NLS_V,"Relations vector ", ((DATA*)data)->simulationInfo->relations, ((DATA*)data)->modelData->nRelations);
2345 if (isNotEqualVectorInt(((DATA*)data)->modelData->nRelations, ((DATA*)data)->simulationInfo->relations, relationsPreBackup)>0)
2346 {
2347 /* re-run the solution process, since relations in the system have changed */
2348 success = 0;
2349 giveUp = 0;
2350 runHomotopy = 0;
2351 alreadyTested = 1;
2352 vecCopy(solverData->n, solverData->x0, solverData->x);
2353 vecCopy(solverData->n, solverData->fx0, solverData->f1);
2354 vecCopy(solverData->n*solverData->m, solverData->fJacx0, solverData->fJac);
2355
2356 /* calculate scaling factor of residuals */
2357 matVecMultAbsBB(solverData->n, solverData->fJac, solverData->ones, solverData->resScaling);
2358 scaleMatrixRows(solverData->n, solverData->m, solverData->fJac);
2359
2360 pos = solverData->n;
2361 solveSystemWithTotalPivotSearch(solverData->n, solverData->dy0, solverData->fJac, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet);
2362 debugDouble(LOG_NLS_V,"solve mixed system at time : ", solverData->timeValue);
2363 continue;
2364 }
2365 }
2366 if (success)
2367 {
2368 debugString(LOG_NLS_V,"SYSTEM SOLVED");
2369 debugInt(LOG_NLS_V, "homotopy method: ",runHomotopy);
2370 debugInt(LOG_NLS_V, "number of function calls: ",solverData->numberOfFunctionEvaluations-numberOfFunctionEvaluationsOld);
2371 printUnknowns(LOG_NLS_V, solverData);
2372 debugString(LOG_NLS_V, "------------------------------------------------------");
2373 /* take the solution */
2374 vecCopy(solverData->n, solverData->x, systemData->nlsx);
2375 debugVectorDouble(LOG_NLS_V,"Solution", solverData->x, solverData->n);
2376 /* reset continous flag */
2377 ((DATA*)data)->simulationInfo->solveContinuous = 0;
2378 break;
2379 }
2380 }
2381 if (!success && runHomotopy>=3) break;
2382 /* Start homotopy search for new start values */
2383 vecCopy(solverData->n, solverData->x0, solverData->x);
2384 runHomotopy++;
2385 /* debug output */
2386 debugString(LOG_NLS_HOMOTOPY, "======================================================");
2387
2388 if (solverData->initHomotopy) {
2389 if (runHomotopy == 1) {
2390 solverData->h_function = wrapper_fvec;
2391 solverData->hJac_dh = wrapper_fvec_der;
2392 solverData->startDirection = omc_flag[FLAG_HOMOTOPY_NEG_START_DIR] ? -1.0 : 1.0;
2393 debugInt(LOG_INIT_HOMOTOPY, "Homotopy run: ", runHomotopy);
2394 debugDouble(LOG_INIT_HOMOTOPY,"startDirection = ", solverData->startDirection);
2395 }
2396
2397 if (runHomotopy == 2) {
2398 solverData->h_function = wrapper_fvec;
2399 solverData->hJac_dh = wrapper_fvec_der;
2400 solverData->startDirection = omc_flag[FLAG_HOMOTOPY_NEG_START_DIR] ? 1.0 : -1.0;
2401 infoStreamPrint(LOG_ASSERT, 0, "The homotopy algorithm is started again with opposing start direction.");
2402 debugInt(LOG_INIT_HOMOTOPY, "Homotopy run: ", runHomotopy);
2403 debugDouble(LOG_INIT_HOMOTOPY,"Try again with startDirection = ", solverData->startDirection);
2404 }
2405
2406 if (runHomotopy == 3) {
2407 success = 0;
2408 break;
2409 }
2410 }
2411 else {
2412 debugInt(LOG_NLS_HOMOTOPY, "Homotopy run: ", runHomotopy);
2413 if (runHomotopy == 1)
2414 {
2415 /* store x0 and calculate f(x0) -> newton homotopy, fJac(x0) -> taylor, affin homotopy */
2416 solverData->h_function = wrapper_fvec_homotopy_newton;
2417 solverData->hJac_dh = wrapper_fvec_homotopy_newton_der;
2418 solverData->startDirection = 1.0;
2419 debugDouble(LOG_NLS_HOMOTOPY,"STARTING NEWTON HOMOTOPY METHOD; startDirection = ", solverData->startDirection);
2420 }
2421 if (runHomotopy == 2)
2422 {
2423 /* store x0 and calculate f(x0) -> newton homotopy, fJac(x0) -> taylor, affin homotopy */
2424 solverData->h_function = wrapper_fvec_homotopy_newton;
2425 solverData->hJac_dh = wrapper_fvec_homotopy_newton_der;
2426 solverData->startDirection = -1.0;
2427 debugDouble(LOG_NLS_HOMOTOPY,"STARTING NEWTON HOMOTOPY METHOD; startDirection = ", solverData->startDirection);
2428 }
2429 if (runHomotopy == 3)
2430 {
2431 solverData->h_function = wrapper_fvec_homotopy_fixpoint;
2432 solverData->hJac_dh = wrapper_fvec_homotopy_fixpoint_der;
2433 solverData->startDirection = 1.0;
2434 debugDouble(LOG_NLS_HOMOTOPY,"STARTING FIXPOINT HOMOTOPY METHOD = ", solverData->startDirection);
2435 }
2436 }
2437
2438 homotopyAlgorithm(solverData, solverData->x);
2439
2440 if (solverData->info<1)
2441 {
2442 skipNewton = 1;
2443 giveUp = runHomotopy>=3;
2444
2445 } else if (solverData->initHomotopy && solverData->info==1) {
2446 /* take the solution */
2447 vecCopy(solverData->n, solverData->x, systemData->nlsx);
2448 debugVectorDouble(LOG_NLS_V,"Solution", solverData->x, solverData->n);
2449 success = 1;
2450 }
2451
2452 else {
2453 assert = 1;
2454#ifndef OMC_EMCC
2455 MMC_TRY_INTERNAL(simulationJumpBuffer){ jmp_buf new_mmc_jumper, *old_jumper __attribute__((unused))
= threadData->simulationJumpBuffer; threadData->simulationJumpBuffer
= &new_mmc_jumper; if (_setjmp (new_mmc_jumper) == 0) {
2456#endif
2457 if (solverData->casualTearingSet){
2458 constraintViolated = solverData->f_con(solverData, solverData->x, solverData->f1);
2459 if (constraintViolated){
2460 success = 0;
2461 break;
2462 }
2463 }
2464 else
2465 solverData->f(solverData, solverData->x, solverData->f1);
2466
2467 solverData->fJac_f(solverData, solverData->x, solverData->fJac);
2468 vecCopy(solverData->n, solverData->f1, solverData->fJac + solverData->n*solverData->n);
2469 /* calculate scaling factor of residuals */
2470 matVecMultAbsBB(solverData->n, solverData->fJac, solverData->ones, solverData->resScaling);
2471 debugVectorDouble(LOG_NLS_JAC, "residuum scaling:", solverData->resScaling, solverData->n);
2472 scaleMatrixRows(solverData->n, solverData->m, solverData->fJac);
2473
2474 pos = solverData->n;
2475 assert = (solveSystemWithTotalPivotSearch(solverData->n, solverData->dy0, solverData->fJac, solverData->indRow, solverData->indCol, &pos, &rank, solverData->casualTearingSet) == -1);
2476 if (!assert)
2477 debugString(LOG_NLS_V, "regular initial point!!!");
2478#ifndef OMC_EMCC
2479 MMC_CATCH_INTERNAL(simulationJumpBuffer)} threadData->simulationJumpBuffer = old_jumper;mmc_catch_dummy_fn
();}
2480 #endif
2481 if (assert)
2482 {
2483 giveUp = 1;
2484 } else
2485 {
2486 giveUp = 0;
2487 skipNewton = 0;
2488 }
2489 }
2490 }
2491 if (!success)
2492 {
2493 debugString(LOG_NLS_V,"Homotopy solver did not converge!");
2494 }
2495 free(relationsPreBackup);
2496
2497 /* write statistics */
2498 systemData->numberOfFEval = solverData->numberOfFunctionEvaluations;
2499 systemData->numberOfIterations = solverData->numberOfIterations;
2500
2501 return success;
2502}
2503
2504#endif