Hello, Thanks for your report. The problem is that on the first iteration, the lambda value jumps to -12.8, and the value of exp(12.8 * 60) is too large to represent in double precision, so infs and nans appear in the residual vector and Jacobian matrix.
For this particular problem, the value of lambda should be > 0, however there is no constraint placed to assure that, so when the first iteration jumps to -12.8, bad numerical things happen. I think the best solution for now is to adjust the model parameters a little to prevent this type of behavior, and bring TMAX closer to 1 (instead of 40). A better solution would be to implement a constraint so that lambda > 0, but this would require a lot of development effort. I will modify the test program to make it less likely to run into this problem in the future. Thanks, Patrick On 12/7/20 3:18 PM, Zhoulai Fu@Gmail wrote: > Running the attached program nlfit_bug.c triggers undefined behavior, > calculating unreasonable > results. The program "nlfit_debug.c" is adapted from the example program > provided > by GSL's doc/examples/nlfit.c. Their difference lies in only one parameter > 40.0 > versus 60.0: > > 1d0 > < > 12c11,12 > < #define TMAX (60.0) /* time variable in [0,TMAX] */ //40.0 instead of > 60.0 in nlfit.c > --- >> //#define TMAX (40.0) /* time variable in [0,TMAX] */ >> #define TMAX (60.0) /* time variable in [0,TMAX] */ > > The undefined behavior can be better visualized if GCC's Undefined > Behavior Sanitizer is enabled (-fsanitize=undefined). Below is the > execution results: > > data: 59.3939 1.03153 0.101317 > data: 60 1.14181 0.101239 > iter 0: A = 1.0000, lambda = 1.0000, b = 0.0000, cond(J) = inf, > |f(x)| = 101.0200 > iter 1: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = 92.8216, > |f(x)| = -nan > iter 2: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 3: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 4: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > > *nielsen.c:98:7: runtime error: left shift of 4611686018427387904 by 1 > places cannot be represented in type 'long int'*iter 5: A = 3.5110, lambda > = -12.8820, b = 1.2364, cond(J) = nan, |f(x)| = -nan > iter 6: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 7: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 8: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 9: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 10: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 11: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 12: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 13: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 14: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 15: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 16: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 17: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 18: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 19: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 20: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 21: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 22: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 23: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 24: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 25: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 26: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 27: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 28: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 29: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 30: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 31: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 32: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 33: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 34: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 35: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 36: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 37: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 38: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 39: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 40: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 41: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 42: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 43: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 44: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 45: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 46: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 47: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 48: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 49: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 50: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 51: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 52: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 53: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 54: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 55: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 56: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 57: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 58: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 59: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 60: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 61: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 62: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 63: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 64: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 65: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 66: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 67: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 68: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 69: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 70: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 71: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 72: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 73: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 74: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 75: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 76: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 77: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 78: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 79: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 80: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 81: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 82: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 83: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 84: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 85: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 86: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 87: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 88: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 89: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 90: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 91: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 92: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 93: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 94: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 95: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 96: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 97: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 98: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 99: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > iter 100: A = 3.5110, lambda = -12.8820, b = 1.2364, cond(J) = nan, > |f(x)| = -nan > summary from method 'trust-region/levenberg-marquardt' > number of iterations: 100 > function evaluations: 1586 > Jacobian evaluations: 2 > reason for stopping: small gradient > initial |f(x)| = 101.020000 > final |f(x)| = inf > chisq/dof = inf > A = 3.51103 +/- nan > lambda = -12.88201 +/- nan > b = 1.23642 +/- nan > status = exceeded max number of iterations > > Preliminary debugging: The undefined behavior occurs in > multifit_nlinear/nielsen.c, line 98, triggered by overflow of a left shift > of a long int. > > *nu <<= 1; > > BR, > Zhoulai
