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 / 181
Functions: 0.0% 0 / 0 / 10
Branches: 0.0% 0 / 0 / 130
OMCompiler/SimulationRuntime/c/simulation/solver/gbode_sparse.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 gbode_sparse.c
29 */
30
31 #include "gbode_main.h"
32 #include "gbode_util.h"
33
34 #include "model_help.h"
35 #include "simulation_data.h"
36 #include "solver_main.h"
37
38 #include <limits.h>
39
40 /**
41 * Greedily pack structurally independent columns into colors.
42 *
43 * This works directly on CSC data in O(colors * (nnz + columns)) time and
44 * O(rows) workspace.
45 */
46 ✗ static void colorSparsePattern(SPARSE_PATTERN *pattern, int sizeRows, int sizeCols, int nStages, unsigned int *work)
47 {
48 ✗ assertStreamPrint(NULL, sizeRows >= 0 && sizeCols >= 0 && nStages > 0 && sizeCols % nStages == 0,
49 "Invalid GBODE sparse coloring dimensions %d x %d with %d stages.", sizeRows, sizeCols, nStages);
50 ✗ if (sizeCols == 0)
51 {
52 ✗ pattern->maxColors = 0;
53 ✗ return;
54 }
55 modelica_boolean ownsWork = work == NULL;
56 ✗ if (ownsWork)
57 {
58 ✗ work = malloc(sizeRows * sizeof(unsigned int));
59 }
60
61 unsigned int color = 0;
62 ✗ const int stageSize = sizeCols / nStages;
63 ✗ memset(pattern->colorCols, 0, sizeCols * sizeof(unsigned int));
64 ✗ memset(work, 0, sizeRows * sizeof(unsigned int));
65 ✗ for (int stage = 0; stage < nStages; stage++)
66 {
67 ✗ const int firstCol = stage * stageSize;
68 ✗ const int endCol = firstCol + stageSize;
69 int remaining = stageSize;
70 ✗ while (remaining > 0)
71 {
72 ✗ color++;
73
74 ✗ for (int col = firstCol; col < endCol; col++)
75 {
76 ✗ if (pattern->colorCols[col])
77 {
78 ✗ continue;
79 }
80
81 modelica_boolean conflict = FALSE;
82 ✗ for (unsigned int nz = pattern->leadindex[col]; nz < pattern->leadindex[col + 1]; nz++)
83 {
84 ✗ if (work[pattern->index[nz]] == color)
85 {
86 conflict = TRUE;
87 break;
88 }
89 }
90 ✗ if (conflict)
91 {
92 ✗ continue;
93 }
94
95 ✗ pattern->colorCols[col] = color;
96 ✗ remaining--;
97 ✗ for (unsigned int nz = pattern->leadindex[col]; nz < pattern->leadindex[col + 1]; nz++)
98 {
99 ✗ work[pattern->index[nz]] = color;
100 }
101 }
102 }
103 }
104
105 ✗ pattern->maxColors = color;
106 ✗ if (ownsWork)
107 {
108 ✗ free(work);
109 }
110 }
111
112 ✗ static SPARSE_PATTERN *copySparsePattern(const SPARSE_PATTERN *source, int size, modelica_boolean copyColoring)
113 {
114 ✗ SPARSE_PATTERN *copy = allocSparsePattern(size, source->nnz, size);
115 ✗ memcpy(copy->leadindex, source->leadindex, (size + 1) * sizeof(unsigned int));
116 ✗ memcpy(copy->index, source->index, source->nnz * sizeof(unsigned int));
117 ✗ if (copyColoring)
118 {
119 ✗ memcpy(copy->colorCols, source->colorCols, size * sizeof(unsigned int));
120 ✗ copy->maxColors = source->maxColors;
121 }
122 else
123 {
124 ✗ memset(copy->colorCols, 0, size * sizeof(unsigned int));
125 ✗ copy->maxColors = 0;
126 }
127 ✗ return copy;
128 }
129
130 // Build struct(I + J), optionally coloring it, allocating a pattern when target is NULL
131 ✗ static SPARSE_PATTERN *sparsePatternWithDiagonal(const SPARSE_PATTERN *source, int size, SPARSE_PATTERN *target, unsigned int *work,
132 modelica_boolean colorPattern, modelica_boolean reuseSourceColoring)
133 {
134 int diagonalCount = 0;
135
136 ✗ if (target == NULL)
137 {
138 ✗ for (int col = 0; col < size; col++)
139 {
140 ✗ for (unsigned int nz = source->leadindex[col]; nz < source->leadindex[col + 1]; nz++)
141 {
142 ✗ if (source->index[nz] == col)
143 {
144 ✗ diagonalCount++;
145 ✗ break;
146 }
147 }
148 }
149 ✗ target = allocSparsePattern(size, source->nnz + size - diagonalCount, size);
150 }
151
152 unsigned int targetNz = 0;
153 diagonalCount = 0;
154 ✗ target->leadindex[0] = 0;
155 ✗ for (int col = 0; col < size; col++)
156 {
157 modelica_boolean diagonalPresent = FALSE;
158 ✗ for (unsigned int nz = source->leadindex[col]; nz < source->leadindex[col + 1]; nz++)
159 {
160 ✗ unsigned int row = source->index[nz];
161 ✗ if (!diagonalPresent && row > col)
162 {
163 ✗ target->index[targetNz++] = col;
164 diagonalPresent = TRUE;
165 }
166 ✗ if (row == col)
167 {
168 diagonalPresent = TRUE;
169 ✗ diagonalCount++;
170 }
171 ✗ target->index[targetNz++] = row;
172 }
173 ✗ if (!diagonalPresent)
174 {
175 ✗ target->index[targetNz++] = col;
176 }
177 ✗ target->leadindex[col + 1] = targetNz;
178 }
179 ✗ target->nnz = targetNz;
180
181 ✗ if (!colorPattern)
182 {
183 ✗ memset(target->colorCols, 0, size * sizeof(unsigned int));
184 ✗ target->maxColors = 0;
185 }
186 ✗ else if (reuseSourceColoring && diagonalCount == size && source->maxColors > 0)
187 {
188 ✗ memcpy(target->colorCols, source->colorCols, size * sizeof(unsigned int));
189 ✗ target->maxColors = source->maxColors;
190 }
191 else
192 {
193 ✗ colorSparsePattern(target, size, size, 1, work);
194 }
195
196 ✗ return target;
197 }
198
199 // Reduce a square CSC pattern to the principal submatrix selected by indices
200 ✗ static void reduceSparsePattern(const SPARSE_PATTERN *source, int sourceSize, SPARSE_PATTERN *target,
201 const int *indices, int targetSize, unsigned int *work)
202 {
203 ✗ for (int i = 0; i < sourceSize; i++)
204 {
205 ✗ work[i] = UINT_MAX;
206 }
207 ✗ for (int i = 0; i < targetSize; i++)
208 {
209 ✗ work[indices[i]] = i;
210 }
211
212 unsigned int targetNz = 0;
213 ✗ target->leadindex[0] = 0;
214 ✗ for (int col = 0; col < targetSize; col++)
215 {
216 ✗ int sourceCol = indices[col];
217 ✗ for (unsigned int nz = source->leadindex[sourceCol]; nz < source->leadindex[sourceCol + 1]; nz++)
218 {
219 ✗ unsigned int row = work[source->index[nz]];
220 ✗ if (row != UINT_MAX)
221 {
222 ✗ target->index[targetNz++] = row;
223 }
224 }
225 ✗ target->leadindex[col + 1] = targetNz;
226 }
227 ✗ target->nnz = targetNz;
228 ✗ memset(target->colorCols, 0, targetSize * sizeof(unsigned int));
229 ✗ target->maxColors = 0;
230 ✗ }
231
232 ✗ void gbodeMapSparsePattern(const SPARSE_PATTERN *source, const SPARSE_PATTERN *target,
233 int size, int *sourceToTarget, int *targetDiagonal)
234 {
235 ✗ for (int col = 0; col < size; col++)
236 {
237 ✗ unsigned int sourceNz = source->leadindex[col];
238 ✗ unsigned int sourceEnd = source->leadindex[col + 1];
239 ✗ targetDiagonal[col] = -1;
240
241 ✗ for (unsigned int targetNz = target->leadindex[col]; targetNz < target->leadindex[col + 1]; targetNz++)
242 {
243 ✗ unsigned int targetRow = target->index[targetNz];
244 ✗ if (targetRow == col)
245 {
246 ✗ targetDiagonal[col] = targetNz;
247 }
248 ✗ if (sourceNz < sourceEnd && source->index[sourceNz] == targetRow)
249 {
250 ✗ sourceToTarget[sourceNz++] = targetNz;
251 }
252 }
253
254 ✗ assertStreamPrint(NULL, targetDiagonal[col] >= 0 && sourceNz == sourceEnd, "GBODE sparse pattern mapping failed in column %d.", col);
255 }
256 ✗ }
257
258 // Create the struct(I + J) pattern used by block solves
259 ✗ static SPARSE_PATTERN* initializeSparsePatternBlock(DATA* data, modelica_boolean colorPattern)
260 {
261 ✗ JACOBIAN *jacobian = &data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A];
262 ✗ return sparsePatternWithDiagonal(jacobian->sparsePattern, jacobian->sizeRows, NULL, NULL, colorPattern, TRUE);
263 }
264
265 static SPARSE_PATTERN* initializeSparsePattern_IRK(DATA* data);
266
267 ✗ void initializeSparsePattern_GBODE(DATA* data, DATA_GBODE* gbData)
268 {
269 ✗ assertStreamPrint(NULL, !gbData->isExplicit,
270 "Cannot initialize GBODE sparsity for an explicit method.");
271 ✗ if (gbData->type == GM_TYPE_IMPLICIT && gbData->nlsSolverMethod != GB_NLS_INTERNAL)
272 {
273 ✗ gbData->sparsePattern_NLS = initializeSparsePattern_IRK(data);
274 }
275 else
276 {
277 ✗ gbData->sparsePattern_NLS = initializeSparsePatternBlock(data, gbData->nlsSolverMethod != GB_NLS_INTERNAL);
278 }
279 ✗ }
280
281 ✗ void initializeSparsePattern_GBODEF(DATA* data, DATA_GBODEF* gbfData)
282 {
283 ✗ assertStreamPrint(NULL, !gbfData->isExplicit,
284 "Cannot initialize GBODEF sparsity for an explicit method.");
285 ✗ JACOBIAN *jacobian = &data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A];
286 ✗ const modelica_boolean useInternal = gbfData->nlsSolverMethod == GB_NLS_INTERNAL;
287 ✗ gbfData->sparseWork = malloc(jacobian->sizeRows * sizeof(unsigned int));
288 ✗ gbfData->sparsePattern_ODE = copySparsePattern(jacobian->sparsePattern, jacobian->sizeRows, useInternal);
289 ✗ gbfData->sparsePattern_NLS = sparsePatternWithDiagonal(jacobian->sparsePattern, jacobian->sizeRows,
290 NULL, gbfData->sparseWork, !useInternal, TRUE);
291 ✗ }
292
293 ✗ void updateSparsePattern_GBODEF(DATA* data, DATA_GBODE* gbData)
294 {
295 ✗ DATA_GBODEF *gbfData = gbData->gbfData;
296 ✗ SPARSE_PATTERN *fullPattern = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern;
297
298 ✗ reduceSparsePattern(fullPattern, gbData->nStates, gbfData->sparsePattern_ODE, gbData->fastStatesIdx, gbData->nFastStates, gbfData->sparseWork);
299 ✗ if (gbfData->nlsSolverMethod == GB_NLS_INTERNAL)
300 {
301 ✗ colorSparsePattern(gbfData->sparsePattern_ODE, gbData->nFastStates, gbData->nFastStates, 1, gbfData->sparseWork);
302 }
303
304 ✗ sparsePatternWithDiagonal(gbfData->sparsePattern_ODE, gbData->nFastStates, gbfData->sparsePattern_NLS,
305 ✗ gbfData->sparseWork, gbfData->nlsSolverMethod != GB_NLS_INTERNAL, FALSE);
306
307 ✗ printSparseStructure(gbfData->sparsePattern_NLS, gbData->nFastStates, gbData->nFastStates, OMC_LOG_GBODE_V, "sparsePattern_MR");
308 ✗ }
309
310 /**
311 * @brief Initialize sparsity pattern for non-linear system of full implicit Runge-Kutta methods.
312 *
313 * Get sparsity pattern of ODE Jacobian and map it on the different stages taking into account
314 * the non-zero elements of the A matrix in the Butcher-tableau
315 * Coloring will be calculated, whereby different stages will have different colors, due to the
316 * column-wise calculation of the Jacobian
317 *
318 * @param data Runtime data struct.
319 * @return SPARSE_PATTERN* Pointer to sparsity pattern of non-linear system.
320 */
321 ✗ static SPARSE_PATTERN* initializeSparsePattern_IRK(DATA* data)
322 {
323 unsigned int i,j,k,l;
324 unsigned int row, col;
325 unsigned int missingZeros = 0;
326 unsigned int nDiags = 0, nDiags_A, nnz_A;
327 unsigned int shift = 0;
328 modelica_boolean diagElemNonZero;
329 SPARSE_PATTERN* sparsePattern_IRK;
330 ✗ DATA_GBODE* gbData = (DATA_GBODE*) data->simulationInfo->backupSolverData;
331
332 /* Get Sparsity of ODE Jacobian */
333 ✗ JACOBIAN* jacobian = getSymbolicOdeJacobian(data);
334 ✗ SPARSE_PATTERN* sparsePattern_ODE = getJacobianCscPattern(jacobian);
335
336 ✗ int sizeRows = jacobian->sizeRows;
337 ✗ int sizeCols = jacobian->sizeCols;
338 ✗ int nStages = gbData->tableau->nStages;
339 ✗ int nStates = gbData->nStates;
340 ✗ double* A = gbData->tableau->A;
341
342 ✗ printSparseStructure(sparsePattern_ODE,
343 sizeRows,
344 sizeCols,
345 OMC_LOG_GBODE_V,
346 "sparsePatternODE");
347
348 nnz_A = 0;
349 nDiags_A = 0;
350 ✗ for (i=0; i<nStages; i++) {
351 ✗ if (A[i*nStages + i] != 0) nDiags_A++;
352 ✗ for (j=0; j<nStages; j++) {
353 ✗ if (A[i*nStages + j] != 0) nnz_A++;
354 }
355 }
356
357 i = 0;
358 ✗ for(col=0; col < sizeRows; col++) {
359 ✗ for(; i < sparsePattern_ODE->leadindex[col+1];) {
360 ✗ if(sparsePattern_ODE->index[i++] == col) {
361 ✗ nDiags++;
362 }
363 }
364 }
365 ✗ int missingDiags = jacobian->sizeRows - nDiags;
366 ✗ int nnz = nnz_A*sparsePattern_ODE->nnz + nDiags_A*missingDiags + (nStages-nDiags_A)*nStates;
367
368 // first generated a coordinate format and transform this later to Column pressed format
369 ✗ int *coo_col = (int*) malloc(nnz*sizeof(int));
370 ✗ int *coo_row = (int*) malloc(nnz*sizeof(int));
371
372 i = 0;
373 ✗ for (k=0; k<nStages; k++)
374 {
375 ✗ for (col=0; col < nStates; col++)
376 {
377 diagElemNonZero = FALSE;
378 ✗ for (l=0; l<nStages; l++)
379 {
380 ✗ for (j=sparsePattern_ODE->leadindex[col]; j<sparsePattern_ODE->leadindex[col+1]; j++)
381 {
382 ✗ if (((col + k*nStates) < (sparsePattern_ODE->index[j] + l*nStates)) && !diagElemNonZero)
383 {
384 ✗ coo_col[i] = col + k*nStates;
385 ✗ coo_row[i] = col + k*nStates;
386 ✗ i++;
387 diagElemNonZero = TRUE;
388 }
389 // if the entry in A is non-zero, the sparsity pattern of the ODE-Jacobian will be inserted,
390 // respectively
391 ✗ if (A[l*nStages + k] != 0)
392 {
393 ✗ if ((col + k*nStates) == (sparsePattern_ODE->index[j] + l*nStates))
394 diagElemNonZero = TRUE;
395 ✗ coo_col[i] = col + k*nStates;
396 ✗ coo_row[i] = sparsePattern_ODE->index[j] + l*nStates;
397 ✗ i++;
398 }
399 }
400 }
401 ✗ if (!diagElemNonZero) {
402 ✗ coo_col[i] = col + k*nStates;
403 ✗ coo_row[i] = col + k*nStates;
404 ✗ i++;
405 diagElemNonZero = TRUE;
406 }
407 }
408 }
409
410 ✗ nnz = i;
411
412 ✗ if (OMC_ACTIVE_STREAM(OMC_LOG_GBODE_V)){
413 ✗ printIntVector_gb(OMC_LOG_GBODE_V, "rows", coo_row, nnz, 0.0);
414 ✗ printIntVector_gb(OMC_LOG_GBODE_V, "cols", coo_col, nnz, 0.0);
415 }
416
417 ✗ int length_row_indices = jacobian->sizeCols*nStages+1;
418
419 // Allocate memory for new sparsity pattern
420 ✗ sparsePattern_IRK = allocSparsePattern(jacobian->sizeCols*nStages, nnz, jacobian->sizeCols*nStages);
421
422 /* Set diagonal elements of sparsitiy pattern to non-zero */
423 ✗ for (i=0; i<length_row_indices; i++)
424 ✗ sparsePattern_IRK->leadindex[i] = 0;
425
426 ✗ for (int i = 0; i < nnz; i++)
427 {
428 ✗ sparsePattern_IRK->index[i] = coo_row[i];
429 ✗ sparsePattern_IRK->leadindex[coo_col[i] + 1]++;
430 }
431 ✗ for (int i = 0; i < sizeCols*nStages; i++)
432 {
433 ✗ sparsePattern_IRK->leadindex[i + 1] += sparsePattern_IRK->leadindex[i];
434 }
435
436 ✗ free(coo_col);
437 ✗ free(coo_row);
438
439 ✗ colorSparsePattern(sparsePattern_IRK, sizeRows*nStages, sizeCols*nStages, nStages, NULL);
440
441 ✗ return sparsePattern_IRK;
442 }
443