Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/nonlinearSolverHomotopy.c
Warning:line 531, column 15
Use of zero-allocated memory

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]))
87
Taking true branch
526 {
527 int i;
528 char *buffer = (char*)malloc(sizeof(char)*n*20);
88
Memory is allocated
529
530 infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
531 buffer[0] = 0;
89
Use of zero-allocated memory
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++) {
77
Loop condition is false. Execution continues on line 1117
1115 indRow[i] = i;
1116 }
1117 for (i=0; i<m; i++) {
78
Assuming 'i' is < 'm'
79
Loop condition is true. Entering loop body
80
Loop condition is false. Execution continues on line 1120
1118 indCol[i] = i;
1119 }
1120 if (*pos>=0) {
81
Taking true branch
1121 indCol[n] = *pos;
1122 indCol[*pos] = n;
1123 } else {
1124 nPivot = n+1;
1125 }
1126
1127 for (i = 0; i < n; i++) {
82
Loop condition is false. Execution continues on line 1159
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];
83
Loop condition is false. Execution continues on line 1161
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))
){
84
Taking false branch
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--) {
85
Loop condition is false. Execution continues on line 1191
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);
86
Calling 'debugVectorInt'
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;
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);
1
Assuming the condition is false
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)
2
Assuming the condition is false
3
Taking false branch
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)
4
Assuming the condition is false
5
Taking false branch
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) {
6
Assuming the condition is false
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++){
7
Assuming the condition is false
8
Loop condition is false. Execution continues on line 2177
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) {
9
Taking false branch
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)
10
Loop condition is true. Entering loop body
27
Loop condition is true. Entering loop body
46
Loop condition is true. Entering loop body
59
Loop condition is true. Entering loop body
2291 {
2292 giveUp = 1;
2293
2294 if (!skipNewton && !solverData->initHomotopy){
11
Taking false branch
28
Assuming the condition is false
29
Taking false branch
47
Assuming the condition is false
48
Taking false branch
60
Assuming the condition is false
61
Taking false branch
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)
12
Taking false branch
30
Assuming the condition is false
31
Taking false branch
49
Assuming the condition is true
50
Taking true branch
62
Assuming the condition is false
63
Taking false branch
2325 {
2326 success = 1;
2327 /* This case may be switched off, because of event chattering!!!*/
2328 if(mixedSystem && data->simulationInfo->discreteCall && (alreadyTested<1))
51
Assuming 'mixedSystem' is not equal to 0
52
Assuming the condition is true
53
Taking true branch
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){
54
Assuming the condition is false
55
Taking false branch
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)
56
Assuming the condition is true
57
Taking true branch
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;
58
Execution continues on line 2290
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;
13
Taking false branch
32
Taking false branch
64
Taking false branch
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) {
14
Taking true branch
33
Taking true branch
65
Taking true branch
2389 if (runHomotopy == 1) {
15
Taking true branch
34
Taking false branch
66
Taking true branch
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;
16
Assuming the condition is false
17
'?' condition is false
67
Assuming the condition is true
68
'?' condition is true
2393 debugInt(LOG_INIT_HOMOTOPY, "Homotopy run: ", runHomotopy);
2394 debugDouble(LOG_INIT_HOMOTOPY,"startDirection = ", solverData->startDirection);
2395 }
2396
2397 if (runHomotopy == 2) {
18
Taking false branch
35
Taking true branch
69
Taking false branch
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;
36
Assuming the condition is false
37
'?' condition is false
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) {
19
Taking false branch
38
Taking false branch
70
Taking false branch
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)
20
Assuming the condition is false
21
Taking false branch
39
Assuming the condition is false
40
Taking false branch
71
Assuming the condition is false
72
Taking false branch
2441 {
2442 skipNewton = 1;
2443 giveUp = runHomotopy>=3;
2444
2445 } else if (solverData->initHomotopy && solverData->info==1) {
22
Assuming the condition is false
41
Assuming the condition is false
73
Assuming the condition is false
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){
23
Assuming the condition is false
24
Taking false branch
42
Assuming the condition is false
43
Taking false branch
74
Assuming the condition is false
75
Taking false branch
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);
76
Calling 'solveSystemWithTotalPivotSearch'
2476 if (!assert)
25
Taking true branch
44
Taking true branch
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)
26
Taking false branch
45
Taking false branch
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