Bug Summary

File:OMCompiler/SimulationRuntime/c/./simulation/solver/linearSolverUmfpack.c
Warning:line 140, 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 linearSolverUmfpack.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 "linearSolverUmfpack.h"
51
52void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n);
53void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n);
54int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x);
55
56/*! \fn allocate memory for linear system solver UmfPack
57 *
58 */
59int
60allocateUmfPackData(int n_row, int n_col, int nz, void** voiddata)
61{
62 DATA_UMFPACK* data = (DATA_UMFPACK*) malloc(sizeof(DATA_UMFPACK));
63 assertStreamPrint(NULL, 0 != data, "Could not allocate data for linear solver UmfPack.")if (!(0 != data)) {throwStreamPrint((((void*)0)), "Could not allocate data for linear solver UmfPack."
); ((void) sizeof ((0) ? 1 : 0), __extension__ ({ if (0) ; else
__assert_fail ("0", "./simulation/solver/linearSolverUmfpack.c"
, 63, __extension__ __PRETTY_FUNCTION__); }));}
;
64
65 data->symbolic = NULL((void*)0);
66 data->numeric = NULL((void*)0);
67
68 data->n_col = n_col;
69 data->n_row = n_row;
70 data->nnz = nz;
71
72
73 data->Ap = (int*) calloc((n_row+1),sizeof(int));
74
75 data->Ai = (int*) calloc(nz,sizeof(int));
76 data->Ax = (double*) calloc(nz,sizeof(double));
77 data->work = (double*) calloc(n_col,sizeof(double));
78
79 data->Wi = (int*) malloc(n_row * sizeof(int));
80 data->W = (double*) malloc(5*n_row * sizeof(double));
81
82 data->numberSolving=0;
83 umfpack_di_defaults(data->control);
84
85 data->control[UMFPACK_PIVOT_TOLERANCE3] = 0.1;
86 data->control[UMFPACK_IRSTEP7] = 2;
87 data->control[UMFPACK_SCALE16] = 1;
88 data->control[UMFPACK_STRATEGY5] = 5;
89
90
91
92 *voiddata = (void*)data;
93
94 return 0;
95}
96
97
98/*! \fn free memory for linear system solver UmfPack
99 *
100 */
101int
102freeUmfPackData(void **voiddata)
103{
104 TRACE_PUSH
105
106 DATA_UMFPACK* data = (DATA_UMFPACK*) *voiddata;
107
108 free(data->Ap);
109 free(data->Ai);
110 free(data->Ax);
111 free(data->work);
112
113 free(data->Wi);
114 free(data->W);
115
116 if(data->symbolic)
117 umfpack_di_free_symbolic (&data->symbolic);
118 if(data->numeric)
119 umfpack_di_free_numeric (&data->numeric);
120
121 TRACE_POP
122 return 0;
123}
124
125/*! \fn getAnalyticalJacobian
126 *
127 * function calculates analytical jacobian
128 *
129 * \param [ref] [data]
130 * \param [in] [sysNumber]
131 *
132 * \author wbraun
133 *
134 */
135int getAnalyticalJacobianUmfPack(DATA* data, threadData_t *threadData, int sysNumber)
136{
137 int i,ii,j,k,l;
138 LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]);
139
140 const int index = systemData->jacobianIndex;
Value stored to 'index' during its initialization is never read
141 ANALYTIC_JACOBIAN* jacobian = systemData->parDynamicData[omc_get_thread_num()].jacobian;
142 ANALYTIC_JACOBIAN* parentJacobian = systemData->parDynamicData[omc_get_thread_num()].parentJacobian;
143
144 int nth = 0;
145 int nnz = jacobian->sparsePattern->numberOfNoneZeros;
146
147 for(i=0; i < jacobian->sizeRows; i++)
148 {
149 jacobian->seedVars[i] = 1;
150
151 ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian);
152
153 for(j = 0; j < jacobian->sizeCols; j++)
154 {
155 if(jacobian->seedVars[j] == 1)
156 {
157 ii = jacobian->sparsePattern->leadindex[j];
158 while(ii < jacobian->sparsePattern->leadindex[j+1])
159 {
160 l = jacobian->sparsePattern->index[ii];
161 /* infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -jacobian->resultVars[l]); */
162 systemData->setAElement(i, l, -jacobian->resultVars[l], nth, (void*) systemData, threadData);
163 nth++;
164 ii++;
165 };
166 }
167 };
168
169 /* de-activate seed variable for the corresponding color */
170 jacobian->seedVars[i] = 0;
171 }
172
173 return 0;
174}
175
176/*! \fn wrapper_fvec_umfpack for the residual function
177 *
178 */
179static int wrapper_fvec_umfpack(double* x, double* f, void** data, int sysNumber)
180{
181 int iflag = 0;
182
183 (*((DATA*)data[0])->simulationInfo->linearSystemData[sysNumber].residualFunc)(data, x, f, &iflag);
184 return 0;
185}
186
187/*! \fn solve linear system with UmfPack method
188 *
189 * \param [in] [data]
190 * [sysNumber] index of the corresponding linear system
191 *
192 *
193 * author: kbalzereit, wbraun
194 */
195int
196solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
197{
198 void *dataAndThreadData[2] = {data, threadData};
199 LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
200 DATA_UMFPACK* solverData = (DATA_UMFPACK*)systemData->parDynamicData[omc_get_thread_num()].solverData[0];
201 _omc_scalar residualNorm = 0;
202
203 int i, j, status = UMFPACK_OK(0), success = 0, ni=0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
204 int casualTearingSet = systemData->strictTearingFunctionCall != NULL((void*)0);
205 double tmpJacEvalTime;
206 int reuseMatrixJac = (data->simulationInfo->currentContext == CONTEXT_SYM_JACOBIAN && data->simulationInfo->currentJacobianEval > 0);
207
208 infoStreamPrintWithEquationIndexes(LOG_LS, 0, indexes, "Start solving Linear System %d (size %d) at time %g with UMFPACK Solver",
209 eqSystemNumber, (int) systemData->size,
210 data->localData[0]->timeValue);
211
212 rt_ext_tp_tick(&(solverData->timeClock));
213 if (0 == systemData->method)
214 {
215 if (!reuseMatrixJac){
216 /* set A matrix */
217 solverData->Ap[0] = 0;
218 systemData->setA(data, threadData, systemData);
219 solverData->Ap[solverData->n_row] = solverData->nnz;
220 }
221
222 /* set b vector */
223 systemData->setb(data, threadData, systemData);
224 } else {
225
226 if (!reuseMatrixJac){
227 solverData->Ap[0] = 0;
228 /* calculate jacobian -> matrix A*/
229 if(systemData->jacobianIndex != -1){
230 getAnalyticalJacobianUmfPack(data, threadData, sysNumber);
231 } else {
232 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/linearSolverUmfpack.c"
, 232, __extension__ __PRETTY_FUNCTION__); }));}
;
233 }
234 solverData->Ap[solverData->n_row] = solverData->nnz;
235 }
236
237 /* calculate vector b (rhs) */
238 memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row);
239 wrapper_fvec_umfpack(solverData->work, systemData->parDynamicData[omc_get_thread_num()].b, dataAndThreadData, sysNumber);
240 }
241 tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
242 systemData->jacobianTime += tmpJacEvalTime;
243 infoStreamPrint(LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
244
245 if (ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V]))
246 {
247 infoStreamPrint(LOG_LS_V, 1, "Old solution x:");
248 for(i = 0; i < solverData->n_row; ++i)
249 infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
250 messageClose(LOG_LS_V);
251
252 infoStreamPrint(LOG_LS_V, 1, "Matrix A n_rows = %d", solverData->n_row);
253 for (i=0; i<solverData->n_row; i++){
254 infoStreamPrint(LOG_LS_V, 0, "%d. Ap => %d -> %d", i, solverData->Ap[i], solverData->Ap[i+1]);
255 for (j=solverData->Ap[i]; j<solverData->Ap[i+1]; j++){
256 infoStreamPrint(LOG_LS_V, 0, "A[%d,%d] = %f", i, solverData->Ai[j], solverData->Ax[j]);
257 }
258 }
259 messageClose(LOG_LS_V);
260
261 for (i=0; i<solverData->n_row; i++) {
262 // ToDo Rework stream prints like this one to work in parallel regions
263 infoStreamPrint(LOG_LS_V, 0, "b[%d] = %e", i, systemData->parDynamicData[omc_get_thread_num()].b[i]);
264 }
265 }
266 rt_ext_tp_tick(&(solverData->timeClock));
267
268 /* symbolic pre-ordering of A to reduce fill-in of L and U */
269 if (0 == solverData->numberSolving) {
270 status = umfpack_di_symbolic(solverData->n_col, solverData->n_row, solverData->Ap, solverData->Ai, solverData->Ax, &(solverData->symbolic), solverData->control, solverData->info);
271 }
272
273 /* compute the LU factorization of A */
274 /* if reuseMatrixJac use also previous factorization */
275 if (!reuseMatrixJac)
276 {
277 if (0 == status){
278 umfpack_di_free_numeric(&(solverData->numeric));
279 status = umfpack_di_numeric(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, &(solverData->numeric), solverData->control, solverData->info);
280 }
281 }
282
283 if (0 == status){
284 if (1 == systemData->method){
285 status = umfpack_di_wsolve(UMFPACK_A(0), solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->parDynamicData[omc_get_thread_num()].b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W);
286 } else {
287 status = umfpack_di_wsolve(UMFPACK_Aat(2), solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->parDynamicData[omc_get_thread_num()].b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W);
288 }
289 }
290
291 if (status == UMFPACK_OK(0)){
292 success = 1;
293 }
294 else if ((status == UMFPACK_WARNING_singular_matrix(1)) && (casualTearingSet==0))
295 {
296 if (!solveSingularSystem(systemData, aux_x))
297 {
298 success = 1;
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 if (1 == systemData->method){
306 /* take the solution */
307 for(i = 0; i < solverData->n_row; ++i)
308 aux_x[i] += solverData->work[i];
309
310 /* update inner equations */
311 wrapper_fvec_umfpack(aux_x, solverData->work, dataAndThreadData, sysNumber);
312 residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
313
314 if ((isnan(residualNorm)(sizeof ((residualNorm)) == sizeof (float) ? __isnanf (residualNorm
) : sizeof ((residualNorm)) == sizeof (double) ? __isnan (residualNorm
) : __isnanl (residualNorm))
) || (residualNorm>1e-4)){
315 warningStreamPrint(LOG_LS, 0,
316 "Failed to solve linear system of equations (no. %d) at time %f. Residual norm is %.15g.",
317 (int)systemData->equationIndex, data->localData[0]->timeValue, residualNorm);
318 success = 0;
319 }
320 } else {
321 /* the solution is automatically in x */
322 }
323
324 if (ACTIVE_STREAM(LOG_LS_V)(useStream[LOG_LS_V]))
325 {
326 if (1 == systemData->method) {
327 infoStreamPrint(LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
328 } else {
329 infoStreamPrint(LOG_LS_V, 1, "Solution x:");
330 }
331 infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar);
332
333 for(i = 0; i < systemData->size; ++i)
334 infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
335
336 messageClose(LOG_LS_V);
337 }
338 }
339 else
340 {
341 warningStreamPrint(LOG_STDOUT, 0,
342 "Failed to solve linear system of equations (no. %d) at time %f, system status %d.",
343 (int)systemData->equationIndex, data->localData[0]->timeValue, status);
344 }
345 solverData->numberSolving += 1;
346
347 return success;
348}
349
350/*! \fn solve a singular linear system with UmfPack methods
351 *
352 * \param [in/out] [systemData]
353 *
354 *
355 * solve even singular system
356 * (note that due to initialization A is given in its transposed form A^T)
357 *
358 * A * x = b
359 * <=> P * R * A * Q * Q * x = P * R * b | P * R * A * Q = L * U
360 * <=> L * U * Q * x = P * R * b
361 *
362 * note that P and Q are orthogonal permutation matrices, so P^(-1) = P^T and Q^(-1) = Q^T
363 *
364 * (1) L * y = P * R * b <=> P^T * L * y = R * b (L is always regular so this can be solved by umfpack)
365 *
366 * (2) U * z = y (U is singular, this cannot be solved by umfpack)
367 *
368 * (3) Q * x = z <=> x = Q^T * z
369 *
370 *
371 * author: kbalzereit, wbraun
372 */
373int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x)
374{
375 DATA_UMFPACK* solverData = (DATA_UMFPACK*) systemData->parDynamicData[omc_get_thread_num()].solverData[0];
376 double *Ux, *Rs, r_ii, *b, sum, *y, *z;
377 int *Up, *Ui, *Q, do_recip, rank = 0, current_rank, current_unz, i, j, k, l,
378 success = 0, status, stop = 0;
379
380 int unz = solverData->info[UMFPACK_UNZ44];
381
382 Up = (int*) malloc((solverData->n_row + 1) * sizeof(int));
383 Ui = (int*) malloc(unz * sizeof(int));
384 Ux = (double*) malloc(unz * sizeof(double));
385
386 Q = (int*) malloc(solverData->n_col * sizeof(int));
387 Rs = (double*) malloc(solverData->n_row * sizeof(double));
388
389 b = (double*) malloc(solverData->n_col * sizeof(double));
390 y = (double*) malloc(solverData->n_col * sizeof(double));
391 z = (double*) malloc(solverData->n_col * sizeof(double));
392
393 infoStreamPrint(LOG_LS_V, 0, "Solve singular system");
394
395 status = umfpack_di_get_numeric((int*) NULL((void*)0), (int*) NULL((void*)0), (double*) NULL((void*)0), Up,
396 Ui, Ux, (int*) NULL((void*)0), Q, (double*) NULL((void*)0), &do_recip, Rs,
397 solverData->numeric);
398
399 switch (status)
400 {
401 case UMFPACK_WARNING_singular_matrix(1):
402 case UMFPACK_ERROR_out_of_memory(-1):
403 case UMFPACK_ERROR_argument_missing(-5):
404 case UMFPACK_ERROR_invalid_system(-13):
405 case UMFPACK_ERROR_invalid_Numeric_object(-3):
406 infoStreamPrint(LOG_LS_V, 0, "error: %d", status);
407 }
408
409 /* calculate R*b */
410 if (do_recip == 0)
411 {
412 for (i = 0; i < solverData->n_row; i++)
413 {
414 b[i] = systemData->parDynamicData[omc_get_thread_num()].b[i] / Rs[i];
415 }
416 }
417 else
418 {
419 for (i = 0; i < solverData->n_row; i++) {
420 b[i] = systemData->parDynamicData[omc_get_thread_num()].b[i] * Rs[i];
421 }
422 }
423
424 /* solve L * y = P * R * b <=> P^T * L * y = R * b */
425 status = umfpack_di_wsolve(UMFPACK_Pt_L(3), solverData->Ap, solverData->Ai,
426 solverData->Ax, y, b, solverData->numeric, solverData->control,
427 solverData->info, solverData->Wi, solverData->W);
428
429 switch (status)
430 {
431 case UMFPACK_WARNING_singular_matrix(1):
432 case UMFPACK_ERROR_out_of_memory(-1):
433 case UMFPACK_ERROR_argument_missing(-5):
434 case UMFPACK_ERROR_invalid_system(-13):
435 case UMFPACK_ERROR_invalid_Numeric_object(-3):
436 infoStreamPrint(LOG_LS_V, 0, "error: %d", status);
437 }
438
439 /* rank is at most as high as the maximum in Ui */
440 for (i = 0; i < unz; i++)
441 {
442 if (rank < Ui[i])
443 rank = Ui[i];
444 }
445
446 /* if rank is already smaller than n set last component of result zero */
447 for (i = rank + 1; i < solverData->n_col; i++)
448 {
449 if (y[i] < 1e-12)
450 {
451 z[i] = 0.0;
452 }
453 else
454 {
455 infoStreamPrint(LOG_LS_V, 0, "error: system is not solvable*");
456 /* free all used memory */
457 free(Up);
458 free(Ui);
459 free(Ux);
460
461 free(Q);
462 free(Rs);
463
464 free(b);
465 free(y);
466 free(z);
467 return -1;
468 }
469 }
470
471 current_rank = rank;
472 current_unz = unz;
473
474 while ((stop == 0) && (current_rank > 1))
475 {
476 /* check if last two rows of U are the same */
477 if ((Ux[current_unz] == Ux[current_unz - 1])
478 && (Ui[current_unz] == Ui[current_unz - 1])
479 && (Up[current_rank] - Up[current_rank - 1] > 1))
480 {
481 /* if diagonal entry on second to last row is nonzero, remaining matrix is regular */
482 if (Ui[Up[current_rank] - 1] == current_rank - 1)
483 {
484 stop = 1;
485 }
486 /* last two rows are the same -> under-determined system, calculate one value and set the other one zero */
487 else
488 {
489 z[current_rank] = y[current_rank] / Ux[current_unz];
490
491 /* reduce system */
492 for (i = Up[current_rank]; i < current_unz; i++)
493 {
494 y[Ui[i]] -= z[current_rank] * Ux[i];
495 }
496
497 current_unz = Up[current_rank] - 1;
498 current_rank--;
499
500 /* now last row has only zero entries */
501 if (y[current_rank] < 1e-12)
502 {
503 z[current_rank] = 0.0;
504 }
505 else
506 {
507 infoStreamPrint(LOG_LS_V, 0, "error: system is not solvable");
508 /* free all used memory */
509 free(Up);
510 free(Ui);
511 free(Ux);
512
513 free(Q);
514 free(Rs);
515
516 free(b);
517 free(y);
518 free(z);
519 return -1;
520 }
521
522 current_rank--;
523 }
524 }
525 else
526 {
527 stop = 1;
528 }
529 }
530
531 /* remaining system is regular so solve system by back substitution */
532 z[current_rank] = Ux[current_unz] * y[current_rank];
533
534 for (i = current_rank - 1; i >= 0; i--)
535 {
536 /* get diagonal element r_ii, j shows where the element is in vector Ux, Ui */
537 j = Up[i];
538 while (Ui[j] != i)
539 {
540 j++;
541 }
542 r_ii = Ux[j];
543 sum = 0.0;
544 for (k = i + 1; k < current_rank; k++)
545 {
546 for (l = Up[k]; l < Up[k + 1]; l++)
547 {
548 if (Ui[l] == Ui[i])
549 {
550 sum += Ux[i] * z[k];
551 }
552 }
553 }
554 z[i] = (y[i] - sum) / r_ii;
555 }
556
557 /* x = Q^T * z */
558 for (i = 0; i < solverData->n_col; i++)
559 {
560 aux_x[Q[i]] = z[i];
561 }
562
563 /* free all used memory */
564 free(Up);
565 free(Ui);
566 free(Ux);
567
568 free(Q);
569 free(Rs);
570
571 free(b);
572 free(y);
573 free(z);
574
575 return success;
576}
577
578void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n)
579{
580 int i, j, k, l;
581
582 char **buffer = (char**)malloc(sizeof(char*)*n);
583 for (l=0; l<n; l++)
584 {
585 buffer[l] = (char*)malloc(sizeof(char)*n*20);
586 buffer[l][0] = 0;
587 }
588
589 k = 0;
590 for (i = 0; i < n; i++)
591 {
592 for (j = 0; j < n; j++)
593 {
594 if ((k < Ap[i + 1]) && (Ai[k] == j))
595 {
596 sprintf(buffer[j], "%s %5g ", buffer[j], Ax[k]);
597 k++;
598 }
599 else
600 {
601 sprintf(buffer[j], "%s %5g ", buffer[j], 0.0);
602 }
603 }
604 }
605 for (l=0; l<n; l++)
606 {
607 infoStreamPrint(LOG_LS_V, 0, "%s", buffer[l]);
608 free(buffer[l]);
609 }
610 free(buffer);
611}
612
613void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n)
614{
615 int i, j, k;
616 char *buffer = (char*)malloc(sizeof(char)*n*20);
617 k = 0;
618 for (i = 0; i < n; i++)
619 {
620 buffer[0] = 0;
621 for (j = 0; j < n; j++)
622 {
623 if ((k < Ap[i + 1]) && (Ai[k] == j))
624 {
625 sprintf(buffer, "%s %5.2g ", buffer, Ax[k]);
626 k++;
627 }
628 else
629 {
630 sprintf(buffer, "%s %5.2g ", buffer, 0.0);
631 }
632 }
633 infoStreamPrint(LOG_LS_V, 0, "%s", buffer);
634 }
635 free(buffer);
636}
637
638#endif