diff --git a/unsupported/Eigen/NonLinear b/unsupported/Eigen/NonLinear index 5606b7b42..a01d705cd 100644 --- a/unsupported/Eigen/NonLinear +++ b/unsupported/Eigen/NonLinear @@ -116,9 +116,16 @@ int lmstr ( minpack_funcderstr_mn fcn, void *p, int m, void chkder ( int m, int n, const double *x, double *fvec, double *fjac, int ldfjac, double *xp, double *fvecp, int mode, double *err ); + +void covar(int n, double *r__, int ldr, const int *ipvt, double tol, double *wa); #endif +template +Scalar ei_enorm ( int n, const Scalar *x ){ + return Map< Matrix< Scalar, Dynamic, 1 > >(x,n).stableNorm(); +// return Map< Matrix< Scalar, Dynamic, 1 > >(x,n).blueNorm(); +} #include "src/NonLinear/lmder1.h" #include "src/NonLinear/lmder.h" diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h index 53019a003..6a28b63e4 100644 --- a/unsupported/Eigen/src/NonLinear/hybrd.h +++ b/unsupported/Eigen/src/NonLinear/hybrd.h @@ -77,7 +77,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(n, &fvec[1]); + fnorm = ei_enorm(n, &fvec[1]); /* determine the number of calls to fcn needed to compute */ /* the jacobian matrix. */ @@ -140,7 +140,7 @@ L50: wa3[j] = diag[j] * x[j]; /* L60: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -250,7 +250,7 @@ L190: wa3[j] = diag[j] * wa1[j]; /* L200: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -265,7 +265,7 @@ L190: if (iflag < 0) { goto L300; } - fnorm1 = enorm(n, &wa4[1]); + fnorm1 = ei_enorm(n, &wa4[1]); /* compute the scaled actual reduction. */ @@ -291,7 +291,7 @@ L190: wa3[i__] = qtf[i__] + sum; /* L220: */ } - temp = enorm(n, &wa3[1]); + temp = ei_enorm(n, &wa3[1]); prered = 0.; if (temp < fnorm) { /* Computing 2nd power */ @@ -344,7 +344,7 @@ L240: fvec[j] = wa4[j]; /* L250: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L260: diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h index 30712e5db..013757ae7 100644 --- a/unsupported/Eigen/src/NonLinear/hybrj.h +++ b/unsupported/Eigen/src/NonLinear/hybrj.h @@ -78,7 +78,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(n, &fvec[1]); + fnorm = ei_enorm(n, &fvec[1]); /* initialize iteration counter and monitors. */ @@ -133,7 +133,7 @@ L50: wa3[j] = diag[j] * x[j]; /* L60: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -243,7 +243,7 @@ L190: wa3[j] = diag[j] * wa1[j]; /* L200: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -258,7 +258,7 @@ L190: if (iflag < 0) { goto L300; } - fnorm1 = enorm(n, &wa4[1]); + fnorm1 = ei_enorm(n, &wa4[1]); /* compute the scaled actual reduction. */ @@ -284,7 +284,7 @@ L190: wa3[i__] = qtf[i__] + sum; /* L220: */ } - temp = enorm(n, &wa3[1]); + temp = ei_enorm(n, &wa3[1]); prered = 0.; if (temp < fnorm) { /* Computing 2nd power */ @@ -337,7 +337,7 @@ L240: fvec[j] = wa4[j]; /* L250: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L260: diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h index 6df4dcbb2..99b912f4b 100644 --- a/unsupported/Eigen/src/NonLinear/lmder.h +++ b/unsupported/Eigen/src/NonLinear/lmder.h @@ -71,7 +71,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(m, &fvec[1]); + fnorm = ei_enorm(m, &fvec[1]); /* initialize levenberg-marquardt parameter and iteration counter. */ @@ -136,7 +136,7 @@ L60: wa3[j] = diag[j] * x[j]; /* L70: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -242,7 +242,7 @@ L200: wa3[j] = diag[j] * wa1[j]; /* L210: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -257,7 +257,7 @@ L200: if (iflag < 0) { goto L300; } - fnorm1 = enorm(m, &wa4[1]); + fnorm1 = ei_enorm(m, &wa4[1]); /* compute the scaled actual reduction. */ @@ -283,7 +283,7 @@ L200: } /* L230: */ } - temp1 = enorm(n, &wa3[1]) / fnorm; + temp1 = ei_enorm(n, &wa3[1]) / fnorm; temp2 = sqrt(par) * pnorm / fnorm; /* Computing 2nd power */ d__1 = temp1; @@ -351,7 +351,7 @@ L260: fvec[i__] = wa4[i__]; /* L280: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L290: diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h index 136451a09..da04caacf 100644 --- a/unsupported/Eigen/src/NonLinear/lmdif.h +++ b/unsupported/Eigen/src/NonLinear/lmdif.h @@ -73,7 +73,7 @@ L20: if (iflag < 0) { goto L300; } - fnorm = enorm(m, &fvec[1]); + fnorm = ei_enorm(m, &fvec[1]); /* initialize levenberg-marquardt parameter and iteration counter. */ @@ -139,7 +139,7 @@ L60: wa3[j] = diag[j] * x[j]; /* L70: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -245,7 +245,7 @@ L200: wa3[j] = diag[j] * wa1[j]; /* L210: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -260,7 +260,7 @@ L200: if (iflag < 0) { goto L300; } - fnorm1 = enorm(m, &wa4[1]); + fnorm1 = ei_enorm(m, &wa4[1]); /* compute the scaled actual reduction. */ @@ -286,7 +286,7 @@ L200: } /* L230: */ } - temp1 = enorm(n, &wa3[1]) / fnorm; + temp1 = ei_enorm(n, &wa3[1]) / fnorm; temp2 = sqrt(par) * pnorm / fnorm; /* Computing 2nd power */ d__1 = temp1; @@ -354,7 +354,7 @@ L260: fvec[i__] = wa4[i__]; /* L280: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L290: diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h index 17aba6865..015e0bd13 100644 --- a/unsupported/Eigen/src/NonLinear/lmstr.h +++ b/unsupported/Eigen/src/NonLinear/lmstr.h @@ -72,7 +72,7 @@ L20: if (iflag < 0) { goto L340; } - fnorm = enorm(m, &fvec[1]); + fnorm = ei_enorm(m, &fvec[1]); /* initialize levenberg-marquardt parameter and iteration counter. */ @@ -136,7 +136,7 @@ L40: sing = TRUE_; } ipvt[j] = j; - wa2[j] = enorm(j, &fjac[j * fjac_dim1 + 1]); + wa2[j] = ei_enorm(j, &fjac[j * fjac_dim1 + 1]); /* L80: */ } if (! sing) { @@ -194,7 +194,7 @@ L150: wa3[j] = diag[j] * x[j]; /* L160: */ } - xnorm = enorm(n, &wa3[1]); + xnorm = ei_enorm(n, &wa3[1]); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -269,7 +269,7 @@ L240: wa3[j] = diag[j] * wa1[j]; /* L250: */ } - pnorm = enorm(n, &wa3[1]); + pnorm = ei_enorm(n, &wa3[1]); /* on the first iteration, adjust the initial step bound. */ @@ -284,7 +284,7 @@ L240: if (iflag < 0) { goto L340; } - fnorm1 = enorm(m, &wa4[1]); + fnorm1 = ei_enorm(m, &wa4[1]); /* compute the scaled actual reduction. */ @@ -310,7 +310,7 @@ L240: } /* L270: */ } - temp1 = enorm(n, &wa3[1]) / fnorm; + temp1 = ei_enorm(n, &wa3[1]) / fnorm; temp2 = sqrt(par) * pnorm / fnorm; /* Computing 2nd power */ d__1 = temp1; @@ -378,7 +378,7 @@ L300: fvec[i__] = wa4[i__]; /* L320: */ } - xnorm = enorm(n, &wa2[1]); + xnorm = ei_enorm(n, &wa2[1]); fnorm = fnorm1; ++iter; L330: diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 00ac38795..ff94198c2 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -1142,7 +1142,7 @@ void testNistLanczos1(void) VERIFY( 79 == nfev); VERIFY( 72 == njev); // check norm^2 - VERIFY_IS_APPROX(fvec.squaredNorm(), 1.4292388868910E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats + VERIFY_IS_APPROX(fvec.squaredNorm(), 1.429604433690E-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats // check x VERIFY_IS_APPROX(x[0], 9.5100000027E-02 ); VERIFY_IS_APPROX(x[1], 1.0000000001E+00 ); @@ -1294,7 +1294,7 @@ void testNistMGH10(void) info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); // check return value - VERIFY( 3 == info); + VERIFY( 2 == info); VERIFY( 285 == nfev); VERIFY( 250 == njev); // check norm^2 @@ -1312,7 +1312,7 @@ void testNistMGH10(void) info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); // check return value - VERIFY( 3 == info); + VERIFY( 2 == info); VERIFY( 126 == nfev); VERIFY( 116 == njev); // check norm^2 @@ -1388,8 +1388,8 @@ void testNistBoxBOD(void) // check return value VERIFY( 1 == info); - VERIFY( 16 == nfev); - VERIFY( 15 == njev); + VERIFY( 15 == nfev); + VERIFY( 14 == njev); // check norm^2 VERIFY_IS_APPROX(fvec.squaredNorm(), 1.1680088766E+03); // check x @@ -1531,15 +1531,15 @@ void testNistMGH09(void) // check return value VERIFY( 1 == info); - VERIFY( 487== nfev); - VERIFY( 378 == njev); + VERIFY( 503== nfev); + VERIFY( 385 == njev); // check norm^2 VERIFY_IS_APPROX(fvec.squaredNorm(), 3.0750560385E-04); // check x - VERIFY_IS_APPROX(x[0], 0.19280590); // should be 1.9280693458E-01 - VERIFY_IS_APPROX(x[1], 0.19130543); // should be 1.9128232873E-01 - VERIFY_IS_APPROX(x[2], 0.12306085); // should be 1.2305650693E-01 - VERIFY_IS_APPROX(x[3], 0.13607303); // should be 1.3606233068E-01 + VERIFY_IS_APPROX(x[0], 0.19280624); // should be 1.9280693458E-01 + VERIFY_IS_APPROX(x[1], 0.19129774); // should be 1.9128232873E-01 + VERIFY_IS_APPROX(x[2], 0.12305940); // should be 1.2305650693E-01 + VERIFY_IS_APPROX(x[3], 0.13606946); // should be 1.3606233068E-01 /* * Second try @@ -1924,3 +1924,13 @@ void test_NonLinear() CALL_SUBTEST(testNistEckerle4()); } +/* + * Can be useful for debugging... + printf("info, nfev, jev : %d, %d, %d\n", info, nfev, njev); + printf("x[0] : %.32g\n", x[0]); + printf("x[1] : %.32g\n", x[1]); + printf("x[2] : %.32g\n", x[2]); + printf("x[3] : %.32g\n", x[3]); + printf("fvec.squaredNorm() : %.32g\n", fvec.squaredNorm()); +*/ +