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 / 107
Functions: 0.0% 0 / 0 / 5
Branches: 0.0% 0 / 0 / 88
OMCompiler/SimulationRuntime/c/simulation/solver/sundials_util.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 sundials_util.c
29 */
30
31 #include "util/omc_error.h"
32
33 #include "sundials_util.h"
34 #include <sunmatrix/sunmatrix_sparse.h>
35 #include <nvector/nvector_serial.h>
36
37
38 #ifdef WITH_SUNDIALS
39
40 #define UNUSED(x) (void)(x) /* Surpress compiler warnings for unused function input */
41
42 /**
43 * @brief Set element of sparse Sundials matrix.
44 *
45 * Jac(row, column) = val. The column pointers are not touched, set them with
46 * setSundialsSparseColPtrs.
47 *
48 * @param row Row of matrix element.
49 * @param column Column of matrix element, unused.
50 * @param nth Sparsity pattern lead index.
51 * @param value Value to set in position (i,j).
52 * @param Jac Pointer to double array storing matrix.
53 * @param nRows Number of rows of Jacobian matrix, unused.
54 */
55 ✗ void setJacElementSundialsSparse(int row, int column, int nth, double value, void* Jac, int nRows)
56 {
57 UNUSED(column);
58 UNUSED(nRows); /* Disables compiler warning */
59
60 SUNMatrix A = (SUNMatrix) Jac;
61 /* TODO: Remove this check for performance reasons? */
62 ✗ if (SM_SPARSETYPE_S(A) != SUN_CSC_MAT) {
63 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0,
64 "In function setJacElementSundialsSparse: Wrong sparse format "
65 "of SUNMatrix A.");
66 }
67
68 ✗ SM_INDEXVALS_S(A)[nth] = row;
69 ✗ SM_DATA_S(A)[nth] = value;
70 ✗ }
71
72 /**
73 * @brief Set the column pointers of a CSC Sundials matrix from a sparsity pattern.
74 *
75 * @param sp Column oriented (CSC) sparsity pattern of the matrix.
76 * @param Jac Sundials Matrix with the same number of columns as sp.
77 */
78 ✗ void setSundialsSparseColPtrs(const SPARSE_PATTERN* sp, SUNMatrix Jac)
79 {
80 sunindextype column;
81
82 ✗ for (column = 0; column <= SM_COLUMNS_S(Jac); column++) {
83 ✗ SM_INDEXPTRS_S(Jac)[column] = sp->leadindex[column];
84 }
85 ✗ }
86
87 /**
88 * @brief Set Sundials sparse pattern from SimRuntime SPARSE_PATTERN
89 *
90 * @param jacobian Jacobian whose column oriented (CSC) pattern is used, see
91 * getJacobianCscPattern().
92 * @param Jac Sundials Matrix
93 */
94 ✗ void setSundialsSparsePattern(JACOBIAN* jacobian, SUNMatrix Jac) {
95 ✗ const SPARSE_PATTERN *sp = getJacobianCscPattern(jacobian);
96 long int column, row, nz;
97
98 ✗ for (column = 0; column < jacobian->sizeCols; column++) {
99 ✗ for (nz = sp->leadindex[column]; nz < sp->leadindex[column + 1]; nz++) {
100 ✗ row = sp->index[nz];
101 ✗ SM_INDEXVALS_S(Jac)[nz] = row;
102 }
103 }
104 ✗ setSundialsSparseColPtrs(sp, Jac);
105 ✗ }
106
107 /**
108 * @brief Scaling of a sparse matrix column-wise by a vector
109 *
110 * Might get fixed in a future version of SUDNIALS: https://github.com/llnl/sundials/issues/590
111 *
112 * @param A Sparse matrix in CSC. Will be scaled on return.
113 * @param vScale Vector for scaling.
114 * @return Return `SUN_SUCCESS` on success.
115 */
116 ✗ SUNErrCode _omc_SUNSparseMatrixVecScaling(SUNMatrix A, N_Vector vScale)
117 {
118
119 /* should not be called unless A is a sparse matrix in CSC format;
120 otherwise return immediately */
121 ✗ if (SUNMatGetID(A) != SUNMATRIX_SPARSE || SM_SPARSETYPE_S(A) == SUN_CSR_MAT) {
122 return SUN_ERR_ARG_INCOMPATIBLE;
123 }
124
125 sunindextype i, j;
126 char *matrixtype;
127 char *indexname;
128 ✗ sunrealtype *vScaling = N_VGetArrayPointer(vScale);
129
130 ✗ for (j=0; j<SM_NP_S(A); j++) {
131 ✗ for (i=(SM_INDEXPTRS_S(A))[j]; i<(SM_INDEXPTRS_S(A))[j+1]; i++) {
132 ✗ (SM_DATA_S(A))[i] = (SM_DATA_S(A))[i]/vScaling[j];
133 }
134 }
135
136 return SUN_SUCCESS;
137 }
138
139 /**
140 * @brief Calculates A+c*I and stores the result in A.
141 *
142 * SUNDIALS has no A+c*I operation. Its SUNMatScaleAddI computes c*A+I, which
143 * scales the matrix and adds an unscaled identity - the other way around. The
144 * exact substitute would be SUNMatScaleAdd(1, A, B) with B holding c*I, but
145 * SUNMatScaleAdd_Sparse still zeroes an M-sized work array once per column and
146 * so runs in O(M*N); see https://github.com/LLNL/sundials/issues/253 and
147 * https://github.com/llnl/sundials/issues/590.
148 *
149 * This implementation is O(NNZ): the diagonal scan runs over the stored entries
150 * of each column, and the reallocating path fills its work arrays per column
151 * from that column's own entries rather than clearing them over all M rows.
152 *
153 * TODO: Drop this if SUNDIALS ever grows an A+c*I operation, or if
154 * SUNMatScaleAdd_Sparse becomes linear in the number of nonzeros.
155 *
156 * @param c Constant to scale identity matrix I.
157 * @param A Sparse matrix in CSC or CSR format.
158 * @return Returns SUN_SUCCESS on success
159 * or SUN_ERR_MEM_FAIL if failed to allocate memory.
160 */
161 ✗ SUNErrCode _omc_SUNMatScaleIAdd_Sparse(sunrealtype c, SUNMatrix A)
162 {
163 sunindextype j, i, p, nz, newvals, M, N, cend, nw;
164 sunbooleantype newmat, found;
165 sunindextype *w, *Ap, *Ai, *Cp, *Ci;
166 sunrealtype *x, *Ax, *Cx;
167 SUNMatrix C;
168
169 /* store shortcuts to matrix dimensions (M is inner dimension, N is outer) */
170 ✗ if (SM_SPARSETYPE_S(A) == SUN_CSC_MAT) {
171 ✗ M = SM_ROWS_S(A);
172 ✗ N = SM_COLUMNS_S(A);
173 }
174 else {
175 ✗ M = SM_COLUMNS_S(A);
176 ✗ N = SM_ROWS_S(A);
177 }
178
179 /* access data arrays from A (return if failure) */
180 Ap = Ai = NULL;
181 Ax = NULL;
182 ✗ if (SM_INDEXPTRS_S(A)) Ap = SM_INDEXPTRS_S(A);
183 else return (SUN_ERR_MEM_FAIL);
184 ✗ if (SM_INDEXVALS_S(A)) Ai = SM_INDEXVALS_S(A);
185 else return (SUN_ERR_MEM_FAIL);
186 ✗ if (SM_DATA_S(A)) Ax = SM_DATA_S(A);
187 else return (SUN_ERR_MEM_FAIL);
188
189
190 /* determine if A: contains values on the diagonal (so c*I can just be added in);
191 if not, then increment counter for extra storage that should be required. */
192 newvals = 0;
193 ✗ for (j = 0; j < (M < N ? M : N); j++) {
194 /* scan column (row if CSR) of A, searching for diagonal value */
195 found = SUNFALSE;
196 ✗ for (i = Ap[j]; i < Ap[j+1]; i++) {
197 ✗ if (Ai[i] == j) {
198 found = SUNTRUE;
199 break;
200 }
201 }
202 /* if no diagonal found, increment necessary storage counter */
203 ✗ if (!found) newvals += 1;
204 }
205
206 /* If extra nonzeros required, check whether matrix has sufficient storage space
207 for new nonzero entries (so c*I can be inserted into existing storage) */
208 newmat = SUNFALSE; /* no reallocation needed */
209 ✗ if (newvals > (SM_NNZ_S(A) - Ap[N]))
210 newmat = SUNTRUE;
211
212
213 /* perform operation based on existing/necessary structure */
214
215 /* case 1: A already contains the diagonal */
216 ✗ if (newvals == 0) {
217
218 /* iterate through diagonal, adding c */
219 ✗ for (j = 0; j < (M < N ? M : N); j++)
220 ✗ for (i = Ap[j]; i < Ap[j+1]; i++)
221 ✗ if (Ai[i] == j) {
222 ✗ Ax[i] += c;
223 ✗ break;
224 }
225
226
227 /* case 2: A has sufficient storage, but does not already contain a diagonal */
228 ✗ } else if (!newmat) {
229
230 /* create work arrays for nonzero row (column) indices and values in a single column (row) */
231 ✗ w = (sunindextype *) malloc(M * sizeof(sunindextype));
232 ✗ x = (sunrealtype *) malloc(M * sizeof(sunrealtype));
233
234 /* determine storage location where last column (row) should end */
235 ✗ nz = Ap[N] + newvals;
236
237 /* store pointer past last column (row) from original A,
238 and store updated value in revised A */
239 cend = Ap[N];
240 ✗ Ap[N] = nz;
241
242 /* iterate through columns (rows) backwards */
243 ✗ for (j = N-1; j >= 0; j--) {
244
245 /* reset diagonal entry, in case it's not in A */
246 ✗ x[j] = 0.0;
247
248 /* iterate down column (row) of A, collecting nonzeros */
249 ✗ for (p = Ap[j], i = 0; p < cend; p++, i++) {
250 ✗ w[i] = Ai[p]; /* collect index */
251 ✗ x[Ai[p]] = Ax[p]; /* collect value */
252 }
253
254 /* store nnz of this column (row) */
255 ✗ nw = cend - Ap[j];
256
257 /* add identity to this column (row) */
258 ✗ if (j < M) {
259 ✗ x[j] += c; /* update value */
260 }
261
262 /* fill entries of A with this column's (row's) data */
263 /* fill entries past diagonal */
264 ✗ for (i = nw-1; i >= 0 && w[i] > j; i--) {
265 ✗ Ai[--nz] = w[i];
266 ✗ Ax[nz] = x[w[i]];
267 }
268 /* fill diagonal if applicable (i < 0 when all entries were above diagonal) */
269 ✗ if ((i < 0) || (w[i] != j)) {
270 ✗ Ai[--nz] = j;
271 ✗ Ax[nz] = x[j];
272 }
273 /* fill entries before diagonal */
274 ✗ for (; i >= 0; i--) {
275 ✗ Ai[--nz] = w[i];
276 ✗ Ax[nz] = x[w[i]];
277 }
278
279 /* store ptr past this col (row) from orig A, update value for new A */
280 ✗ cend = Ap[j];
281 ✗ Ap[j] = nz;
282
283 }
284
285 /* clean up */
286 ✗ free(w);
287 ✗ free(x);
288
289
290 /* case 3: A must be reallocated with sufficient storage */
291 } else {
292
293 /* create work array for nonzero values in a single column (row) */
294 ✗ x = (sunrealtype *) malloc(M * sizeof(sunrealtype));
295
296 /* create new matrix for sum, reusing A's SUNDIALS context */
297 ✗ C = SUNSparseMatrix(SM_ROWS_S(A), SM_COLUMNS_S(A),
298 Ap[N] + newvals,
299 SM_SPARSETYPE_S(A),
300 A->sunctx);
301
302 /* access data from CSR structures (return if failure) */
303 Cp = Ci = NULL;
304 Cx = NULL;
305 ✗ if (SM_INDEXPTRS_S(C)) Cp = SM_INDEXPTRS_S(C);
306 else return (SUN_ERR_MEM_FAIL);
307 ✗ if (SM_INDEXVALS_S(C)) Ci = SM_INDEXVALS_S(C);
308 else return (SUN_ERR_MEM_FAIL);
309 ✗ if (SM_DATA_S(C)) Cx = SM_DATA_S(C);
310 else return (SUN_ERR_MEM_FAIL);
311
312 /* initialize total nonzero count */
313 nz = 0;
314
315 /* iterate through columns (rows for CSR) */
316 ✗ for (j = 0; j < N; j++) {
317
318 /* set current column (row) pointer to current # nonzeros */
319 ✗ Cp[j] = nz;
320
321 /* reset diagonal entry, in case it's not in A */
322 ✗ x[j] = 0.0;
323
324 /* iterate down column (along row) of A, collecting nonzeros */
325 ✗ for (p = Ap[j]; p < Ap[j+1]; p++) {
326 ✗ x[Ai[p]] = Ax[p]; /* collect value */
327 }
328
329 /* add identity to this column (row) */
330 ✗ if (j < M) {
331 ✗ x[j] += c; /* update value */
332 }
333
334 /* fill entries of C with this column's (row's) data */
335 /* fill entries before diagonal */
336 ✗ for (p = Ap[j]; p < Ap[j+1] && Ai[p] < j; p++) {
337 ✗ Ci[nz] = Ai[p];
338 ✗ Cx[nz++] = x[Ai[p]];
339 }
340 /* fill diagonal if applicable */
341 ✗ if (Ai[p] != j || Ap[j] == Ap[j+1]) {
342 ✗ Ci[nz] = j;
343 ✗ Cx[nz++] = x[j];
344 }
345 /* fill entries past diagonal */
346 ✗ for (; p < Ap[j+1]; p++) {
347 ✗ Ci[nz] = Ai[p];
348 ✗ Cx[nz++] = x[Ai[p]];
349 }
350 }
351
352 /* indicate end of data */
353 ✗ Cp[N] = nz;
354
355 /* update A's structure with C's values; nullify C's pointers */
356 ✗ SM_NNZ_S(A) = SM_NNZ_S(C);
357
358 ✗ if (SM_DATA_S(A))
359 ✗ free(SM_DATA_S(A));
360 ✗ SM_DATA_S(A) = SM_DATA_S(C);
361 ✗ SM_DATA_S(C) = NULL;
362
363 ✗ if (SM_INDEXVALS_S(A))
364 ✗ free(SM_INDEXVALS_S(A));
365 ✗ SM_INDEXVALS_S(A) = SM_INDEXVALS_S(C);
366 ✗ SM_INDEXVALS_S(C) = NULL;
367
368 ✗ if (SM_INDEXPTRS_S(A))
369 ✗ free(SM_INDEXPTRS_S(A));
370 ✗ SM_INDEXPTRS_S(A) = SM_INDEXPTRS_S(C);
371 ✗ SM_INDEXPTRS_S(C) = NULL;
372
373 /* clean up */
374 ✗ SUNMatDestroy_Sparse(C);
375 ✗ free(x);
376
377 }
378 return SUN_SUCCESS;
379
380 }
381
382 #endif /* WITH_SUNDIALS */
383