Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 240
Functions: 0.0% 0 / 0 / 7
Branches: 0.0% 0 / 0 / 144
OMCompiler/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c
Line Branch Exec Source
1 /*
2 * This file belongs to the OpenModelica Run-Time System
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings
5 * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights
6 * reserved.
7 *
8 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
9 * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY
10 * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
11 * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL
12 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
13 *
14 * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License
15 * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs:
16 * http://www.openmodelica.org or https://github.com/OpenModelica/ or
17 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU
18 * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW
19 * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY
23 * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF
24 * OSMC-PL.
25 *
26 */
27
28 /*! \file linearSolverUmfpack.c
29 */
30
31 #include "omc_config.h"
32
33 #ifdef WITH_SUITESPARSE
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 "omc_math.h"
42 #include "util/varinfo.h"
43 #include "model_help.h"
44
45 #include "linearSystem.h"
46 #include "linearSolverUmfpack.h"
47
48 void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n);
49 void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n);
50 int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x);
51
52 /*! \fn allocate memory for linear system solver UmfPack
53 *
54 */
55 int
56 ✗ allocateUmfPackData(int n_row, int n_col, int nz, void** voiddata)
57 {
58 ✗ DATA_UMFPACK* data = (DATA_UMFPACK*) malloc(sizeof(DATA_UMFPACK));
59 ✗ assertStreamPrint(NULL, 0 != data, "Could not allocate data for linear solver UmfPack.");
60
61 ✗ data->symbolic = NULL;
62 ✗ data->numeric = NULL;
63
64 ✗ data->n_col = n_col;
65 ✗ data->n_row = n_row;
66 ✗ data->nnz = nz;
67
68
69 ✗ data->Ap = (int*) calloc((n_row+1),sizeof(int));
70
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->Wi = (int*) malloc(n_row * sizeof(int));
76 ✗ data->W = (double*) malloc(5*n_row * sizeof(double));
77
78 ✗ data->numberSolving=0;
79 ✗ umfpack_di_defaults(data->control);
80
81 ✗ data->control[UMFPACK_PIVOT_TOLERANCE] = 0.1;
82 ✗ data->control[UMFPACK_IRSTEP] = 2;
83 ✗ data->control[UMFPACK_SCALE] = 1;
84 ✗ data->control[UMFPACK_STRATEGY] = 5;
85
86
87
88 ✗ *voiddata = (void*)data;
89
90 ✗ return 0;
91 }
92
93
94 /*! \fn free memory for linear system solver UmfPack
95 *
96 */
97 int
98 ✗ freeUmfPackData(void **voiddata)
99 {
100 ✗ DATA_UMFPACK* data = (DATA_UMFPACK*) *voiddata;
101
102 ✗ free(data->Ap);
103 ✗ free(data->Ai);
104 ✗ free(data->Ax);
105 ✗ free(data->work);
106
107 ✗ free(data->Wi);
108 ✗ free(data->W);
109
110 ✗ if(data->symbolic)
111 ✗ umfpack_di_free_symbolic (&data->symbolic);
112 ✗ if(data->numeric)
113 ✗ umfpack_di_free_numeric (&data->numeric);
114
115 ✗ return 0;
116 }
117
118 /*! \fn getAnalyticalJacobian
119 *
120 * function calculates analytical jacobian
121 *
122 * \param [ref] [data]
123 * \param [in] [sysNumber]
124 *
125 * \author wbraun
126 *
127 */
128 ✗ void getAnalyticalJacobianUmfPack(DATA* data, threadData_t *threadData, LINEAR_SYSTEM_DATA* systemData)
129 {
130 int i,j,l,nth;
131 ✗ JACOBIAN* jacobian = systemData->jacobian;
132 ✗ JACOBIAN* parentJacobian = systemData->parentJacobian;
133 ✗ const SPARSE_PATTERN* sp = jacobian->sparsePattern;
134
135 /* evaluate constant equations of Jacobian */
136 ✗ if (jacobian->constantEqns != NULL) {
137 ✗ jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
138 }
139
140 /* evaluate Jacobian */
141 ✗ for (i = 0; i < sp->maxColors; i++) {
142 /* activate seed variable for the corresponding color */
143 ✗ for (j = 0; j < jacobian->sizeCols; j++)
144 ✗ if (sp->colorCols[j]-1 == i)
145 ✗ jacobian->seedVars[j] = 1.0;
146
147 /* evaluate Jacobian column */
148 ✗ jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
149
150 ✗ for (j = 0; j < jacobian->sizeCols; j++) {
151 ✗ if (sp->colorCols[j]-1 == i) {
152 ✗ for (nth = sp->leadindex[j]; nth < sp->leadindex[j+1]; nth++) {
153 ✗ l = sp->index[nth];
154 ✗ systemData->setAElement(j, l, -jacobian->resultVars[l], nth, systemData, threadData);
155 }
156 /* de-activate seed variable for the corresponding color */
157 ✗ jacobian->seedVars[j] = 0.0;
158 }
159 }
160 }
161 ✗ }
162
163 /*! \fn wrapper_fvec_umfpack for the residual function
164 *
165 */
166 static int wrapper_fvec_umfpack(double* x, double* f, RESIDUAL_USERDATA* resUserData, int sysNumber)
167 {
168 ✗ int iflag = 0;
169
170 ✗ resUserData->data->simulationInfo->linearSystemData[sysNumber].residualFunc(resUserData, x, f, &iflag);
171 ✗ return 0;
172 }
173
174 /*! \fn solve linear system with UmfPack method
175 *
176 * \param [in] [data]
177 * [sysNumber] index of the corresponding linear system
178 *
179 *
180 * author: kbalzereit, wbraun
181 */
182 int
183 ✗ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x)
184 {
185 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
186 ✗ LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]);
187 ✗ DATA_UMFPACK* solverData = (DATA_UMFPACK*)systemData->solverData[0];
188 _omc_scalar residualNorm = 0;
189
190 ✗ int i, j, status = UMFPACK_OK, success = 0, ni=0, n = systemData->size, eqSystemNumber = systemData->equationIndex, indexes[2] = {1,eqSystemNumber};
191 ✗ int casualTearingSet = systemData->strictTearingFunctionCall != NULL;
192 double tmpJacEvalTime;
193 ✗ int reuseMatrixJac = (data->simulationInfo->currentContext == CONTEXT_SYM_JACOBIAN && data->simulationInfo->currentJacobianEval > 0);
194
195 ✗ infoStreamPrintWithEquationIndexes(OMC_LOG_LS, omc_dummyFileInfo, 0, indexes,
196 "Start solving Linear System %d (size %d) at time %g with UMFPACK Solver",
197 ✗ eqSystemNumber, (int) systemData->size, data->localData[0]->timeValue);
198
199 ✗ rt_ext_tp_tick(&(solverData->timeClock));
200 ✗ if (0 == systemData->method)
201 {
202 ✗ if (!reuseMatrixJac){
203 /* set A matrix */
204 ✗ solverData->Ap[0] = 0;
205 ✗ systemData->setA(data, threadData, systemData);
206 ✗ solverData->Ap[solverData->n_row] = solverData->nnz;
207 }
208
209 /* set b vector */
210 ✗ systemData->setb(data, threadData, systemData);
211 } else {
212
213 ✗ if (!reuseMatrixJac){
214 ✗ solverData->Ap[0] = 0;
215 /* calculate jacobian -> matrix A*/
216 ✗ if(systemData->jacobianIndex != -1){
217 ✗ getAnalyticalJacobianUmfPack(data, threadData, systemData);
218 } else {
219 assertStreamPrint(threadData, 1, "jacobian function pointer is invalid" );
220 }
221 ✗ solverData->Ap[solverData->n_row] = solverData->nnz;
222 }
223
224 /* calculate vector b (rhs) */
225 ✗ memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row);
226 ✗ wrapper_fvec_umfpack(solverData->work, systemData->b, &resUserData, sysNumber);
227 }
228 ✗ tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock));
229 ✗ systemData->jacobianTime += tmpJacEvalTime;
230 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "### %f time to set Matrix A and vector b.", tmpJacEvalTime);
231
232 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_LS_V))
233 {
234 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Old solution x:");
235 ✗ for(i = 0; i < solverData->n_row; ++i)
236 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
237 ✗ messageClose(OMC_LOG_LS_V);
238
239 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Matrix A n_rows = %d", solverData->n_row);
240 ✗ for (i=0; i<solverData->n_row; i++){
241 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "%d. Ap => %d -> %d", i, solverData->Ap[i], solverData->Ap[i+1]);
242 ✗ for (j=solverData->Ap[i]; j<solverData->Ap[i+1]; j++){
243 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "A[%d,%d] = %f", i, solverData->Ai[j], solverData->Ax[j]);
244 }
245 }
246 ✗ messageClose(OMC_LOG_LS_V);
247
248 ✗ for (i=0; i<solverData->n_row; i++) {
249 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "b[%d] = %e", i, systemData->b[i]);
250 }
251 }
252 ✗ rt_ext_tp_tick(&(solverData->timeClock));
253
254 /* symbolic pre-ordering of A to reduce fill-in of L and U */
255 ✗ if (0 == solverData->numberSolving) {
256 ✗ status = umfpack_di_symbolic(solverData->n_col, solverData->n_row, solverData->Ap, solverData->Ai, solverData->Ax, &(solverData->symbolic), solverData->control, solverData->info);
257 }
258
259 /* compute the LU factorization of A */
260 /* if reuseMatrixJac use also previous factorization */
261 ✗ if (!reuseMatrixJac)
262 {
263 ✗ if (0 == status){
264 ✗ umfpack_di_free_numeric(&(solverData->numeric));
265 ✗ status = umfpack_di_numeric(solverData->Ap, solverData->Ai, solverData->Ax, solverData->symbolic, &(solverData->numeric), solverData->control, solverData->info);
266 }
267 }
268
269 ✗ if (0 == status){
270 ✗ if (1 == systemData->method){
271 ✗ status = umfpack_di_wsolve(UMFPACK_A, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W);
272 } else {
273 ✗ status = umfpack_di_wsolve(UMFPACK_Aat, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W);
274 }
275 }
276
277 ✗ if (status == UMFPACK_OK){
278 success = 1;
279 }
280 ✗ else if ((status == UMFPACK_WARNING_singular_matrix) && (casualTearingSet==0))
281 {
282 ✗ if (!solveSingularSystem(systemData, aux_x))
283 {
284 success = 1;
285 }
286 }
287 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Solve System: %f", rt_ext_tp_tock(&(solverData->timeClock)));
288
289 /* print solution */
290 ✗ if (1 == success){
291 ✗ if (1 == systemData->method){
292 /* take the solution */
293 ✗ for(i = 0; i < solverData->n_row; ++i)
294 ✗ aux_x[i] += solverData->work[i];
295
296 /* update inner equations */
297 ✗ wrapper_fvec_umfpack(aux_x, solverData->work, &resUserData, sysNumber);
298 ✗ residualNorm = _omc_gen_euclideanVectorNorm(solverData->work, solverData->n_row);
299
300 ✗ if ((isnan(residualNorm)) || (residualNorm>1e-4)){
301 ✗ warningStreamPrintWithLimit(OMC_LOG_LS, 0, ++(systemData->numberOfFailures) /* Update counter */, data->simulationInfo->maxWarnDisplays,
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 /* the solution is automatically in x */
308 }
309
310 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_LS_V))
311 {
312 ✗ if (1 == systemData->method) {
313 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Residual Norm %.15g of solution x:", residualNorm);
314 } else {
315 ✗ infoStreamPrint(OMC_LOG_LS_V, 1, "Solution x:");
316 }
317 ✗ infoStreamPrint(OMC_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(OMC_LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]);
321
322 ✗ messageClose(OMC_LOG_LS_V);
323 }
324 }
325 else
326 {
327 ✗ warningStreamPrintWithLimit(OMC_LOG_LS, 0, ++(systemData->numberOfFailures) /* Update counter */, data->simulationInfo->maxWarnDisplays,
328 "Failed to solve linear system of equations (no. %d) at time %f, system status %d.",
329 ✗ (int)systemData->equationIndex, data->localData[0]->timeValue, status);
330 }
331 ✗ solverData->numberSolving += 1;
332
333 ✗ return success;
334 }
335
336 /*! \fn solve a singular linear system with UmfPack methods
337 *
338 * \param [in/out] [systemData]
339 *
340 *
341 * solve even singular system
342 * (note that due to initialization A is given in its transposed form A^T)
343 *
344 * A * x = b
345 * <=> P * R * A * Q * Q * x = P * R * b | P * R * A * Q = L * U
346 * <=> L * U * Q * x = P * R * b
347 *
348 * note that P and Q are orthogonal permutation matrices, so P^(-1) = P^T and Q^(-1) = Q^T
349 *
350 * (1) L * y = P * R * b <=> P^T * L * y = R * b (L is always regular so this can be solved by umfpack)
351 *
352 * (2) U * z = y (U is singular, this cannot be solved by umfpack)
353 *
354 * (3) Q * x = z <=> x = Q^T * z
355 *
356 *
357 * author: kbalzereit, wbraun
358 */
359 ✗ int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x)
360 {
361 ✗ DATA_UMFPACK* solverData = (DATA_UMFPACK*) systemData->solverData[0];
362 double *Ux, *Rs, r_ii, *b, sum, *y, *z;
363 int *Up, *Ui, *Q, do_recip, rank = 0, current_rank, current_unz, i, j, k, l,
364 success = 0, status, stop = 0;
365
366 ✗ int unz = solverData->info[UMFPACK_UNZ];
367
368 /* umfpack_di_get_numeric writes Up[n_col+1], Ui[unz] and Ux[unz] */
369 ✗ Up = (int*) malloc((solverData->n_col + 1) * sizeof(int));
370 ✗ Ui = (int*) malloc(unz * sizeof(int));
371 ✗ Ux = (double*) malloc(unz * sizeof(double));
372
373 ✗ Q = (int*) malloc(solverData->n_col * sizeof(int));
374 ✗ Rs = (double*) malloc(solverData->n_row * sizeof(double));
375
376 ✗ b = (double*) malloc(solverData->n_col * sizeof(double));
377 ✗ y = (double*) malloc(solverData->n_col * sizeof(double));
378 ✗ z = (double*) malloc(solverData->n_col * sizeof(double));
379
380 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "Solve singular system");
381
382 ✗ status = umfpack_di_get_numeric((int*) NULL, (int*) NULL, (double*) NULL, Up,
383 Ui, Ux, (int*) NULL, Q, (double*) NULL, &do_recip, Rs,
384 solverData->numeric);
385
386 ✗ switch (status)
387 {
388 ✗ case UMFPACK_WARNING_singular_matrix:
389 case UMFPACK_ERROR_out_of_memory:
390 case UMFPACK_ERROR_argument_missing:
391 case UMFPACK_ERROR_invalid_system:
392 case UMFPACK_ERROR_invalid_Numeric_object:
393 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "error: %d", status);
394 }
395
396 /* calculate R*b */
397 ✗ if (do_recip == 0)
398 {
399 ✗ for (i = 0; i < solverData->n_row; i++)
400 {
401 ✗ b[i] = systemData->b[i] / Rs[i];
402 }
403 }
404 else
405 {
406 ✗ for (i = 0; i < solverData->n_row; i++) {
407 ✗ b[i] = systemData->b[i] * Rs[i];
408 }
409 }
410
411 /* solve L * y = P * R * b <=> P^T * L * y = R * b */
412 ✗ status = umfpack_di_wsolve(UMFPACK_Pt_L, solverData->Ap, solverData->Ai,
413 ✗ solverData->Ax, y, b, solverData->numeric, solverData->control,
414 ✗ solverData->info, solverData->Wi, solverData->W);
415
416 ✗ switch (status)
417 {
418 ✗ case UMFPACK_WARNING_singular_matrix:
419 case UMFPACK_ERROR_out_of_memory:
420 case UMFPACK_ERROR_argument_missing:
421 case UMFPACK_ERROR_invalid_system:
422 case UMFPACK_ERROR_invalid_Numeric_object:
423 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "error: %d", status);
424 }
425
426 /* rank is at most as high as the maximum in Ui */
427 ✗ for (i = 0; i < unz; i++)
428 {
429 ✗ if (rank < Ui[i])
430 rank = Ui[i];
431 }
432
433 /* if rank is already smaller than n set last component of result zero */
434 ✗ for (i = rank + 1; i < solverData->n_col; i++)
435 {
436 ✗ if (y[i] < 1e-12)
437 {
438 ✗ z[i] = 0.0;
439 }
440 else
441 {
442 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "error: system is not solvable*");
443 success = -1;
444 ✗ goto cleanup;
445 }
446 }
447
448 current_rank = rank;
449 /* U is column-stored, so column j owns Ui/Ux[Up[j] .. Up[j+1]-1] and its last
450 * entry is the diagonal; current_unz is that entry's index, as every use below
451 * assumes. unz is one past the end of the whole array. */
452 ✗ if (Up[current_rank + 1] <= Up[current_rank])
453 {
454 /* no pivot in this column - nothing to back-substitute with */
455 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "error: system is not solvable*");
456 success = -1;
457 ✗ goto cleanup;
458 }
459 ✗ current_unz = Up[current_rank + 1] - 1;
460
461 ✗ while ((stop == 0) && (current_rank > 1))
462 {
463 /* check if last two rows of U are the same */
464 ✗ if ((Ux[current_unz] == Ux[current_unz - 1])
465 ✗ && (Ui[current_unz] == Ui[current_unz - 1])
466 ✗ && (Up[current_rank] - Up[current_rank - 1] > 1))
467 {
468 /* if diagonal entry on second to last row is nonzero, remaining matrix is regular */
469 ✗ if (Ui[Up[current_rank] - 1] == current_rank - 1)
470 {
471 stop = 1;
472 }
473 /* last two rows are the same -> under-determined system, calculate one value and set the other one zero */
474 else
475 {
476 ✗ z[current_rank] = y[current_rank] / Ux[current_unz];
477
478 /* reduce system */
479 ✗ for (i = Up[current_rank]; i < current_unz; i++)
480 {
481 ✗ y[Ui[i]] -= z[current_rank] * Ux[i];
482 }
483
484 ✗ current_unz = Up[current_rank] - 1;
485 current_rank--;
486
487 /* now last row has only zero entries */
488 ✗ if (y[current_rank] < 1e-12)
489 {
490 ✗ z[current_rank] = 0.0;
491 }
492 else
493 {
494 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "error: system is not solvable");
495 success = -1;
496 ✗ goto cleanup;
497 }
498
499 ✗ current_rank--;
500 }
501 }
502 else
503 {
504 stop = 1;
505 }
506 }
507
508 /* remaining system is regular so solve system by back substitution */
509 ✗ z[current_rank] = Ux[current_unz] * y[current_rank];
510
511 ✗ for (i = current_rank - 1; i >= 0; i--)
512 {
513 /* get diagonal element r_ii, j shows where the element is in vector Ux, Ui */
514 ✗ j = Up[i];
515 ✗ while ((j < Up[i + 1]) && (Ui[j] != i))
516 {
517 ✗ j++;
518 }
519 ✗ if (j >= Up[i + 1])
520 {
521 /* a singular U can miss a diagonal; searching on would run off Ui */
522 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "error: system is not solvable*");
523 success = -1;
524 ✗ goto cleanup;
525 }
526 ✗ r_ii = Ux[j];
527 sum = 0.0;
528 ✗ for (k = i + 1; k < current_rank; k++)
529 {
530 ✗ for (l = Up[k]; l < Up[k + 1]; l++)
531 {
532 ✗ if (Ui[l] == Ui[i])
533 {
534 ✗ sum += Ux[i] * z[k];
535 }
536 }
537 }
538 ✗ z[i] = (y[i] - sum) / r_ii;
539 }
540
541 /* x = Q^T * z */
542 ✗ for (i = 0; i < solverData->n_col; i++)
543 {
544 ✗ aux_x[Q[i]] = z[i];
545 }
546
547 ✗ cleanup:
548 /* free all used memory */
549 ✗ free(Up);
550 ✗ free(Ui);
551 ✗ free(Ux);
552
553 ✗ free(Q);
554 ✗ free(Rs);
555
556 ✗ free(b);
557 ✗ free(y);
558 ✗ free(z);
559
560 ✗ return success;
561 }
562
563 ✗ void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n)
564 {
565 int i, j, k, l;
566
567 ✗ char **buffer = (char**)malloc(sizeof(char*)*n);
568 ✗ for (l=0; l<n; l++)
569 {
570 ✗ buffer[l] = (char*)malloc(sizeof(char)*n*20);
571 ✗ buffer[l][0] = 0;
572 }
573
574 ✗ char **p = (char**)malloc(sizeof(char*)*n);
575 ✗ for (l=0; l<n; l++)
576 ✗ p[l] = buffer[l];
577
578 k = 0;
579 ✗ for (i = 0; i < n; i++)
580 {
581 ✗ for (j = 0; j < n; j++)
582 {
583 ✗ if ((k < Ap[i + 1]) && (Ai[k] == j))
584 {
585 ✗ p[j] += sprintf(p[j], " %5g ", Ax[k]);
586 ✗ k++;
587 }
588 else
589 {
590 ✗ p[j] += sprintf(p[j], " %5g ", 0.0);
591 }
592 }
593 }
594 ✗ for (l=0; l<n; l++)
595 {
596 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "%s", buffer[l]);
597 ✗ free(buffer[l]);
598 }
599 ✗ free(p);
600 ✗ free(buffer);
601 ✗ }
602
603 ✗ void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n)
604 {
605 int i, j, k;
606 ✗ char *buffer = (char*)malloc(sizeof(char)*n*20);
607 char *q;
608 k = 0;
609 ✗ for (i = 0; i < n; i++)
610 {
611 q = buffer;
612 ✗ for (j = 0; j < n; j++)
613 {
614 ✗ if ((k < Ap[i + 1]) && (Ai[k] == j))
615 {
616 ✗ q += sprintf(q, " %5.2g ", Ax[k]);
617 ✗ k++;
618 }
619 else
620 {
621 ✗ q += sprintf(q, " %5.2g ", 0.0);
622 }
623 }
624 ✗ infoStreamPrint(OMC_LOG_LS_V, 0, "%s", buffer);
625 }
626 ✗ free(buffer);
627 ✗ }
628
629 #endif
630