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 / 643
Functions: 0.0% 0 / 0 / 15
Branches: 0.0% 0 / 0 / 510
OMCompiler/SimulationRuntime/c/simulation/solver/newton_diagnostics.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 newton_diagnostics.c
29 * Containing all functions to run newton diagnostics on non-linear loops.
30 * Improve start values for non-linear loops.
31 */
32
33
34 /**
35 * @brief Start point for newton diagnostics.
36 *
37 * Calculation of:
38 * 1) alpha : alpha coefficients
39 * 2) Gamma_ijk : curvature factors
40 * 3) sigma_ij : solution sensitivities
41 *
42 * According to
43 * F. Casella and B. Bachman
44 * On the choice of initial guesses for the Newton-Raphson algorithm
45 * Applied Mathematics and Computation 398 (2021) 125991
46 *
47 * By Teus van der Stelt, Asimptote bv, the Netherlands
48 * Carried out on behalf of the Delft University of Technology, 2023
49 *
50 * @param data Pointer to all simulation data.
51 * @param threadData Pointer to thread data for error handling mainly.
52 */
53
54 #include "newton_diagnostics.h"
55 #include "../simulation_info_json.h"
56 #include "../jacobian_util.h"
57
58 extern int dgesv_(int *n, int *nrhs, double *a, int *lda,
59 int *ipiv, double *b, int *ldb, int *info);
60
61 extern int dgetrf_(int *n, int *nrhs, double *a, int *lda,
62 int *ipiv, int *info);
63
64 extern int dgetri_(int *n, double *a, int *lda,
65 int *ipiv, double *work, int *lwork, int *info);
66
67 // --------------------------------------------------------------------------------------------------------------------------------
68
69 ✗ unsigned var_id( unsigned idx, DATA* data, NONLINEAR_SYSTEM_DATA* systemData)
70 {
71 // Returns index of "modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[idx])"
72 // in "data->modelData->realVarsData[i]"
73
74 ✗ const char *name = modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[idx];
75 ✗ for (unsigned int i = 0; i < data->modelData->nVariablesReal; ++i) {
76 ✗ if (!strcmp(data->modelData->realVarsData[i].info.name, name)) {
77 ✗ return i;
78 }
79 }
80 return -1;
81 }
82
83 ✗ double** MatMult( unsigned rA, unsigned cArB, unsigned cB, double** A, double** B)
84 {
85 // Matrix multiplication A[rA][cArB] * B[cArB][cB] = C[rA][cB]
86
87 ✗ double** C = (double**)malloc(rA * sizeof(double*));
88 ✗ assertStreamPrint(NULL, NULL != C, "out of memory");
89 ✗ for (unsigned i = 0; i < rA; i++) {
90 ✗ C[i] = (double*)malloc(cB * sizeof(double));
91 ✗ assertStreamPrint(NULL, NULL != C[i], "out of memory");
92 }
93
94 ✗ for (unsigned i = 0; i < rA; i++) {
95 ✗ for (unsigned j = 0; j < cB; j++) {
96 ✗ C[i][j] = 0;
97 ✗ for (unsigned k = 0; k < cArB; k++)
98 ✗ C[i][j] += A[i][k] * B[k][j];
99 }
100 }
101
102 ✗ return C;
103 }
104
105 // --------------------------------------------------------------------------------------------------------------------------------
106
107 ✗ double** getJacobian( DATA* data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* systemData)
108 {
109 unsigned i, j;
110 ✗ size_t m = systemData->size;
111 JACOBIAN* jacobian = NULL;
112
113 modelica_real* jac = NULL;
114
115 // Allocate memory for fx (m * m matrix)
116 ✗ double** fx = (double**)malloc(m * sizeof(double*)); // freed by the caller
117 ✗ assertStreamPrint(threadData, NULL != fx, "out of memory");
118 ✗ for (i = 0; i < m; i++) {
119 ✗ fx[i] = (double*)malloc(m * sizeof(double)); // freed by the caller
120 ✗ assertStreamPrint(threadData, NULL != fx[i], "out of memory");
121 }
122
123 // Order of Jacobian elements:
124 // variable 1: df_1/dv_1, df_1/dv_2, .... df_1/dv_n
125 // variable 2: df_2/dv_1, df_2/dv_2, .... df_2/dv_n
126 // ...
127 // variable n: df_n/dv_1, df_2/dv_2, .... df_n/dv_n
128
129 ✗ if (systemData->jacobianIndex != -1) {
130 ✗ jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]);
131
132 ✗ jac = (modelica_real*) calloc(jacobian->sizeRows * jacobian->sizeCols, sizeof(modelica_real));
133 ✗ assertStreamPrint(threadData, NULL != jac, "out of memory");
134
135 /* call generic dense Jacobian */
136 ✗ evalJacobian(data, threadData, jacobian, NULL, jac, TRUE);
137
138 /* copy jacobian from column-major to row-major */
139 ✗ for (i = 0; i < jacobian->sizeRows; ++i)
140 ✗ for (j = 0; j < jacobian->sizeCols; ++j)
141 ✗ fx[i][j] = jac[j*jacobian->sizeRows + i];
142
143 ✗ free(jac);
144
145 } else {
146 ✗ assertStreamPrint(threadData, FALSE, "NEWTON_DIAGNOSTICS: numeric jacobian not yet supported.");
147 }
148
149 ✗ return fx;
150 }
151
152 // --------------------------------------------------------------------------------------------------------------------------------
153
154 ✗ double* getFirstNewtonStep( unsigned m, double* f, double** fx)
155 {
156 // Function values iteration 0: vector f(x0)
157 // Values Jacobian iteration 0: vector fx(x0)
158 // Newton step: dx = -f(x0)/fx(x0)
159
160 // Allocate memory for Newton steps
161 ✗ double* dx = (double*)malloc(m * sizeof(double));
162 ✗ assertStreamPrint(NULL, NULL != dx, "out of memory");
163
164 // Variables for Lapack routines
165 ✗ int N = m; // number of rows and columns of Jacobian
166 ✗ int NRHS = 1; // number of columns of b, i.e. f(x)
167 ✗ int LDA = N;
168 ✗ int LDB = N;
169 ✗ int* ipiv = (int*)malloc(N* sizeof(int));
170 ✗ assertStreamPrint(NULL, NULL != ipiv, "out of memory");
171 int info;
172
173 ✗ double* a = (double*)malloc( LDA * N * sizeof(double));
174 ✗ assertStreamPrint(NULL, NULL != a, "out of memory");
175 ✗ double* b = (double*)malloc( LDB * NRHS * sizeof(double));
176 ✗ assertStreamPrint(NULL, NULL != b, "out of memory");
177
178 unsigned i, j;
179
180 // Store Jacobian values J(x0) in a
181 ✗ for (i = 0; i < m; i++)
182 ✗ for (j = 0; j < m; j++)
183 ✗ a[m*i+j] = fx[j][i];
184
185 // Store function values f(x0) in b
186 ✗ for (i = 0; i < m; i++)
187 ✗ b[i] = f[i];
188
189 // Call Lapack function dgesv; after return, b contains the Newton steps
190 ✗ dgesv_(&N, &NRHS, a, &LDA, ipiv, b, &LDB, &info);
191
192 ✗ if (info > 0)
193 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "getFirstNewtonStep: the first Newton step could not be computed; the info satus is : %d", info);
194 else {
195 // Store Newton steps in dx
196 ✗ for (j = 0; j < m; j++)
197 ✗ dx[j] = -b[j];
198 }
199
200 ✗ free(ipiv);
201 ✗ free(a);
202 ✗ free(b);
203
204 ✗ return dx;
205 }
206
207 // --------------------------------------------------------------------------------------------------------------------------------
208
209 ✗ double maxNonLinearResiduals( unsigned m, unsigned l, unsigned* z_idx,
210 double* f, double** fx, double* dx)
211 {
212 // Calculate the absolute maximum value of the non-linear residuals r_x0 = f_x0 + fz * (z1 - z0)
213 // at iteration point x0, where z1 - z0 = dx and fz = J for the linear values and equations.
214
215 // l = m - q: number of linear unknowns
216 // z_idx : index of linear dependable in f, fx, dx
217
218 double r_x0, fz_dz;
219 double maxRes = 0; // Initialize to 0 for maximum search
220 unsigned i, j;
221
222 ✗ for (i = 0; i < m; i++) {
223 fz_dz = 0;
224 ✗ if (z_idx)
225 ✗ for (j = 0; j < l; j++) // iteration point x0 ==> j = 1 as r_x(j-1) = f_x(j-1) + fz * (z(j) - z(j-1)) = f_x(j-1) + fz * dz(j-1) ????
226 ✗ fz_dz += fx[i][z_idx[j]] * dx[z_idx[j]];
227
228 ✗ r_x0 = fabs(f[i] + fz_dz);
229 ✗ if (r_x0 > maxRes)
230 maxRes = r_x0;
231 }
232
233 ✗ return maxRes;
234 }
235
236 // --------------------------------------------------------------------------------------------------------------------------------
237
238 ✗ double*** getHessian( DATA* data, threadData_t *threadData, unsigned sysNumber, unsigned m)
239 {
240 ✗ NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
241 ✗ JACOBIAN* jac = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]);
242
243 unsigned i, j, k;
244 const modelica_real eps = 1.e-7;
245 const modelica_real nominal_x = 1.e-4;
246 ✗ SIMULATION_DATA *sData = data->localData[0];
247
248 // Allocate memory for Hessian fxx (m * m * m doubles)
249 ✗ double*** fxx = (double***)malloc(m * sizeof(double**));
250 ✗ assertStreamPrint(NULL, NULL != fxx, "out of memory");
251 ✗ for (i = 0; i < m; i++) {
252 ✗ fxx[i] = (double**)malloc(m * sizeof(double*));
253 ✗ assertStreamPrint(NULL, NULL != fxx[i], "out of memory");
254 ✗ for (j = 0; j < m; j++) {
255 ✗ fxx[i][j] = (double*)malloc(m * sizeof(double));
256 ✗ assertStreamPrint(NULL, NULL != fxx[i][j], "out of memory");
257 }
258 }
259
260 // Allocate memory for Jacobians
261 ✗ double** fxPls = (double**)malloc(m * sizeof(double*));
262 ✗ assertStreamPrint(NULL, NULL != fxPls, "out of memory");
263 ✗ double** fxMin = (double**)malloc(m * sizeof(double*));
264 ✗ assertStreamPrint(NULL, NULL != fxMin, "out of memory");
265 ✗ for (i = 0; i < m; i++) {
266 ✗ fxPls[i] = (double*)malloc(m * sizeof(double));
267 ✗ assertStreamPrint(NULL, NULL != fxPls[i], "out of memory");
268 ✗ fxMin[i] = (double*)malloc(m * sizeof(double));
269 ✗ assertStreamPrint(NULL, NULL != fxMin[i], "out of memory");
270 }
271
272 // ----------------------------------------------- Debug -------------------------------------------------
273 /*printf( "\n");
274 for ( k = 0; k < m; k++) {
275 unsigned id = var_id(k, data, systemData);
276 printf( " k = %d: id = %d (%s)\n", k, id, data->modelData->realVarsData[id].info.name);
277 }*/
278 // -------------------------------------------- end of Debug ---------------------------------------------
279
280 ✗ for (k = 0; k < m; k++) {
281 ✗ unsigned id = var_id(k, data, systemData);
282
283 ✗ double tmp_x = sData->realVars[id];
284 ✗ const modelica_real delta_x = eps * fmax( fabs(tmp_x), nominal_x);
285
286 ✗ sData->realVars[id] = tmp_x + delta_x;
287 ✗ for (j = 0; j < m; j++) {
288 ✗ jac->seedVars[j] = 1.0;
289 ✗ systemData->analyticalJacobianColumn(data, threadData, jac, NULL);
290 ✗ for (i = 0; i < m; i++)
291 ✗ fxPls[i][j] = jac->resultVars[i];
292 ✗ jac->seedVars[j] = 0.0;
293 }
294
295 ✗ sData->realVars[id] = tmp_x - delta_x;
296 ✗ for (j = 0; j < m; j++) {
297 ✗ jac->seedVars[j] = 1.0;
298 ✗ systemData->analyticalJacobianColumn(data, threadData, jac, NULL);
299 ✗ for (i = 0; i < m; i++)
300 ✗ fxMin[i][j] = jac->resultVars[i];
301 ✗ jac->seedVars[j] = 0.0;
302 }
303
304 ✗ sData->realVars[id] = tmp_x;
305
306 ✗ for (j = 0; j < m; j++)
307 ✗ for (i = 0; i < m; i++) {
308 ✗ fxx[i][k][j] = (fxPls[i][j] - fxMin[i][j]) / (2 * delta_x);
309 ✗ if (isnan(fxx[i][k][j])) {
310 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "NaN detected: fxx[%d][%d][%d]: fxPls[%d][%d] = %f, fxMin[%d][%d] = %f, delta_x = %f\n",
311 i+1,j+1,k+1, i+1,j+1,fxPls[i][j], i+1,j+1,fxMin[i][j], delta_x);
312 ✗ return fxx;
313 }
314 }
315 }
316
317 ✗ for (i = 0; i < m; i++) {
318 ✗ free(fxPls[i]);
319 ✗ free(fxMin[i]);
320 }
321 ✗ free(fxPls);
322 ✗ free(fxMin);
323
324 // ----------------------------------------------- Debug -------------------------------------------------
325 /*printf( "\n");
326 for (k = 0; k < m; k++) {
327 // For each eqn. k print m x m matrix
328 for (i = 0; i < m; i++) {
329 if (i == 0)
330 printf( "\n\neqn k = %2d: ", k);
331 else
332 printf( "\n ");
333
334 for (j = 0; j < m; j++)
335 printf( "%7.2f ", fxx[k][i][j]);
336 }
337 }
338 printf( "\n\n");*/
339 // -------------------------------------------- end of Debug ---------------------------------------------
340
341 ✗ return fxx;
342 }
343
344 // --------------------------------------------------------------------------------------------------------------------------------
345
346 ✗ double*** calcGamma( unsigned m, unsigned p, unsigned q, unsigned* n_idx,
347 unsigned* w_idx, double* dx, double*** fxx, double maxRes)
348 {
349 // Calculation of curvature factors Gamma_ijk
350 // ------------------------------------------
351 //
352 // m : total number of equations/independents
353 // p : number of non-linear equations
354 // q : number of non-linear independents
355 // n_idx : index of non-linear equation, i.e. of i in fxx[i][j][k]
356 // w_idx : index of non-linear dependent, i.e. of j in dx[j] and of j and k in fxx[i][j][k]
357 // dx : Newton step first iteration = x1 - x0, and dx[w_idx] = w1 - w0
358 // fxx : Hessian as function of x0
359 // maxRes: absolute maximum value of the non-linear residuals
360
361 unsigned i, j, k;
362
363 // Allocate memory for Gamma_ijk (p * q * q doubles)
364 ✗ double*** Gamma_ijk = (double***)malloc(p * sizeof(double**));
365 ✗ assertStreamPrint(NULL, NULL != Gamma_ijk, "out of memory");
366 ✗ for (i = 0; i < p; i++) {
367 ✗ Gamma_ijk[i] = (double**)malloc(q * sizeof(double*));
368 ✗ assertStreamPrint(NULL, NULL != Gamma_ijk[i], "out of memory");
369 ✗ for (j = 0; j < q; j++) {
370 ✗ Gamma_ijk[i][j] = (double*)malloc(q * sizeof(double));
371 ✗ assertStreamPrint(NULL, NULL != Gamma_ijk[i][j], "out of memory");
372 }
373 }
374
375 // Calculate Gamma_ijk
376 ✗ for (i = 0; i < p; i++)
377 ✗ for (j = 0; j < q; j++)
378 ✗ for (k = 0; k < q; k++)
379 ✗ if (!isnan(fxx[n_idx[i]][w_idx[j]][w_idx[k]]) && fxx[n_idx[i]][w_idx[j]][w_idx[k]] != 0)
380 ✗ Gamma_ijk[i][j][k] = fabs(0.5 * fxx[n_idx[i]][w_idx[j]][w_idx[k]] * (dx[w_idx[j]] * dx[w_idx[k]]) / maxRes);
381 else
382 ✗ Gamma_ijk[i][j][k] = 0;
383
384 ✗ return Gamma_ijk;
385 }
386
387 // --------------------------------------------------------------------------------------------------------------------------------
388
389 ✗ double* calcAlpha( DATA* data, threadData_t* threadData, unsigned sysNumber, unsigned m, unsigned p,
390 unsigned q, unsigned* n_idx, unsigned* w_idx, double* x, double* dx,
391 double* f, double*** fxx, double lambda, double maxRes)
392 {
393 // Calculation of alpha coefficients for all non-linear equations
394 // --------------------------------------------------------------
395 //
396 // m : total number of equations/independents
397 // p : number of non-linear equations
398 // q : number of non-linear independents
399 // n_idx : index of non-linear equation, ie of i in f[i] and fxx[i][j][k]
400 // w_idx : index of non-linear dependent, ie of j in x[j] and dx[j] and of j and k in fxx[i][j][k]
401 // x : all independents (non-linear & linear)
402 // dx : Newton step first iteration = x1 - x0, and dx[w_idx] = w1 - w0
403 // f : Function values (ie residuals) as function of x0
404 // fxx : Hessian as function of x0
405 // lambda: damping factor
406 // maxRes: absolute maximum value of the non-linear residuals of iteration 0
407
408 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
409 ✗ NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
410
411 unsigned i, j, k;
412
413 // Allocate memory for alpha (p doubles)
414 ✗ double* alpha = (double*)malloc(p * sizeof(double));
415 ✗ assertStreamPrint(NULL, NULL != alpha, "out of memory");
416
417 // Get damped guess x1_star for second iteration step
418 ✗ double* x1_star = (double*)malloc(m * sizeof(double));
419 ✗ assertStreamPrint(NULL, NULL != x1_star, "out of memory");
420 ✗ for (j = 0; j < m; j++)
421 ✗ x1_star[j] = x[j] + lambda * dx[j];
422
423 // Calculate residuals f_x1_star for damped guess x1_star
424 ✗ double* f_x1_star = (double*)malloc(m * sizeof(double));
425 ✗ assertStreamPrint(NULL, NULL != f_x1_star, "out of memory");
426 ✗ systemData->residualFunc(&resUserData, x1_star, f_x1_star, (int*)&systemData->size);
427
428 // For each non-linear independent get w1_star - w0
429 ✗ double* w1_star_w0 = (double*)malloc(q * sizeof(double));
430 ✗ assertStreamPrint(NULL, NULL != w1_star_w0, "out of memory");
431 ✗ for (j = 0; j < q; j++)
432 ✗ w1_star_w0[j] = lambda * dx[w_idx[j]];
433
434 // Calculate alpha for each non-linear equation i
435 ✗ double* w_times_fww_w0 = (double*)malloc(q * sizeof(double));
436 ✗ assertStreamPrint(NULL, NULL != w_times_fww_w0, "out of memory");
437
438 ✗ for (i = 0; i < p; i++) {
439 // Vector w_times_fww_w0 = (w1_star - w0)' * fww_w0 (1 x q * q x q --> 1 x q vector)
440 ✗ for (j = 0; j < q; j++) {
441 // For each independent
442 ✗ w_times_fww_w0[j] = 0;
443 ✗ for (k = 0; k < q; k++) {
444 ✗ if (!isnan(fxx[n_idx[i]][w_idx[k]][w_idx[j]]) && fabs(fxx[n_idx[i]][w_idx[k]][w_idx[j]]) != 0)
445 ✗ w_times_fww_w0[j] += w1_star_w0[k] * fxx[n_idx[i]][w_idx[k]][w_idx[j]];
446 }
447 }
448
449 // Scalar w_times_f_times_w = w_times_f_i_ww_w0 * (w1_star - w0) = (w1_star - w0)' * f_i_ww_w0 * (w1_star - w0)
450 // (1 x q * q x 1 vector --> scalar)
451 double w_times_fww_times_w = 0;
452 ✗ for (k = 0; k < q; k++)
453 ✗ w_times_fww_times_w += w_times_fww_w0[k] * w1_star_w0[k];
454
455 // Calculate alpha for the non-linear equations
456 ✗ alpha[i] = fabs(f_x1_star[n_idx[i]] - (1 - lambda) * f[n_idx[i]] - 0.5 * w_times_fww_times_w) / (pow(lambda,3) * maxRes);
457 }
458
459 ✗ free(w_times_fww_w0);
460 ✗ free(w1_star_w0);
461 ✗ free(f_x1_star);
462 ✗ free(x1_star);
463
464 ✗ return alpha;
465 }
466
467 // --------------------------------------------------------------------------------------------------------------------------------
468
469 ✗ double** getInvJacobian( unsigned m, double** fx)
470 {
471 // Calculates inverse matrix of Jacobian fx as function of x0 (m x m matrix)
472 // -------------------------------------------------------------------------
473 //
474 // m : total number of equations/independents
475 // fx : Jacobian as function of x0
476
477 unsigned i, j;
478
479 // Intialize inverse a with fx
480 ✗ double* a = (double*)malloc(m * m * sizeof(double));
481 ✗ assertStreamPrint(NULL, NULL != a, "out of memory");
482 ✗ for (i = 0; i < m; i++)
483 ✗ for (j = 0; j < m; j++)
484 ✗ a[m*i+j] = fx[j][i];
485
486 // Variables for Lapack routines
487 ✗ int N = m;
488 ✗ int LWORK = N * N;
489 ✗ int* ipiv = (int*)malloc(N * sizeof(int));
490 ✗ assertStreamPrint(NULL, NULL != ipiv, "out of memory");
491 int info;
492 ✗ double* WORK = (double*)malloc(LWORK * sizeof(double));
493 ✗ assertStreamPrint(NULL, NULL != WORK, "out of memory");
494
495 // Call Lapack function dgetrf to compute the LU factorization of fx
496 ✗ dgetrf_(&N, &N, a, &N, ipiv, &info);
497 ✗ if (info > 0)
498 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "getInvJacobian: LU factorization could not be computed; the info status is : %d", info);
499
500 // Call Lapack function dgetri to compute the inverse of fx
501 ✗ dgetri_(&N, a, &N, ipiv, WORK, &LWORK, &info);
502 ✗ if (info > 0)
503 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "getInvJacobian: inverse Jacobian could not be computed; the info satus is : %d", info);
504
505 // Return two dimensional array
506 ✗ double** inv_fx = (double**)malloc(m * sizeof(double*));
507 ✗ assertStreamPrint(NULL, NULL != inv_fx, "out of memory");
508 ✗ for (i = 0; i < m; i++) {
509 ✗ inv_fx[i] = (double*)malloc(m * sizeof(double));
510 ✗ assertStreamPrint(NULL, NULL != inv_fx[i], "out of memory");
511 }
512 ✗ for (i = 0; i < m; i++)
513 ✗ for (j = 0; j < m; j++)
514 ✗ inv_fx[j][i] = a[m*i+j];
515
516 ✗ free(ipiv);
517 ✗ free(WORK);
518 ✗ free(a);
519
520 ✗ return inv_fx;
521 }
522
523 // --------------------------------------------------------------------------------------------------------------------------------
524
525 ✗ double** calcSigma( unsigned m, unsigned q, unsigned* w_idx,
526 double* dx, double** fx, double*** fxx)
527 {
528 // Calculation of solution sensitivities Sigma_ij
529 // ----------------------------------------------
530 //
531 // m : total number of equations/independents
532 // q : number of non-linear variables
533 // w_idx : index of non-linear dependent, i.e. of j in dx[j] and of j and k in fxx[i][j][k]
534 // dx : Newton step first iteration = x1 - x0, and dx[w_idx] = w1 - w0
535 // fx : Jacobian as function of x0
536 // fxx : Hessian as function of x0
537
538 unsigned i, j, k;
539
540 // Calculate inverse Jacobian, i.e. inverse matrix of fx
541 ✗ double** inv_fx = getInvJacobian( m, fx);
542
543 // Get matrix H[i] = (x1 - x0)' * fxx = dx' * fxx (1 x m * m x m matrix --> m vector)
544 ✗ double** H_i = (double**)malloc(m * sizeof(double*)); // m functions * m vectors --> m x m matrix
545 ✗ assertStreamPrint(NULL, NULL != H_i, "out of memory");
546 ✗ for (i = 0; i < m; i++) {
547 ✗ H_i[i] = (double*)malloc(m * sizeof(double));
548 ✗ assertStreamPrint(NULL, NULL != H_i[i], "out of memory");
549 }
550 ✗ for (i = 0; i < m; i++) {
551 ✗ for (j = 0; j < m; j++) {
552 ✗ H_i[i][j] = 0;
553 ✗ for (k = 0; k < m; k++)
554 ✗ H_i[i][j] += dx[k] * fxx[i][k][j];
555 }
556 }
557
558 // Calculate tmp1 = -inv_fx * H_i
559 // (m x m matrix) * (m x m matrix) --> m x m matrix
560 ✗ for (i = 0; i < m; i++)
561 ✗ for (j = 0; j < m; j++)
562 ✗ inv_fx[i][j] = -inv_fx[i][j];
563
564 ✗ double** tmp1 = MatMult( m, m, m, inv_fx, H_i);
565
566 // Extract matrix tmp2 from tmp1 for only non-linears (q x q matrix)
567 ✗ double** tmp2 = (double**)malloc(q * sizeof(double*));
568 ✗ assertStreamPrint(NULL, NULL != tmp2, "out of memory");
569 ✗ for (i = 0; i < q; i++) {
570 ✗ tmp2[i] = (double*)malloc(q * sizeof(double));
571 ✗ assertStreamPrint(NULL, NULL != tmp2[i], "out of memory");
572 }
573 ✗ for (i = 0; i < q; i++)
574 ✗ for (j = 0; j < q; j++)
575 ✗ tmp2[i][j] = tmp1[w_idx[i]][w_idx[j]];
576
577
578 // Create a q x q matrix wDiag with w1 - w0 = dx[w_idx] on diagonal
579 ✗ double** wDiag = (double**)malloc(q * sizeof(double*));
580 ✗ assertStreamPrint(NULL, NULL != wDiag, "out of memory");
581 ✗ for (i = 0; i < q; i++) {
582 ✗ wDiag[i] = (double*)malloc(q * sizeof(double));
583 ✗ assertStreamPrint(NULL, NULL != wDiag[i], "out of memory");
584 }
585 ✗ for (i = 0; i < q; i++) {
586 ✗ for (j = 0; j < q; j++)
587 ✗ if (i == j)
588 ✗ wDiag[i][j] = dx[w_idx[i]];
589 else
590 ✗ wDiag[i][j] = 0;
591 }
592
593 // Get inverse matrix inv_wDiag of wDiag
594 ✗ double** inv_wDiag = getInvJacobian( q, wDiag);
595
596 // Calculate tmp3 = | inv_wDiag | * tmp2
597 // (q x q matrix) * (q x q matrix) --> q x q matrix
598 ✗ for (i = 0; i < q; i++)
599 ✗ for (j = 0; j < q; j++)
600 ✗ inv_wDiag[i][j] = fabs(inv_wDiag[i][j]);
601
602 ✗ double** tmp3 = MatMult( q, q, q, inv_wDiag, tmp2);
603
604 // Calculate Sigma = tmp3 * wDiag = | inv_wDiag | * tmp2 * wDiag = | inv_wDiag | * -inv_fx * H_i * wDiag
605 ✗ double** Sigma = MatMult( q, q, q, tmp3, wDiag);
606
607 // Free dynamically allocated memory
608 ✗ for (i = 0; i < m; i++) {
609 ✗ free(inv_fx[i]);
610 ✗ free(H_i[i]);
611 ✗ free(tmp1[i]);
612 }
613 ✗ free(inv_fx);
614 ✗ free(H_i);
615 ✗ free(tmp1);
616
617 ✗ for (i = 0; i < q; i++) {
618 ✗ free(wDiag[i]);
619 ✗ free(inv_wDiag[i]);
620 ✗ free(tmp2[i]);
621 ✗ free(tmp3[i]);
622 }
623 ✗ free(wDiag);
624 ✗ free(inv_wDiag);
625 ✗ free(tmp2);
626 ✗ free(tmp3);
627
628 ✗ return Sigma;
629 }
630
631 // --------------------------------------------------------------------------------------------------------------------------------
632
633 ✗ void PrintResults( DATA* data, unsigned sysNumber, unsigned m, unsigned p, unsigned q, unsigned* n_idx, unsigned* w_idx,
634 double* x0, double* alpha, double*** Gamma_ijk, double** Sigma_ij)
635 {
636 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Values of relevant indicators");
637
638 ✗ NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
639 modelica_integer sizeOfTorns;
640 unsigned i, j, k;
641 double eps = 1e-2;
642
643 // ----------------------------------------------- Debug -------------------------------------------------
644 /*printf(" Equations\n");
645 for (i = 0; i < p; i++)
646 printf("\n i =%2d: %s", n_idx[i]+1, "???"); //, modelInfoGetFunction(&data->modelData->modelDataXml, i).name);
647 printf("\n\n");
648
649 printf(" Variables, initial guesses\n");
650 for (j = 0; j < q; j++)
651 printf("\n j =%2d: %8s = %10.8g", w_idx[j]+1, data->modelData->realVarsData[var_id(w_idx[j], data, systemData)].info.name, x0[w_idx[j]]);
652 printf("\n\n");*/
653 // ------------------------------------------- end of Debug ----------------------------------------------
654
655 // Print alpha, Gamma, and Sigma if value > eps
656 // --------------------------------------------
657
658 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "alpha_i > %5.3f", eps);
659 ✗ for (i = 0; i < p; ++i)
660 ✗ if (alpha[i] > eps)
661 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "alpha_%-3d = %5.2f", n_idx[i]+1, alpha[i]);
662 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
663
664 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Gamma_ijk > %5.3f", eps);
665 ✗ for (i = 0; i < p; i++)
666 ✗ for (j = 0; j < q; j++)
667 ✗ for (k = j; k < q; k++)
668 ✗ if (Gamma_ijk[i][j][k] > eps)
669 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Gamma_%-4d_%-4d_%-4d = %5.2f", n_idx[i]+1, w_idx[j]+1, w_idx[k]+1, Gamma_ijk[i][j][k]);
670 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
671
672 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "sigma_jj > %5.3f", eps);
673 ✗ for (i = 0; i < q; i++)
674 ✗ if (fabs(Sigma_ij[i][i]) > eps)
675 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "sigma_%-4d_%-4d = %5.2f", w_idx[i]+1, w_idx[i]+1, fabs(Sigma_ij[i][i]));
676 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
677
678 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS); // This closes the "Values of relevant indicators" section
679
680
681 // Select values of Gamma and Sigma > eps and store them in descending order
682 // -------------------------------------------------------------------------
683
684 double val_largest_alpha, val_largest_Sigma, val_largest_Gamma;
685
686 unsigned idx_largest_alpha, idx_largest_Sigma, l, n_gt_eps = 0,
687 idx_largest_G_i, idx_largest_G_j, idx_largest_G_k;
688
689 ✗ unsigned* alpha_checked = (unsigned*)malloc(p * sizeof(unsigned));
690 ✗ assertStreamPrint(NULL, NULL != alpha_checked, "out of memory");
691 ✗ unsigned* Sigma_checked = (unsigned*)malloc(q * sizeof(unsigned));
692 ✗ assertStreamPrint(NULL, NULL != Sigma_checked, "out of memory");
693 ✗ unsigned*** Gamma_checked = (unsigned***)malloc(p * sizeof(unsigned**));
694 ✗ assertStreamPrint(NULL, NULL != Gamma_checked, "out of memory");
695 ✗ for (i = 0; i < p; i++) {
696 ✗ Gamma_checked[i] = (unsigned**)malloc(q * sizeof(unsigned*));
697 ✗ assertStreamPrint(NULL, NULL != Gamma_checked[i], "out of memory");
698 ✗ for (j = 0; j < q; j++) {
699 ✗ Gamma_checked[i][j] = (unsigned*)malloc(q * sizeof(unsigned));
700 ✗ assertStreamPrint(NULL, NULL != Gamma_checked[i][j], "out of memory");
701 }
702 }
703 ✗ unsigned* index_alpha = (unsigned*)malloc((p * q * q + m) * sizeof(unsigned));
704 ✗ assertStreamPrint(NULL, NULL != index_alpha, "out of memory");
705 ✗ unsigned* index_Sigma = (unsigned*)malloc((p * q * q + m) * sizeof(unsigned));
706 ✗ assertStreamPrint(NULL, NULL != index_Sigma, "out of memory");
707 ✗ unsigned* index_Gamma_i = (unsigned*)malloc((p * q * q + m) * sizeof(unsigned));
708 ✗ assertStreamPrint(NULL, NULL != index_Gamma_i, "out of memory");
709 ✗ unsigned* index_Gamma_j = (unsigned*)malloc((p * q * q + m) * sizeof(unsigned));
710 ✗ assertStreamPrint(NULL, NULL != index_Gamma_j, "out of memory");
711 ✗ unsigned* index_Gamma_k = (unsigned*)malloc((p * q * q + m) * sizeof(unsigned));
712 ✗ assertStreamPrint(NULL, NULL != index_Gamma_k, "out of memory");
713
714 // Initialize tmp arrays for sorting
715 ✗ for (i = 0; i < p; i++) {
716 ✗ for (j = 0; j < q; j++)
717 ✗ for (k = 0; k < q; k++)
718 ✗ Gamma_checked[i][j][k] = 0;
719 }
720 ✗ for (j = 0; j < q; j++)
721 ✗ Sigma_checked[j] = 0;
722
723 ✗ for (l = 0; l < p * q * q + m; l++) {
724 // Select largest Gamma variable and its value
725 val_largest_Gamma = -1.e10;
726 idx_largest_G_i = 0;
727 idx_largest_G_j = 0;
728 idx_largest_G_k = 0;
729 ✗ for (i = 0; i < p; i++)
730 ✗ for (j = 0; j < q; j++)
731 ✗ for (k = j; k < q; k++)
732 ✗ if (Gamma_ijk[i][j][k] > val_largest_Gamma && !Gamma_checked[i][j][k]) {
733 val_largest_Gamma = Gamma_ijk[i][j][k];
734 idx_largest_G_i = i;
735 idx_largest_G_j = j;
736 idx_largest_G_k = k;
737 }
738
739 // Select largest Sigma variable and its value
740 val_largest_Sigma = -1.e10;
741 idx_largest_Sigma = 0;
742 ✗ for (i = 0; i < q; i++)
743 ✗ if (fabs(Sigma_ij[i][i]) > val_largest_Sigma && !Sigma_checked[i]) {
744 val_largest_Sigma = fabs(Sigma_ij[i][i]);
745 idx_largest_Sigma = i;
746 }
747
748 // Values < 0 , i.e. less than eps are not considered
749 ✗ if (val_largest_Gamma < eps && val_largest_Sigma < eps) break;
750
751 // Checkmark and store indices of largest value
752 ✗ if (val_largest_Gamma > val_largest_Sigma) {
753 ✗ index_Gamma_i[n_gt_eps] = idx_largest_G_i;
754 ✗ index_Gamma_j[n_gt_eps] = idx_largest_G_j;
755 ✗ index_Gamma_k[n_gt_eps] = idx_largest_G_k;
756 ✗ index_Sigma[n_gt_eps] = -1;
757 ✗ Gamma_checked[idx_largest_G_i][idx_largest_G_j][idx_largest_G_k] = 1;
758
759 // -------------------------------------------- Debug ----------------------------------------------
760 //printf("\n Gamma_%d_%d_%d = %8.3f", n_idx[idx_largest_G_i]+1, w_idx[idx_largest_G_j]+1,
761 // w_idx[idx_largest_G_k]+1, val_largest_Gamma);
762 // ---------------------------------------- end of Debug -------------------------------------------
763 } else {
764 ✗ index_Sigma[n_gt_eps] = idx_largest_Sigma;
765 ✗ index_Gamma_i[n_gt_eps] = -1;
766 ✗ index_Gamma_j[n_gt_eps] = -1;
767 ✗ index_Gamma_k[n_gt_eps] = -1;
768 ✗ Sigma_checked[idx_largest_Sigma] = 1;
769
770 // -------------------------------------------- Debug ---------------------------------------------
771 //printf("\n Sigma_%d_%d = %8.3f", w_idx[idx_largest_Sigma]+1, w_idx[idx_largest_Sigma]+1,
772 // fabs(Sigma_ij[idx_largest_Sigma][idx_largest_Sigma]));
773 // ---------------------------------------- end of Debug -------------------------------------------
774 }
775
776 // Increment number of values found
777 ✗ n_gt_eps++;
778 }
779
780 // Print ranked indicators
781 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Ranked indicators");
782
783
784 // Print variables referenced by Sigma and Gamma values > eps and concerned Sigma or Gamma value
785 // ---------------------------------------------------------------------------------------------
786 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "By variable");
787
788 ✗ unsigned* printedIdx = (unsigned*)malloc(2 * n_gt_eps * sizeof(unsigned));
789 ✗ assertStreamPrint(NULL, NULL != printedIdx, "out of memory");
790 unsigned nPrinted = 0;
791 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Var no. Var name Initial guess max(Gamma,sigma)");
792 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "------- ---------------------------------------- ------------- ----------------");
793 ✗ for (l = 0; l < n_gt_eps; l++) {
794 ✗ printedIdx[nPrinted] = -1;
795 ✗ if (0 <= index_Sigma[l] && index_Sigma[l] < q) {
796 // Check if variable l referenced by Sigma has already been printed for Gamma
797 unsigned alreadyPrinted = 0;
798 ✗ for (unsigned nP = 0; nP < nPrinted && !alreadyPrinted; nP++)
799 ✗ alreadyPrinted = index_Sigma[l] == printedIdx[nP];
800
801 ✗ if (!alreadyPrinted) {
802 // Print variable referenced l by Sigma, its init value and the max value between Gamma and Sigma
803 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "%7d %40s %13.7g %5.2f",
804 w_idx[index_Sigma[l]]+1,
805 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[index_Sigma[l]]],
806 ✗ x0[w_idx[index_Sigma[l]]],
807 ✗ fabs(Sigma_ij[index_Sigma[l]][index_Sigma[l]]));
808 ✗ printedIdx[nPrinted++] = index_Sigma[l];
809 }
810 ✗ } else if (0 <= index_Gamma_i[l] && index_Gamma_i[l] < p &&
811 ✗ 0 <= index_Gamma_j[l] && index_Gamma_j[l] < q &&
812 ✗ 0 <= index_Gamma_k[l] && index_Gamma_k[l] < q)
813 {
814 // Check if variable l referenced by Gamma has already been printed for Sigma
815 unsigned alreadyPrinted_j = 0;
816 unsigned alreadyPrinted_k = 0;
817 ✗ for (unsigned nP = 0; nP < nPrinted; nP++) {
818 ✗ alreadyPrinted_j = alreadyPrinted_j || index_Gamma_j[l] == printedIdx[nP];
819 ✗ alreadyPrinted_k = alreadyPrinted_k || index_Gamma_k[l] == printedIdx[nP];
820 }
821
822 ✗ if (!alreadyPrinted_j) {
823 // Print variable referenced l by Gamma, its init value and the value of Gamma_ilk
824 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "%7d %40s %13.7g %5.2f",
825 w_idx[index_Gamma_j[l]]+1,
826 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[index_Gamma_j[l]]],
827 ✗ x0[w_idx[index_Gamma_j[l]]],
828 ✗ Gamma_ijk[index_Gamma_i[l]][index_Gamma_j[l]][index_Gamma_k[l]]);
829 ✗ printedIdx[nPrinted++] = index_Gamma_j[l];
830 }
831 ✗ if (!alreadyPrinted_k) {
832 // Print variable referenced l by Gamma, its init value and the value of Gamma_ijl
833 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "%7d %40s %13.7g %5.2f",
834 w_idx[index_Gamma_k[l]]+1,
835 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[index_Gamma_k[l]]],
836 ✗ x0[w_idx[index_Gamma_k[l]]],
837 ✗ Gamma_ijk[index_Gamma_i[l]][index_Gamma_j[l]][index_Gamma_k[l]]);
838 ✗ printedIdx[nPrinted++] = index_Gamma_k[l];
839 }
840 }
841 }
842 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
843
844 // Select values of alpha and Gamma > eps and store them in descending order
845 // -------------------------------------------------------------------------
846
847 n_gt_eps = 0;
848
849 // Initialize tmp arrays for sorting
850 ✗ for (i = 0; i < p; i++) {
851 ✗ for (j = 0; j < q; j++)
852 ✗ for (k = 0; k < q; k++)
853 ✗ Gamma_checked[i][j][k] = 0;
854 ✗ alpha_checked[i] = 0;
855 }
856
857 ✗ for (l = 0; l < p * q * q + m; l++) {
858 // Select largest Gamma variable and its value
859 val_largest_Gamma = -1.e10;
860 idx_largest_G_i = 0;
861 idx_largest_G_j = 0;
862 idx_largest_G_k = 0;
863 ✗ for (i = 0; i < p; i++)
864 ✗ for (j = 0; j < q; j++)
865 ✗ for (k = j; k < q; k++)
866 ✗ if ( Gamma_ijk[i][j][k] > val_largest_Gamma && !Gamma_checked[i][j][k]) {
867 val_largest_Gamma = Gamma_ijk[i][j][k];
868 idx_largest_G_i = i;
869 idx_largest_G_j = j;
870 idx_largest_G_k = k;
871 }
872
873 // Select largest alpha variable and its value
874 val_largest_alpha = -1.e10;
875 idx_largest_alpha = 0;
876 ✗ for (i = 0; i < p; i++)
877 ✗ if (alpha[i] > val_largest_alpha && !alpha_checked[i]) {
878 val_largest_alpha = alpha[i];
879 idx_largest_alpha = i;
880 }
881
882 // Values < 0 , i.e. less than eps, are not considered
883 ✗ if (val_largest_Gamma < eps && val_largest_alpha < eps) break;
884
885 // Checkmark and store indices of largest value
886 ✗ if (val_largest_Gamma > val_largest_alpha) {
887 ✗ index_Gamma_i[n_gt_eps] = idx_largest_G_i;
888 ✗ index_Gamma_j[n_gt_eps] = idx_largest_G_j;
889 ✗ index_Gamma_k[n_gt_eps] = idx_largest_G_k;
890 ✗ index_alpha[n_gt_eps] = -1;
891 ✗ Gamma_checked[idx_largest_G_i][idx_largest_G_j][idx_largest_G_k] = 1;
892
893 // -------------------------------------------- Debug ----------------------------------------------
894 //printf("\n Gamma_%d_%d_%d = %8.3f", n_idx[idx_largest_G_i]+1, w_idx[idx_largest_G_j]+1,
895 // w_idx[idx_largest_G_k]+1, val_largest_Gamma);
896 // ---------------------------------------- end of Debug -------------------------------------------
897 } else {
898 ✗ index_alpha[n_gt_eps] = idx_largest_alpha;
899 ✗ index_Gamma_i[n_gt_eps] = -1;
900 ✗ index_Gamma_j[n_gt_eps] = -1;
901 ✗ index_Gamma_k[n_gt_eps] = -1;
902 ✗ alpha_checked[idx_largest_alpha] = 1;
903
904 // -------------------------------------------- Debug ---------------------------------------------
905 //printf("\n alpha_%d = %8.3f", w_idx[idx_largest_alpha]+1, alpha[idx_largest_alpha]);
906 // ---------------------------------------- end of Debug -------------------------------------------
907 }
908
909 // Increment number of values found
910 ✗ n_gt_eps++;
911 }
912 //printf("\n\n");
913
914 // Print equations referenced by alpha and Gamma values > eps and concerned alpha or Gamma value
915 // ---------------------------------------------------------------------------------------------
916
917 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "By equation");
918
919 ✗ printedIdx = (unsigned*)realloc(printedIdx, n_gt_eps * sizeof(unsigned));
920 nPrinted = 0;
921 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Eq no. Eq idx max(alpha,Gamma)\n");
922 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "------ ------ ----------------");
923 ✗ for (l = 0; l < n_gt_eps; l++) {
924 ✗ printedIdx[nPrinted] = -1;
925 ✗ sizeOfTorns = systemData->torn_plus_residual_size - systemData->size;
926 ✗ if (0 <= index_alpha[l] && index_alpha[l] < p) {
927 // Check if equation l referenced by alpha has already been printed for Gamma
928 unsigned alreadyPrinted = 0;
929 ✗ for (unsigned nP = 0; nP < nPrinted && !alreadyPrinted; nP++)
930 ✗ alreadyPrinted = index_alpha[l] == printedIdx[nP];
931
932 ✗ if (!alreadyPrinted) {
933 // Print equation l referenced by alpha and the value of alpha_i
934 ✗ if (alpha[index_alpha[l]] < 1.e3)
935 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "%6d %6d %5.2f", n_idx[index_alpha[l]]+1,
936 ✗ systemData->eqn_simcode_indices[sizeOfTorns + n_idx[index_alpha[l]]], alpha[index_alpha[l]]);
937 else
938 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "%6d %6d %5.2e", n_idx[index_alpha[l]]+1,
939 ✗ systemData->eqn_simcode_indices[sizeOfTorns + n_idx[index_alpha[l]]], alpha[index_alpha[l]]);
940 ✗ printedIdx[nPrinted++] = index_alpha[l];
941 }
942 ✗ } else if (0 <= index_Gamma_i[l] && index_Gamma_i[l] < p &&
943 ✗ 0 <= index_Gamma_j[l] && index_Gamma_j[l] < q &&
944 ✗ 0 <= index_Gamma_k[l] && index_Gamma_k[l] < q)
945 {
946 // Check if equation l referenced by Gamma has already been printed for alpha
947 unsigned alreadyPrinted = 0;
948 ✗ for (unsigned nP = 0; nP < nPrinted && !alreadyPrinted; nP++)
949 ✗ alreadyPrinted = index_Gamma_i[l] == printedIdx[nP];
950
951 ✗ if (!alreadyPrinted) {
952 // Print equation l referenced by Gamma and the value of Gamma_ljk
953 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "%6d %6d %5.2f", n_idx[index_Gamma_i[l]]+1,
954 ✗ systemData->eqn_simcode_indices[sizeOfTorns + n_idx[index_Gamma_i[l]]],
955 ✗ Gamma_ijk[index_Gamma_i[l]][index_Gamma_j[l]][index_Gamma_k[l]]);
956 ✗ printedIdx[nPrinted++] = index_Gamma_i[l];
957 }
958 }
959 }
960 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
961
962 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
963
964 ✗ free(printedIdx);
965 ✗ free(alpha_checked);
966 ✗ free(Sigma_checked);
967 ✗ for (i = 0; i < p; i++) {
968 ✗ for (j = 0; j < q; j++)
969 ✗ free(Gamma_checked[i][j]);
970 ✗ free(Gamma_checked[i]);
971 }
972 ✗ free(Gamma_checked);
973 ✗ free(index_alpha);
974 ✗ free(index_Sigma);
975 ✗ free(index_Gamma_i);
976 ✗ free(index_Gamma_j);
977 ✗ free(index_Gamma_k);
978 ✗ }
979
980 // --------------------------------------------------------------------------------------------------------------------------------
981
982 ✗ unsigned* getNonlinearEqns( DATA* data, threadData_t* threadData, unsigned sysNumber,
983 unsigned m, double* f_x0, double* x0, double* dx, double* lambda, unsigned* p)
984 {
985 // If |f^i(x1)| > 0, then f^i is a nonlinear function
986
987 ✗ RESIDUAL_USERDATA resUserData = {.data=data, .threadData=threadData, .solverData=NULL};
988 ✗ NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
989
990 unsigned i;
991 double eps = 1.e-9;
992
993 // Calculate x1 from NewtonFirstStep data dx
994 ✗ double* x1 = (double*)malloc(m * sizeof(double));
995 ✗ assertStreamPrint(NULL, NULL != x1, "out of memory");
996 ✗ for (i = 0; i < m; ++i)
997 ✗ x1[i] = x0[i] + *lambda * dx[i];
998
999 ✗ modelica_boolean failed = TRUE;
1000 ✗ double* f_x1 = (double*)malloc(m * sizeof(double));
1001 ✗ assertStreamPrint(NULL, NULL != f_x1, "out of memory");
1002
1003 // Try
1004 #if !defined(OMC_EMCC)
1005 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1006 #endif
1007
1008 // Calculate residuals f_x1 for x1
1009 ✗ systemData->residualFunc(&resUserData, x1, f_x1, (int*)&systemData->size);
1010
1011 // Catch
1012 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { failed = FALSE; }
1013 #if !defined(OMC_EMCC)
1014 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1015 #endif
1016
1017 // Lower the dampening factor until the function call succeeds
1018 ✗ while (failed) {
1019 double d_lambda = 0.7;
1020 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Dampening factor lowered from %7.3f to %7.3f", *lambda, *lambda * d_lambda);
1021
1022 // Handle failure
1023 ✗ *lambda *= d_lambda;
1024
1025 // Update x1 based on new lambda
1026 ✗ for (i = 0; i < m; ++i)
1027 ✗ x1[i] = x0[i] + *lambda * dx[i];
1028
1029 // Retry the function call
1030 #if !defined(OMC_EMCC)
1031 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
1032 #endif
1033
1034 ✗ systemData->residualFunc(&resUserData, x1, f_x1, (int*)&systemData->size);
1035
1036 ✗ if (OMC_ERROR_RAISED()) { OMC_ERROR_CLEAR(); } else { failed = FALSE; }
1037 #if !defined(OMC_EMCC)
1038 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
1039 #endif
1040 }
1041
1042 // Count number of nonlinear functions, i.e. all functions satifying: |f(x1)| > eps
1043 ✗ *p = 0;
1044 ✗ for (i = 0; i < m; ++i)
1045 ✗ if (fabs(f_x1[i] + ((*lambda) - 1)*f_x0[i]) > eps)
1046 ✗ (*p)++;
1047
1048 // Get indices of nonlinear functions of f^i
1049 unsigned* n_idx = NULL;
1050 ✗ if (*p > 0) {
1051 ✗ n_idx = (unsigned*)malloc(*p * sizeof(unsigned));
1052 ✗ assertStreamPrint(NULL, NULL != n_idx, "out of memory");
1053 unsigned n = 0;
1054 ✗ for (i = 0; i < m; ++i)
1055 ✗ if (fabs(f_x1[i] + ((*lambda) - 1)*f_x0[i]) > eps)
1056 ✗ n_idx[n++] = i;
1057 }
1058
1059 // Free allocated memory
1060 ✗ free(x1);
1061 ✗ free(f_x1);
1062
1063 ✗ return n_idx;
1064 }
1065
1066 ✗ unsigned* getNonlinearVars( unsigned m, double*** fxx, unsigned* q)
1067 {
1068 // If at least one value in the entire column j of f_xx[k][i][j] > eps, then x[j] is a nonlinear variable
1069
1070 unsigned i, j, k;
1071 double eps = 1.e-9;
1072
1073 // Allocate memory for indicator of value of column j != 0
1074 ✗ unsigned* aValueOfColumn_gt_0 = (unsigned*)malloc(m * sizeof(unsigned));
1075 ✗ assertStreamPrint(NULL, NULL != aValueOfColumn_gt_0, "out of memory");
1076
1077 // Initialize indicator of value of column j != 0
1078 ✗ for (j = 0; j < m; j++)
1079 ✗ aValueOfColumn_gt_0[j] = 0;
1080
1081 // Retrieve indicator of value of column j != 0
1082 ✗ for (k = 0; k < m; k++)
1083 ✗ for (i = 0; i < m; i++)
1084 ✗ for (j = 0; j < m; j++)
1085 ✗ if (fabs(fxx[k][i][j]) > eps)
1086 ✗ aValueOfColumn_gt_0[j] = 1;
1087
1088 // Count number of columns where a value > 0 <==> number of nonlinear variables
1089 ✗ *q = 0;
1090 ✗ for (j = 0; j < m; j++)
1091 ✗ *q += aValueOfColumn_gt_0[j];
1092
1093 // Get indices of nonlinear variables of x[j]
1094 unsigned *w_idx = NULL;
1095 ✗ if (*q > 0) {
1096 ✗ w_idx = (unsigned*)malloc(*q * sizeof(unsigned));
1097 ✗ assertStreamPrint(NULL, NULL != w_idx, "out of memory");
1098 unsigned n = 0;
1099 ✗ for (j = 0; j < m; j++)
1100 ✗ if (aValueOfColumn_gt_0[j] == 1)
1101 ✗ w_idx[n++] = j;
1102 }
1103
1104 ✗ free(aValueOfColumn_gt_0);
1105
1106 ✗ return w_idx;
1107 }
1108
1109 ✗ unsigned* getLinearVars( unsigned m, unsigned q, unsigned *w_idx )
1110 {
1111 // Linear dependables "z": store the remaining ones (those not being in w_idx) in z_idx
1112
1113 unsigned i, j, k, i_in_w;
1114 unsigned* z_idx = NULL;
1115
1116 ✗ if (m > q) {
1117 ✗ z_idx = (unsigned*)malloc((m - q) * sizeof(unsigned));
1118 ✗ assertStreamPrint(NULL, NULL != z_idx, "out of memory");
1119 j = 0;
1120 ✗ for (i = 0; i < m; i++) {
1121 i_in_w = 0;
1122 ✗ for (k = 0; k < q; k++) {
1123 ✗ if (w_idx[k] == i) {
1124 i_in_w = 1;
1125 break;
1126 }
1127 }
1128 ✗ if (!i_in_w) {
1129 ✗ z_idx[j] = i;
1130 ✗ j++;
1131 }
1132 }
1133 }
1134 ✗ return z_idx;
1135 }
1136
1137 // --------------------------------------------------------------------------------------------------------------------------------
1138
1139 ✗ void newtonDiagnostics(DATA* data, threadData_t *threadData, int sysNumber)
1140 {
1141 // infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Newton diagnostics starting ...."); THIS IS NOT REALLY NECESSARY
1142
1143 /**** This section is not really required for the Newton diagnostics. It could be useful if only parameters and variables
1144 **** relevant for the specific system being analyzed were printed out. Otherwise, for large systems this section would be
1145 **** huge and just annoying
1146
1147 printf("\n ****** Model name: %s\n", data->modelData->modelName);
1148 printf(" ****** Initial : %d\n" , data->simulationInfo->initial);
1149
1150 printf(" ****** Number of integer parameters : %ld\n", data->modelData->nParametersInteger);
1151 for( unsigned int i = 0; i < data->modelData->nParametersInteger; ++i)
1152 printf(" ****** %2d: id=%d, name=%10s, value=%10ld\n", i+1, (data->modelData->integerParameterData[i].info.id),
1153 (data->modelData->integerParameterData[i].info.name),
1154 (data->modelData->integerParameterData[i].attribute.start));
1155
1156 printf(" ****** Number of discrete real params : %ld\n", data->modelData->nDiscreteRealArray);
1157 printf(" ****** Number of real parameters : %ld\n", data->modelData->nParametersReal);
1158 for( unsigned int i = 0; i < data->modelData->nParametersReal; ++i)
1159 printf(" ****** %2d: id=%d, name=%10s, value=%10f\n", i+1, (data->modelData->realParameterData[i].info.id),
1160 (data->modelData->realParameterData[i].info.name),
1161 (data->modelData->realParameterData[i].attribute.start));
1162
1163 printf(" ****** Number of integer variables : %ld\n", data->modelData->nVariablesInteger);
1164 for( unsigned int i = 0; i < data->modelData->nVariablesInteger; ++i)
1165 printf(" ****** %2d: id=%d, name=%10s, value=%10ld\n", i+1, (data->modelData->integerVarsData[i].info.id),
1166 (data->modelData->integerVarsData[i].info.name),
1167 (data->modelData->integerVarsData[i].attribute.start));
1168
1169 printf(" ****** Number of real variables : %ld\n", data->modelData->nVariablesReal);
1170 for( unsigned int i = 0; i < data->modelData->nVariablesReal; ++i)
1171 printf(" ****** %2d: id=%d, name=%10s, value=%10f\n", i+1, (data->modelData->realVarsData[i].info.id),
1172 (data->modelData->realVarsData[i].info.name),
1173 (data->modelData->realVarsData[i].attribute.start));
1174
1175 */
1176
1177 // --------------------------------------------------------------------------------------------------------------------------------
1178
1179 // Damping factor
1180 ✗ double lambda = 1.0;
1181
1182 // m: total number of equations f(x)
1183 // p: number of non-linear equations n(x)
1184 // q: number of variables on which non-linear equations n(x) just depend
1185
1186 ✗ NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[sysNumber]);
1187 ✗ unsigned m = systemData->size;
1188 unsigned i, j, k, p, q;
1189
1190 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Running newton diagnostics for system " OMC_INT_FORMAT, systemData->equationIndex);
1191
1192 // Store all dependents in "x0" and function values as function of x0 in f
1193 ✗ double* x0 = (double*)malloc(m * sizeof(double));
1194 ✗ assertStreamPrint(NULL, NULL != x0, "out of memory");
1195 ✗ double* f = (double*)malloc(m * sizeof(double));
1196 ✗ assertStreamPrint(NULL, NULL != f, "out of memory");
1197 ✗ for( i = 0; i < m; i++) {
1198 ✗ x0[i] = systemData->nlsx[i];
1199 ✗ f[i] = systemData->resValues[i];
1200 }
1201
1202 // Get Jacobian fx from system data
1203 ✗ double** fx = getJacobian(data, threadData, systemData);
1204
1205 // Obtain Newton steps dx = -f(x0)/fx(x0)
1206 ✗ double* dx = getFirstNewtonStep(m, f, fx);
1207
1208 // Get Hessian fxx from numerical differentiation of fx
1209 ✗ double*** fxx = getHessian(data, threadData, sysNumber, m);
1210
1211 // Obtain indices of non-linear functions "n" (p is the number of non-linear functions)
1212 ✗ unsigned* n_idx = getNonlinearEqns(data, threadData, sysNumber, m, f, x0, dx, &lambda, &p);
1213
1214 ✗ if (p == 0) {
1215 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Newton diagnostics terminated: no non-linear equations!");
1216 ✗ free(x0);
1217 ✗ free(f);
1218 ✗ free(dx);
1219 ✗ for (i = 0; i < m; i++)
1220 ✗ free(fx[i]);
1221 ✗ free(fx);
1222 ✗ for (i = 0; i < m; i++) {
1223 ✗ for (j = 0; j < m; j++)
1224 ✗ free(fxx[i][j]);
1225 ✗ free(fxx[i]);
1226 }
1227 ✗ free(fxx);
1228 ✗ free(n_idx);
1229 ✗ return;
1230 }
1231
1232 // Obtain vector "w0": initial guesses of vars where Jacobian matrix J(w) of f(x) only depends on
1233 ✗ unsigned* w_idx = getNonlinearVars( m, fxx, &q);
1234
1235 // Obtain vector "z": linear dependents
1236 ✗ unsigned* z_idx = getLinearVars( m, q, w_idx);
1237
1238 // --------------------------------------------------------------------------------------------------------------------------------
1239
1240 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Information about the system from non-linear pattern");
1241 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Total number of equations = %d", systemData->nonlinearPattern->numberOfEqns);
1242 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Number of unknowns = %d", systemData->nonlinearPattern->numberOfVars);
1243 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Number of non-linear entries = %d", systemData->nonlinearPattern->numberOfNonlinear);
1244 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
1245
1246 // --------------------------------------------------------------------------------------------------------------------------------
1247
1248 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Information about the initial guess");
1249
1250 // Prints values of unknown vector x0 - printed indeces range from 1 to m (as in mathematics, not in C)
1251 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Vector x0 of unknowns");
1252 ✗ for (i = 0; i < m; i++) {
1253 ✗ if(m < 10)
1254 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "x0[%1d] = %14.10f (%s)", i+1, x0[i],
1255 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[i]);
1256 ✗ else if(m < 100)
1257 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "x0[%2d] = %14.10f (%s)", i+1, x0[i],
1258 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[i]);
1259 ✗ else if(m < 1000)
1260 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "x0[%3d] = %14.10f (%s)", i+1, x0[i],
1261 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[i]);
1262 else if(m < 100)
1263 infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "x0[%4d] = %14.10f (%s)", i+1, x0[i],
1264 modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[i]);
1265 else
1266 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "x0[%5d] = %14.10f (%s)", i+1, x0[i],
1267 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[i]);
1268 }
1269 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
1270
1271 // Prints residual function values at x0: vector f(x0)
1272 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Residual function values of all equations f(x0)");
1273 ✗ for (i = 0; i < m; i++) {
1274 ✗ if (fabs(f[i]) > 1.e-9) {
1275 ✗ if (m < 10)
1276 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "f[%1d] = %14.10f", i+1, f[i]);
1277 ✗ else if (m < 100)
1278 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "f[%2d] = %14.10f", i+1, f[i]);
1279 ✗ else if (m < 1000)
1280 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "f[%3d] = %14.10f", i+1, f[i]);
1281 ✗ else if (m < 10000)
1282 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "f[%4d] = %14.10f", i+1, f[i]);
1283 else
1284 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "f[%5d] = %14.10f", i+1, f[i]);
1285 }
1286 }
1287 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
1288
1289 // Prints values of nonlinear unknown vector w0 - printed indeces range from 1 to m (as in mathematics, not in C)
1290 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Vector w0 of nonlinear unknowns");
1291 ✗ for (i = 0; i < q; i++) {
1292 ✗ if (m < 10)
1293 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "w0[%1d] = x0[%1d] = %14.10f (%s)", i+1, w_idx[i] + 1, x0[w_idx[i]],
1294 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[i]]);
1295 ✗ else if (m < 100)
1296 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "w0[%2d] = x0[%2d] = %14.10f (%s)", i+1, w_idx[i] + 1, x0[w_idx[i]],
1297 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[i]]);
1298 ✗ else if (m < 1000)
1299 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "w0[%3d] = x0[%3d] = %14.10f (%s)", i+1, w_idx[i] + 1, x0[w_idx[i]],
1300 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[i]]);
1301 ✗ else if (m < 10000)
1302 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "w0[%4d] = x0[%4d] = %14.10f (%s)", i+q+1, w_idx[i] + 1, x0[w_idx[i]],
1303 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[i]]);
1304 else
1305 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "w0[%5d] = x0[%5d] = %14.10f (%s)", i+q+1, w_idx[i] + 1, x0[w_idx[i]],
1306 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[w_idx[i]]);
1307 }
1308 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
1309
1310 // Prints valuse of linear unknown vector z0 - printed indeces range from 1 to m (as in mathematics, not in C)
1311 ✗ if (m > q) {
1312 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Vector z0 of nonlinear unknowns");
1313 ✗ for (i = 0; i < m-q; i++) {
1314 ✗ if (m - q < 10)
1315 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "z0[%1d] = %14.10f (%s)", i + 1, x0[z_idx[i]],
1316 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[z_idx[i]]);
1317 ✗ else if (m - q < 100)
1318 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "z0[%2d] = %14.10f (%s)", i + 1, x0[z_idx[i]],
1319 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[z_idx[i]]);
1320 ✗ else if (m - q < 1000)
1321 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "z0[%3d] = %14.10f (%s)", i + 1, x0[z_idx[i]],
1322 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[z_idx[i]]);
1323 ✗ else if (m - q < 10000)
1324 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "z0[%4d] = %14.10f (%s)", i + 1, x0[z_idx[i]],
1325 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[z_idx[i]]);
1326 else
1327 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "z0[%5d] = %14.10f (%s)", i + 1, x0[z_idx[i]],
1328 ✗ modelInfoGetEquation(&data->modelData->modelDataXml, systemData->equationIndex).vars[z_idx[i]]);
1329 }
1330 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
1331 }
1332
1333 // Prints nonlinear residual function values at x0: vector n(x0)
1334 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 1, "Residual function values of all nonlinear equations n(w0)");
1335 ✗ for (i = 0; i < p; ++i) {
1336 ✗ if (m < 10)
1337 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "n[%1d] = f[%1d] = %14.10f", i+1, n_idx[i]+1, f[n_idx[i]]);
1338 ✗ else if (m < 100)
1339 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "n[%2d] = f[%2d] = %14.10f", i+1, n_idx[i]+1, f[n_idx[i]]);
1340 ✗ else if (m < 1000)
1341 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "n[%3d] = f[%3d] = %14.10f", i+1, n_idx[i]+1, f[n_idx[i]]);
1342 ✗ else if (m < 10000)
1343 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "n[%4d] = f[%4d] = %14.10f", i+1, n_idx[i]+1, f[n_idx[i]]);
1344 else
1345 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "n[%5d] = f[%5d] = %14.10f", i+1, n_idx[i]+1, f[n_idx[i]]);
1346 }
1347 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS);
1348
1349
1350 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Final damping factor lambda = %.3g", lambda);
1351
1352 ✗ messageClose(OMC_LOG_NLS_NEWTON_DIAGNOSTICS); // End of information about the initial guess
1353
1354 // --------------------------------------------------------------------------------------------------------------------------------
1355
1356 ✗ double maxRes = maxNonLinearResiduals(m, m - q, z_idx, f, fx, dx);
1357
1358 ✗ double* alpha = calcAlpha(data, threadData, sysNumber, m, p, q, n_idx, w_idx, x0, dx, f, fxx, lambda, maxRes);
1359
1360 ✗ double*** Gamma_ijk = calcGamma(m, p, q, n_idx, w_idx, dx, fxx, maxRes);
1361
1362 ✗ double** Sigma = calcSigma(m, q, w_idx, dx, fx, fxx);
1363
1364 ✗ PrintResults(data, sysNumber, m, p, q, n_idx, w_idx, x0, alpha, Gamma_ijk, Sigma);
1365
1366 // --------------------------------------------------------------------------------------------------------------------------------
1367
1368 // Free dynamically allocated memory
1369 ✗ free(x0);
1370 ✗ free(f);
1371 ✗ free(dx);
1372
1373 ✗ for (i = 0; i < m; i++)
1374 ✗ free(fx[i]);
1375 ✗ free(fx);
1376
1377 ✗ for (i = 0; i < m; i++) {
1378 ✗ for (j = 0; j < m; j++)
1379 ✗ free(fxx[i][j]);
1380 ✗ free(fxx[i]);
1381 }
1382 ✗ free(fxx);
1383
1384 ✗ free(n_idx);
1385 ✗ free(w_idx);
1386 ✗ if (z_idx)
1387 ✗ free(z_idx);
1388
1389 ✗ free(alpha);
1390
1391 ✗ for (i = 0; i < p; i++) {
1392 ✗ for (j = 0; j < q; j++)
1393 ✗ free(Gamma_ijk[i][j]);
1394 ✗ free(Gamma_ijk[i]);
1395 }
1396 ✗ free(Gamma_ijk);
1397
1398 ✗ for (i = 0; i < q; i++)
1399 ✗ free(Sigma[i]);
1400 ✗ free(Sigma);
1401
1402 ✗ infoStreamPrint(OMC_LOG_NLS_NEWTON_DIAGNOSTICS, 0, "Newton diagnostics complete!");
1403
1404 }
1405