Bug Summary

File:OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverLapack.c
Warning:line 114, 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 "linearSolverLapack.h"
48
49#ifdef USE_PARJAC
50 #include <omp.h>
51#endif
52
53extern int dgesv_(int *n, int *nrhs, double *a, int *lda,
54 int *ipiv, double *b, int *ldb, int *info);
55
56extern int dgetrs_(char* tran, int *n, int *nrhs, double *a, int *lda,
57 int *ipiv, double *b, int *ldb, int *info);
58/*! \fn allocate memory for linear system solver lapack
59 *
60 */
61int allocateLapackData(int size, void** voiddata)
62{
63 DATA_LAPACK* data = (DATA_LAPACK*) calloc(1, sizeof(DATA_LAPACK));
64
65 data->ipiv = (int*) calloc(size, sizeof(int));
66 assertStreamPrint(NULL, 0 != data->ipiv, "Could not allocate data for linear solver lapack.")if (!(0 != data->ipiv)) {throwStreamPrint((((void*)0)), "Could not allocate data for linear solver lapack."
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "simulation/solver/linearSolverLapack.c"
, 66, __extension__ __PRETTY_FUNCTION__); }));}
;
67 data->nrhs = 1;
68 data->info = 0;
69 data->work = _omc_allocateVectorData(size);
70
71 data->x = _omc_createVector(size, NULL((void*)0));
72 data->b = _omc_createVector(size, NULL((void*)0));
73 data->A = _omc_createMatrix(size, size, NULL((void*)0));
74
75 *voiddata = (void*)data;
76 return 0;
77}
78
79/*! \fn free memory of lapack
80 *
81 */
82int freeLapackData(void **voiddata)
83{
84 DATA_LAPACK* data = (DATA_LAPACK*) *voiddata;
85
86 free(data->ipiv);
87 _omc_deallocateVectorData(data->work);
88
89 _omc_destroyVector(data->x);
90 _omc_destroyVector(data->b);
91 _omc_destroyMatrix(data->A);
92
93 free(data);
94 voiddata[0] = NULL((void*)0);
95
96 return 0;
97}
98
99/*! \fn getAnalyticalJacobian
100 *
101 * function calculates analytical jacobian
102 *
103 * \param [ref] [data]
104 * \param [out] [jac]
105 *
106 * \author wbraun
107 *
108 */
109int getAnalyticalJacobianLapack(DATA* data, threadData_t *threadData, double* jac, int sysNumber)
110{
111 int i,j,k,l,ii,currentSys = sysNumber;
112 LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[currentSys]);
113
114 const int index = systemData->jacobianIndex;
Value stored to 'index' during its initialization is never read
115 ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian;
116 ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian;
117
118 memset(jac, 0, (systemData->size)*(systemData->size)*sizeof(double));
119
120 if (jacobian->constantEqns != NULL((void*)0)) {
121 jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
122 }
123
124 for(i=0; i < jacobian->sparsePattern->maxColors; i++)
125 {
126 /* activate seed variable for the corresponding color */
127 for(ii=0; ii < jacobian->sizeCols; ii++)
128 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
129 jacobian->seedVars[ii] = 1;
130
131 ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian);
132
133 for(j = 0; j < jacobian->sizeCols; j++) {
134 if(jacobian->seedVars[j] == 1) {
135 ii = jacobian->sparsePattern->leadindex[j];
136 while(ii < jacobian->sparsePattern->leadindex[j+1]) {
137 l = jacobian->sparsePattern->index[ii];
138 k = j*jacobian->sizeRows + l;
139 jac[k] = -jacobian->resultVars[l];
140 ii++;
141 };
142 }
143 /* de-activate seed variable for the corresponding color */
144 if(jacobian->sparsePattern->colorCols[j]-1 == i) {
145 jacobian->seedVars[j] = 0;
146 }
147 }
148 }
149
150 return 0;
151}
152
153/*! \fn wrapper_fvec_lapack for the residual function
154 *
155 */
156static int wrapper_fvec_lapack(_omc_vector* x, _omc_vector* f, int* iflag, void** data, int sysNumber)
157{
158 int currentSys = sysNumber;
159
160 (*((DATA*)data[0])->simulationInfo->linearSystemData[currentSys].residualFunc)(data, x->data, f->data, iflag);
161 return 0;
162}
163
164/*! \fn solve linear system with lapack method
165 *
166 * \param [in] [data]
167 * [sysNumber] index of the corresponding linear system
168 *
169 * \author wbraun
170 */
171int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
172{
173 void *dataAndThreadData[2] = {data, threadData};
174 int i, iflag = 1;
175 LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
176
177 DATA_LAPACK* solverData = (DATA_LAPACK*) systemData->parDynamicData[omc_get_thread_num()].solverData[0];
178 int success = 1;
179
180 /* We are given the number of the linear system.
181 * We want to look it up among all equations. */
182 int eqSystemNumber = systemData->equationIndex;
183 int indexes[2] = {1,eqSystemNumber};
184 _omc_scalar residualNorm = 0;
185 double tmpJacEvalTime;
186 int reuseMatrixJac = (data->simulationInfo->currentContext == CONTEXT_SYM_JACOBIAN && data->simulationInfo->currentJacobianEval > 0);
187
188 infoStreamPrintWithEquationIndexes(LOG_LS, 0, indexes, "Start solving Linear System %d (size %d) at time %g with Lapack Solver",
189 eqSystemNumber, (int) systemData->size,
190 data->localData[0]->timeValue);
191
192 /* set data */
193 _omc_setVectorData(solverData->x, aux_x);
194 _omc_setVectorData(solverData->b, systemData->parDynamicData[omc_get_thread_num()].b);
195 _omc_setMatrixData(solverData->A, systemData->parDynamicData[omc_get_thread_num()].A);
196
197 // ToDo: Make time variables thread safe as this can be called in a parallel region
198 rt_ext_tp_tick(&(solverData->timeClock));
199 if (0 == systemData->method) {
200
201 if (!reuseMatrixJac) {
202 /* reset matrix A */
203 memset(systemData->parDynamicData[omc_get_thread_num()].A, 0, (systemData->size)*(systemData->size)*sizeof(double));
204 /* update matrix A */
205 systemData->setA(data, threadData, systemData);
206 }
207
208 /* update vector b (rhs) */
209 systemData->setb(data, threadData, systemData);
210 } else {
211 if (!reuseMatrixJac) {
212 /* calculate jacobian -> matrix A*/
213 if(systemData->jacobianIndex != -1) {
214 getAnalyticalJacobianLapack(data, threadData, solverData->A->data, sysNumber);
215 } else {
216 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/linearSolverLapack.c"
, 216, __extension__ __PRETTY_FUNCTION__); }));}
;
217 }
218 }
219 /* calculate vector b (rhs) */
220 _omc_copyVector(solverData->work, solverData->x);
221
222 wrapper_fvec_lapack(solverData->work, solverData->b, &iflag, dataAndThreadData, sysNumber);
223 }
224 tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
225 systemData->jacobianTime += tmpJacEvalTime;
226 infoStreamPrint(LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
227
228 /* Log A*x=b */
229 if(ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V])){
230 _omc_printVector(solverData->x, "Vector old x", LOG_LS_V);
231 _omc_printMatrix(solverData->A, "Matrix A", LOG_LS_V);
232 _omc_printVector(solverData->b, "Vector b", LOG_LS_V);
233 }
234
235 rt_ext_tp_tick(&(solverData->timeClock));
236
237 /* if reuseMatrixJac use also previous factorization */
238 if (!reuseMatrixJac)
239 {
240 /* Solve system */
241 dgesv_((int*) &systemData->size,
242 (int*) &solverData->nrhs,
243 solverData->A->data,
244 (int*) &systemData->size,
245 solverData->ipiv,
246 solverData->b->data,
247 (int*) &systemData->size,
248 &solverData->info);
249
250 } /* further Jacobian evaluations */
251 else
252 {
253 char trans = 'N';
254 /* Solve system */
255 dgetrs_(&trans,
256 (int*) &systemData->size,
257 (int*) &solverData->nrhs,
258 solverData->A->data,
259 (int*) &systemData->size,
260 solverData->ipiv,
261 solverData->b->data,
262 (int*) &systemData->size,
263 &solverData->info);
264 }
265
266
267 infoStreamPrint(LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
268
269 if(solverData->info < 0)
270 {
271 warningStreamPrint(LOG_LS, 0, "Error solving linear system of equations (no. %d) at time %f. Argument %d illegal.", (int)systemData->equationIndex, data->localData[0]->timeValue, (int)solverData->info);
272 success = 0;
273 }
274 else if(solverData->info > 0)
275 {
276 warningStreamPrint(LOG_LS, 0,
277 "Failed to solve linear system of equations (no. %d) at time %f, system is singular for U[%d, %d].",
278 (int)systemData->equationIndex, data->localData[0]->timeValue, (int)solverData->info+1, (int)solverData->info+1);
279
280 success = 0;
281
282 /* debug output */
283 if (ACTIVE_STREAM(LOG_LS)(useStream[LOG_LS])){
284 _omc_printMatrix(solverData->A, "Matrix U", LOG_LS);
285
286 _omc_printVector(solverData->b, "Output vector x", LOG_LS);
287 }
288 }
289
290 if (1 == success){
291
292 if (1 == systemData->method){
293 /* take the solution */
294 solverData->x = _omc_addVectorVector(solverData->x, solverData->work, solverData->b); // x = xold(work) + xnew(b)
295
296 /* update inner equations */
297 wrapper_fvec_lapack(solverData->x, solverData->work, &iflag, dataAndThreadData, sysNumber);
298 residualNorm = _omc_euclideanVectorNorm(solverData->work);
299
300 if ((isnan(residualNorm)(sizeof ((residualNorm)) == sizeof (float) ? __isnanf (residualNorm
) : sizeof ((residualNorm)) == sizeof (double) ? __isnan (residualNorm
) : __isnanl (residualNorm))
) || (residualNorm>1e-4)){
301 warningStreamPrint(LOG_LS, 0,
302 "Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
303 (int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
304 success = 0;
305 }
306 } else {
307 /* take the solution */
308 _omc_copyVector(solverData->x, solverData->b);
309 }
310
311 if (ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V])){
312 if (1 == systemData->method) {
313 infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
314 } else {
315 infoStreamPrint(LOG_LS_V, 1, "Solution x:");
316 }
317 infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
318
319 for(i = 0; i < systemData->size; ++i) {
320 infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %.15g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
321 }
322
323 messageClose(LOG_LS_V);
324 }
325 }
326
327 return success;
328}