Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/linearSolverKlu.c
Warning:line 130, column 7
Value stored to 'nnz' 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 linearSolverKlu.c
32 */
33
34#include "omc_config.h"
35
36#ifdef WITH_UMFPACK
37#include <math.h>
38#include <stdlib.h>
39#include <string.h>
40
41#include "simulation_data.h"
42#include "simulation/simulation_info_json.h"
43#include "util/omc_error.h"
44#include "util/parallel_helper.h"
45#include "omc_math.h"
46#include "util/varinfo.h"
47#include "model_help.h"
48
49#include "linearSystem.h"
50#include "linearSolverKlu.h"
51
52static void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n);
53static void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n);
54
55/*! \fn allocate memory for linear system solver Klu
56 *
57 */
58int allocateKluData(int n_row, int n_col, int nz, void** voiddata)
59{
60 DATA_KLU* data = (DATA_KLU*) malloc(sizeof(DATA_KLU));
61 assertStreamPrint(NULL, 0 != data, "Could not allocate data for linear solver Klu.")if (!(0 != data)) {throwStreamPrint((((void*)0)), "Could not allocate data for linear solver Klu."
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "./simulation/solver/linearSolverKlu.c",
61, __extension__ __PRETTY_FUNCTION__); }));}
;
62
63 data->symbolic = NULL((void*)0);
64 data->numeric = NULL((void*)0);
65
66 data->n_col = n_col;
67 data->n_row = n_row;
68 data->nnz = nz;
69
70 data->Ap = (int*) calloc((n_row+1),sizeof(int));
71 data->Ai = (int*) calloc(nz,sizeof(int));
72 data->Ax = (double*) calloc(nz,sizeof(double));
73 data->work = (double*) calloc(n_col,sizeof(double));
74
75 data->numberSolving = 0;
76 klu_defaults(&(data->common));
77
78 *voiddata = (void*)data;
79
80 return 0;
81}
82
83
84/*! \fn free memory for linear system solver Klu
85 *
86 */
87int freeKluData(void **voiddata)
88{
89 TRACE_PUSH
90
91 DATA_KLU* data = (DATA_KLU*) *voiddata;
92
93 free(data->Ap);
94 free(data->Ai);
95 free(data->Ax);
96 free(data->work);
97
98
99 if(data->symbolic)
100 klu_free_symbolic(&data->symbolic, &data->common);
101 if(data->numeric)
102 klu_free_numeric(&data->numeric, &data->common);
103
104 TRACE_POP
105 return 0;
106}
107
108/*! \fn getAnalyticalJacobian
109 *
110 * function calculates analytical jacobian
111 *
112 * \param [ref] [data]
113 * \param [in] [sysNumber]
114 *
115 * \author wbraun
116 *
117 */
118static int getAnalyticalJacobian(DATA* data, threadData_t *threadData,
119 int sysNumber)
120{
121 int i,ii,j,k,l;
122
123 LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]);
124
125 const int index = systemData->jacobianIndex;
126 ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian;
127 ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian;
128
129 int nth = 0;
130 int nnz = jacobian->sparsePattern->numberOfNoneZeros;
Value stored to 'nnz' during its initialization is never read
131
132 if (jacobian->constantEqns != NULL((void*)0)) {
133 jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
134 }
135
136 for(i=0; i < jacobian->sparsePattern->maxColors; i++)
137 {
138 /* activate seed variable for the corresponding color */
139 for(ii=0; ii < jacobian->sizeCols; ii++)
140 {
141 if(jacobian->sparsePattern->colorCols[ii]-1 == i)
142 {
143 jacobian->seedVars[ii] = 1;
144 }
145 }
146
147 ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian);
148
149 for(j = 0; j < jacobian->sizeCols; j++)
150 {
151 if(jacobian->seedVars[j] == 1)
152 {
153 nth = jacobian->sparsePattern->leadindex[j];
154 while(nth < jacobian->sparsePattern->leadindex[j+1])
155 {
156 l = jacobian->sparsePattern->index[nth];
157 systemData->setAElement(j, l, -jacobian->resultVars[l], nth, (void*) systemData, threadData);
158 nth++;
159 }
160 /* de-activate seed variable for the corresponding color */
161 jacobian->seedVars[j] = 0;
162 }
163 }
164 }
165 return 0;
166}
167
168/*! \fn residual_wrapper for the residual function
169 *
170 */
171static int residual_wrapper(double* x, double* f, void** data, int sysNumber)
172{
173 int iflag = 0;
174
175 (*((DATA*)data[0])->simulationInfo->linearSystemData[sysNumber].residualFunc)(data, x, f, &iflag);
176 return 0;
177}
178
179/*! \fn solve linear system with Klu method
180 *
181 * \param [in] [data]
182 * [sysNumber] index of the corresponding linear system
183 *
184 *
185 * author: wbraun
186 */
187int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
188{
189 void *dataAndThreadData[2] = {data, threadData};
190 LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
191 DATA_KLU* solverData = (DATA_KLU*)systemData->parDynamicData[omc_get_thread_num()].solverData[0];
192 _omc_scalar residualNorm = 0;
193
194 int i, j, status = 0, success = 0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
195 double tmpJacEvalTime;
196 int reuseMatrixJac = (data->simulationInfo->currentContext == CONTEXT_SYM_JACOBIAN && data->simulationInfo->currentJacobianEval > 0);
197
198 infoStreamPrintWithEquationIndexes(LOG_LS, 0, indexes, "Start solving Linear System %d (size %d) at time %g with Klu Solver",
199 eqSystemNumber, (int) systemData->size,
200 data->localData[0]->timeValue);
201
202 rt_ext_tp_tick(&(solverData->timeClock));
203 if (0 == systemData->method)
204 {
205 if (!reuseMatrixJac){
206 /* set A matrix */
207 solverData->Ap[0] = 0;
208 systemData->setA(data, threadData, systemData);
209 solverData->Ap[solverData->n_row] = solverData->nnz;
210 }
211
212 /* set b vector */
213 systemData->setb(data, threadData, systemData);
214 } else {
215
216 if (!reuseMatrixJac){
217 solverData->Ap[0] = 0;
218 /* calculate jacobian -> matrix A*/
219 if(systemData->jacobianIndex != -1){
220 getAnalyticalJacobian(data, threadData, sysNumber);
221 } else {
222 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/linearSolverKlu.c",
222, __extension__ __PRETTY_FUNCTION__); }));}
;
223 }
224 solverData->Ap[solverData->n_row] = solverData->nnz;
225 }
226
227 /* calculate vector b (rhs) */
228 memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row);
229
230 residual_wrapper(solverData->work, systemData->parDynamicData[omc_get_thread_num()].b, dataAndThreadData, sysNumber);
231 }
232 tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
233 systemData->jacobianTime += tmpJacEvalTime;
234 infoStreamPrint(LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
235
236 if (ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V]))
237 {
238 infoStreamPrint(LOG_LS_V, 1, "Old solution x:");
239 for(i = 0; i < solverData->n_row; ++i)
240 infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
241 messageClose(LOG_LS_V);
242
243 infoStreamPrint(LOG_LS_V, 1, "Matrix A n_rows = %d", solverData->n_row);
244 for (i=0; i<solverData->n_row; i++){
245 infoStreamPrint(LOG_LS_V, 0, "%d. Ap => %d -> %d", i, solverData->Ap[i], solverData->Ap[i+1]);
246 for (j=solverData->Ap[i]; j<solverData->Ap[i+1]; j++){
247 infoStreamPrint(LOG_LS_V, 0, "A[%d,%d] = %f", i, solverData->Ai[j], solverData->Ax[j]);
248 }
249 }
250 messageClose(LOG_LS_V);
251
252 for (i=0; i<solverData->n_row; i++)
253 {
254 // ToDo Rework stream prints like this one to work in parallel regions
255 infoStreamPrint(LOG_LS_V, 0, "b[%d] = %e", i, systemData->parDynamicData[omc_get_thread_num()].b[i]);
256 }
257 }
258 rt_ext_tp_tick(&(solverData->timeClock));
259
260 /* symbolic pre-ordering of A to reduce fill-in of L and U */
261 if (0 == solverData->numberSolving)
262 {
263 infoStreamPrint(LOG_LS_V, 0, "Perform analyze settings:\n - ordering used: %d\n - current status: %d", solverData->common.ordering, solverData->common.status);
264 solverData->symbolic = klu_analyze(solverData->n_col, solverData->Ap, solverData->Ai, &solverData->common);
265 }
266
267 /* if reuseMatrixJac use also previous factorization */
268 if (!reuseMatrixJac)
269 {
270 /* compute the LU factorization of A */
271 if (0 == solverData->common.status){
272 if(solverData->numeric){
273 /* Just refactor using the same pivots, but check that the refactor is still accurate */
274 klu_refactor(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, solverData->numeric, &solverData->common);
275 klu_rgrowth(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, solverData->numeric, &solverData->common);
276 infoStreamPrint(LOG_LS_V, 0, "Klu rgrowth after refactor: %f", solverData->common.rgrowth);
277 /* If rgrowth is small then do a whole factorization with new pivots (What should this tolerance be?) */
278 if (solverData->common.rgrowth < 1e-3){
279 klu_free_numeric(&solverData->numeric, &solverData->common);
280 solverData->numeric = klu_factor(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, &solverData->common);
281 infoStreamPrint(LOG_LS_V, 0, "Klu new factorization performed.");
282 }
283 } else {
284 solverData->numeric = klu_factor(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, &solverData->common);
285 }
286 }
287 }
288
289 if (0 == solverData->common.status){
290 if (1 == systemData->method){
291 if (klu_solve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->parDynamicData[omc_get_thread_num()].b, &solverData->common)){
292 success = 1;
293 }
294 } else {
295 if (klu_tsolve(solverData->symbolic, solverData->numeric, solverData->n_col, 1, systemData->parDynamicData[omc_get_thread_num()].b, &solverData->common)){
296 success = 1;
297 }
298 }
299 }
300
301 infoStreamPrint(LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
302
303 /* print solution */
304 if (1 == success){
305
306 if (1 == systemData->method){
307 /* take the solution */
308 for(i = 0; i < solverData->n_row; ++i)
309 aux_x[i] += systemData->parDynamicData[omc_get_thread_num()].b[i];
310
311 /* update inner equations */
312 residual_wrapper(aux_x, solverData->work, dataAndThreadData, sysNumber);
313 residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
314
315 if ((isnan(residualNorm)(sizeof ((residualNorm)) == sizeof (float) ? __isnanf (residualNorm
) : sizeof ((residualNorm)) == sizeof (double) ? __isnan (residualNorm
) : __isnanl (residualNorm))
) || (residualNorm>1e-4)){
316 warningStreamPrint(LOG_LS, 0,
317 "Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
318 (int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
319 success = 0;
320 }
321 } else {
322 /* the solution is automatically in x */
323 memcpy(aux_x, systemData->parDynamicData[omc_get_thread_num()].b, sizeof(double)*systemData->size);
324 }
325
326 if (ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V]))
327 {
328 if (1 == systemData->method) {
329 infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
330 } else {
331 infoStreamPrint(LOG_LS_V, 1, "Solution x:");
332 }
333 infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
334
335 for(i = 0; i < systemData->size; ++i)
336 infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
337
338 messageClose(LOG_LS_V);
339 }
340 }
341 else
342 {
343 warningStreamPrint(LOG_STDOUT, 0,
344 "Failed to solve linear system of equations (no. %d) at time %f, system status %d.",
345 (int)systemData->equationIndex, data->localData[0]->timeValue, status);
346 }
347 solverData->numberSolving += 1;
348
349 return success;
350}
351
352static
353void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n)
354{
355 int i, j, k, l;
356
357 char **buffer = (char**)malloc(sizeof(char*)*n);
358 for (l=0; l<n; l++)
359 {
360 buffer[l] = (char*)malloc(sizeof(char)*n*20);
361 buffer[l][0] = 0;
362 }
363
364 k = 0;
365 for (i = 0; i < n; i++)
366 {
367 for (j = 0; j < n; j++)
368 {
369 if ((k < Ap[i + 1]) && (Ai[k] == j))
370 {
371 sprintf(buffer[j], "%s %5g ", buffer[j], Ax[k]);
372 k++;
373 }
374 else
375 {
376 sprintf(buffer[j], "%s %5g ", buffer[j], 0.0);
377 }
378 }
379 }
380 for (l = 0; l < n; l++)
381 {
382 infoStreamPrint(LOG_LS_V, 0, "%s", buffer[l]);
383 free(buffer[l]);
384 }
385 free(buffer);
386}
387
388static
389void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n)
390{
391 int i, j, k;
392 char *buffer = (char*)malloc(sizeof(char)*n*15);
393 k = 0;
394 for (i = 0; i < n; i++)
395 {
396 buffer[0] = 0;
397 for (j = 0; j < n; j++)
398 {
399 if ((k < Ap[i + 1]) && (Ai[k] == j))
400 {
401 sprintf(buffer, "%s %5.2g ", buffer, Ax[k]);
402 k++;
403 }
404 else
405 {
406 sprintf(buffer, "%s %5.2g ", buffer, 0.0);
407 }
408 }
409 infoStreamPrint(LOG_LS_V, 0, "%s", buffer);
410 }
411 free(buffer);
412}
413
414#endif