Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/nonlinearSolverNewton.c
Warning:line 84, column 13
Value stored to 'index' during its initialization is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2014, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
11 * GPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
14 * ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the OSMC (Open Source Modelica Consortium)
17 * Public License (OSMC-PL) are obtained from OSMC, either from the above
18 * address, from the URLs: http://www.openmodelica.org or
19 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica
20 * distribution. GNU version 3 is obtained from:
21 * http://www.gnu.org/copyleft/gpl.html. The New BSD License is obtained from:
22 * http://www.opensource.org/licenses/BSD-3-Clause.
23 *
24 * This program is distributed WITHOUT ANY WARRANTY; without even the implied
25 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS
26 * EXPRESSLY SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE
27 * CONDITIONS OF OSMC-PL.
28 *
29 */
30
31/*! \file nonlinearSolverNewton.c
32 */
33
34#ifdef __cplusplus
35extern "C" {
36#endif
37
38#include <math.h>
39#include <stdlib.h>
40#include <string.h> /* memcpy */
41
42#include "simulation/simulation_info_json.h"
43#include "util/omc_error.h"
44#include "util/varinfo.h"
45#include "model_help.h"
46
47#include "nonlinearSystem.h"
48#include "nonlinearSolverNewton.h"
49#include "newtonIteration.h"
50
51#include "external_input.h"
52
53
54extern double enorm_(int *n, double *x);
55int wrapper_fvec_newton(int* n, double* x, double* fvec, void* userdata, int fj);
56
57
58#ifdef __cplusplus
59extern "C" {
60#endif
61
62extern int dgesv_(int *n, int *nrhs, doublereal *a, int *lda, int *ipiv, doublereal *b, int *ldb, int *info);
63
64#ifdef __cplusplus
65}
66#endif
67
68
69/*! \fn getAnalyticalJacobian
70 *
71 * function calculates analytical jacobian
72 *
73 * \param [ref] [data]
74 * \param [out] [jac]
75 *
76 * \author wbraun
77 *
78 */
79int getAnalyticalJacobianNewton(DATA* data, threadData_t *threadData, double* jac, int sysNumber)
80{
81 int i,j,k,l,ii,currentSys = sysNumber;
82 NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->nonlinearSystemData[currentSys]);
83 DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData);
84 const int index = systemData->jacobianIndex;
Value stored to 'index' during its initialization is never read
85 ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]);
86
87 memset(jac, 0, (solverData->n)*(solverData->n)*sizeof(double));
88
89 for(i=0; i < jacobian->sparsePattern->maxColors; i++)
90 {
91 /* activate seed variable for the corresponding color */
92 for(ii=0; ii < jacobian->sizeCols; ii++)
93 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
94 jacobian->seedVars[ii] = 1;
95
96 systemData->analyticalJacobianColumn(data, threadData, jacobian, NULL((void*)0));
97
98 for(j = 0; j < jacobian->sizeCols; j++)
99 {
100 if(jacobian->seedVars[j] == 1)
101 {
102 ii = jacobian->sparsePattern->leadindex[j];
103 while(ii < jacobian->sparsePattern->leadindex[j+1])
104 {
105 l = jacobian->sparsePattern->index[ii];
106 k = j*jacobian->sizeRows + l;
107 jac[k] = jacobian->resultVars[l];
108 ii++;
109 };
110 }
111 /* de-activate seed variable for the corresponding color */
112 if(jacobian->sparsePattern->colorCols[j]-1 == i)
113 jacobian->seedVars[j] = 0;
114 }
115
116 }
117
118 return 0;
119}
120
121
122/*! \fn wrapper_fvec_newton for the residual Function
123 * tensolve calls for the subroutine fcn(n, x, fvec, iflag, data)
124 *
125 * fj decides whether the function values or the jacobian matrix shall be calculated
126 * fj = 1 ==> calculate function values
127 * fj = 0 ==> calculate jacobian matrix
128 */
129int wrapper_fvec_newton(int* n, double* x, double* fvec, void* userdata, int fj)
130{
131 DATA_USER* uData = (DATA_USER*) userdata;
132 DATA* data = (DATA*)(uData->data);
133 void *dataAndThreadData[2] = {data, uData->threadData};
134 int currentSys = ((DATA_USER*)userdata)->sysNumber;
135 NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[currentSys]);
136 DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData);
137 int flag = 1;
138 int *iflag=&flag;
139
140 if (fj) {
141 (data->simulationInfo->nonlinearSystemData[currentSys].residualFunc)(dataAndThreadData, x, fvec, iflag);
142 } else {
143 /* performance measurement */
144 rt_ext_tp_tick(&systemData->jacobianTimeClock);
145
146 if(systemData->jacobianIndex != -1) {
147 getAnalyticalJacobianNewton(data, uData->threadData, solverData->fjac, currentSys);
148 } else {
149 double delta_h = sqrt(solverData->epsfcn);
150 double delta_hh;
151 double xsave;
152
153 int i,j,l, linear=0;
154
155 for(i = 0; i < *n; i++) {
156 delta_hh = fmax(delta_h * fmax(fabs(x[i]), fabs(fvec[i])), delta_h);
157 delta_hh = ((fvec[i] >= 0) ? delta_hh : -delta_hh);
158 delta_hh = x[i] + delta_hh - x[i];
159 xsave = x[i];
160 x[i] += delta_hh;
161 delta_hh = 1. / delta_hh;
162
163 wrapper_fvec_newton(n, x, solverData->rwork, userdata, 1);
164 solverData->nfev++;
165
166 for(j = 0; j < *n; j++) {
167 l = i * *n + j;
168 solverData->fjac[l] = (solverData->rwork[j] - fvec[j]) * delta_hh;
169 }
170 x[i] = xsave;
171 }
172 }
173 /* performance measurement and statistics */
174 systemData->jacobianTime += rt_ext_tp_tock(&(systemData->jacobianTimeClock));
175 systemData->numberOfJEval++;
176 }
177 return *iflag;
178}
179
180/*! \fn solve non-linear system with newton method
181 *
182 * \param [in] [data]
183 * [sysNumber] index of the corresponding non-linear system
184 *
185 * \author wbraun
186 */
187int solveNewton(DATA *data, threadData_t *threadData, int sysNumber)
188{
189 NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
190 DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData);
191 int eqSystemNumber = 0;
192 int i;
193 double xerror = -1, xerror_scaled = -1;
194 int success = 0;
195 int nfunc_evals = 0;
196 int continuous = 1;
197 double local_tol = solverData->ftol;
198
199 int giveUp = 0;
200 int retries = 0;
201 int retries2 = 0;
202 int nonContinuousCase = 0;
203 modelica_boolean *relationsPreBackup = NULL((void*)0);
204 int casualTearingSet = data->simulationInfo->nonlinearSystemData[sysNumber].strictTearingFunctionCall != NULL((void*)0);
205
206 DATA_USER* userdata = (DATA_USER*)malloc(sizeof(DATA_USER));
207 assert(userdata != NULL)((void) sizeof ((userdata != ((void*)0)) ? 1 : 0), __extension__
({ if (userdata != ((void*)0)) ; else __assert_fail ("userdata != NULL"
, "./simulation/solver/nonlinearSolverNewton.c", 207, __extension__
__PRETTY_FUNCTION__); }))
;
208
209 userdata->data = (void*)data;
210 userdata->threadData = threadData;
211 userdata->sysNumber = sysNumber;
212
213 /*
214 * We are given the number of the non-linear system.
215 * We want to look it up among all equations.
216 */
217 eqSystemNumber = systemData->equationIndex;
218
219 local_tol = solverData->ftol;
220
221 relationsPreBackup = (modelica_boolean*) malloc(data->modelData->nRelations*sizeof(modelica_boolean));
222
223 solverData->nfev = 0;
224
225 /* try to calculate jacobian only once at the beginning of the iteration */
226 solverData->calculate_jacobian = 0;
227
228 // Initialize lambda variable
229 if (data->simulationInfo->nonlinearSystemData[sysNumber].homotopySupport) {
230 solverData->x[solverData->n] = 1.0;
231 solverData->x_new[solverData->n] = 1.0;
232 }
233 else {
234 solverData->x[solverData->n] = 0.0;
235 solverData->x_new[solverData->n] = 0.0;
236 }
237
238 /* debug output */
239 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
240 {
241 int indexes[2] = {1,eqSystemNumber};
242 infoStreamPrintWithEquationIndexes(LOG_NLS_V, 1, indexes, "Start solving Non-Linear System %d at time %g with Newton Solver",
243 eqSystemNumber, data->localData[0]->timeValue);
244
245 for(i = 0; i < solverData->n; i++) {
246 infoStreamPrint(LOG_NLS_V, 1, "x[%d] = %.15e", i, data->simulationInfo->discreteCall ? systemData->nlsx[i] : systemData->nlsxExtrapolation[i]);
247 infoStreamPrint(LOG_NLS_V, 0, "nominal = %g +++ nlsx = %g +++ old = %g +++ extrapolated = %g",
248 systemData->nominal[i], systemData->nlsx[i], systemData->nlsxOld[i], systemData->nlsxExtrapolation[i]);
249 messageClose(LOG_NLS_V);
250 }
251 messageClose(LOG_NLS_V);
252 }
253
254 /* set x vector */
255 if(data->simulationInfo->discreteCall) {
256 memcpy(solverData->x, systemData->nlsx, solverData->n*(sizeof(double)));
257 } else {
258 memcpy(solverData->x, systemData->nlsxExtrapolation, solverData->n*(sizeof(double)));
259 }
260
261 /* start solving loop */
262 while(!giveUp && !success)
263 {
264
265 giveUp = 1;
266 solverData->newtonStrategy = data->simulationInfo->newtonStrategy;
267 _omc_newton(wrapper_fvec_newton, solverData, (void*)userdata);
268
269 /* check for proper inputs */
270 if(solverData->info == 0)
271 printErrorEqSyst(IMPROPER_INPUT, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber), data->localData[0]->timeValue);
272
273 /* reset non-contunuousCase */
274 if(nonContinuousCase && xerror > local_tol && xerror_scaled > local_tol)
275 {
276 memcpy(data->simulationInfo->relationsPre, relationsPreBackup, sizeof(modelica_boolean)*data->modelData->nRelations);
277 nonContinuousCase = 0;
278 }
279
280 /* check for error */
281 xerror_scaled = enorm_(&solverData->n, solverData->fvecScaled);
282 xerror = enorm_(&solverData->n, solverData->fvec);
283
284 /* solution found */
285 if((xerror <= local_tol || xerror_scaled <= local_tol) && solverData->info > 0)
286 {
287 success = 1;
288 nfunc_evals += solverData->nfev;
289 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
290 {
291 infoStreamPrint(LOG_NLS_V, 0, "*** System solved ***\n%d restarts", retries);
292 infoStreamPrint(LOG_NLS_V, 0, "nfunc = %d +++ error = %.15e +++ error_scaled = %.15e", nfunc_evals, xerror, xerror_scaled);
293 for(i = 0; i < solverData->n; i++)
294 infoStreamPrint(LOG_NLS_V, 0, "x[%d] = %.15e\n\tresidual = %e", i, solverData->x[i], solverData->fvec[i]);
295 }
296
297 /* take the solution */
298 memcpy(systemData->nlsx, solverData->x, solverData->n*(sizeof(double)));
299
300 /* Then try with old values (instead of extrapolating )*/
301 }
302 // If this is the casual tearing set (only exists for dynamic tearing), break after first try
303 else if(retries < 1 && casualTearingSet)
304 {
305 giveUp = 1;
306 infoStreamPrint(LOG_NLS_V, 0, "### No Solution for the casual tearing set at the first try! ###");
307 }
308 else if(retries < 1)
309 {
310 memcpy(solverData->x, systemData->nlsxOld, solverData->n*(sizeof(double)));
311
312 retries++;
313 giveUp = 0;
314 nfunc_evals += solverData->nfev;
315 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t try old values.");
316 /* try to vary the initial values */
317
318 /* evaluate jacobian in every step now */
319 solverData->calculate_jacobian = 1;
320 }
321 else if(retries < 2)
322 {
323 for(i = 0; i < solverData->n; i++)
324 solverData->x[i] += systemData->nominal[i] * 0.01;
325 retries++;
326 giveUp = 0;
327 nfunc_evals += solverData->nfev;
328 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t vary solution point by 1%%.");
329 /* try to vary the initial values */
330 }
331 else if(retries < 3)
332 {
333 for(i = 0; i < solverData->n; i++)
334 solverData->x[i] = systemData->nominal[i];
335 retries++;
336 giveUp = 0;
337 nfunc_evals += solverData->nfev;
338 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t try nominal values as initial solution.");
339 }
340 else if(retries < 4 && data->simulationInfo->discreteCall)
341 {
342 /* try to solve non-continuous
343 * work-a-round: since other wise some model does
344 * stuck in event iteration. e.g.: Modelica.Mechanics.Rotational.Examples.HeatLosses
345 */
346
347 memcpy(solverData->x, systemData->nlsxOld, solverData->n*(sizeof(double)));
348 retries++;
349
350 /* try to solve a discontinuous system */
351 continuous = 0;
352
353 nonContinuousCase = 1;
354 memcpy(relationsPreBackup, data->simulationInfo->relationsPre, sizeof(modelica_boolean)*data->modelData->nRelations);
355
356 giveUp = 0;
357 nfunc_evals += solverData->nfev;
358 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t try to solve a discontinuous system.");
359 }
360 else if(retries2 < 4)
361 {
362 memcpy(solverData->x, systemData->nlsxOld, solverData->n*(sizeof(double)));
363 /* reduce tolarance */
364 local_tol = local_tol*10;
365
366 retries = 0;
367 retries2++;
368 giveUp = 0;
369 nfunc_evals += solverData->nfev;
370 infoStreamPrint(LOG_NLS_V, 0, " - iteration making no progress:\t reduce the tolerance slightly to %e.", local_tol);
371 }
372 else
373 {
374 printErrorEqSyst(ERROR_AT_TIME, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber), data->localData[0]->timeValue);
375 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
376 {
377 infoStreamPrint(LOG_NLS_V, 0, "### No Solution! ###\n after %d restarts", retries);
378 infoStreamPrint(LOG_NLS_V, 0, "nfunc = %d +++ error = %.15e +++ error_scaled = %.15e", nfunc_evals, xerror, xerror_scaled);
379 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
380 for(i = 0; i < solverData->n; i++)
381 infoStreamPrint(LOG_NLS_V, 0, "x[%d] = %.15e\n\tresidual = %e", i, solverData->x[i], solverData->fvec[i]);
382 }
383 }
384 }
385 if(ACTIVE_STREAM(LOG_NLS_V)(useStream[LOG_NLS_V]))
386 messageClose(LOG_NLS_V);
387
388 free(relationsPreBackup);
389
390 /* write statistics */
391 systemData->numberOfFEval = solverData->numberOfFunctionEvaluations;
392 systemData->numberOfIterations = solverData->numberOfIterations;
393
394 return success;
395}