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