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 / 469
Functions: 0.0% 0 / 0 / 24
Branches: 0.0% 0 / 0 / 200
OMCompiler/SimulationRuntime/c/simulation/solver/jacobian_analysis.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 #include "jacobian_analysis.h"
29
30 // LAPACK dense SVD routine
31 extern void dgesvd_(char *jobu, char *jobvt, int *m, int *n,
32 modelica_real *a, int *lda, modelica_real *s,
33 modelica_real *u, int *ldu, modelica_real *vt, int *ldvt,
34 modelica_real *work, int *lwork, int *info);
35
36 // cmp for sorting singular vectors by magnitude
37 ✗ static int cmp_fabs_desc(const void *a, const void *b)
38 {
39 ✗ modelica_real abs_a = fabs(((SVD_Component*)a)->value);
40 ✗ modelica_real abs_b = fabs(((SVD_Component*)b)->value);
41 ✗ if (abs_a < abs_b) return 1;
42 ✗ if (abs_a > abs_b) return -1;
43 return 0;
44 }
45
46 /**
47 * @brief Create and initialize SVD data structure for a given nonlinear system.
48 *
49 * Builds a dense matrix from a sparse pattern (if provided), applies scaling,
50 * and allocates buffers for the SVD results (singular values and vectors).
51 *
52 * @param data Pointer to the global simulation DATA structure.
53 * @param nls_data Pointer to the nonlinear system data.
54 * @param values Non-zero values of the sparse Jacobian (CSC format).
55 * @param x_scale Optional scaling factors for variables (NULL if not used).
56 * @param f_scale Optional scaling factors for residuals (NULL if not used).
57 *
58 * @return Pointer to an allocated SVD_DATA structure, or NULL on allocation failure.
59 */
60 ✗ static SVD_DATA *svd_dense_create(DATA *data, NONLINEAR_SYSTEM_DATA *nls_data, modelica_real *values, modelica_boolean scaled, SolverCaller caller)
61 {
62 ✗ SVD_DATA *svd_data = calloc(1, sizeof(SVD_DATA));
63 ✗ if (!svd_data) return NULL;
64
65 ✗ int rows = nls_data->size;
66 int cols = nls_data->size;
67 ✗ SPARSE_PATTERN *sparse_pattern = nls_data->sparsePattern;
68 unsigned int *lead, *index;
69 unsigned int row, column, nz;
70
71 ✗ svd_data->data = data;
72 ✗ svd_data->nls_data = nls_data;
73 ✗ svd_data->rows = rows;
74 ✗ svd_data->cols = cols;
75 ✗ svd_data->sparse_pattern = sparse_pattern;
76 ✗ svd_data->sp_values = values;
77 ✗ svd_data->min_rows_cols = rows < cols ? rows : cols;
78 ✗ svd_data->scaled = scaled;
79 ✗ svd_data->caller = caller;
80
81 ✗ svd_data->A_dense = calloc(rows * cols, sizeof(modelica_real));
82
83 // for now, create dense matrix from sparse CSC
84 ✗ if (sparse_pattern)
85 {
86 ✗ lead = sparse_pattern->leadindex;
87 ✗ index = sparse_pattern->index;
88
89 ✗ for (column = 0; column < cols; column++)
90 {
91 ✗ for (nz = lead[column]; nz < lead[column + 1]; nz++)
92 {
93 ✗ row = index[nz];
94 ✗ svd_data->A_dense[column * rows + row] = values[nz];
95 }
96 }
97 }
98 else
99 {
100 ✗ memcpy(svd_data->A_dense, values, rows * cols * sizeof(modelica_real));
101 }
102
103 // allocate SVD result buffers
104 ✗ svd_data->S = malloc(svd_data->min_rows_cols * sizeof(modelica_real));
105 ✗ svd_data->U = malloc(rows * rows * sizeof(modelica_real));
106 ✗ svd_data->VT = malloc(cols * cols * sizeof(modelica_real));
107
108 ✗ return svd_data;
109 }
110
111 ✗ static void svd_dense_free(SVD_DATA* svd_data)
112 {
113 ✗ if (!svd_data) return;
114 ✗ free(svd_data->A_dense);
115 ✗ free(svd_data->S);
116 ✗ free(svd_data->U);
117 ✗ free(svd_data->VT);
118 ✗ free(svd_data);
119 }
120
121 /**
122 * @brief Computes the singular value decomposition (SVD) of a matrix using LAPACK's DGESVD.
123 *
124 * This function performs an SVD on the matrix stored in svd_data->A_dense,
125 * producing singular values in svd_data->S and singular vectors in svd_data->U and svd_data->VT.
126 *
127 * @param svd_data Pointer to the structure containing SVD results and statistics.
128 * @return LAPACK info code
129 */
130 ✗ static int svd_dense_compute_lapack(SVD_DATA* svd_data)
131 {
132 ✗ int rows = svd_data->rows;
133 ✗ int cols = svd_data->cols;
134 ✗ int lda = rows;
135 ✗ int ldu = rows;
136 ✗ int ldvt = cols;
137 int info;
138 ✗ char jobu = 'A';
139 ✗ char jobvt = 'A';
140 modelica_real *work;
141
142 // workspace query
143 ✗ int lwork = -1;
144 modelica_real wkopt;
145 ✗ dgesvd_(&jobu, &jobvt, &rows, &cols,
146 svd_data->A_dense, &lda,
147 svd_data->S, svd_data->U, &ldu, svd_data->VT, &ldvt,
148 &wkopt, &lwork, &info);
149
150 ✗ if (info != 0) return info;
151
152 ✗ lwork = (int)wkopt;
153 ✗ work = malloc(sizeof(modelica_real) * lwork);
154
155 // actual SVD, O(n^3)
156 ✗ dgesvd_(&jobu, &jobvt, &rows, &cols,
157 svd_data->A_dense, &lda,
158 svd_data->S, svd_data->U, &ldu, svd_data->VT, &ldvt,
159 work, &lwork, &info);
160
161 // U = U - column major
162 // S = diag(S)
163 // VT = V^T - column major
164
165 ✗ return info;
166 }
167
168 /**
169 * @brief Calculates statistics from the computed singular values.
170 *
171 * Updates condition number, estimated rank, and identifies the index of the first singular value
172 * below 1% of the maximum singular value.
173 *
174 * @param svd_data Pointer to the structure containing SVD results and statistics.
175 */
176 ✗ static void svd_dense_calculate_statistics(SVD_DATA* svd_data)
177 {
178 int dim, low, mid, high, first_below;
179 modelica_real sigma_max, threshold;
180
181 // condition statistics
182 ✗ svd_data->sigma_max = svd_data->S[0];
183 ✗ svd_data->sigma_min = svd_data->S[svd_data->min_rows_cols - 1];
184 ✗ svd_data->cond = svd_data->sigma_min > 0.0 ? svd_data->sigma_max / svd_data->sigma_min : INFINITY;
185
186 // rank estimation
187 ✗ svd_data->estimated_rank = 0;
188 ✗ svd_data->rank_est_tol = _svd_max2(svd_data->rows, svd_data->cols) * DBL_EPSILON * svd_data->sigma_max;
189 ✗ for (dim = 0; dim < svd_data->min_rows_cols; dim++)
190 {
191 ✗ if (svd_data->S[dim] > svd_data->rank_est_tol)
192 {
193 ✗ svd_data->estimated_rank++;
194 }
195 }
196
197 // binary search to find first singular value < threshold, O(log(n))
198 ✗ sigma_max = svd_data->S[0];
199 ✗ threshold = 0.01 * sigma_max;
200
201 low = 0;
202 ✗ high = svd_data->min_rows_cols - 1;
203 first_below = svd_data->min_rows_cols;
204
205 ✗ while (low <= high)
206 {
207 ✗ mid = (low + high) / 2;
208 ✗ if (svd_data->S[mid] < threshold)
209 {
210 first_below = mid;
211 ✗ high = mid - 1;
212 }
213 else
214 {
215 ✗ low = mid + 1;
216 }
217 }
218 ✗ svd_data->least_one_percent = first_below;
219 ✗ }
220
221 ✗ static void svd_general_matrix_print_info(DATA *data, NONLINEAR_SYSTEM_DATA *nls_data)
222 {
223 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Matrix Info");
224 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "NLS eq index = " OMC_INT_FORMAT, nls_data->equationIndex);
225 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Columns = " OMC_INT_FORMAT, nls_data->size);
226 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Rows = " OMC_INT_FORMAT, nls_data->size);
227 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "NNZ = %u", nls_data->sparsePattern->nnz);
228 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Curr Time = %-11.5e", data->localData[0]->timeValue);
229 ✗ messageClose(OMC_LOG_NLS_SVD);
230 ✗ }
231
232 ✗ static void svd_general_matrix_print_cond(modelica_real cond)
233 {
234 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Matrix condition");
235 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Cond(M) = %.8e", cond);
236 ✗ if (cond > 1e12)
237 {
238 ✗ warningStreamPrint(OMC_LOG_NLS_SVD, 0, "Matrix is very ill-conditioned: 1e12 < Cond(M) = %.8e", cond);
239 }
240 ✗ else if (cond > 1e8)
241 {
242 ✗ warningStreamPrint(OMC_LOG_NLS_SVD, 0, "Matrix is fairly ill-conditioned: 1e8 < Cond(M) = %.8e < 1e12", cond);
243 }
244 ✗ else if (cond > 1e4)
245 {
246 ✗ warningStreamPrint(OMC_LOG_NLS_SVD, 0, "Matrix is moderately ill-conditioned: 1e4 < Cond(M) = %.8e < 1e8", cond);
247 }
248 else
249 {
250 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Matrix is well conditioned: Cond(M) = %.8e < 1e4", cond);
251 }
252 ✗ messageClose(OMC_LOG_NLS_SVD);
253 ✗ }
254
255 /**
256 * @brief Logs computed SVD statistics.
257 *
258 * Outputs singular value statistics such as condition number, estimated rank, and others.
259 *
260 * @param svd_data Pointer to the structure containing SVD results and statistics.
261 * @param scaled If true: statistics are marked as scaled.
262 */
263 ✗ static void svd_dense_dump_statistics(const SVD_DATA *svd_data)
264 {
265 int i, u, v, var_idx, eq_idx, start, end, count;
266 modelica_real val;
267 modelica_integer size_of_torns;
268 SVD_Component *entries = NULL;
269 ✗ NONLINEAR_SYSTEM_DATA *nls_data = svd_data->nls_data;
270 NONLINEAR_SOLVER solver = nls_data->nlsMethod;
271
272 ✗ if (!svd_data || !svd_data->S) {
273 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "No SVD data available.");
274 ✗ return;
275 }
276
277 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "%s: dense SVD analysis (scaled = %s, Caller: %s).",
278 ✗ SolverCaller_callerString(svd_data->caller), svd_data->scaled ? "true" : "false", SolverCaller_toString(svd_data->caller));
279 ✗ svd_general_matrix_print_info(svd_data->data, nls_data);
280 ✗ svd_general_matrix_print_cond(svd_data->cond);
281
282 // singular values
283 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Singular values");
284 ✗ for (i = 0; i < svd_data->min_rows_cols; i++)
285 {
286 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "sigma_%-3d = %.8e", i + 1, svd_data->S[i]);
287 }
288 ✗ messageClose(OMC_LOG_NLS_SVD);
289
290 // rank estimation
291 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Rank estimation");
292 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "estimated = %d", svd_data->estimated_rank);
293 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "actual = %d", svd_data->min_rows_cols);
294 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "estimation tolerance = %.8e (= sigma_max * max(rows, cols) * DBL_EPSILON)", svd_data->rank_est_tol);
295 ✗ if (svd_data->estimated_rank < svd_data->min_rows_cols)
296 {
297 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Matrix may be rank-deficient.");
298 }
299 else
300 {
301 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Matrix should have full rank.");
302 }
303 ✗ messageClose(OMC_LOG_NLS_SVD);
304
305 // print right singular vectors for singular values below 1% of sigma_max
306 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Smallest right singular vectors (variable space)");
307
308 ✗ entries = (SVD_Component*)malloc(svd_data->rows * sizeof(SVD_Component));
309
310 ✗ if (svd_data->least_one_percent == svd_data->min_rows_cols)
311 {
312 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "No singular values below %.8e (1%% of max)", 0.01 * svd_data->sigma_max);
313 }
314 else
315 {
316 ✗ start = svd_data->min_rows_cols - 1;
317 end = svd_data->least_one_percent;
318 ✗ count = start - end + 1;
319
320 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0,
321 ✗ "Found %d singular %s below %.8e (1%% of sigma_max)", count, count > 1 ? "values" : "value", 0.01 * svd_data->sigma_max);
322
323 ✗ for (v = start; v >= end; v--)
324 {
325 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "V[:,%d] (singular value %.8e)", v + 1, svd_data->S[v]);
326
327 ✗ for (i = 0; i < svd_data->cols; i++)
328 {
329 // V[i][v] = VT[v][i]
330 ✗ entries[i].index = i;
331 ✗ entries[i].value = svd_data->VT[v + i * svd_data->rows]; // VT = V^T when reading column-wise
332 }
333
334 /* Same gauge as the sparse dump. */
335 {
336 int lead = 0;
337 ✗ for (i = 1; i < svd_data->cols; i++)
338 {
339 ✗ if (fabs(entries[i].value) > fabs(entries[lead].value)) lead = i;
340 }
341 ✗ if (entries[lead].value < 0.0)
342 {
343 ✗ for (i = 0; i < svd_data->cols; i++) entries[i].value = -entries[i].value;
344 }
345 }
346
347 // sort by abs value descending O(n * log(n))
348 ✗ qsort(entries, svd_data->cols, sizeof(SVD_Component), cmp_fabs_desc);
349
350 ✗ for (i = 0; i < svd_data->cols; i++)
351 {
352 ✗ var_idx = entries[i].index;
353 ✗ val = entries[i].value;
354 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "V[%d][%d] = %+.8e for NLS Var: %d with Name: %s", var_idx + 1, v + 1, val, var_idx + 1,
355 ✗ modelInfoGetEquation(&svd_data->data->modelData->modelDataXml, nls_data->equationIndex).vars[var_idx]);
356 }
357
358 ✗ messageClose(OMC_LOG_NLS_SVD);
359 }
360 }
361 ✗ messageClose(OMC_LOG_NLS_SVD);
362
363 // print left singular vectors for singular values below 1% of sigma_max
364 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Smallest left singular vectors (function space)");
365
366 ✗ if (svd_data->least_one_percent == svd_data->min_rows_cols)
367 {
368 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "No singular values below %.8e (1%% of max)", 0.01 * svd_data->sigma_max);
369 }
370 else
371 {
372 ✗ start = svd_data->min_rows_cols - 1;
373 end = svd_data->least_one_percent;
374 ✗ count = start - end + 1;
375
376 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0,
377 ✗ "Found %d singular %s below %.8e (1%% of sigma_max)", count, count > 1 ? "values" : "value", 0.01 * svd_data->sigma_max);
378
379 ✗ for (u = start; u >= end; u--)
380 {
381 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "U[:,%d] (singular value %.8e)", u + 1, svd_data->S[u]);
382
383 ✗ for (i = 0; i < svd_data->rows; i++)
384 {
385 ✗ entries[i].index = i;
386 ✗ entries[i].value = svd_data->U[i + svd_data->min_rows_cols * u];
387 }
388
389 // sort by abs value descending O(n * log(n))
390 ✗ qsort(entries, svd_data->rows, sizeof(SVD_Component), cmp_fabs_desc);
391
392 ✗ size_of_torns = nls_data->torn_plus_residual_size - nls_data->size;
393 ✗ for (i = 0; i < svd_data->rows; i++)
394 {
395 ✗ eq_idx = entries[i].index;
396 ✗ val = entries[i].value;
397
398 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "U[%d][%d] = %+.8e for NLS Eqn: %d with transformational debugger Idx: %d", eq_idx + 1, u + 1, val, eq_idx + 1,
399 ✗ nls_data->eqn_simcode_indices[size_of_torns + entries[i].index]);
400 }
401
402 ✗ messageClose(OMC_LOG_NLS_SVD);
403 }
404 }
405 ✗ free(entries);
406
407 ✗ messageClose(OMC_LOG_NLS_SVD);
408 ✗ messageClose(OMC_LOG_NLS_SVD);
409 }
410
411 ✗ static int svd_dense_main(DATA *data, NONLINEAR_SYSTEM_DATA *nls_data, modelica_real *values, modelica_boolean scaled, SolverCaller caller)
412 {
413 int ret = 0;
414 ✗ SVD_DATA *svd_data = svd_dense_create(data, nls_data, values, scaled, caller);
415 ✗ ret = svd_dense_compute_lapack(svd_data);
416 ✗ if (ret != 0) return ret;
417 ✗ svd_dense_calculate_statistics(svd_data);
418 ✗ svd_dense_dump_statistics(svd_data);
419 ✗ svd_dense_free(svd_data);
420 ✗ return ret;
421 }
422
423 #ifdef OMC_HAVE_PRIMME
424
425 #include <primme_svds.h>
426
427 /**
428 * @brief Function pointer for the matrix-vector product (transpose and standard).
429 */
430 typedef void (*primme_mvp_fn_t)(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize,
431 int *transpose, primme_svds_params *primme_svds, int *ierr);
432
433
434 typedef void (*primme_prec_fn_t)(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize,
435 int *mode, primme_svds_params *primme_svds, int *ierr);
436 /**
437 * @brief Public struct containing the results of the SVD computation.
438 * @attention The user is responsible for freeing this struct using `svd_sparse_free`.
439 */
440 typedef struct primme_result_t
441 {
442 int rows; /* Number of rows of the original matrix. */
443 int cols; /* Number of columns of the original matrix. */
444 int target_size; /* Number of singular values found. */
445 double *svals; /* Array with the computed singular values. */
446 double *svecs; /* Array with the computed singular vectors. */
447 double *rnorms; /* Array with the computed residual norms. */
448 } primme_result_t;
449
450 /**
451 * @brief Encapsulates SVD problem data and results.
452 * Contains matrix dimensions, config, and result pointers for SVD computations.
453 */
454 typedef struct primme_handle_t
455 {
456 primme_result_t result; /* The results of the SVD. */
457 primme_svds_params primme_svds; /* PRIMME's internal state. */
458 } primme_handle_t;
459
460 /**
461 * @brief Context for callback Matrix-vector products, stored in primme_svds->matrix field.
462 */
463 typedef struct primme_callback_ctx_t
464 {
465 DATA *data;
466 NONLINEAR_SYSTEM_DATA *nls_data;
467 modelica_real *values;
468 modelica_boolean scaled;
469 SolverCaller caller;
470 int svd_count;
471
472 // must be freed in svd_sparse_free_ctx
473 double *inv_diag_AtA;
474 double *inv_diag_AAt;
475 } primme_callback_ctx_t;
476
477 static void svd_sparse_free_ctx(primme_callback_ctx_t *ctx)
478 {
479 ✗ free(ctx->inv_diag_AAt);
480 ✗ free(ctx->inv_diag_AtA);
481 }
482
483 /**
484 * @brief Computes both Jacobi scaling vectors for preconditioning:
485 * inv_diag_AtA = 1 / diag(A^T * A)
486 * inv_diag_AAt = 1 / diag(A * A^T)
487 */
488 ✗ static void compute_jacobi_diags(const NONLINEAR_SYSTEM_DATA *nls_data,
489 const double *values,
490 primme_callback_ctx_t *ctx)
491 {
492 ✗ const SPARSE_PATTERN *sp = nls_data->sparsePattern;
493 ✗ const modelica_integer size = nls_data->size;
494 double sigma = 1e-8;
495
496 ✗ if(omc_flag[FLAG_SVD_SPARSE_SIGMA])
497 {
498 ✗ sigma = fabs(atof(omc_flagValue[FLAG_SVD_SPARSE_SIGMA]));
499 }
500
501 ✗ const double reg = sigma * sigma;
502
503 ✗ for (modelica_integer j = 0; j < size; j++)
504 {
505 ✗ for (modelica_integer nz = sp->leadindex[j]; nz < sp->leadindex[j + 1]; nz++)
506 {
507 ✗ modelica_integer i = sp->index[nz];
508 ✗ double val = values[nz];
509
510 ✗ ctx->inv_diag_AtA[j] += val * val;
511 ✗ ctx->inv_diag_AAt[i] += val * val;
512 }
513 }
514
515 ✗ for (modelica_integer j = 0; j < size; j++)
516 {
517 ✗ double a_AtA = ctx->inv_diag_AtA[j] + reg;
518 ✗ double a_AAt = ctx->inv_diag_AAt[j] + reg;
519
520 ✗ ctx->inv_diag_AtA[j] = 1.0 / a_AtA;
521 ✗ ctx->inv_diag_AAt[j] = 1.0 / a_AAt;
522 }
523 ✗ }
524
525 /**
526 * @brief Allocates and initializes an SVD computation handle.
527 * @param rows [in] The number of rows of the matrix.
528 * @param cols [in] The number of columns of the matrix.
529 * @param target_size [in] The number of singular values to compute.
530 * @param linear_operator [in] The callback function for the matrix-vector product.
531 * @return A handle to the internal SVD state.
532 */
533 ✗ static primme_handle_t* svd_sparse_allocate(primme_callback_ctx_t *ctx, primme_mvp_fn_t linear_operator, primme_prec_fn_t precond)
534 {
535 ✗ primme_handle_t *handle = (primme_handle_t*)malloc(sizeof(primme_handle_t));
536 ✗ primme_svds_initialize(&handle->primme_svds);
537 ✗ handle->primme_svds.m = ctx->nls_data->size;
538 ✗ handle->primme_svds.n = ctx->nls_data->size;
539 ✗ handle->primme_svds.numSvals = ctx->svd_count < ctx->nls_data->size ? ctx->svd_count : ctx->nls_data->size;
540 ✗ handle->primme_svds.matrixMatvec = linear_operator;
541 ✗ handle->primme_svds.matrix = ctx;
542
543 ✗ if (ctx->inv_diag_AAt == NULL && ctx->inv_diag_AtA == NULL)
544 {
545 ✗ ctx->inv_diag_AAt = (double *)malloc(handle->primme_svds.n * sizeof(double));
546 ✗ ctx->inv_diag_AtA = (double *)malloc(handle->primme_svds.n * sizeof(double));
547 ✗ compute_jacobi_diags(ctx->nls_data, ctx->values, ctx);
548 }
549
550 ✗ handle->primme_svds.applyPreconditioner = precond;
551 ✗ handle->primme_svds.preconditioner = ctx;
552
553 ✗ handle->result.rows = handle->primme_svds.m;
554 ✗ handle->result.cols = handle->primme_svds.n;;
555 ✗ handle->result.target_size = handle->primme_svds.numSvals;
556 ✗ handle->result.svals = (double *) malloc(handle->primme_svds.numSvals * sizeof(double));
557 ✗ handle->result.svecs = (double *) malloc((handle->primme_svds.n + handle->primme_svds.m) * handle->primme_svds.numSvals * sizeof(double));
558 ✗ handle->result.rnorms = (double *) malloc(handle->primme_svds.numSvals * sizeof(double));
559
560 ✗ return handle;
561 }
562
563 /**
564 * @brief Deallocates all memory associated with the SVD computation handle.
565 *
566 * @param handle [in] The handle returned by `svd_sparse_allocate`.
567 */
568 ✗ static void svd_sparse_free(primme_handle_t* handle)
569 {
570 ✗ primme_svds_free(&handle->primme_svds);
571 ✗ free(handle->result.svals);
572 ✗ free(handle->result.svecs);
573 ✗ free(handle->result.rnorms);
574 ✗ free(handle);
575 ✗ }
576
577 /**
578 * @brief Performs the singular value decomposition.
579 *
580 * @param handle [in] The handle returned by `svd_sparse_allocate`.
581 * @param target [in] Specifies whether to find the `TOP` or `LEAST` singular values.
582 * @return A reference pointer to a `primme_result_t` struct on success, or `NULL` on error (owned by handle).
583 */
584 ✗ static primme_result_t* svd_sparse_compute(primme_handle_t* handle, primme_svds_target target)
585 {
586 primme_callback_ctx_t * ctx = (primme_callback_ctx_t *)(handle->primme_svds.matrix);
587
588 double eps = 1e-8;
589
590 ✗ if (omc_flag[FLAG_SVD_SPARSE_TOL])
591 {
592 ✗ eps = fabs(atof(omc_flagValue[FLAG_SVD_SPARSE_TOL]));
593 }
594
595 /* ||r|| <= eps * ||matrix|| */
596 ✗ handle->primme_svds.eps = eps;
597 ✗ handle->primme_svds.target = target;
598
599 // we only need the largest for the condition
600 ✗ handle->primme_svds.numSvals = (target == primme_svds_largest) ? 1 : handle->primme_svds.numSvals;
601
602 /* Normal equations resolve nothing below sqrt(DBL_EPSILON)*||A||, however tight
603 eps is. Do not "fix" that with primme_svds_hybrid (reports sigma_max as the
604 smallest when its first stage cannot resolve sigma_min) or with
605 primme_svds_augmented targeted at zero (does not terminate). */
606 ✗ primme_svds_set_method(primme_svds_normalequations, PRIMME_DEFAULT_MIN_TIME,
607 PRIMME_DEFAULT_MIN_MATVECS, &handle->primme_svds);
608
609 ✗ if (omc_useStream[OMC_LOG_NLS_SVD_V])
610 {
611 // we write these to stdout, since we cant really redirect them
612 ✗ handle->primme_svds.printLevel = 2;
613 ✗ primme_svds_display_params(handle->primme_svds);
614 }
615 else
616 {
617 ✗ handle->primme_svds.printLevel = 0;
618 }
619
620 ✗ handle->primme_svds.precondition = (target == primme_svds_smallest) ? 1 : 0;
621
622
623 ✗ int ret = dprimme_svds(handle->result.svals, handle->result.svecs, handle->result.rnorms, &handle->primme_svds);
624
625 ✗ if (ret != 0)
626 {
627 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Error: primme_svds returned with nonzero exit status: %d\n", ret);
628 ✗ return NULL;
629 }
630
631 ✗ return &handle->result;
632 }
633
634 ✗ static void matrix_vector(const NONLINEAR_SYSTEM_DATA *nls_data, const double *values, const double *x, double *y)
635 {
636 modelica_integer row, column, nz;
637 ✗ const SPARSE_PATTERN *sparsity = nls_data->sparsePattern;
638 ✗ memset(y, 0, nls_data->size * sizeof(double));
639
640 ✗ for (column = 0; column < nls_data->size; column++)
641 {
642 ✗ for (nz = sparsity->leadindex[column]; nz < sparsity->leadindex[column + 1]; nz++)
643 {
644 ✗ row = sparsity->index[nz];
645 ✗ y[row] += values[nz] * x[column];
646 }
647 }
648 ✗ }
649
650 ✗ static void matrix_vector_transpose(const NONLINEAR_SYSTEM_DATA *nls_data, const double *values, const double *x, double *y)
651 {
652 modelica_integer row, column, nz;
653 ✗ const SPARSE_PATTERN *sparsity = nls_data->sparsePattern;
654 ✗ memset(y, 0, nls_data->size * sizeof(double));
655
656 ✗ for (column = 0; column < nls_data->size; column++)
657 {
658 ✗ for (nz = sparsity->leadindex[column]; nz < sparsity->leadindex[column + 1]; nz++)
659 {
660 ✗ row = sparsity->index[nz];
661 ✗ y[column] += values[nz] * x[row];
662 }
663 }
664 ✗ }
665
666 /**
667 * @brief Implements the matrix-vector products for the given matrix.
668 * It operates on blocks of vectors for improved performance and
669 * in general looks like this (depending on input):
670 * Y := A * X, for transpose = 0
671 * or Y := A^T * X, for transpose = 1
672 *
673 * @attention get column i of x: (double *)x + (*ldx) * i;
674 * @attention get column i of y: (double *)y + (*ldy) * i;
675 *
676 * @param x [in] Input dense matrix of vectors `X`.
677 * @param ldx [in] Leading dimension of the input matrix `X`.
678 * @param y [out] Output dense matrix of vectors `Y`.
679 * @param ldy [in] Leading dimension of the output matrix `Y`.
680 * @param blockSize [in] Number of vectors in the current block, number of columns of the X matrix.
681 * @param transpose [in] Flag indicating if the transpose is applied (0 for A*x, 1 for A^T*x).
682 * @param primme_svds [in] PRIMME configuration struct.
683 * @param err [out] Error status; must be set to 0 on success.
684 */
685 ✗ static void LinearOperator(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize,
686 int *transpose, primme_svds_params *primme_svds, int *err)
687 {
688 int i, j; /* vector index, from 0 to *blockSize-1 */
689 double *xvec; /* pointer to i-th input vector x */
690 double *yvec; /* pointer to i-th output vector y */
691
692 ✗ primme_callback_ctx_t *ctx = (void*) primme_svds->matrix;
693 ✗ NONLINEAR_SYSTEM_DATA *nls_data = ctx->nls_data;
694 ✗ double *values = ctx->values;
695
696 ✗ if (*transpose == 0)
697 {
698 /* Do y <- A * x */
699 ✗ for (i = 0; i < *blockSize; i++)
700 {
701 ✗ xvec = (double *)x + (*ldx) * i;
702 ✗ yvec = (double *)y + (*ldy) * i;
703 ✗ matrix_vector(nls_data, values, xvec, yvec);
704 }
705 }
706 else
707 {
708 /* Do y <- A^t * x */
709 ✗ for (i = 0; i < *blockSize; i++)
710 {
711 ✗ xvec = (double *)x + (*ldx) * i;
712 ✗ yvec = (double *)y + (*ldy) * i;
713 ✗ matrix_vector_transpose(nls_data, values, xvec, yvec);
714 }
715 }
716 ✗ *err = 0;
717 ✗ }
718
719 ✗ void GenericJacobiPreconditioner(void *x, PRIMME_INT *ldx, void *y, PRIMME_INT *ldy, int *blockSize,
720 int *mode, primme_svds_params *primme_svds, int *ierr)
721 {
722 int i, j; /* vector index, from 0 to *blockSize-1 */
723 double *xvec; /* pointer to i-th input vector x */
724 double *yvec; /* pointer to i-th output vector y */
725
726 ✗ primme_callback_ctx_t *ctx = (primme_callback_ctx_t*)primme_svds->matrix;
727 ✗ int size = ctx->nls_data->size;
728 ✗ const double *d_AtA = ctx->inv_diag_AtA;
729 ✗ const double *d_AAt = ctx->inv_diag_AAt;
730
731 ✗ int modeAtA = primme_svds_op_AtA;
732 ✗ int modeAAt = primme_svds_op_AAt;
733 int modeAug = primme_svds_op_augmented;
734 ✗ PRIMME_INT ldaux = 2 * size;
735 ✗ int notrans = 0;
736 ✗ int trans = 1;
737 double *aux;
738
739 ✗ if (*mode == modeAtA)
740 {
741 /* Preconditioner for A^t * A, diag(A^t * A + sigma_est * I)^{-1} */
742 ✗ for (i = 0; i < *blockSize; i++)
743 {
744 ✗ xvec = (double *)x + (*ldx) * i;
745 ✗ yvec = (double *)y + (*ldy) * i;
746 ✗ for (j = 0; j < size; j++)
747 {
748 ✗ yvec[j] = xvec[j] * d_AtA[j];
749 }
750 }
751 ✗ *ierr = 0;
752 }
753 ✗ else if (*mode == modeAAt)
754 {
755 /* Preconditioner for A * A^t, diag(A * A^t + sigma_est * I)^{-1} */
756 ✗ for (i = 0; i<*blockSize; i++)
757 {
758 ✗ xvec = (double *)x + (*ldx) * i;
759 ✗ yvec = (double *)y + (*ldy) * i;
760 ✗ for (j = 0; j < size; j++)
761 {
762 ✗ yvec[j] = xvec[j] * d_AAt[j];
763 }
764 }
765 ✗ *ierr = 0;
766 }
767 ✗ else if (*mode == modeAug)
768 {
769 /* Preconditioner for [0 A^t; A 0],
770 [diag(A^t * A + sigma_est * I) 0 ]^{-1} * [0 A^t]
771 [ 0 diag(A * A^t + sigma_est * I)] [A 0 ]
772 */
773
774 // [y0; y1] <- [0 A^t; A 0] * [x0; x1]
775 ✗ aux = (double*)malloc((*blockSize) * ldaux * sizeof(double));
776 ✗ primme_svds->matrixMatvec(x, ldx, &aux[size], &ldaux, blockSize, &notrans, primme_svds, ierr);
777
778 ✗ xvec = (double *)x + size;
779 ✗ primme_svds->matrixMatvec(xvec, ldx, aux, &ldaux, blockSize, &trans, primme_svds, ierr);
780
781 /* y0 <- preconditioner for A^t*A * y0 */
782 ✗ GenericJacobiPreconditioner(aux, &ldaux, y, ldy, blockSize, &modeAtA, primme_svds, ierr);
783
784 /* y1 <- preconditioner for A*A^t * y1 */
785 ✗ yvec = (double *)y + size;
786 ✗ GenericJacobiPreconditioner(&aux[size], &ldaux, yvec, ldy, blockSize, &modeAAt, primme_svds, ierr);
787 ✗ free(aux);
788 }
789 ✗ }
790
791 ✗ static void svd_sparse_print_singular_values(primme_callback_ctx_t *ctx, primme_handle_t *handle_top, primme_handle_t *handle_least,
792 primme_result_t *res_top, primme_result_t *res_least)
793 {
794 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Smallest Singular values");
795 ✗ for (int i = 0; i < handle_least->primme_svds.numSvals; i++)
796 {
797 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "sigma_%-3d = %.8e, rnorm_%-3d = %.8e", i + 1, res_least->svals[i], i + 1, res_least->rnorms[i]);
798 }
799 ✗ messageClose(OMC_LOG_NLS_SVD);
800 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "Largest Singular values");
801 ✗ for (int i = 0; i < handle_top->primme_svds.numSvals; i++)
802 {
803 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "sigma_%-3d = %.8e, rnorm_%-3d = %.8e", i + 1, res_top->svals[i], i + 1, res_top->rnorms[i]);
804 }
805 ✗ messageClose(OMC_LOG_NLS_SVD);
806 ✗ }
807
808 /* A triplet is defined up to a common sign, which the iteration picks by rounding.
809 * Gauge: largest entry of v positive, u follows it. */
810 static modelica_real svd_sparse_vector_sign(const modelica_real *v, int size)
811 {
812 int i, lead = 0;
813 ✗ for (i = 1; i < size; i++)
814 {
815 ✗ if (fabs(v[i]) > fabs(v[lead]))
816 {
817 lead = i;
818 }
819 }
820 ✗ return v[lead] < 0.0 ? -1.0 : 1.0;
821 }
822
823 ✗ static void svd_sparse_print_vectors(primme_callback_ctx_t *ctx, primme_handle_t *handle, primme_result_t *res, modelica_boolean smallest)
824 {
825 modelica_real sign;
826 int i, u, v, var_idx, eq_idx, sing_value_idx;
827 modelica_real val;
828 modelica_integer size_of_torns;
829 ✗ int size = res->rows;
830 ✗ SVD_Component *entries = (SVD_Component*)malloc(size * sizeof(SVD_Component));
831 ✗ NONLINEAR_SYSTEM_DATA *nls_data = ctx->nls_data;
832 NONLINEAR_SOLVER solver = nls_data->nlsMethod;
833
834 ✗ const char* target_string = smallest ? "Smallest" : "Largest";
835
836 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "%s right singular vectors (variable space)", target_string);
837 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Found %d singular vectors.", handle->primme_svds.numSvals);
838
839 ✗ for (v = 0; v < handle->primme_svds.numSvals; v++)
840 {
841 ✗ sing_value_idx = smallest ? size - v : v + 1;
842 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "V[:,%d] (singular value %.8e)", sing_value_idx, res->svals[v]);
843
844 ✗ sign = svd_sparse_vector_sign(&res->svecs[size * (res->target_size + v)], size);
845 ✗ for (i = 0; i < size; i++)
846 {
847 // V[i][v] = VT[v][i]
848 ✗ entries[i].index = i;
849 ✗ entries[i].value = sign * res->svecs[size * (res->target_size + v) + i];
850 }
851
852 // sort by abs value descending O(n * log(n))
853 ✗ qsort(entries, size, sizeof(SVD_Component), cmp_fabs_desc);
854
855 ✗ for (i = 0; i < size; i++)
856 {
857 ✗ var_idx = entries[i].index;
858 ✗ val = entries[i].value;
859 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "V[%d][%d] = %+.8e for NLS Var: %d with Name: %s", var_idx + 1, sing_value_idx, val, var_idx + 1,
860 ✗ modelInfoGetEquation(&ctx->data->modelData->modelDataXml, nls_data->equationIndex).vars[var_idx]);
861 }
862 ✗ messageClose(OMC_LOG_NLS_SVD);
863 }
864 ✗ messageClose(OMC_LOG_NLS_SVD);
865
866 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "%s left singular vectors (function space)", target_string);
867 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "Found %d singular vectors.", handle->primme_svds.numSvals);
868
869 ✗ for (u = 0; u < handle->primme_svds.numSvals; u++)
870 {
871 ✗ sing_value_idx = smallest ? size - u : u + 1;
872 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "U[:,%d] (singular value %.8e)", sing_value_idx, res->svals[u]);
873
874 /* Its right vector's flip, so the pair stays a triplet of A. */
875 ✗ sign = svd_sparse_vector_sign(&res->svecs[size * (res->target_size + u)], size);
876 ✗ for (i = 0; i < size; i++)
877 {
878 ✗ entries[i].index = i;
879 ✗ entries[i].value = sign * res->svecs[size * u + i];
880 }
881
882 // sort by abs value descending O(n * log(n))
883 ✗ qsort(entries, size, sizeof(SVD_Component), cmp_fabs_desc);
884
885 ✗ size_of_torns = nls_data->torn_plus_residual_size - nls_data->size;
886 ✗ for (i = 0; i < size; i++)
887 {
888 ✗ eq_idx = entries[i].index;
889 ✗ val = entries[i].value;
890
891 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 0, "U[%d][%d] = %+.8e for NLS Eqn: %d with transformational debugger Idx: %d", eq_idx + 1, sing_value_idx, val, eq_idx + 1,
892 ✗ nls_data->eqn_simcode_indices[size_of_torns + entries[i].index]);
893 }
894 ✗ messageClose(OMC_LOG_NLS_SVD);
895 }
896 ✗ messageClose(OMC_LOG_NLS_SVD);
897 ✗ free(entries);
898 ✗ }
899
900 ✗ static void svd_sparse_dump_statistics(primme_callback_ctx_t *ctx, primme_handle_t *handle_top, primme_handle_t *handle_least,
901 primme_result_t *res_top, primme_result_t *res_least)
902 {
903 ✗ if (res_top == NULL || res_least == NULL)
904 {
905 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Error: primme_result_t* is NULL, no statistics available.\n");
906 ✗ return;
907 }
908
909 ✗ infoStreamPrint(OMC_LOG_NLS_SVD, 1, "%s: sparse SVD analysis (scaled = %s, Caller: %s).",
910 ✗ SolverCaller_callerString(ctx->caller), ctx->scaled ? "true" : "false", SolverCaller_toString(ctx->caller));
911 ✗ svd_general_matrix_print_info(ctx->data, ctx->nls_data);
912
913 ✗ modelica_real sigma_max = res_top->svals[0];
914 ✗ modelica_real sigma_min = res_least->svals[0];
915 ✗ modelica_real cond = sigma_min != 0.0 ? sigma_max / sigma_min : INFINITY;
916
917 ✗ svd_general_matrix_print_cond(cond);
918 ✗ svd_sparse_print_singular_values(ctx, handle_top, handle_least, res_top, res_least);
919 ✗ svd_sparse_print_vectors(ctx, handle_least, res_least, TRUE);
920
921 ✗ messageClose(OMC_LOG_NLS_SVD);
922 }
923
924 ✗ static int svd_sparse_main(DATA *data, NONLINEAR_SYSTEM_DATA *nls_data, modelica_real *values, modelica_boolean scaled, SolverCaller caller, int svd_count) {
925 ✗ primme_callback_ctx_t ctx = { .data = data, .nls_data = nls_data, .values = values, .scaled = scaled, .caller = caller,
926 .svd_count = svd_count, .inv_diag_AtA = NULL, .inv_diag_AAt = NULL};
927
928 ✗ primme_handle_t *handle_top = svd_sparse_allocate(&ctx, LinearOperator, GenericJacobiPreconditioner);
929 ✗ primme_handle_t *handle_least = svd_sparse_allocate(&ctx, LinearOperator, GenericJacobiPreconditioner);
930
931 ✗ primme_result_t* res_top = svd_sparse_compute(handle_top, primme_svds_largest);
932 ✗ primme_result_t* res_least = svd_sparse_compute(handle_least, primme_svds_smallest);
933
934 ✗ svd_sparse_dump_statistics(&ctx, handle_top, handle_least, res_top, res_least);
935
936 ✗ svd_sparse_free(handle_top);
937 ✗ svd_sparse_free(handle_least);
938 svd_sparse_free_ctx(&ctx);
939
940 ✗ return 0;
941 }
942
943 #endif // OMC_HAVE_PRIMME
944
945 /**
946 * @brief Main routine to compute the SVD of the Jacobian matrix.
947 *
948 * Creates the SVD data structure, performs the SVD, calculates statistics,
949 * and outputs the results. Currently computes the unscaled SVD.
950 *
951 * @param data Pointer to simulation data.
952 * @param nls_data Pointer to the nonlinear system data.
953 * @param values Pointer to the matrix values to decompose.
954 * @param scaled Boolean if matrix is scaled (only for printout)
955 * @param caller Caller of the routine (only for printout)
956 * @return return code: 0 = success
957 */
958 ✗ int svd_compute(DATA *data, NONLINEAR_SYSTEM_DATA *nls_data, modelica_real *values, modelica_boolean scaled, SolverCaller caller)
959 {
960 ✗ const char* cflags = omc_flagValue[FLAG_SVD_SPARSE_COUNT];
961 ✗ int sparse_svd_count = (cflags ? atoi(cflags) : 0);
962
963 ✗ if (sparse_svd_count > 0)
964 {
965 #ifdef OMC_HAVE_PRIMME
966 ✗ return svd_sparse_main(data, nls_data, values, scaled, caller, sparse_svd_count);
967 #else
968 errorStreamPrint(OMC_LOG_STDOUT, 0, "Cannot call sparse SVD analysis, because OpenModelica was not build with PRIMME. "
969 "Set FLAG_SVD_SPARSE_COUNT=0 to perform dense SVD or build OpenModelica with "
970 "PRIMME via -DOM_OMC_ENABLE_PRIMME=ON.");
971 return -1;
972 #endif
973 }
974 ✗ else if (sparse_svd_count < 0)
975 {
976 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Invalid argument specified for SVD_SPARSE_COUNT (must be >= 0).");
977 ✗ return -1;
978 }
979 else
980 {
981 ✗ return svd_dense_main(data, nls_data, values, scaled, caller);
982 }
983 }
984
985 // ================================ Sums of absolute values of Jacobian Columns and Rows ================================ //
986
987 // quick struct + cmp operator, to sort the arrays of col / row sums and keep their respective index
988 typedef struct {
989 modelica_real value;
990 int index;
991 } IndexedValue;
992
993 ✗ static int compare_desc(const void *a, const void *b) {
994 ✗ modelica_real diff = ((IndexedValue*)b)->value - ((IndexedValue*)a)->value;
995 ✗ return (diff > 0) - (diff < 0); // returns 1 if b > a, -1 if a > b
996 }
997
998 /**
999 * @brief analyze absolute row and column sums of a sparse KINSOL Jacobian matrix
1000 *
1001 * computes the absolute row and column sums of a sparse Jacobian (CSC format)
1002 * and prints them sorted in descending order. This is useful for diagnosing
1003 * scaling issues, structural sparsity, or ill-conditioning in nonlinear systems.
1004 *
1005 * @param data
1006 * @param nlsData pointer to nonlinear system data
1007 * @param J sparse Jacobian matrix in CSC format
1008 * @param caller caller of the method (solver + where in the code it was called)
1009 * @param scaled boolean indicating if the passed Jacobian is scaled (only used for printout)
1010 */
1011 ✗ void nlsJacobianRowColSums(DATA *data, NONLINEAR_SYSTEM_DATA *nlsData, SUNMatrix J,
1012 SolverCaller caller, modelica_boolean scaled)
1013 {
1014 int i, row, col, nz, count;
1015 modelica_real value;
1016 ✗ const int size = (int)nlsData->size;
1017 ✗ const int size_of_torns = (int)nlsData->torn_plus_residual_size - size;
1018
1019 ✗ sunindextype nnz = SUNSparseMatrix_NNZ(J);
1020
1021 ✗ sunindextype *colPointers = SM_INDEXPTRS_S(J);
1022 ✗ sunindextype *rowIndices = SM_INDEXVALS_S(J);
1023 ✗ sunrealtype *values = SM_DATA_S(J);
1024
1025 ✗ modelica_real *rowSumsRaw = (modelica_real*)calloc(size, sizeof(modelica_real));
1026 ✗ modelica_real *colSumsRaw = (modelica_real*)calloc(size, sizeof(modelica_real));
1027 ✗ IndexedValue *rowSums = (IndexedValue*)malloc(size * sizeof(IndexedValue));
1028 ✗ IndexedValue *colSums = (IndexedValue*)malloc(size * sizeof(IndexedValue));
1029
1030 ✗ for (col = 0; col < size; col++)
1031 {
1032 ✗ for (nz = colPointers[col]; nz < colPointers[col + 1]; nz++)
1033 {
1034 ✗ row = rowIndices[nz];
1035 ✗ value = values[nz];
1036
1037 ✗ rowSumsRaw[row] += fabs(value);
1038 ✗ colSumsRaw[col] += fabs(value);
1039 }
1040 }
1041
1042 ✗ for (int i = 0; i < size; i++)
1043 {
1044 ✗ rowSums[i].value = rowSumsRaw[i];
1045 ✗ rowSums[i].index = i;
1046
1047 ✗ colSums[i].value = colSumsRaw[i];
1048 ✗ colSums[i].index = i;
1049 }
1050
1051 ✗ qsort(rowSums, size, sizeof(IndexedValue), compare_desc);
1052 ✗ qsort(colSums, size, sizeof(IndexedValue), compare_desc);
1053
1054 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "%s: Jacobian absolute row & col sum analysis (scaled = %s, Caller: %s).",
1055 SolverCaller_callerString(caller), scaled ? "true" : "false", SolverCaller_toString(caller));
1056
1057 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "Matrix Info");
1058 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "NLS eq index = " OMC_INT_FORMAT, nlsData->equationIndex);
1059 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "Columns = %d", size);
1060 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "Rows = %d", size);
1061 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "NNZ = %u", nlsData->sparsePattern->nnz);
1062 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "Curr Time = %-11.5e", data->localData[0]->timeValue);
1063 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1064
1065 ✗ int print_count = (size < 5) ? size : 5;
1066
1067 // top row sums
1068 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "Top %d Jacobian row abs sums (sorted by descending value):", print_count);
1069 ✗ for (i = 0; i < print_count; i++)
1070 {
1071 ✗ row = rowSums[i].index;
1072 ✗ modelica_integer eq_debug_idx = nlsData->eqn_simcode_indices[size_of_torns + row];
1073 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "fabs(Row[%d]) = %+.5e for NLS Eq ID (debugger): " OMC_INT_FORMAT, row + 1, rowSums[i].value, eq_debug_idx);
1074 }
1075 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1076
1077 // bottom row sums
1078 ✗ if (size > 5)
1079 {
1080 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "Bottom %d Jacobian row abs sums (sorted by descending value):", print_count);
1081 ✗ for (i = size - print_count; i < size; i++)
1082 {
1083 ✗ row = rowSums[i].index;
1084 ✗ modelica_integer eq_debug_idx = nlsData->eqn_simcode_indices[size_of_torns + row];
1085 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "fabs(Row[%d]) = %+.5e for NLS Eq ID (debugger): " OMC_INT_FORMAT, row + 1, rowSums[i].value, eq_debug_idx);
1086 }
1087 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1088 }
1089
1090
1091 // top column sums
1092 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "Top %d Jacobian column abs sums (sorted by descending value):", print_count);
1093 ✗ for (i = 0; i < print_count; i++)
1094 {
1095 ✗ col = colSums[i].index;
1096 ✗ const char *var_name = modelInfoGetEquation(&data->modelData->modelDataXml, nlsData->equationIndex).vars[col];
1097 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "fabs(Col[%d]) = %+.5e for Variable %d: %s", col + 1, colSums[i].value, col + 1, var_name);
1098 }
1099 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1100
1101 // bottom column sums
1102 ✗ if (size > 5)
1103 {
1104 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "Bottom %d Jacobian column abs sums (sorted by descending value):", print_count);
1105 ✗ for (i = size - print_count; i < size; i++)
1106 {
1107 ✗ col = colSums[i].index;
1108 ✗ const char *var_name = modelInfoGetEquation(&data->modelData->modelDataXml, nlsData->equationIndex).vars[col];
1109 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "fabs(Col[%d]) = %+.5e for Variable %d: %s", col + 1, colSums[i].value, col + 1, var_name);
1110 }
1111 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1112 }
1113
1114 // row sums
1115 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "All Jacobian row abs sums (sorted by descending value):");
1116 ✗ for (i = 0; i < size; i++)
1117 {
1118 ✗ row = rowSums[i].index;
1119 ✗ modelica_integer eq_debug_idx = nlsData->eqn_simcode_indices[size_of_torns + row];
1120 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "fabs(Row[%d]) = %+.5e for NLS Eq ID (debugger): " OMC_INT_FORMAT, row + 1, rowSums[i].value, eq_debug_idx);
1121 }
1122 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1123
1124 // column sums
1125 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 1, "All Jacobian column abs sums (sorted by descending value):");
1126 ✗ for (i = 0; i < size; i++)
1127 {
1128 ✗ col = colSums[i].index;
1129 ✗ const char *var_name = modelInfoGetEquation(&data->modelData->modelDataXml, nlsData->equationIndex).vars[col];
1130 ✗ infoStreamPrint(OMC_LOG_NLS_JAC_SUMS, 0, "fabs(Col[%d]) = %+.5e for Variable %d: %s", col + 1, colSums[i].value, col + 1, var_name);
1131 }
1132 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1133
1134 ✗ messageClose(OMC_LOG_NLS_JAC_SUMS);
1135
1136 ✗ free(rowSumsRaw);
1137 ✗ free(colSumsRaw);
1138 ✗ free(rowSums);
1139 ✗ free(colSums);
1140 ✗ }
1141