Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 24.3% 119 / 0 / 490
Functions: 43.3% 13 / 0 / 30
Branches: 10.6% 34 / 0 / 322
OMCompiler/SimulationRuntime/c/simulation/jacobian_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 jacobian_util.c
29 */
30
31 #include "jacobian_util.h"
32 #include "options.h"
33 #include "../util/omc_file.h"
34 #include "eval_dep.h"
35 #include "jacobian_colpack.h"
36
37 /**
38 * @brief Initialize analytic jacobian.
39 *
40 * Jacobian has to be allocatd already.
41 *
42 * @param jacobian Jacobian to initialized.
43 * @param sizeCols Number of columns of Jacobian
44 * @param sizeRows Number of rows of Jacobian
45 * @param sizeTmpVars Size of tmp vars array.
46 * @param constantEqns Function pointer for constant equations of Jacobian.
47 * NULL if not available.
48 * @param sparsePattern Pointer to sparsity pattern of Jacobian.
49 */
50 1 void initJacobian(JACOBIAN* jacobian, unsigned int sizeCols, unsigned int sizeRows, unsigned int sizeTmpVars, EVAL_DAG* dag, jacobianColumn_func_ptr evalColumn, jacobianColumn_func_ptr constantEqns, SPARSE_PATTERN* sparsePattern)
51 {
52 /* isRowEval is only known after this call, so make both vectors large enough for
53 * either orientation. For square Jacobians (the common case) this is exact. */
54 1 const unsigned int sizeDirection = sizeCols > sizeRows ? sizeCols : sizeRows;
55
56 1 jacobian->sizeCols = sizeCols;
57 1 jacobian->sizeRows = sizeRows;
58 1 jacobian->sizeTmpVars = sizeTmpVars;
59 1 jacobian->seedVars = (modelica_real*) calloc(sizeDirection, sizeof(modelica_real));
60 1 jacobian->resultVars = (modelica_real*) calloc(sizeDirection, sizeof(modelica_real));
61 1 jacobian->tmpVars = (modelica_real*) calloc(sizeTmpVars, sizeof(modelica_real));
62 1 jacobian->dag = dag;
63 1 jacobian->evalSelection = NULL;
64 1 jacobian->evalColumn = evalColumn;
65 1 jacobian->constantEqns = constantEqns;
66 1 jacobian->sparsePattern = sparsePattern;
67 1 jacobian->availability = JACOBIAN_UNKNOWN;
68 1 jacobian->dae_cj = 0;
69 1 jacobian->isRowEval = FALSE;
70 1 jacobian->cscPattern = NULL;
71 1 jacobian->isBidirectional = FALSE;
72 1 jacobian->adjointJacobian = NULL;
73 1 jacobian->recoverMask = NULL;
74 1 jacobian->csrToCscMap = NULL;
75 1 }
76
77
78 /**
79 * @brief Copy analytic Jacobian.
80 *
81 * Sparsity pattern and DAG are not copied, only their pointers.
82 *
83 * @param source Jacobian that should be copied.
84 * @return JACOBIAN* Copy of source.
85 */
86 ✗ JACOBIAN* copyJacobian(JACOBIAN* source)
87 {
88 ✗ JACOBIAN* jacobian = (JACOBIAN*) malloc(sizeof(JACOBIAN));
89 ✗ initJacobian(jacobian,
90 ✗ source->sizeCols,
91 ✗ source->sizeRows,
92 ✗ source->sizeTmpVars,
93 source->dag,
94 source->evalColumn,
95 source->constantEqns,
96 source->sparsePattern);
97
98 ✗ jacobian->isRowEval = source->isRowEval;
99 ✗ jacobian->isBidirectional = source->isBidirectional;
100 ✗ jacobian->adjointJacobian = source->adjointJacobian; /* shared pointer, not deep copy */
101 ✗ jacobian->recoverMask = source->recoverMask; /* shared pointer, not deep copy */
102 ✗ jacobian->csrToCscMap = source->csrToCscMap; /* shared pointer, not deep copy */
103 ✗ jacobian->cscPattern = NULL; /* not owned by the copy, rebuild on demand */
104
105 ✗ return jacobian;
106 }
107
108 /**
109 * @brief Free memory of analytic Jacobian.
110 *
111 * Also frees sparse pattern.
112 *
113 * @param jac Pointer to Jacobian.
114 */
115 6 void freeJacobian(JACOBIAN *jac)
116 {
117
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 if (jac) {
118 6 free(jac->seedVars); jac->seedVars = NULL;
119 6 free(jac->tmpVars); jac->tmpVars = NULL;
120 6 free(jac->resultVars); jac->resultVars = NULL;
121 6 freeSparsePattern(jac->sparsePattern); jac->sparsePattern = NULL;
122 6 freeSparsePattern(jac->cscPattern); jac->cscPattern = NULL;
123 6 freeEvalDAG(jac->dag); jac->dag = NULL;
124 6 freeEvalSelection(jac->evalSelection); jac->evalSelection = NULL;
125 6 free(jac->recoverMask); jac->recoverMask = NULL;
126 6 free(jac->csrToCscMap); jac->csrToCscMap = NULL;
127 /* adjointJacobian is not owned; do not free */
128 6 jac->adjointJacobian = NULL;
129 6 jac->availability = JACOBIAN_UNKNOWN;
130 }
131 6 }
132
133 /**
134 * @brief Free memory of analytic Jacobian.
135 *
136 * Does not free sparsity pattern and DAG.
137 * Call this for Jacobians that were copied from another Jacobian.
138 *
139 * @param jac Pointer to Jacobian.
140 */
141 ✗ void freeJacobianCopy(JACOBIAN *jac)
142 {
143 ✗ if (jac) {
144 ✗ free(jac->seedVars);
145 ✗ free(jac->tmpVars);
146 ✗ free(jac->resultVars);
147 ✗ freeSparsePattern(jac->cscPattern);
148 ✗ freeEvalSelection(jac->evalSelection);
149 ✗ free(jac);
150 }
151 ✗ }
152
153
154 /*!
155 * \brief Row-wise (adjoint / reverse mode) Jacobian evaluation.
156 *
157 * Assumptions (see JACOBIAN::isRowEval):
158 * - jacobian->evalColumn evaluates a row-direction seed, i.e. it computes s^T * J
159 * - the struct describes J^T: sizeCols == number of rows of J, sizeRows == number of columns of J
160 * - sparsePattern is CSC of J^T (== CSR of J) with row coloring in colorCols
161 *
162 * Output:
163 * - If isDense == false: jac is an nnz-sized buffer. If jacobian->csrToCscMap is set
164 * (see getJacobianCscPattern) the values are written in the CSC
165 * order of J, otherwise in the pattern's own (CSR) order.
166 * - If isDense == true: jac is a dense column-major buffer of J with
167 * J(row, col) stored at jac[col * nRowsJ + row].
168 */
169 ✗ void evalJacobianRow(DATA* data, threadData_t *threadData,
170 JACOBIAN* jacobian, JACOBIAN* parentJacobian,
171 modelica_real* jac, modelica_boolean isDense)
172 {
173 int color, row, col, nz;
174 ✗ const SPARSE_PATTERN* sp = jacobian->sparsePattern;
175 ✗ const unsigned int nRowsJ = jacobian->sizeCols;
176 ✗ const unsigned int nColsJ = jacobian->sizeRows;
177 ✗ const unsigned int* csrToCsc = jacobian->csrToCscMap;
178
179 ✗ if (!jacobian->isRowEval) {
180 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "cant perform row-wise evaluation on column-evaluation Jacobian\n");
181 ✗ return;
182 }
183
184 /* evaluate constant equations of Jacobian (if any) */
185 ✗ if (jacobian->constantEqns != NULL) {
186 ✗ jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
187 }
188
189 /* memset to zero for dense, since solvers might destroy "hard zeros" */
190 ✗ if (isDense) {
191 ✗ memset(jac, 0, (size_t)nRowsJ * (size_t)nColsJ * sizeof(modelica_real));
192 }
193
194 /* Ensure seeds are zeroed before use (one seed per row of J) */
195 ✗ memset(jacobian->seedVars, 0, nRowsJ * sizeof(modelica_real));
196
197 /* evaluate Jacobian row-wise using row-coloring */
198 ✗ for (color = 0; color < (int)sp->maxColors; color++) {
199 /* activate seed variable(s) for the corresponding color (rows) */
200 ✗ for (row = 0; row < (int)nRowsJ; row++) {
201 ✗ if ((int)sp->colorCols[row] - 1 == color) {
202 ✗ jacobian->seedVars[row] = 1.0;
203 }
204 }
205
206 /* evaluate all active rows at once (evalColumn acts as evalRow here) */
207 ✗ jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
208
209 /* scatter results */
210 ✗ for (row = 0; row < (int)nRowsJ; row++) {
211 ✗ if ((int)sp->colorCols[row] - 1 == color) {
212 ✗ for (nz = sp->leadindex[row]; nz < (int)sp->leadindex[row + 1]; nz++) {
213 ✗ col = sp->index[nz];
214 ✗ if (!isDense) {
215 ✗ jac[csrToCsc ? (int)csrToCsc[nz] : nz] = jacobian->resultVars[col];
216 } else {
217 /* dense case, column-major of J */
218 ✗ jac[col * nRowsJ + row] = jacobian->resultVars[col];
219 }
220 }
221 /* de-activate seed variable for the corresponding color (row) */
222 ✗ jacobian->seedVars[row] = 0.0;
223 }
224 }
225
226 /* Row evaluators accumulate adjoints; reset between colors. */
227 ✗ memset(jacobian->resultVars, 0, nColsJ * sizeof(modelica_real));
228 ✗ memset(jacobian->tmpVars, 0, jacobian->sizeTmpVars * sizeof(modelica_real));
229 }
230 }
231
232 /*! \fn evalJacobian
233 *
234 * compute entries of Jacobian in sparse CSC or dense format
235 * uses coloring (sparsePattern non NULL)
236 *
237 * \param [ref] [data]
238 * \param [ref] [threadData]
239 * \param [ref] [jacobian] Pointer to Jacobian
240 * \param [ref] [parentJacobian] Pointer to parent Jacobian
241 * \param [out] [jac] Output buffer, size nnz (sparse) or #rows * #cols (dense), non zero-initialized
242 * \param [ref] [isDense] Flag to set dense / sparse output
243 */
244 ✗ void evalJacobian(DATA* data, threadData_t *threadData, JACOBIAN* jacobian, JACOBIAN* parentJacobian, modelica_real* jac, modelica_boolean isDense)
245 {
246 int color, column, row, nz;
247 ✗ const SPARSE_PATTERN* sp = jacobian->sparsePattern;
248
249 /* Dispatch to bidirectional evaluation if applicable */
250 ✗ if (jacobian->isBidirectional && jacobian->adjointJacobian) {
251 ✗ evalJacobianBidirectional(data, threadData, jacobian, parentJacobian, jac, isDense);
252 ✗ return;
253 }
254
255 ✗ if (jacobian->isRowEval) {
256 ✗ evalJacobianRow(data, threadData, jacobian, parentJacobian, jac, isDense);
257 ✗ return;
258 }
259
260 /* evaluate constant equations of Jacobian */
261 ✗ if (jacobian->constantEqns != NULL) {
262 ✗ jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
263 }
264
265 /* Dense buffer callers use two different conventions:
266 * - NLS/torn-system solvers (nonlinearSolverNewton.c etc.) allocate a square
267 * sizeCols x sizeCols buffer, since sizeRows can exceed sizeCols with
268 * auxiliary residual rows beyond the NLS size that are not part of the
269 * square system to solve.
270 * - Genuinely rectangular Jacobians (e.g. state-selection candidate
271 * matrices in stateset.c, where sizeCols > sizeRows is normal, not an
272 * "auxiliary rows" case) allocate exactly sizeRows * sizeCols.
273 * min(sizeRows, sizeCols) as the stride is correct for both: it equals
274 * sizeCols in the NLS case (matching its square buffer) and sizeRows in the
275 * rectangular case (matching its exact buffer, with every column still
276 * written since column is only ever bounded by sizeCols, never clamped). */
277 ✗ const int denseRows = jacobian->sizeRows < jacobian->sizeCols ? jacobian->sizeRows : jacobian->sizeCols;
278
279 ✗ if (isDense) {
280 ✗ memset(jac, 0, (size_t)denseRows * jacobian->sizeCols * sizeof(modelica_real));
281 }
282
283 ✗ if (!sp) return; /* no sparsity pattern; Jacobian entries cannot be filled */
284
285 /* evaluate Jacobian */
286 ✗ for (color = 0; color < sp->maxColors; color++) {
287 /* activate seed variable for the corresponding color */
288 ✗ for (column = 0; column < jacobian->sizeCols; column++)
289 ✗ if (sp->colorCols[column]-1 == color)
290 ✗ jacobian->seedVars[column] = 1.0;
291
292 /* evaluate Jacobian column */
293 ✗ jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
294 // increaseJacContext(data); // should this be added as is done in ida_solver.c?
295
296 ✗ for (column = 0; column < jacobian->sizeCols; column++) {
297 ✗ if (sp->colorCols[column]-1 == color) {
298 ✗ for (nz = sp->leadindex[column]; nz < sp->leadindex[column+1]; nz++) {
299 ✗ row = sp->index[nz];
300 ✗ if (!isDense) {
301 /* sparse case */
302 ✗ jac[nz] = jacobian->resultVars[row]; //* solverData->xScaling[j];
303 }
304 ✗ else if (row < denseRows) {
305 /* dense case: column-major, denseRows rows per column.
306 * Skip auxiliary rows (row >= denseRows) that lie outside the NLS matrix
307 * (only relevant when sizeRows > sizeCols; never true for rectangular
308 * Jacobians where denseRows == sizeRows). */
309 ✗ jac[column * denseRows + row] = jacobian->resultVars[row];
310 }
311 }
312 /* de-activate seed variable for the corresponding color */
313 ✗ jacobian->seedVars[column] = 0.0;
314 }
315 }
316 }
317 }
318
319 /**
320 * @brief Initialize bidirectional recovery masks for star bicoloring.
321 *
322 * For each nonzero in forward (CSC) and adjoint (CSR) patterns, determines
323 * whether the entry is recoverable from the respective direction.
324 * Also computes CSR-to-CSC index mapping for sparse output.
325 *
326 * Must be called after both jacobians are fully initialized (patterns + colors)
327 * and linked (fwd->adjointJacobian != NULL).
328 *
329 * @param fwd Forward jacobian with CSC pattern + column coloring.
330 */
331 ✗ void initBidirectionalRecovery(JACOBIAN* fwd)
332 {
333 ✗ JACOBIAN* adj = fwd->adjointJacobian;
334 ✗ if (!adj) return;
335
336 ✗ const SPARSE_PATTERN* fwdsp = fwd->sparsePattern;
337 ✗ const SPARSE_PATTERN* adjsp = adj->sparsePattern;
338 ✗ const unsigned int nCols = fwd->sizeCols;
339 ✗ const unsigned int nRows = fwd->sizeRows;
340 ✗ const unsigned int nnz = fwdsp->nnz;
341 unsigned int j, i, nz, k, j2, i2;
342
343 ✗ fwd->recoverMask = (unsigned char*) calloc(nnz, sizeof(unsigned char));
344 ✗ adj->recoverMask = (unsigned char*) calloc(nnz, sizeof(unsigned char));
345 ✗ adj->csrToCscMap = (unsigned int*) calloc(nnz, sizeof(unsigned int));
346
347 /* Forward recoverMask: entry (i,j) is column-recoverable if j is the ONLY
348 * column with its column color among all columns having a nonzero in row i. */
349 // iterate over all columns
350 ✗ for (j = 0; j < nCols; j++) {
351 ✗ unsigned int cj = fwdsp->colorCols[j];
352 // iterate over nonzeros (rows with nonzero) in this column via forward CSC pattern
353 ✗ for (nz = fwdsp->leadindex[j]; nz < fwdsp->leadindex[j+1]; nz++) {
354 ✗ i = fwdsp->index[nz]; // row index of current nonzero
355 int unique = 1; // assume current column is unique for this nonzero until we find otherwise
356 // check all other columns with nonzero in the same row i via adjoint CSR pattern
357 ✗ for (k = adjsp->leadindex[i]; k < adjsp->leadindex[i+1]; k++) {
358 ✗ j2 = adjsp->index[k]; // column index of nonzero in same row
359 // check its a different column and has the same color, if so current column is not unique for this nonzero
360 ✗ if (j2 != j && fwdsp->colorCols[j2] == cj) {
361 unique = 0;
362 break;
363 }
364 }
365 // mark as unique (column-recoverable) or not
366 // if unique, this nonzero can be recovered from forward evaluation when column j is seeded, otherwise it cannot and must be recovered from adjoint evaluation
367 // if not unique, it gives a wrong value when recovered from forward evaluation and thus can not be written into result vector
368 ✗ fwd->recoverMask[nz] = (unsigned char)unique;
369 }
370 }
371
372 /* Adjoint recoverMask: entry (i,j) is row-recoverable if i is the ONLY
373 * row with its row color among all rows having a nonzero in column j. */
374 // same logic as forward, but now iterate over rows and check uniqueness of row color among rows with nonzero in same column via forward pattern
375 ✗ for (i = 0; i < nRows; i++) {
376 ✗ unsigned int ri = adjsp->colorCols[i];
377 ✗ for (nz = adjsp->leadindex[i]; nz < adjsp->leadindex[i+1]; nz++) {
378 ✗ j = adjsp->index[nz];
379 int unique = 1;
380 ✗ for (k = fwdsp->leadindex[j]; k < fwdsp->leadindex[j+1]; k++) {
381 ✗ i2 = fwdsp->index[k];
382 ✗ if (i2 != i && adjsp->colorCols[i2] == ri) {
383 unique = 0;
384 break;
385 }
386 }
387 ✗ adj->recoverMask[nz] = (unsigned char)unique;
388 }
389 }
390
391 /* CSR-to-CSC mapping: for each adjoint CSR position (nonzero), find forward CSC position */
392 // iterate over all rows
393 ✗ for (i = 0; i < nRows; i++) {
394 // iterate over all nonzeros in this row via adjoint CSR pattern
395 ✗ for (nz = adjsp->leadindex[i]; nz < adjsp->leadindex[i+1]; nz++) {
396 ✗ j = adjsp->index[nz]; // get column index of current nonzero
397 ✗ adj->csrToCscMap[nz] = 0;
398 // iterate over all nonzeros in this column via forward CSC pattern, so the nonzero rows
399 ✗ for (k = fwdsp->leadindex[j]; k < fwdsp->leadindex[j+1]; k++) {
400 // if row index matches, we found the same nonzero in forward pattern
401 // and can record its position k for later indexing into forward result vector when recovering this nonzero from adjoint evaluation
402 ✗ if (fwdsp->index[k] == i) {
403 ✗ adj->csrToCscMap[nz] = k;
404 ✗ break;
405 }
406 }
407 }
408 }
409 }
410
411 /**
412 * @brief Evaluate Jacobian using bidirectional (star bicoloring) approach.
413 *
414 * Uses both forward (column) and adjoint (row) evaluations to recover all
415 * nonzero entries with fewer total colors than unidirectional coloring.
416 *
417 * Dense output: column-major jac[col * nRows + row].
418 * Sparse output: CSC-indexed jac[nz] matching forward sparse pattern.
419 *
420 * @param data Runtime data struct.
421 * @param threadData Thread data for error handling.
422 * @param fwd Forward jacobian (isBidirectional=TRUE, adjointJacobian set).
423 * @param parentJacobian Parent Jacobian for nested use (can be NULL).
424 * @param jac Output buffer.
425 * @param isDense TRUE for dense, FALSE for sparse CSC.
426 */
427 ✗ void evalJacobianBidirectional(DATA* data, threadData_t *threadData,
428 JACOBIAN* fwd, JACOBIAN* parentJacobian,
429 modelica_real* jac, modelica_boolean isDense)
430 {
431 ✗ JACOBIAN* adj = fwd->adjointJacobian;
432 ✗ const SPARSE_PATTERN* fwdsp = fwd->sparsePattern;
433 ✗ const SPARSE_PATTERN* adjsp = adj->sparsePattern;
434 ✗ const int nRows = (int)fwd->sizeRows;
435 ✗ const int nCols = (int)fwd->sizeCols;
436 int color, column, row, nz, j;
437
438 ✗ if (fwd->constantEqns) fwd->constantEqns(data, threadData, fwd, parentJacobian);
439 ✗ if (adj->constantEqns) adj->constantEqns(data, threadData, adj, parentJacobian);
440
441 ✗ if (isDense) {
442 ✗ memset(jac, 0, (size_t)nRows * (size_t)nCols * sizeof(modelica_real));
443 }
444
445 /* Column phase (forward mode, CSC + column coloring) */
446 ✗ for (color = 0; color < (int)fwdsp->maxColors; color++) {
447 ✗ for (column = 0; column < nCols; column++)
448 ✗ if ((int)fwdsp->colorCols[column] - 1 == color)
449 ✗ fwd->seedVars[column] = 1.0;
450
451 ✗ fwd->evalColumn(data, threadData, fwd, parentJacobian);
452
453 ✗ for (column = 0; column < nCols; column++) {
454 ✗ if ((int)fwdsp->colorCols[column] - 1 == color) {
455 ✗ for (nz = (int)fwdsp->leadindex[column]; nz < (int)fwdsp->leadindex[column + 1]; nz++) {
456 ✗ if (fwd->recoverMask[nz]) {
457 ✗ row = (int)fwdsp->index[nz];
458 ✗ if (isDense)
459 ✗ jac[column * nRows + row] = fwd->resultVars[row];
460 else
461 ✗ jac[nz] = fwd->resultVars[row];
462 }
463 }
464 ✗ fwd->seedVars[column] = 0.0;
465 }
466 }
467 }
468
469 /* Row phase (adjoint mode, CSR + row coloring) */
470 ✗ for (color = 0; color < (int)adjsp->maxColors; color++) {
471 ✗ for (row = 0; row < nRows; row++)
472 ✗ if ((int)adjsp->colorCols[row] - 1 == color)
473 ✗ adj->seedVars[row] = 1.0;
474
475 ✗ adj->evalColumn(data, threadData, adj, parentJacobian);
476
477 ✗ for (row = 0; row < nRows; row++) {
478 ✗ if ((int)adjsp->colorCols[row] - 1 == color) {
479 ✗ for (nz = (int)adjsp->leadindex[row]; nz < (int)adjsp->leadindex[row + 1]; nz++) {
480 ✗ if (adj->recoverMask[nz]) {
481 ✗ column = (int)adjsp->index[nz];
482 ✗ if (isDense)
483 ✗ jac[column * nRows + row] = adj->resultVars[column];
484 else
485 ✗ jac[adj->csrToCscMap[nz]] = adj->resultVars[column];
486 }
487 }
488 ✗ adj->seedVars[row] = 0.0;
489 }
490 }
491 /* Reset adjoint result vars to zero after reading to prevent accumulation across colors */
492 ✗ memset(adj->resultVars, 0, (size_t)nCols * sizeof(modelica_real));
493 // also for tmp vars
494 ✗ memset(adj->tmpVars, 0, (size_t)adj->sizeTmpVars * sizeof(modelica_real));
495 }
496 ✗ }
497
498 /**
499 * @brief Compute Jacobian-Vector product y = J * s.
500 *
501 * @param data Runtime data struct.
502 * @param threadData Thread data for error handling.
503 * @param jacobian Jacobian object (must have evalColumn and sparsePattern set).
504 * @param parentJacobian Parent Jacobian (if nested), can be NULL.
505 * @param seed Input seed vector s, length = jacobian->sizeCols.
506 * @param out Output vector y, length = jacobian->sizeRows.
507 * @param zero_out If true, zero-initialize out before accumulation.
508 */
509 ✗ void jvp(DATA* data, threadData_t *threadData,
510 JACOBIAN* jacobian, JACOBIAN* parentJacobian,
511 const modelica_real* seed, modelica_real* out,
512 modelica_boolean zero_out)
513 {
514 ✗ if (jacobian->isRowEval) {
515 /* Error: jvp called on row-evaluation Jacobian */
516 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "cant perform jvp on row-evaluation Jacobian\n");
517 ✗ return;
518 }
519 ✗ const unsigned int nCols = jacobian->sizeCols;
520 ✗ const unsigned int nRows = jacobian->sizeRows;
521
522 /* Optional: zero output before accumulation */
523 ✗ if (zero_out) {
524 ✗ memset(out, 0, nRows * sizeof(modelica_real));
525 }
526
527 /* Ensure seeds are zeroed before use */
528 ✗ memset(jacobian->seedVars, 0, nCols * sizeof(modelica_real));
529
530 /* Evaluate constant equations (if any) */
531 ✗ if (jacobian->constantEqns != NULL) {
532 ✗ jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
533 }
534
535 /* Set all seeds */
536 ✗ for (unsigned int col = 0; col < nCols; col++) {
537 ✗ jacobian->seedVars[col] = seed[col];
538 }
539
540 /* Evaluate J * s into resultVars */
541 ✗ jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
542
543 /* Accumulate results into out */
544 ✗ for (unsigned int row = 0; row < nRows; row++) {
545 ✗ out[row] += jacobian->resultVars[row];
546 }
547 }
548
549
550 /**
551 * @brief Compute Vector-Jacobian product y = J^T * s.
552 *
553 * @param data Runtime data struct.
554 * @param threadData Thread data for error handling.
555 * @param jacobian Jacobian object (must have evalColumn and sparsePattern set).
556 * @param parentJacobian Parent Jacobian (if nested), can be NULL.
557 * @param seed Input seed vector s, length = jacobian->sizeRows.
558 * @param out Output vector y, length = jacobian->sizeCols.
559 * @param zero_out If true, zero-initialize out before accumulation.
560 */
561 ✗ void vjp(DATA* data, threadData_t *threadData,
562 JACOBIAN* jacobian, JACOBIAN* parentJacobian,
563 const modelica_real* seed, modelica_real* out,
564 modelica_boolean zero_out)
565 {
566 ✗ if (!jacobian->isRowEval) {
567 /* Error: vjp called on column-evaluation Jacobian */
568 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "cant perform vjp on column-evaluation Jacobian\n");
569 ✗ return;
570 }
571 ✗ const unsigned int nCols = jacobian->sizeCols;
572 ✗ const unsigned int nRows = jacobian->sizeRows;
573
574 /* Optional: zero output before accumulation */
575 ✗ if (zero_out) {
576 ✗ memset(out, 0, nCols * sizeof(modelica_real));
577 }
578
579 /* Ensure seeds are zeroed before use */
580 ✗ memset(jacobian->seedVars, 0, nRows * sizeof(modelica_real));
581
582 /* Evaluate constant equations (if any) */
583 ✗ if (jacobian->constantEqns != NULL) {
584 ✗ jacobian->constantEqns(data, threadData, jacobian, parentJacobian);
585 }
586
587 /* Set all seeds */
588 ✗ for (unsigned int row = 0; row < nRows; row++) {
589 ✗ jacobian->seedVars[row] = seed[row];
590 }
591
592 /* Evaluate J * s into resultVars */
593 // this is actually evalRow
594 ✗ jacobian->evalColumn(data, threadData, jacobian, parentJacobian);
595
596 /* Accumulate results into out */
597 ✗ for (unsigned int col = 0; col < nCols; col++) {
598 ✗ out[col] += jacobian->resultVars[col];
599 }
600 }
601
602 /**
603 * @brief Allocate memory for sparsity pattern.
604 *
605 * @param n_leadIndex Number of rows or columns of Matrix.
606 * Depending on compression type CSR (-->rows) or CSC (-->columns).
607 * @param nnz Number of non-zero elements in Matrix.
608 * @param maxColors Maximum number of colors of Matrix.
609 * @return SPARSE_PATTERN* Pointer to allocated sparsity pattern of Matrix.
610 */
611 1 SPARSE_PATTERN* allocSparsePattern(unsigned int n_leadIndex, unsigned int nnz, unsigned int maxColors)
612 {
613 1 SPARSE_PATTERN* sparsePattern = (SPARSE_PATTERN*) malloc(sizeof(SPARSE_PATTERN));
614 1 sparsePattern->nnz = nnz;
615 1 sparsePattern->leadindex = (unsigned int*) malloc((n_leadIndex+1)*sizeof(unsigned int));
616 1 sparsePattern->index = (unsigned int*) malloc(nnz*sizeof(unsigned int));
617 1 sparsePattern->colorCols = (unsigned int*) malloc(n_leadIndex*sizeof(unsigned int));
618 1 sparsePattern->maxColors = maxColors;
619 1 sparsePattern->sizeCols = n_leadIndex;
620
621 1 return sparsePattern;
622 }
623
624
625 /**
626 * @brief Transpose a compressed sparse pattern.
627 *
628 * Works in both directions, since CSC of A is CSR of A^T:
629 * - CSC -> CSR: transposeSparsePattern(csc, nRows, nCols, ...)
630 * - CSR -> CSC: transposeSparsePattern(csr, nCols, nRows, ...)
631 *
632 * Input (A):
633 * - Ap = in->leadindex (size nLeadIn+1), lead pointers
634 * - Ai = in->index (size nnz), secondary indices in [0, nLeadOut)
635 *
636 * Output (B):
637 * - Bp = out->leadindex (size nLeadOut+1)
638 * - Bj = out->index (size nnz), sorted ascending within each lead slot
639 *
640 * @param in Input pattern.
641 * @param nLeadOut Number of lead slots of the result (== range of in->index). For CSC its nRows, for CSR its nCols
642 * @param nLeadIn Number of lead slots of the input. For CSC its nCols, for CSR its nRows
643 * @param nzMap If non-NULL, receives a newly allocated array of size nnz mapping
644 * each nonzero position of `in` to its position in the result.
645 * @return SPARSE_PATTERN* Transposed pattern without coloring, or NULL on error.
646 *
647 * Complexity: O(nnz + max(nLeadIn, nLeadOut))
648 */
649 ✗ SPARSE_PATTERN* transposeSparsePattern(const SPARSE_PATTERN* in,
650 unsigned int nLeadOut,
651 unsigned int nLeadIn,
652 unsigned int** nzMap)
653 {
654 ✗ if (!in) return NULL;
655
656 ✗ const unsigned int nnz = in->nnz;
657
658 /* Allocate result pattern: leadindex size = nLeadOut+1, index size = nnz */
659 ✗ SPARSE_PATTERN* out = allocSparsePattern(nLeadOut, nnz, /*maxColors*/ 0);
660 ✗ if (!out) return NULL;
661
662 // only fill map if nzMap is not NULL
663 unsigned int* map = NULL;
664 ✗ if (nzMap) {
665 ✗ map = (unsigned int*) malloc((nnz ? nnz : 1) * sizeof(unsigned int));
666 ✗ if (!map) {
667 ✗ freeSparsePattern(out);
668 ✗ return NULL;
669 }
670 }
671
672 /* Aliases for conciseness */
673 ✗ const unsigned int* Ap = in->leadindex;
674 ✗ const unsigned int* Ai = in->index;
675 ✗ unsigned int* Bp = out->leadindex;
676 ✗ unsigned int* Bj = out->index;
677
678 /* 1) Count nnz per output lead slot */
679 ✗ memset(Bp, 0, (nLeadOut+1) * sizeof(unsigned int));
680 ✗ for (unsigned int k = 0; k < nnz; k++) {
681 ✗ if (Ai[k] >= nLeadOut) {
682 /* Out of bounds. Clean up and abort. */
683 ✗ freeSparsePattern(out);
684 ✗ free(map);
685 ✗ return NULL;
686 }
687 ✗ Bp[Ai[k]]++;
688 }
689
690 /* 2) Exclusive prefix sum over Bp to get lead pointers; set Bp[nLeadOut] = nnz */
691 {
692 unsigned int presum = 0;
693 ✗ for (unsigned int r = 0; r < nLeadOut; r++) {
694 ✗ const unsigned int tmp = Bp[r];
695 ✗ Bp[r] = presum;
696 ✗ presum += tmp;
697 }
698 ✗ Bp[nLeadOut] = nnz;
699 }
700
701 /* 3) Fill result indices using running heads in Bp.
702 * Iterating the input lead slots in ascending order yields ascending
703 * secondary indices in the result, which KLU and printSparseStructure require. */
704 ✗ for (unsigned int lead = 0; lead < nLeadIn; lead++) {
705 ✗ const unsigned int start = Ap[lead];
706 ✗ const unsigned int stop = Ap[lead + 1];
707 ✗ if (stop < start || stop > nnz) {
708 /* Corrupt pointers. Clean up and abort. */
709 ✗ freeSparsePattern(out);
710 ✗ free(map);
711 ✗ return NULL;
712 }
713 ✗ for (unsigned int jj = start; jj < stop; jj++) {
714 ✗ const unsigned int dest = Bp[Ai[jj]]; /* next free slot */
715 ✗ Bj[dest] = lead;
716 ✗ if (map) map[jj] = dest;
717 ✗ Bp[Ai[jj]]++; /* advance head */
718 }
719 }
720
721 /* 4) Restore Bp to lead pointers by shifting heads back */
722 {
723 unsigned int last = 0;
724 ✗ for (unsigned int r = 0; r <= nLeadOut; r++) {
725 ✗ const unsigned int tmp = Bp[r];
726 ✗ Bp[r] = last;
727 last = tmp;
728 }
729 }
730
731 /* We don't have a coloring for the transposed pattern; keep defaults. */
732 ✗ out->maxColors = 0;
733 ✗ memset(out->colorCols, 0, nLeadOut * sizeof(unsigned int));
734
735 ✗ if (nzMap) *nzMap = map;
736 return out;
737 }
738
739 /**
740 * @brief Convert a CSC-format sparsity pattern to CSR-format.
741 *
742 * @param csc Pattern with column pointers and row indices.
743 * @param nRows Number of rows of the matrix.
744 * @param nCols Number of columns of the matrix.
745 * @return SPARSE_PATTERN* CSR pattern with row pointers and column indices.
746 */
747 ✗ SPARSE_PATTERN* cscToCsr(const SPARSE_PATTERN* csc,
748 unsigned int nRows,
749 unsigned int nCols)
750 {
751 ✗ return transposeSparsePattern(csc, nRows, nCols, NULL);
752 }
753
754 /**
755 * @brief Get the column oriented (CSC of J) sparsity pattern of a Jacobian.
756 *
757 * Forward Jacobians already store their pattern in CSC, so the pattern is returned
758 * as is. Row evaluated (adjoint) Jacobians store CSR of J; for those the transposed
759 * pattern and the CSR->CSC nonzero mapping are built once and cached, so that
760 * evalJacobian() can emit sparse values directly in the CSC order expected by the
761 * solvers (KLU / SUNDIALS / GBODE).
762 *
763 * @param jac Jacobian with an initialized sparsity pattern.
764 * @return SPARSE_PATTERN* Column oriented pattern (not owned by the caller).
765 */
766 2 SPARSE_PATTERN* getJacobianCscPattern(JACOBIAN* jac)
767 {
768
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 if (jac == NULL || jac->sparsePattern == NULL) {
769 return NULL;
770 }
771 // If the Jacobian is not row-evaluated, it already has a CSC pattern.
772
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (!jac->isRowEval) {
773 return jac->sparsePattern;
774 }
775 // If the Jacobian is row-evaluated, we need to transpose the CSR pattern to get CSC and store it in jac->cscPattern.
776 ✗ if (jac->cscPattern == NULL) {
777 ✗ unsigned int* map = NULL;
778 ✗ jac->cscPattern = transposeSparsePattern(jac->sparsePattern,
779 ✗ (unsigned int) jac->sizeRows,
780 ✗ (unsigned int) jac->sizeCols,
781 &map);
782 ✗ if (jac->cscPattern == NULL) {
783 ✗ free(map);
784 ✗ return jac->sparsePattern;
785 }
786 ✗ if (jac->csrToCscMap == NULL) {
787 ✗ jac->csrToCscMap = map;
788 } else {
789 /* Already set up by initBidirectionalRecovery, keep it. */
790 ✗ free(map);
791 }
792 /* The transposed pattern has no coloring yet; derive one so that the pattern is
793 * usable wherever a fully featured column oriented pattern is expected. */
794 ✗ computeColumnColoring(jac->cscPattern,
795 ✗ (unsigned int) jac->sizeCols,
796 ✗ (unsigned int) jac->sizeRows);
797 }
798 ✗ return jac->cscPattern;
799 }
800
801
802 /**
803 * @brief Free sparsity pattern
804 *
805 * @param spp Pointer to sparsity pattern
806 */
807 12 void freeSparsePattern(SPARSE_PATTERN *spp)
808 {
809
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 11 times.
 12 if (spp) {
810 1 free(spp->index);
811 1 free(spp->colorCols);
812 1 free(spp->leadindex);
813 1 free(spp);
814 }
815 12 }
816
817 /**
818 * @brief Distance-1 column coloring of a CSC sparse pattern.
819 *
820 * Two columns may share a color only if they have no non-zero row in common.
821 * The rows of every column are sorted and made unique first.
822 * Uses ColPack's partial distance-two column coloring when available, with
823 * a greedy C-only fallback.
824 * The fallback uses the existing cscToCsr helper to build the row→columns map, then
825 * assigns the smallest available color to each column in order.
826 *
827 * Needed for the resizable analytic Jacobian path: the C sparsity pattern
828 * is built at runtime from WHOLEDIM loops that over-approximate array
829 * equations as dense blocks, so the compile-time coloring (derived from the
830 * exact symbolic sparsity) is invalid for the runtime pattern. Recomputing
831 * it here guarantees correctness.
832 *
833 * @param sp CSC sparse pattern (leadindex, index, colorCols already allocated).
834 * @param nRows Number of rows in the Jacobian.
835 * @param nCols Number of columns (== size of sp->colorCols).
836 */
837 ✗ static int compareUnsigned(const void* a, const void* b)
838 {
839 ✗ const unsigned int x = *(const unsigned int*) a, y = *(const unsigned int*) b;
840 ✗ return (x > y) - (x < y);
841 }
842
843 /**
844 * @brief Sorts the rows of every column and removes duplicates in place.
845 *
846 * Runtime built patterns can contain the same entry twice, e.g. an array seed
847 * over a whole dimension and one of its elements in the same row.
848 */
849 1 static void sortUniqueSparsePattern(SPARSE_PATTERN* sp, unsigned int nCols)
850 {
851 unsigned int col, nz, start, end, out = 0;
852
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
 3 for (col = 0; col < nCols; col++) {
853 2 start = sp->leadindex[col];
854 2 end = sp->leadindex[col + 1];
855 2 qsort(sp->index + start, end - start, sizeof(unsigned int), compareUnsigned);
856 2 sp->leadindex[col] = out;
857
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 for (nz = start; nz < end; nz++) {
858
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if (nz == start || sp->index[nz] != sp->index[nz - 1]) {
859 2 sp->index[out++] = sp->index[nz];
860 }
861 }
862 }
863 1 sp->leadindex[nCols] = out;
864 1 sp->nnz = out;
865 1 }
866
867 1 void computeColumnColoring(SPARSE_PATTERN* sp, unsigned int nRows, unsigned int nCols)
868 {
869
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 1 if (!sp || !sp->colorCols) return;
870
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (nCols == 0) {
871 ✗ sp->maxColors = 0;
872 ✗ return;
873 }
874 1 sortUniqueSparsePattern(sp, nCols);
875
876 #if defined(OMC_HAVE_COLPACK)
877
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (computeColPackColumnColoring(
878 1 nRows, nCols, sp->leadindex, sp->index, sp->nnz, sp->colorCols, &sp->maxColors) == 0) {
879 return;
880 }
881 #endif
882
883 ✗ SPARSE_PATTERN* csr = cscToCsr(sp, nRows, nCols);
884 ✗ if (!csr) {
885 /* Fallback: trivial one-column-per-color coloring. */
886 ✗ for (unsigned int c = 0; c < nCols; c++) sp->colorCols[c] = c + 1;
887 ✗ sp->maxColors = nCols;
888 ✗ return;
889 }
890
891 /* forbidden[k] == 1 if color k is already used by an adjacent column.
892 * Index 0 unused; colors are 1-based, max is nCols. */
893 ✗ unsigned char* forbidden = (unsigned char*) calloc(nCols + 2, sizeof(unsigned char));
894 /* Track which forbidden slots were set so we can reset without a full memset. */
895 ✗ unsigned int* setColors = (unsigned int*) malloc(nCols * sizeof(unsigned int));
896
897 ✗ if (!forbidden || !setColors) {
898 ✗ free(forbidden); free(setColors);
899 ✗ freeSparsePattern(csr);
900 ✗ for (unsigned int c = 0; c < nCols; c++) sp->colorCols[c] = c + 1;
901 ✗ sp->maxColors = nCols;
902 ✗ return;
903 }
904
905 unsigned int maxColor = 0;
906
907 ✗ for (unsigned int c = 0; c < nCols; c++) {
908 unsigned int nSet = 0;
909
910 /* Mark colors of already-colored columns that share a row with c. */
911 ✗ for (unsigned int nz = sp->leadindex[c]; nz < sp->leadindex[c + 1]; nz++) {
912 ✗ const unsigned int row = sp->index[nz];
913 ✗ if (row >= nRows) continue;
914 ✗ for (unsigned int nz2 = csr->leadindex[row]; nz2 < csr->leadindex[row + 1]; nz2++) {
915 ✗ const unsigned int c2 = csr->index[nz2];
916 ✗ if (c2 < c) {
917 ✗ const unsigned int used = sp->colorCols[c2];
918 ✗ if (used > 0 && used <= nCols && !forbidden[used]) {
919 ✗ forbidden[used] = 1;
920 ✗ setColors[nSet++] = used;
921 }
922 }
923 }
924 }
925
926 /* Smallest color not forbidden. */
927 unsigned int color = 1;
928 ✗ while (color <= nCols && forbidden[color]) color++;
929 ✗ sp->colorCols[c] = color;
930 if (color > maxColor) maxColor = color;
931
932 /* Reset forbidden markers for next iteration. */
933 ✗ for (unsigned int k = 0; k < nSet; k++) forbidden[setColors[k]] = 0;
934 }
935
936 ✗ sp->maxColors = maxColor;
937
938 ✗ free(setColors);
939 ✗ free(forbidden);
940 ✗ freeSparsePattern(csr);
941 }
942
943 /**
944 * @brief Sort row indices within each column of a CSC sparse pattern.
945 *
946 * KLU and printSparseStructure both require that row indices within each
947 * column are in strictly ascending order. The NBackend-generated
948 * initialResizableAnalyticJacobianA fills entries in equation order which
949 * may not be sorted (e.g. column 0 gets row 10 from one equation and row 0
950 * from another). Call this function once after the pattern is built and
951 * before it is handed to KLU or the print helpers.
952 *
953 * @param sp CSC sparse pattern (leadindex and index already filled).
954 * @param nCols Number of columns (== size of sp->leadindex - 1).
955 */
956 1 void sortSparseColumns(SPARSE_PATTERN* sp, unsigned int nCols)
957 {
958
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (!sp) return;
959
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
 3 for (unsigned int c = 0; c < nCols; c++) {
960 2 unsigned int start = sp->leadindex[c];
961 2 unsigned int end = sp->leadindex[c + 1];
962
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (end <= start + 1) continue;
963 /* Insertion sort — columns typically have very few entries. */
964 ✗ for (unsigned int i = start + 1; i < end; i++) {
965 ✗ unsigned int key = sp->index[i];
966 unsigned int j = i;
967 ✗ while (j > start && sp->index[j - 1] > key) {
968 ✗ sp->index[j] = sp->index[j - 1];
969 j--;
970 }
971 ✗ sp->index[j] = key;
972 }
973 }
974 }
975
976 /**
977 * @brief Opens sparsity pattern file
978 *
979 * @param data Runtime data struct.
980 * @param threadData Thread data for error handling.
981 * @param filename String for the filename.
982 * @return FILE* Pointer to sparsity pattern stream.
983 */
984 ✗ FILE * openSparsePatternFile(DATA* data, threadData_t *threadData, const char* filename)
985 {
986 FILE* pFile;
987 ✗ const char* fullPath = NULL;
988
989 ✗ if (omc_flag[FLAG_INPUT_PATH]) {
990 ✗ GC_asprintf(&fullPath, "%s/%s", omc_flagValue[FLAG_INPUT_PATH], filename);
991 ✗ } else if (data->modelData->resourcesDir) {
992 ✗ GC_asprintf(&fullPath, "%s/%s", data->modelData->resourcesDir, filename);
993 } else {
994 ✗ GC_asprintf(&fullPath, "%s", filename);
995 }
996 ✗ pFile = omc_fopen(fullPath, "rb");
997 ✗ if (pFile == NULL) {
998 ✗ throwStreamPrint(threadData, "Could not open sparsity pattern file %s.", fullPath);
999 }
1000 ✗ omc_rc_release((void*) fullPath);
1001 ✗ return pFile;
1002 }
1003
1004 /**
1005 * @brief Reads one color of sparsity pattern and sets colorCols.
1006 *
1007 * @param threadData Used for error handling.
1008 * @param pFile Pointer to file stream.
1009 * @param colorCols Array of column coloring.
1010 * @param color Current color index.
1011 * @param length Number of columns in color `color`.
1012 */
1013 ✗ void readSparsePatternColor(threadData_t* threadData, FILE * pFile, unsigned int* colorCols, unsigned int color, unsigned int length, unsigned int maxIndex)
1014 {
1015 unsigned int i, index;
1016 size_t count;
1017
1018 ✗ for (i = 0; i < length; i++) {
1019 ✗ count = omc_fread(&index, sizeof(unsigned int), 1, pFile, FALSE);
1020 ✗ if (count != 1) {
1021 ✗ throwStreamPrint(threadData, "Error while reading color %u of sparsity pattern.", color);
1022 }
1023 ✗ if (index < 0 || index >= maxIndex) {
1024 ✗ throwStreamPrint(threadData, "Error while reading color %u of sparsity pattern. Index %d out of bounds", color, index);
1025 }
1026 ✗ colorCols[index] = color;
1027 }
1028 ✗ }
1029
1030 /**
1031 * @brief Read the Jacobian method requested by the user via flag `-jacobian`.
1032 *
1033 * Performs no availability check, this is done in checkJacobianMethod().
1034 *
1035 * @param threadData Used for error handling.
1036 * @return JACOBIAN_METHOD Requested method or JAC_UNKNOWN if the flag is not set.
1037 */
1038 1 JACOBIAN_METHOD getRequestedJacobianMethod(threadData_t* threadData)
1039 {
1040 JACOBIAN_METHOD jacobianMethod = JAC_UNKNOWN;
1041
1042 // Check if the user requested a specific Jacobian method via the `-jacobian` flag
1043 // if not the colored numerical Jacobian is used by default
1044
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (!omc_flag[FLAG_JACOBIAN]) {
1045 return JAC_UNKNOWN;
1046 }
1047 // Set the requested method if it is known
1048 ✗ for (int method=1; method < JAC_MAX; method++) {
1049 ✗ if (!strcmp(omc_flagValue[FLAG_JACOBIAN], JACOBIAN_METHOD_NAME[method])) {
1050 ✗ jacobianMethod = (JACOBIAN_METHOD) method;
1051 ✗ break;
1052 }
1053 }
1054 // Error case if the user requested a method that is not known
1055 ✗ if (jacobianMethod == JAC_UNKNOWN) {
1056 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Unknown value `%s` for flag `-jacobian`", omc_flagValue[FLAG_JACOBIAN]);
1057 ✗ infoStreamPrint(OMC_LOG_STDOUT, 1, "Available options are");
1058 ✗ for (int method=1; method < JAC_MAX; method++) {
1059 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "%s", JACOBIAN_METHOD_NAME[method]);
1060 }
1061 ✗ messageClose(OMC_LOG_STDOUT);
1062 ✗ omc_throw(threadData);
1063 }
1064 return jacobianMethod;
1065 }
1066
1067 /**
1068 * @brief Check that the requested Jacobian method can be used and log it.
1069 *
1070 * @param threadData Used for error handling.
1071 * @param availability Is the symbolic Jacobian available, only the sparsity pattern available or nothing available.
1072 * @param jacobianMethod Requested method, JAC_UNKNOWN selects the default for `availability`.
1073 * @return JACOBIAN_METHOD Jacobian method that will be used.
1074 */
1075 1 JACOBIAN_METHOD checkJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILITY availability, JACOBIAN_METHOD jacobianMethod)
1076 {
1077
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 assertStreamPrint(threadData, availability != JACOBIAN_UNKNOWN, "Jacobian availability status is unknown.");
1078
1079 /* Check if method is available. If all is fine then no case gets triggered. */
1080
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 switch (availability)
1081 {
1082 ✗ case JACOBIAN_NOT_AVAILABLE:
1083 ✗ if (jacobianMethod != INTERNALNUMJAC && jacobianMethod != JAC_UNKNOWN) {
1084 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Jacobian not available, switching to internal numerical Jacobian.");
1085 }
1086 jacobianMethod = INTERNALNUMJAC;
1087 break;
1088 1 case JACOBIAN_ONLY_SPARSITY:
1089
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (jacobianMethod == COLOREDSYMJAC || jacobianMethod == COLOREDSYMJACADJ || jacobianMethod == BICOLOREDSYMJAC) {
1090 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Symbolic Jacobian not available, only sparsity pattern. Switching to colored numerical Jacobian.");
1091 jacobianMethod = COLOREDNUMJAC;
1092
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 } else if(jacobianMethod == SYMJAC) {
1093 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Symbolic Jacobian not available, only sparsity pattern. Switching to uncolored numerical Jacobian.");
1094 jacobianMethod = NUMJAC;
1095 } else if(jacobianMethod == JAC_UNKNOWN) {
1096 jacobianMethod = COLOREDNUMJAC;
1097 }
1098 break;
1099 ✗ case JACOBIAN_AVAILABLE:
1100 ✗ if (jacobianMethod == JAC_UNKNOWN) {
1101 jacobianMethod = COLOREDSYMJAC;
1102 }
1103 break;
1104 ✗ default:
1105 ✗ throwStreamPrint(threadData, "Unhandled case in setJacobianMethod");
1106 break;
1107 }
1108
1109 /* Log Jacobian method */
1110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 assertStreamPrint(threadData, jacobianMethod > JAC_UNKNOWN && jacobianMethod < JAC_MAX, "Unhandled case in setJacobianMethod");
1111 1 infoStreamPrint(OMC_LOG_JAC, 0, "Using Jacobian method: %s.", JACOBIAN_METHOD_NAME[jacobianMethod]);
1112
1113 1 return jacobianMethod;
1114 }
1115
1116 /**
1117 * @brief Set Jacobian method from user flag and available Jacobian.
1118 *
1119 * @param threadData Used for error handling.
1120 * @param availability Is the Jacobian available, only the sparsity pattern available or nothing available.
1121 * @return JACOBIAN_METHOD Returns jacobian method that is availble.
1122 */
1123 ✗ JACOBIAN_METHOD setJacobianMethod(threadData_t* threadData, JACOBIAN_AVAILABILITY availability)
1124 {
1125 ✗ return checkJacobianMethod(threadData, availability, getRequestedJacobianMethod(threadData));
1126 }
1127
1128 /**
1129 * @brief Select and initialize the symbolic ODE Jacobian matching the requested method.
1130 *
1131 * This is the single place where the mapping
1132 *
1133 * COLOREDSYMJAC / other -> forward Jacobian A (column evaluation, CSC)
1134 * COLOREDSYMJACADJ -> adjoint Jacobian ADJ (row evaluation, CSR)
1135 * BICOLOREDSYMJAC -> forward Jacobian A (bidirectional, forward + linked adjoint)
1136 *
1137 * is implemented, so that DASSL, IDA and GBODE cannot diverge. Methods that are not
1138 * backed by the generated code fall back to a method that is, with a warning.
1139 *
1140 * The forward Jacobian A is only initialized if it is actually going to be used. The
1141 * adjoint Jacobian is self contained, so for COLOREDSYMJACADJ the sparsity pattern,
1142 * coloring and evaluation DAG of A are not built at all. Solvers that
1143 * need A regardless of the selected evaluation direction, because they use its
1144 * evaluation DAG or its column function, pass `requireForwardJacobian = TRUE` (GBODE).
1145 *
1146 * The selected Jacobian is stored in `data->simulationInfo->odeJacobian` and can later
1147 * be retrieved with getSymbolicOdeJacobian().
1148 *
1149 * @param data Runtime data struct.
1150 * @param threadData Used for error handling.
1151 * @param jacobianMethod In: requested method (JAC_UNKNOWN for default).
1152 * Out: method that will actually be used.
1153 * @param requireForwardJacobian Initialize the forward Jacobian A even if it is not the
1154 * Jacobian that gets evaluated.
1155 * @return JACOBIAN* Jacobian the solver has to evaluate with evalJacobian().
1156 */
1157 1 JACOBIAN* initSymbolicOdeJacobian(DATA* data, threadData_t* threadData, JACOBIAN_METHOD* jacobianMethod, modelica_boolean requireForwardJacobian)
1158 {
1159 1 JACOBIAN* forwardJacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
1160 1 JACOBIAN* adjointJacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_ADJ]);
1161 JACOBIAN* jacobian;
1162
1163
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 1 if (requireForwardJacobian || *jacobianMethod != COLOREDSYMJACADJ) {
1164 1 data->callback->initialAnalyticJacobianA(data, threadData, forwardJacobian);
1165 }
1166
1167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (*jacobianMethod == COLOREDSYMJACADJ) {
1168 /* If the model was compiled bidirectionally and A was initialized, the adjoint
1169 * Jacobian is already initialized and linked by initialAnalyticJacobianA().
1170 * So this check is true if the adjoint Jacobian was not already initialized but is requested. */
1171 ✗ if (forwardJacobian->adjointJacobian != adjointJacobian) {
1172 ✗ data->callback->initialAnalyticJacobianADJ(data, threadData, adjointJacobian);
1173 }
1174 ✗ if (adjointJacobian->availability == JACOBIAN_AVAILABLE) {
1175 jacobian = adjointJacobian;
1176 } else {
1177 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "No adjoint symbolic Jacobian was generated "
1178 "(compile with --generateDynamicJacobian=symbolicAdjoint or =bidirectional). "
1179 "Switching to the forward symbolic Jacobian.");
1180 ✗ *jacobianMethod = JAC_UNKNOWN;
1181 /* The fallback needs A, which may have been skipped above. */
1182 ✗ if (forwardJacobian->availability == JACOBIAN_UNKNOWN) {
1183 ✗ data->callback->initialAnalyticJacobianA(data, threadData, forwardJacobian);
1184 }
1185 jacobian = forwardJacobian;
1186 }
1187 } else {
1188 jacobian = forwardJacobian;
1189
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (*jacobianMethod == BICOLOREDSYMJAC
1190 ✗ && !(forwardJacobian->adjointJacobian != NULL && forwardJacobian->availability == JACOBIAN_AVAILABLE)) {
1191 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "No bidirectional symbolic Jacobian was generated "
1192 "(compile with --generateDynamicJacobian=bidirectional). "
1193 "Switching to the forward symbolic Jacobian.");
1194 ✗ *jacobianMethod = JAC_UNKNOWN;
1195 }
1196 }
1197 /* Runtime switch for the bidirectional evaluation path in evalJacobian() */
1198 1 forwardJacobian->isBidirectional = (*jacobianMethod == BICOLOREDSYMJAC);
1199
1200
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (jacobian->sparsePattern != NULL) {
1201 /* KLU and the sparse pattern printers require ascending secondary indices. */
1202 1 sortSparseColumns(jacobian->sparsePattern, (unsigned int) jacobian->sizeCols);
1203 /* Build the column oriented view (and the CSR->CSC value mapping) once up front,
1204 * so that evalJacobian() emits sparse values in CSC order for every method. */
1205 1 getJacobianCscPattern(jacobian);
1206 }
1207
1208 // Check that the requested Jacobian method can be used and log it.
1209 1 *jacobianMethod = checkJacobianMethod(threadData, jacobian->availability, *jacobianMethod);
1210
1211
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (jacobian->availability == JACOBIAN_AVAILABLE || jacobian->availability == JACOBIAN_ONLY_SPARSITY) {
1212 1 infoStreamPrint(OMC_LOG_SIMULATION, 1, "Initialized Jacobian:");
1213 1 infoStreamPrint(OMC_LOG_SIMULATION, 0, "columns: %zu rows: %zu", jacobian->sizeCols, jacobian->sizeRows);
1214 1 infoStreamPrint(OMC_LOG_SIMULATION, 0, "NNZ: %u colors: %u", jacobian->sparsePattern->nnz, jacobian->sparsePattern->maxColors);
1215 1 messageClose(OMC_LOG_SIMULATION);
1216 }
1217
1218 // Store the selected Jacobian in the simulation info for later retrieval.
1219 1 data->simulationInfo->odeJacobian = jacobian;
1220 1 return jacobian;
1221 }
1222
1223 /**
1224 * @brief Get the symbolic ODE Jacobian selected by the integrator.
1225 *
1226 * Falls back to the forward Jacobian A if no selection was made yet.
1227 *
1228 * @param data Runtime data struct.
1229 * @return JACOBIAN* Selected ODE Jacobian.
1230 */
1231 1 JACOBIAN* getSymbolicOdeJacobian(DATA* data)
1232 {
1233
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (data->simulationInfo->odeJacobian != NULL) {
1234 return data->simulationInfo->odeJacobian;
1235 }
1236 ✗ return &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
1237 }
1238
1239 /**
1240 * @brief Free the Jacobians initialized by initSymbolicOdeJacobian().
1241 *
1242 * @param data Runtime data struct.
1243 */
1244 1 void freeSymbolicOdeJacobian(DATA* data)
1245 {
1246 1 JACOBIAN* forwardJacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
1247 1 JACOBIAN* adjointJacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_ADJ]);
1248
1249
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (adjointJacobian->availability != JACOBIAN_UNKNOWN) {
1250 ✗ freeJacobian(adjointJacobian);
1251 }
1252
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if (forwardJacobian->availability != JACOBIAN_UNKNOWN) {
1253 1 freeJacobian(forwardJacobian);
1254 }
1255 1 data->simulationInfo->odeJacobian = NULL;
1256 1 }
1257
1258 ✗ void freeNonlinearPattern(NONLINEAR_PATTERN *nlp)
1259 {
1260 ✗ if (nlp != NULL) {
1261 ✗ free(nlp->indexVar);
1262 ✗ free(nlp->indexEqn);
1263 ✗ free(nlp->columns);
1264 ✗ free(nlp->rows);
1265 ✗ free(nlp);
1266 }
1267 ✗ }
1268
1269 ✗ unsigned int* getNonlinearPatternCol(NONLINEAR_PATTERN *nlp, int var_idx)
1270 {
1271 ✗ unsigned int idx_start = nlp->indexVar[var_idx];
1272 unsigned int idx_stop;
1273 ✗ if (var_idx == nlp->numberOfVars) {
1274 ✗ idx_stop = nlp->numberOfNonlinear;
1275 } else {
1276 ✗ idx_stop = nlp->indexVar[var_idx + 1];
1277 }
1278
1279 ✗ unsigned int* col = (unsigned int*) malloc((idx_stop - idx_start + 1)*sizeof(unsigned int));
1280
1281 int index = 0;
1282 ✗ for (int i = idx_start; i < idx_stop + 1; i++) {
1283 ✗ col[index] = nlp->columns[i];
1284 ✗ index++;
1285 }
1286
1287 //for(int j = 0; j < nlp->numberOfNonlinear; j++)
1288 // printf("nlp->columns[%d] = %d\n", j, nlp->columns[j]);
1289 //for(int j = 0; j < nlp->numberOfVars+1; j++)
1290 // printf("nlp->indexVar[%d] = %d\n", j, nlp->indexVar[j]);
1291
1292 ✗ return col;
1293 }
1294
1295 ✗ unsigned int* getNonlinearPatternRow(NONLINEAR_PATTERN *nlp, int eqn_idx)
1296 {
1297 ✗ unsigned int idx_start = nlp->indexEqn[eqn_idx];
1298 unsigned int idx_stop;
1299 ✗ if (eqn_idx == nlp->numberOfEqns) {
1300 ✗ idx_stop = nlp->numberOfNonlinear;
1301 } else {
1302 ✗ idx_stop = nlp->indexEqn[eqn_idx + 1];
1303 }
1304 //printf(" eqn_idx = %d\n", eqn_idx);
1305 //printf(" idx_start = %d\n", idx_start);
1306 //printf(" idx_stop = %d\n", idx_stop);
1307 ✗ unsigned int* row = (unsigned int*) malloc((idx_stop - idx_start + 1)*sizeof(unsigned int));
1308
1309 int index = 0;
1310 ✗ for (int i = idx_start; i < idx_stop + 1; i++) {
1311 ✗ row[index] = nlp->rows[i];
1312 //printf(" row[index] = row[%d] = %d\n", index, row[index]);
1313 ✗ index++;
1314 }
1315
1316 //for(int j = 0; j < nlp->numberOfNonlinear; j++)
1317 // printf("nlp->rows[%d] = %d\n", j, nlp->rows[j]);
1318 //for(int j = 0; j < nlp->numberOfEqns; j++)
1319 // printf("nlp->indexEqn[%d] = %d\n", j, nlp->indexEqn[j]);
1320
1321 ✗ return row;
1322 }
1323
1324