Skip to content

Commit 9059292

Browse files
rocketmarkclaude
andcommitted
mp_fdjac2: guard non-finite numerical derivatives
The only protection against a non-finite finite-difference derivative here was assert(isfinite(fjac[ij])), which a release build (-DNDEBUG) strips out entirely, and the two-sided derivative path had no check at all. If the user function returns a degenerate residual for some parameter value, wa[i]/fvec[i] can be non-finite, writing NaN/Inf straight into fjac. That NaN then poisons fnorm/ratio in mpfit()'s outer Levenberg-Marquardt loop. Since any comparison against NaN is false, the iteration counter (which only advances inside the ratio>=p0001 success branch) never increments, so the maxiter check that's supposed to bound the loop never fires. The result is an unbounded CPU spin instead of a clean failure return -- in our case, confirmed via gdb on a real reproduction: ~465k spin iterations/sec with zero progress, instead of returning MP_MAXITER. Fix: sanitize fjac to 0 on a non-finite derivative in all three non-debug/debug, one-sided/two-sided paths, so the column is treated as having no measurable gradient -- the same fallback mp_qrfac already uses for a zero-norm column. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
1 parent d230bf7 commit 9059292

1 file changed

Lines changed: 6 additions & 2 deletions

File tree

redist/mpfit/mpfit.c

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1228,14 +1228,16 @@ static int mp_fdjac2(mp_func funct, int m, int n, int *ifree, int npar, FLT *x,
12281228
/* Non-debug path for speed */
12291229
for (i = 0; i < m; i++, ij++) {
12301230
fjac[ij] = (wa[i] - fvec[i]) / h; /* fjac[i+m*j] */
1231-
assert(isfinite(fjac[ij]));
1231+
if (!isfinite(fjac[ij]))
1232+
fjac[ij] = 0;
12321233
}
12331234
} else {
12341235
/* Debug path for correctness */
12351236
for (i = 0; i < m; i++, ij++) {
12361237
FLT fjold = fjac[ij];
12371238
fjac[ij] = (wa[i] - fvec[i]) / h; /* fjac[i+m*j] */
1238-
assert(isfinite(fjac[ij]));
1239+
if (!isfinite(fjac[ij]))
1240+
fjac[ij] = 0;
12391241
if ((da == 0 && dr == 0 && (fjold != 0 || fjac[ij] != 0)) ||
12401242
((da != 0 || dr != 0) && (fabs(fjold - fjac[ij]) > da + fabs(fjold) * dr))) {
12411243
fprintf(stderr, " %10d %10.4g %10.4g %10.4g %10.4g %10.4g\n", i, fvec[i], fjold, fjac[ij],
@@ -1264,6 +1266,8 @@ static int mp_fdjac2(mp_func funct, int m, int n, int *ifree, int npar, FLT *x,
12641266
/* Non-debug path for speed */
12651267
for (i = 0; i < m; i++, ij++) {
12661268
fjac[ij] = (wa2[ij] - wa[i]) / (2 * h); /* fjac[i+m*j] */
1269+
if (!isfinite(fjac[ij]))
1270+
fjac[ij] = 0;
12671271
}
12681272
} else {
12691273
/* Debug path for correctness */

0 commit comments

Comments
 (0)