Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/linearSolverLis.c
Warning:line 154, 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 linearSolverLis.c
32 */
33
34#include <math.h>
35#include <stdlib.h>
36#include <string.h>
37
38#include "simulation_data.h"
39#include "simulation/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 "linearSolverLis.h"
48
49/*! \fn allocate memory for linear system solver Lis
50 *
51 */
52int
53allocateLisData(int n_row, int n_col, int nz, void** voiddata)
54{
55 DATA_LIS* data = (DATA_LIS*) malloc(sizeof(DATA_LIS));
56 char buffer[128];
57 assertStreamPrint(NULL, 0 != data, "Could not allocate data for linear solver Lis.")if (!(0 != data)) {throwStreamPrint((((void*)0)), "Could not allocate data for linear solver Lis."
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "./simulation/solver/linearSolverLis.c",
57, __extension__ __PRETTY_FUNCTION__); }));}
;
58
59 data->n_col = n_col;
60 data->n_row = n_row;
61 data->nnz = nz;
62
63 lis_vector_create(LIS_COMM_WORLD((LIS_Comm)0x1), &(data->b));
64 lis_vector_set_size(data->b, data->n_row, 0);
65
66 lis_vector_create(LIS_COMM_WORLD((LIS_Comm)0x1), &(data->x));
67 lis_vector_set_size(data->x, data->n_row, 0);
68
69 lis_matrix_create(LIS_COMM_WORLD((LIS_Comm)0x1), &(data->A));
70 lis_matrix_set_size(data->A, data->n_row, 0);
71 lis_matrix_set_type(data->A, LIS_MATRIX_CSR1);
72
73 lis_solver_create(&(data->solver));
74
75 lis_solver_set_option("-print none", data->solver);
76 sprintf(buffer,"-maxiter %d", n_row*100);
77 lis_solver_set_option(buffer, data->solver);
78 lis_solver_set_option("-scale none", data->solver);
79 lis_solver_set_option("-p none", data->solver);
80 lis_solver_set_option("-initx_zeros 0", data->solver);
81 lis_solver_set_option("-tol 1.0e-12", data->solver);
82
83 data->work = (double*) calloc(n_col,sizeof(double));
84
85 rt_ext_tp_tick(&(data->timeClock));
86
87
88 *voiddata = (void*)data;
89 return 0;
90}
91
92
93/*! \fn free memory for linear system solver Lis
94 *
95 */
96int freeLisData(void **voiddata)
97{
98 DATA_LIS* data = (DATA_LIS*) *voiddata;
99
100 lis_matrix_destroy(data->A);
101 lis_vector_destroy(data->b);
102 lis_vector_destroy(data->x);
103 lis_solver_destroy(data->solver);
104
105 free(data->work);
106
107 return 0;
108}
109
110
111/*!
112 * Print LIS_MATRIX provided in column sparse row format
113 *
114 * \param [in] A Matrix A of linear problem A*x=b in CSR format
115 * \param [in] n Dimension of A
116 */
117void printLisMatrixCSR(LIS_MATRIX A, int n)
118{
119 int i, j;
120 /* A matrix */
121 infoStreamPrint(LOG_LS_V, 1, "A matrix [%dx%d] nnz = %d", n, n, A->nnz);
122 infoStreamPrint(LOG_LS_V, 0, "Column Sparse Row format. Print tuple (index,value) for each row:");
123 for(i=0; i<n; i++)
124 {
125 char *buffer = (char*)malloc(sizeof(char)*A->ptr[i+1]*50);
126 buffer[0] = 0;
127 sprintf(buffer, "column %d: ", i);
128 for(j = A->ptr[i]; j < A->ptr[i+1]; j++)
129 {
130 sprintf(buffer, "%s(%d,%g) ", buffer, A->index[j], A->value[j]);
131 }
132 infoStreamPrint(LOG_LS_V, 0, "%s", buffer);
133 free(buffer);
134 }
135
136 messageClose(LOG_LS_V);
137}
138
139/*! \fn getAnalyticalJacobian
140 *
141 * function calculates analytical jacobian
142 *
143 * \param [ref] [data]
144 * \param [in] [sysNumber]
145 *
146 * \author wbraun
147 *
148 */
149int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber)
150{
151 int i,j,k,l,ii;
152 LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]);
153
154 const int index = systemData->jacobianIndex;
Value stored to 'index' during its initialization is never read
155 ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian;
156 ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian;
157
158 int nth = 0;
159 int nnz = jacobian->sparsePattern->numberOfNoneZeros;
160
161 for(i=0; i < jacobian->sizeRows; i++) {
162 jacobian->seedVars[i] = 1;
163 ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian);
164
165 for(j = 0; j < jacobian->sizeCols; j++) {
166 if(jacobian->seedVars[j] == 1) {
167 ii = jacobian->sparsePattern->leadindex[j];
168 while(ii < jacobian->sparsePattern->leadindex[j+1]) {
169 l = jacobian->sparsePattern->index[ii];
170 /*infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", l, i, nth, -jacobian->resultVars[l]); */
171 systemData->setAElement(l, i, -jacobian->resultVars[l], nth, (void*) systemData, threadData);
172 nth++;
173 ii++;
174 };
175 }
176 }
177 jacobian->seedVars[i] = 0;
178 }
179
180 return 0;
181}
182
183/*! \fn wrapper_fvec_lis for the residual function
184 *
185 */
186static int wrapper_fvec_lis(double* x, double* f, void** data, int sysNumber)
187{
188 int iflag = 0;
189
190 (*((DATA*)data[0])->simulationInfo->linearSystemData[sysNumber].residualFunc)(data, x, f, &iflag);
191 return 0;
192}
193
194
195/*! \fn solve linear system with Lis method
196 *
197 * \param [in] [data]
198 * [sysNumber] index of the corresponding linear system
199 *
200 */
201int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
202{
203 void *dataAndThreadData[2] = {data, threadData};
204 LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
205 DATA_LIS* solverData = (DATA_LIS*)systemData->parDynamicData[omc_get_thread_num()].solverData[0];
206
207 int i, ret, success = 1, ni, iflag = 1, n = systemData->size, eqSystemNumber = systemData->equationIndex;
208 char *lis_returncode[] = {"LIS_SUCCESS", "LIS_ILL_OPTION", "LIS_BREAKDOWN", "LIS_OUT_OF_MEMORY", "LIS_MAXITER", "LIS_NOT_IMPLEMENTED", "LIS_ERR_FILE_IO"};
209 LIS_INT err;
210 _omc_scalar residualNorm = 0;
211
212 int indexes[2] = {1,eqSystemNumber};
213 double tmpJacEvalTime;
214 infoStreamPrintWithEquationIndexes(LOG_LS, 0, indexes, "Start solving Linear System %d (size %d) at time %g with Lis Solver",
215 eqSystemNumber, (int) systemData->size,
216 data->localData[0]->timeValue);
217
218 /* set old values as start value for the iteration */
219 for(i=0; i<n; i++){
220 err = lis_vector_set_value(LIS_INS_VALUE0, i, aux_x[i], solverData->x);
221 }
222
223 rt_ext_tp_tick(&(solverData->timeClock));
224
225 lis_matrix_set_size(solverData->A, solverData->n_row, 0);
226 if (0 == systemData->method)
227 {
228 /* set A matrix */
229 systemData->setA(data, threadData, systemData);
230 lis_matrix_assemble(solverData->A);
231
232 /* set b vector */
233 systemData->setb(data, threadData, systemData);
234
235 } else {
236 /* calculate jacobian -> matrix A*/
237 if(systemData->jacobianIndex != -1){
238 getAnalyticalJacobianLis(data, threadData, sysNumber);
239 } else {
240 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/linearSolverLis.c",
240, __extension__ __PRETTY_FUNCTION__); }));}
;
241 }
242 lis_matrix_assemble(solverData->A);
243
244 /* calculate vector b (rhs) */
245 memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row);
246 wrapper_fvec_lis(solverData->work, systemData->parDynamicData[omc_get_thread_num()].b, dataAndThreadData, sysNumber);
247
248 /* set b vector */
249 for(i=0; i<n; i++) {
250 err = lis_vector_set_value(LIS_INS_VALUE0, i, systemData->parDynamicData[omc_get_thread_num()].b[i], solverData->b);
251 }
252 }
253 tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
254 systemData->jacobianTime += tmpJacEvalTime;
255 infoStreamPrint(LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
256
257 rt_ext_tp_tick(&(solverData->timeClock));
258 err = lis_solve(solverData->A,solverData->b,solverData->x,solverData->solver);
259 infoStreamPrint(LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
260
261 if (err){
262 warningStreamPrint(LOG_LS_V, 0, "lis_solve : %s(code=%d)\n\n ", lis_returncode[err], err);
263 printLisMatrixCSR(solverData->A, solverData->n_row);
264 success = 0;
265 }
266
267
268 /* Log A*x=b */
269 if(ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V]))
270 {
271 char *buffer = (char*)malloc(sizeof(char)*n*25);
272
273 printLisMatrixCSR(solverData->A, n);
274
275 /* b vector */
276 infoStreamPrint(LOG_LS_V, 1, "b vector [%d]", n);
277 for(i=0; i<n; i++)
278 {
279 buffer[0] = 0;
280 sprintf(buffer, "%s%20.12g ", buffer, solverData->b->value[i]);
281 infoStreamPrint(LOG_LS_V, 0, "%s", buffer);
282 }
283 messageClose(LOG_LS_V);
284 free(buffer);
285 }
286
287 /* Log solution */
288 if (1 == success){
289
290 if (1 == systemData->method){ /* Case calculate jacobian -> matrix A*/
291 /* take the solution */
292 lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x);
293 for(i = 0; i < solverData->n_row; ++i)
294 aux_x[i] += solverData->work[i];
295
296 /* update inner equations */
297 wrapper_fvec_lis(aux_x, solverData->work, dataAndThreadData, sysNumber);
298 residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
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 /* write solution */
308 lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x);
309 }
310
311 if (ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V]))
312 {
313 if (1 == systemData->method) {
314 infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
315 } else {
316 infoStreamPrint(LOG_LS_V, 1, "Solution x:");
317 }
318 infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
319
320 for(i = 0; i < systemData->size; ++i)
321 infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
322
323 messageClose(LOG_LS_V);
324 }
325 }
326 else
327 {
328 warningStreamPrint(LOG_STDOUT, 0,
329 "Failed to solve linear system of equations (no. %d) at time %f, system status %d.",
330 (int)systemData->equationIndex, data->localData[0]->timeValue, err);
331 }
332
333 return success;
334}