Bug Summary

File:OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c
Warning:line 327, 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 nonlinear_solver.c
32 */
33
34#include <math.h>
35#include <stdlib.h>
36#include <string.h> /* memcpy */
37
38#include "../../simulation_data.h"
39#include "../simulation_info_json.h"
40#include "../../util/omc_error.h"
41#include "../../util/parallel_helper.h"
42#include "omc_math.h"
43#include "../../util/varinfo.h"
44#include "model_help.h"
45
46#include "linearSystem.h"
47#include "linearSolverTotalPivot.h"
48
49void debugMatrixDoubleLS(int logName, char* matrixName, double* matrix, int n, int m)
50{
51 if(ACTIVE_STREAM(logName)(useStream[logName]))
52 {
53 int i, j;
54 int sparsity = 0;
55 char *buffer = (char*)malloc(sizeof(char)*m*18);
56
57 infoStreamPrint(logName, 1, "%s [%dx%d-dim]", matrixName, n, m);
58 for(i=0; i<n;i++)
59 {
60 buffer[0] = 0;
61 for(j=0; j<m; j++)
62 {
63 if (sparsity)
64 {
65 if (fabs(matrix[i + j*(m-1)])<1e-12)
66 sprintf(buffer, "%s 0", buffer);
67 else
68 sprintf(buffer, "%s *", buffer);
69 }
70 else
71 {
72 sprintf(buffer, "%s%12.4g ", buffer, matrix[i + j*(m-1)]);
73 }
74 }
75 infoStreamPrint(logName, 0, "%s", buffer);
76 }
77 messageClose(logName);
78 free(buffer);
79 }
80}
81
82void debugVectorDoubleLS(int logName, char* vectorName, double* vector, int n)
83{
84 if(ACTIVE_STREAM(logName)(useStream[logName]))
85 {
86 int i;
87 char *buffer = (char*)malloc(sizeof(char)*n*22);
88
89 infoStreamPrint(logName, 1, "%s [%d-dim]", vectorName, n);
90 buffer[0] = 0;
91 for(i=0; i<n;i++)
92 {
93 if (vector[i]<-1e+300)
94 sprintf(buffer, "%s -INF ", buffer);
95 else if (vector[i]>1e+300)
96 sprintf(buffer, "%s +INF ", buffer);
97 else
98 sprintf(buffer, "%s%16.8g ", buffer, vector[i]);
99 }
100 infoStreamPrint(logName, 0, "%s", buffer);
101 free(buffer);
102 messageClose(logName);
103 }
104}
105
106void debugStringLS(int logName, char* message)
107{
108 if(ACTIVE_STREAM(logName)(useStream[logName]))
109 {
110 infoStreamPrint(logName, 1, "%s", message);
111 messageClose(logName);
112 }
113}
114
115void debugIntLS(int logName, char* message, int value)
116{
117 if(ACTIVE_STREAM(logName)(useStream[logName]))
118 {
119 infoStreamPrint(logName, 1, "%s %d", message, value);
120 messageClose(logName);
121 }
122}
123void vecMultScalingLS(int n, double *a, double *b, double *c)
124{
125 int i;
126 for (i=0;i<n;i++)
127 c[i] = a[i]*fabs(b[i]);
128}
129
130void vecAddScalLS(int n, double *a, double *b, double s, double *c)
131{
132 int i;
133 for (i=0;i<n;i++)
134 c[i] = a[i] + s*b[i];
135}
136
137void vecAddLS(int n, double *a, double *b, double *c)
138{
139 int i;
140 for (i=0;i<n;i++)
141 c[i] = a[i] + b[i];
142}
143
144void vecCopyLS(int n, double *a, double *b)
145{
146 memcpy(b, a, n*(sizeof(double)));
147}
148
149void vecConstLS(int n, double value, double *a)
150{
151 int i;
152 for (i=0;i<n;i++)
153 a[i] = value;
154}
155
156void vecScalarMultLS(int n, double *a, double s, double *b)
157{
158 int i;
159 for (i=0;i<n;i++)
160 b[i] = s*a[i];
161}
162
163void getIndicesOfPivotElementLS(int *n, int *m, int *l, double* A, int *indRow, int *indCol, int *pRow, int *pCol, double *absMax)
164{
165 int i, j;
166
167 *absMax = fabs(A[indRow[*l] + indCol[*l]* *n]);
168 *pCol = *l;
169 *pRow = *l;
170 for (i = *l; i < *n; i++) {
171 for (j = *l; j < *m; j++) {
172 if (fabs(A[indRow[i] + indCol[j]* *n]) > *absMax) {
173 *absMax = fabs(A[indRow[i] + indCol[j]* *n]);
174 *pCol = j;
175 *pRow = i;
176 }
177 }
178 }
179}
180
181/*! \ linear solver for A*x = b based on a total pivot search
182 * \ input matrix Ab: first n columns are matrix A
183 * last column is -b
184 * output x: solution dim n+1, last column is 1 for solvable systems!
185 * rank: rank of matrix Ab
186 * det: determinant of matrix A
187 *
188 * \author bbachmann
189 */
190
191int solveSystemWithTotalPivotSearchLS(int n, double* x, double* A, int* indRow, int* indCol, int *rank)
192{
193 int i, k, j, l, m=n+1, nrsh=1, singular=0;
194 int pCol, pRow;
195 double hValue;
196 double hInt;
197 double absMax;
198 int r,s;
199 int permutation = 1;
200
201 /* assume full rank of matrix [n x (n+1)] */
202 *rank = n;
203
204 for (i=0; i<n; i++) {
205 indRow[i] = i;
206 }
207 for (i=0; i<m; i++) {
208 indCol[i] = i;
209 }
210
211 for (i = 0; i < n; i++) {
212 getIndicesOfPivotElementLS(&n, &n, &i, A, indRow, indCol, &pRow, &pCol, &absMax);
213 /* this criteria should be evaluated and may be improved in future */
214 if (absMax<DBL_EPSILON2.2204460492503131e-16) {
215 *rank = i;
216 warningStreamPrint(LOG_LS, 0, "Matrix singular!");
217 debugIntLS(LOG_LS,"rank = ", *rank);
218 break;
219 }
220 /* swap row indices */
221 if (pRow!=i) {
222 hInt = indRow[i];
223 indRow[i] = indRow[pRow];
224 indRow[pRow] = hInt;
225 }
226 /* swap column indices */
227 if (pCol!=i) {
228 hInt = indCol[i];
229 indCol[i] = indCol[pCol];
230 indCol[pCol] = hInt;
231 }
232
233 /* Gauss elimination of row indRow[i] */
234 for (k=i+1; k<n; k++) {
235 hValue = -A[indRow[k] + indCol[i]*n]/A[indRow[i] + indCol[i]*n];
236 for (j=i+1; j<m; j++) {
237 A[indRow[k] + indCol[j]*n] = A[indRow[k] + indCol[j]*n] + hValue*A[indRow[i] + indCol[j]*n];
238 }
239 A[indRow[k] + indCol[i]*n] = 0;
240 }
241 }
242
243 debugMatrixDoubleLS(LOG_LS_V,"LGS: matrix Ab manipulated",A, n, n+1);
244 /* solve even singular matrix !!! */
245 for (i=n-1;i>=0; i--) {
246 if (i>=*rank) {
247 /* this criteria should be evaluated and may be improved in future */
248 if (fabs(A[indRow[i] + n*n])>1e-12) {
249 warningStreamPrint(LOG_LS, 0, "under-determined linear system not solvable!");
250 return -1;
251 } else {
252 x[indCol[i]] = 0.0;
253 }
254 } else {
255 x[indCol[i]] = -A[indRow[i] + n*n];
256 for (j=n-1; j>i; j--) {
257 x[indCol[i]] = x[indCol[i]] - A[indRow[i] + indCol[j]*n]*x[indCol[j]];
258 }
259 x[indCol[i]]=x[indCol[i]]/A[indRow[i] + indCol[i]*n];
260 }
261 }
262 x[n]=1.0;
263 debugVectorDoubleLS(LOG_LS_V,"LGS: solution vector x",x, n+1);
264
265 return 0;
266}
267
268/*! \fn allocate memory for linear system solver totalpivot
269 *
270 * \author bbachmann
271 */
272int allocateTotalPivotData(int size, void** voiddata)
273{
274 DATA_TOTALPIVOT* data = (DATA_TOTALPIVOT*) malloc(sizeof(DATA_TOTALPIVOT));
275
276 /* memory for linear system */
277 data->Ab = (double*) calloc((size*(size+1)),sizeof(double));
278 data->b = (double*) malloc(size*sizeof(double));
279 data->x = (double*) calloc(size+1,sizeof(double));
280
281 /* used for pivot strategy */
282 data->indRow =(int*) calloc(size,sizeof(int));
283 data->indCol =(int*) calloc(size+1,sizeof(int));
284
285 voiddata[1] = (void*)data;
286 return 0;
287}
288
289/*! \fn free memory for nonlinear solver totalpivot
290 *
291 * \author bbachmann
292 */
293int freeTotalPivotData(void** voiddata)
294{
295 DATA_TOTALPIVOT* data = (DATA_TOTALPIVOT*) voiddata[1];
296
297 /* memory for linear system */
298 free(data->Ab);
299 free(data->b);
300 free(data->x);
301
302 /* used for pivot strategy */
303 free(data->indRow);
304 free(data->indCol);
305
306 free(voiddata[1]);
307 voiddata[1] = NULL((void*)0);
308
309 return 0;
310}
311
312/*! \fn getAnalyticalJacobian
313 *
314 * function calculates analytical jacobian
315 *
316 * \param [ref] [data]
317 * \param [out] [jac]
318 *
319 * \author wbraun
320 *
321 */
322int getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, double* jac, int sysNumber)
323{
324 int i,j,k,l,ii;
325 LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]);
326
327 const int index = systemData->jacobianIndex;
Value stored to 'index' during its initialization is never read
328
329 ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian;
330 ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian;
331
332 memset(jac, 0, (systemData->size)*(systemData->size)*sizeof(double));
333
334 for(i=0; i < jacobian->sparsePattern->maxColors; i++)
335 {
336 /* activate seed variable for the corresponding color */
337 for(ii=0; ii < jacobian->sizeCols; ii++)
338 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
339 jacobian->seedVars[ii] = 1;
340
341 ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian);
342
343 for(j = 0; j < jacobian->sizeCols; j++)
344 {
345 if(jacobian->seedVars[j] == 1)
346 {
347 ii = jacobian->sparsePattern->leadindex[j];
348 while(ii < jacobian->sparsePattern->leadindex[j+1]) {
349 l = jacobian->sparsePattern->index[ii];
350 k = j*jacobian->sizeRows + l;
351 jac[k] = jacobian->resultVars[l];
352 ii++;
353 }
354 }
355 /* de-activate seed variable for the corresponding color */
356 if(jacobian->sparsePattern->colorCols[j]-1 == i) {
357 jacobian->seedVars[j] = 0;
358 }
359 }
360
361 }
362 return 0;
363}
364
365/*! \fn wrapper_fvec_hybrd for the residual Function
366 * calls for the subroutine fcn(n, x, fvec, iflag, data)
367 *
368 *
369 */
370static int wrapper_fvec_totalpivot(double* x, double* f, void** data, int sysNumber)
371{
372 int currentSys = sysNumber;
373 int iflag = 0;
374 /* NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->nonlinearSystemData[currentSys]); */
375 /* DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); */
376
377 (*((DATA*)data[0])->simulationInfo->linearSystemData[currentSys].residualFunc)(data, x, f, &iflag);
378 return 0;
379}
380
381/*! \fn solve linear system with totalpivot method
382 *
383 * \param [in] [data]
384 * [sysNumber] index of the corresponing non-linear system
385 *
386 * \author bbachmann
387 */
388int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
389{
390 void *dataAndThreadData[2] = {data, threadData};
391 int i, j;
392 LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
393 DATA_TOTALPIVOT* solverData = (DATA_TOTALPIVOT*) systemData->parDynamicData[omc_get_thread_num()].solverData[1];
394
395 int n = systemData->size, status;
396 double fdeps = 1e-8;
397 double xTol = 1e-8;
398 int eqSystemNumber = systemData->equationIndex;
399 int indexes[2] = {1,eqSystemNumber};
400 int rank;
401 _omc_scalar residualNorm = 0;
402
403 /* We are given the number of the linear system.
404 * We want to look it up among all equations. */
405 /* int eqSystemNumber = systemData->equationIndex; */
406 int success = 1;
407 double tmpJacEvalTime;
408
409 infoStreamPrintWithEquationIndexes(LOG_LS, 0, indexes, "Start solving Linear System %d (size %d) at time %g with Total Pivot Solver",
410 eqSystemNumber, (int) systemData->size,
411 data->localData[0]->timeValue);
412
413 debugVectorDoubleLS(LOG_LS_V,"SCALING",systemData->nominal,n);
414 debugVectorDoubleLS(LOG_LS_V,"Old VALUES",aux_x,n);
415
416 rt_ext_tp_tick(&(solverData->timeClock));
417 if (0 == systemData->method) {
418
419 /* reset matrix A */
420 vecConstLS(n*n, 0.0, systemData->parDynamicData[omc_get_thread_num()].A);
421 /* update matrix A -> first n columns of matrix Ab*/
422 systemData->setA(data, threadData, systemData);
423 vecCopyLS(n*n, systemData->parDynamicData[omc_get_thread_num()].A, solverData->Ab);
424
425 /* update vector b (rhs) -> -b is last column of matrix Ab*/
426 rt_ext_tp_tick(&(solverData->timeClock));
427 systemData->setb(data, threadData, systemData);
428 vecScalarMultLS(n, systemData->parDynamicData[omc_get_thread_num()].b, -1.0, solverData->Ab + n*n);
429 } else {
430
431 /* calculate jacobian -> first n columns of matrix Ab*/
432 if(systemData->jacobianIndex != -1){
433 getAnalyticalJacobianTotalPivot(data, threadData, solverData->Ab, sysNumber);
434 } else {
435 assertStreamPrint(threadData, 1, "jacobian function pointer is invalid" )if (!(1)) {throwStreamPrint((threadData), "jacobian function pointer is invalid"
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "simulation/solver/linearSolverTotalPivot.c"
, 435, __extension__ __PRETTY_FUNCTION__); }));}
;
436 }
437 /* calculate vector b (rhs) -> -b is last column of matrix Ab */
438 wrapper_fvec_totalpivot(aux_x, solverData->Ab + n*n, dataAndThreadData, sysNumber);
439 }
440 tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
441 systemData->jacobianTime += tmpJacEvalTime;
442 infoStreamPrint(LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
443 debugMatrixDoubleLS(LOG_LS_V,"LGS: matrix Ab",solverData->Ab, n, n+1);
444
445 rt_ext_tp_tick(&(solverData->timeClock));
446 status = solveSystemWithTotalPivotSearchLS(n, solverData->x, solverData->Ab, solverData->indRow, solverData->indCol, &rank);
447 infoStreamPrint(LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
448
449 if (status != 0) {
450 // ToDo Rework stream prints like this one to work in parallel regions
451#ifdef USE_PARJAC
452 warningStreamPrint(LOG_STDOUT, 0, "Thread %u: Error solving linear system of equations (no. %d) at time %f.", omc_get_thread_num(), (int)systemData->equationIndex, data->localData[0]->timeValue);
453 success = 0;
454#else
455 warningStreamPrint(LOG_STDOUT, 0, "Error solving linear system of equations (no. %d) at time %f.", (int)systemData->equationIndex, data->localData[0]->timeValue);
456 success = 0;
457#endif
458 } else {
459 debugVectorDoubleLS(LOG_LS_V, "SOLUTION:", solverData->x, n+1);
460 if (1 == systemData->method) {
461 /* add the solution to old solution vector*/
462 vecAddLS(n, aux_x, solverData->x, aux_x);
463 wrapper_fvec_totalpivot(aux_x, solverData->b, dataAndThreadData, sysNumber);
464 } else {
465 /* take the solution */
466 vecCopyLS(n, solverData->x, aux_x);
467 }
468
469 if (ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V])) {
470 if (1 == systemData->method) {
471 infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
472 } else {
473 infoStreamPrint(LOG_LS_V, 1, "Solution x:");
474 }
475 infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
476 for(i=0; i<systemData->size; ++i)
477 {
478 infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
479 }
480 messageClose(LOG_LS_V);
481 }
482 }
483 return success;
484}