From 950eb4a2544eee2b1bcc70c7c5e9e9f945eb574f Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Mon, 24 Aug 2009 08:28:31 +0200 Subject: various cleaning and homogeneization --- unsupported/Eigen/src/NonLinear/hybrd.h | 25 +++------ unsupported/Eigen/src/NonLinear/hybrj.h | 56 +++++++------------- unsupported/Eigen/src/NonLinear/lmder.h | 29 +++++------ unsupported/Eigen/src/NonLinear/lmder1.h | 4 +- unsupported/Eigen/src/NonLinear/lmdif.h | 21 +++----- unsupported/Eigen/src/NonLinear/lmstr.h | 39 ++++---------- unsupported/test/NonLinear.cpp | 88 ++++++++++++++++---------------- 7 files changed, 102 insertions(+), 160 deletions(-) diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h index 987373cb6..183b2a9cb 100644 --- a/unsupported/Eigen/src/NonLinear/hybrd.h +++ b/unsupported/Eigen/src/NonLinear/hybrd.h @@ -192,19 +192,16 @@ L170: /* beginning of the inner loop. */ L180: - - /* if requested, call fcn to enable printing of iterates. */ + /* if requested, call Functor::f to enable printing of iterates. */ if (nprint <= 0) { goto L190; } iflag = 0; - if ((iter - 1) % nprint == 0) { + if ((iter - 1) % nprint == 0) iflag = Functor::debug(x, fvec); - } - if (iflag < 0) { + if (iflag < 0) goto L300; - } L190: /* determine the direction p. */ @@ -220,17 +217,15 @@ L190: /* on the first iteration, adjust the initial step bound. */ - if (iter == 1) { + if (iter == 1) delta = std::min(delta,pnorm); - } /* evaluate the function at x + p and calculate its norm. */ iflag = Functor::f(wa2, wa4); ++nfev; - if (iflag < 0) { + if (iflag < 0) goto L300; - } fnorm1 = wa4.stableNorm(); /* compute the scaled actual reduction. */ @@ -295,7 +290,7 @@ L240: x = wa2; wa2 = diag.cwise() * x; fvec = wa4; - temp = wa2.stableNorm(); + xnorm = wa2.stableNorm(); fnorm = fnorm1; ++iter; L260: @@ -376,12 +371,8 @@ L300: if (iflag < 0) { info = iflag; } - if (nprint > 0) { + if (nprint > 0) iflag = Functor::debug(x, fvec); - } return info; - - /* last card of subroutine hybrd. */ - -} /* hybrd_ */ +} diff --git a/unsupported/Eigen/src/NonLinear/hybrj.h b/unsupported/Eigen/src/NonLinear/hybrj.h index ef4d88c60..587f42099 100644 --- a/unsupported/Eigen/src/NonLinear/hybrj.h +++ b/unsupported/Eigen/src/NonLinear/hybrj.h @@ -57,8 +57,7 @@ int ei_hybrj( } if (mode == 2) for (j = 0; j < n; ++j) - if (diag[j] <= 0.) - goto L300; + if (diag[j] <= 0.) goto L300; /* evaluate the function at the starting point */ /* and calculate its norm. */ @@ -68,7 +67,7 @@ int ei_hybrj( if (iflag < 0) { goto L300; } - fnorm = fvec.stableNorm();; + fnorm = fvec.stableNorm(); /* initialize iteration counter and monitors. */ @@ -93,7 +92,7 @@ L30: /* compute the qr factorization of the jacobian. */ - ei_qrfac(n, n,fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data(), wa3.data()); + ei_qrfac(n, n, fjac.data(), fjac.rows(), false, iwa, 1, wa1.data(), wa2.data(), wa3.data()); /* on the first iteration and if mode is 1, scale according */ /* to the norms of the columns of the initial jacobian. */ @@ -117,7 +116,7 @@ L50: /* and initialize the step bound delta. */ wa3 = diag.cwise() * x; - xnorm = wa3.stableNorm();; + xnorm = wa3.stableNorm(); delta = factor * xnorm; if (delta == 0.) { delta = factor; @@ -175,14 +174,12 @@ L110: goto L170; } /* Computing MAX */ - for (j = 0; j < n; ++j) - diag[j] = std::max(diag[j], wa2[j]); + diag = diag.cwise().max(wa2); L170: /* beginning of the inner loop. */ L180: - /* if requested, call Functor::f to enable printing of iterates. */ if (nprint <= 0) { @@ -191,9 +188,8 @@ L180: iflag = 0; if ((iter - 1) % nprint == 0) iflag = Functor::debug(x, fvec, fjac); - if (iflag < 0) { + if (iflag < 0) goto L300; - } L190: /* determine the direction p. */ @@ -202,26 +198,22 @@ L190: /* store the direction p and x + p. calculate the norm of p. */ - for (j = 0; j < n; ++j) { - wa1[j] = -wa1[j]; - wa2[j] = x[j] + wa1[j]; - wa3[j] = diag[j] * wa1[j]; - } + wa1 = -wa1; + wa2 = x + wa1; + wa3 = diag.cwise() * wa1; pnorm = wa3.stableNorm(); /* on the first iteration, adjust the initial step bound. */ - if (iter == 1) { + if (iter == 1) delta = std::min(delta,pnorm); - } /* evaluate the function at x + p and calculate its norm. */ iflag = Functor::f(wa2, wa4); ++nfev; - if (iflag < 0) { + if (iflag < 0) goto L300; - } fnorm1 = wa4.stableNorm(); /* compute the scaled actual reduction. */ @@ -283,7 +275,7 @@ L240: /* successful iteration. update x, fvec, and their norms. */ - x =wa2; + x = wa2; wa2 = diag.cwise() * x; fvec = wa4; xnorm = wa2.stableNorm(); @@ -319,24 +311,19 @@ L260: info = 2; } /* Computing MAX */ - if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon() * xnorm) { + if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon() * xnorm) info = 3; - } - if (nslow2 == 5) { + if (nslow2 == 5) info = 4; - } - if (nslow1 == 10) { + if (nslow1 == 10) info = 5; - } - if (info != 0) { + if (info != 0) goto L300; - } /* criterion for recalculating jacobian. */ - if (ncfail == 2) { + if (ncfail == 2) goto L290; - } /* calculate the rank one modification to the jacobian */ /* and update qtf if necessary. */ @@ -345,10 +332,8 @@ L260: sum = wa4.dot(fjac.col(j)); wa2[j] = (sum - wa3[j]) / pnorm; wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm); - if (ratio >= Scalar(1e-4)) { + if (ratio >= Scalar(1e-4)) qtf[j] = sum; - } - /* L280: */ } /* compute the qr factorization of the updated jacobian. */ @@ -376,8 +361,5 @@ L300: if (nprint > 0) iflag = Functor::debug(x, fvec, fjac); return info; - - /* last card of subroutine hybrj. */ - -} /* hybrj_ */ +} diff --git a/unsupported/Eigen/src/NonLinear/lmder.h b/unsupported/Eigen/src/NonLinear/lmder.h index 6d7feebac..d279168aa 100644 --- a/unsupported/Eigen/src/NonLinear/lmder.h +++ b/unsupported/Eigen/src/NonLinear/lmder.h @@ -7,6 +7,7 @@ int ei_lmder( int &njev, Matrix< Scalar, Dynamic, Dynamic > &fjac, VectorXi &ipvt, + Matrix< Scalar, Dynamic, 1 > &qtf, Matrix< Scalar, Dynamic, 1 > &diag, int mode=1, Scalar factor = 100., @@ -18,11 +19,12 @@ int ei_lmder( ) { const int m = fvec.size(), n = x.size(); - Matrix< Scalar, Dynamic, 1 > qtf(n), wa1(n), wa2(n), wa3(n), wa4(m); + Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(m); ipvt.resize(n); fjac.resize(m, n); diag.resize(n); + qtf.resize(n); /* Local variables */ int i, j, l; @@ -32,12 +34,11 @@ int ei_lmder( int iflag; Scalar delta; Scalar ratio; - Scalar fnorm, gnorm, pnorm, xnorm, fnorm1, actred, dirder, prered; + Scalar fnorm, gnorm; + Scalar pnorm, xnorm, fnorm1, actred, dirder, prered; int info; - /* Function Body */ - info = 0; iflag = 0; nfev = 0; @@ -51,8 +52,7 @@ int ei_lmder( } if (mode == 2) for (j = 0; j < n; ++j) - if (diag[j] <= 0.) - goto L300; + if (diag[j] <= 0.) goto L300; /* evaluate the function at the starting point */ /* and calculate its norm. */ @@ -188,8 +188,8 @@ L170: if (mode == 2) { goto L190; } - for (j = 0; j < n; ++j) - diag[j] = std::max( diag[j], wa2[j]); + /* Computing MAX */ + diag = diag.cwise().max(wa2); L190: /* beginning of the inner loop. */ @@ -235,7 +235,6 @@ L200: wa3.fill(0.); for (j = 0; j < n; ++j) { - wa3[j] = 0.; l = ipvt[j]; temp = wa1[l]; for (i = 0; i <= j; ++i) { @@ -245,7 +244,7 @@ L200: /* L230: */ } temp1 = ei_abs2(wa3.stableNorm() / fnorm); - temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm); + temp2 = ei_abs2(ei_sqrt(par) * pnorm / fnorm); /* Computing 2nd power */ prered = temp1 + temp2 / Scalar(.5); dirder = -(temp1 + temp2); @@ -269,11 +268,10 @@ L200: if (actred < 0.) { temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred); } - if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) { + if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) temp = Scalar(.1); - } /* Computing MIN */ - delta = temp * std::min(delta, pnorm/Scalar(.1)); + delta = temp * std::min(delta, pnorm / Scalar(.1)); par /= temp; goto L260; L240: @@ -356,8 +354,5 @@ L300: iflag = Functor::debug(x, fvec, fjac); } return info; - - /* last card of subroutine lmder. */ - -} /* lmder_ */ +} diff --git a/unsupported/Eigen/src/NonLinear/lmder1.h b/unsupported/Eigen/src/NonLinear/lmder1.h index f9fb01ce5..a4f968ae2 100644 --- a/unsupported/Eigen/src/NonLinear/lmder1.h +++ b/unsupported/Eigen/src/NonLinear/lmder1.h @@ -10,7 +10,7 @@ int ei_lmder1( const int n = x.size(), m=fvec.size(); int info, nfev=0, njev=0; Matrix< Scalar, Dynamic, Dynamic > fjac(m, n); - Matrix< Scalar, Dynamic, 1> diag; + Matrix< Scalar, Dynamic, 1> diag, qtf; /* check the input parameters for errors. */ if (n <= 0 || m < n || tol < 0.) { @@ -22,7 +22,7 @@ int ei_lmder1( info = ei_lmder( x, fvec, nfev, njev, - fjac, ipvt, diag, + fjac, ipvt, qtf, diag, 1, 100., (n+1)*100, diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h index 64efce600..7a003887f 100644 --- a/unsupported/Eigen/src/NonLinear/lmdif.h +++ b/unsupported/Eigen/src/NonLinear/lmdif.h @@ -19,9 +19,7 @@ int ei_lmdif( ) { const int m = fvec.size(), n = x.size(); - Matrix< Scalar, Dynamic, 1 > - wa1(n), wa2(n), wa3(n), - wa4(m); + Matrix< Scalar, Dynamic, 1 > wa1(n), wa2(n), wa3(n), wa4(m); ipvt.resize(n); fjac.resize(m, n); @@ -53,8 +51,7 @@ int ei_lmdif( } if (mode == 2) for (j = 0; j < n; ++j) - if (diag[j] <= 0.) - goto L300; + if (diag[j] <= 0.) goto L300; /* evaluate the function at the starting point */ /* and calculate its norm. */ @@ -187,8 +184,8 @@ L170: if (mode == 2) { goto L190; } - for (j = 0; j < n; ++j) /* Computing MAX */ - diag[j] = std::max(diag[j], wa2[j]); + /* Computing MAX */ + diag = diag.cwise().max(wa2); L190: /* beginning of the inner loop. */ @@ -232,8 +229,8 @@ L200: /* compute the scaled predicted reduction and */ /* the scaled directional derivative. */ + wa3.fill(0.); for (j = 0; j < n; ++j) { - wa3[j] = 0.; l = ipvt[j]; temp = wa1[l]; for (i = 0; i <= j; ++i) { @@ -267,9 +264,8 @@ L200: if (actred < 0.) { temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred); } - if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) { + if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) temp = Scalar(.1); - } /* Computing MIN */ delta = temp * std::min(delta, pnorm / Scalar(.1)); par /= temp; @@ -354,8 +350,5 @@ L300: iflag = Functor::debug(x, fvec); } return info; - - /* last card of subroutine lmdif. */ - -} /* lmdif_ */ +} diff --git a/unsupported/Eigen/src/NonLinear/lmstr.h b/unsupported/Eigen/src/NonLinear/lmstr.h index 4dd983fbc..45593a505 100644 --- a/unsupported/Eigen/src/NonLinear/lmstr.h +++ b/unsupported/Eigen/src/NonLinear/lmstr.h @@ -54,8 +54,7 @@ int ei_lmstr( if (mode == 2) for (j = 0; j < n; ++j) - if (diag[j] <= 0.) - goto L300; + if (diag[j] <= 0.) goto L300; /* evaluate the function at the starting point */ /* and calculate its norm. */ @@ -95,23 +94,15 @@ L40: /* forming (q transpose)*fvec and storing the first */ /* n components in qtf. */ - for (j = 0; j < n; ++j) { - qtf[j] = 0.; - for (i = 0; i < n; ++i) { - fjac(i,j) = 0.; - /* L50: */ - } - /* L60: */ - } + qtf.fill(0.); + fjac.fill(0.); iflag = 2; for (i = 0; i < m; ++i) { - if (Functor::df(x, wa3, iflag) < 0) { + if (Functor::df(x, wa3, iflag) < 0) goto L340; - } temp = fvec[i]; ei_rwupdt(n, fjac.data(), fjac.rows(), wa3.data(), qtf.data(), &temp, wa1.data(), wa2.data()); ++iflag; - /* L70: */ } ++njev; @@ -126,25 +117,21 @@ L40: ipvt[j] = j; wa2[j] = fjac.col(j).start(j).stableNorm(); } - if (! sing) { + if (! sing) goto L130; - } ipvt.cwise()+=1; ei_qrfac(n, n, fjac.data(), fjac.rows(), true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data()); ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1) for (j = 0; j < n; ++j) { - if (fjac(j,j) == 0.) { + if (fjac(j,j) == 0.) goto L110; - } sum = 0.; for (i = j; i < n; ++i) { sum += fjac(i,j) * qtf[i]; - /* L90: */ } temp = -sum / fjac(j,j); for (i = j; i < n; ++i) { qtf[i] += fjac(i,j) * temp; - /* L100: */ } L110: fjac(j,j) = wa1[j]; @@ -173,10 +160,7 @@ L150: /* on the first iteration, calculate the norm of the scaled x */ /* and initialize the step bound delta. */ - for (j = 0; j < n; ++j) { - wa3[j] = diag[j] * x[j]; - /* L160: */ - } + wa3 = diag.cwise() * x; xnorm = wa3.stableNorm(); delta = factor * xnorm; if (delta == 0.) { @@ -222,8 +206,8 @@ L210: if (mode == 2) { goto L230; } - for (j = 0; j < n; ++j) /* Computing MAX */ - diag[j] = std::max(diag[j], wa2[j]); + /* Computing MAX */ + diag = diag.cwise().max(wa2); L230: /* beginning of the inner loop. */ @@ -390,8 +374,5 @@ L340: iflag = Functor::debug(x, fvec, wa3); } return info; - - /* last card of subroutine lmstr. */ - -} /* lmstr_ */ +} diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp index 354b3bf42..241714e3f 100644 --- a/unsupported/test/NonLinear.cpp +++ b/unsupported/test/NonLinear.cpp @@ -167,7 +167,7 @@ void testLmder() const int m=15, n=3; int info, nfev=0, njev=0; double fnorm, covfac, covar_ftol; - VectorXd x(n), fvec(m), diag(n); + VectorXd x(n), fvec(m), diag(n), qtf; MatrixXd fjac; VectorXi ipvt; @@ -175,7 +175,7 @@ void testLmder() x.setConstant(n, 1.); // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return values VERIFY( 1 == info); @@ -630,7 +630,7 @@ void testNistChwirut2(void) const int m=54, n=3; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -639,7 +639,7 @@ void testNistChwirut2(void) */ x<< 0.1, 0.01, 0.02; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -657,7 +657,7 @@ void testNistChwirut2(void) */ x<< 0.15, 0.008, 0.010; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 400, 1.E6*epsilon(), 1.E6*epsilon()); // check return value @@ -707,7 +707,7 @@ void testNistMisra1a(void) const int m=14, n=2; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -716,7 +716,7 @@ void testNistMisra1a(void) */ x<< 500., 0.0001; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -733,7 +733,7 @@ void testNistMisra1a(void) */ x<< 250., 0.0005; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -792,7 +792,7 @@ void testNistHahn1(void) const int m=236, n=7; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -801,7 +801,7 @@ void testNistHahn1(void) */ x<< 10., -1., .05, -.00001, -.05, .001, -.000001; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -823,7 +823,7 @@ void testNistHahn1(void) */ x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -877,7 +877,7 @@ void testNistMisra1d(void) const int m=14, n=2; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -886,7 +886,7 @@ void testNistMisra1d(void) */ x<< 500., 0.0001; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 3 == info); @@ -903,7 +903,7 @@ void testNistMisra1d(void) */ x<< 450., 0.0003; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -954,7 +954,7 @@ void testNistLanczos1(void) const int m=24, n=6; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -963,7 +963,7 @@ void testNistLanczos1(void) */ x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 2 == info); @@ -984,7 +984,7 @@ void testNistLanczos1(void) */ x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 2 == info); @@ -1039,7 +1039,7 @@ void testNistRat42(void) const int m=9, n=3; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1048,7 +1048,7 @@ void testNistRat42(void) */ x<< 100., 1., 0.1; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -1066,7 +1066,7 @@ void testNistRat42(void) */ x<< 75., 2.5, 0.07; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -1116,7 +1116,7 @@ void testNistMGH10(void) const int m=16, n=3; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1125,7 +1125,7 @@ void testNistMGH10(void) */ x<< 2., 400000., 25000.; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 2 == info); @@ -1143,7 +1143,7 @@ void testNistMGH10(void) */ x<< 0.02, 4000., 250.; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 2 == info); @@ -1191,7 +1191,7 @@ void testNistBoxBOD(void) const int m=6, n=2; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1200,7 +1200,7 @@ void testNistBoxBOD(void) */ x<< 1., 1.; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 10., 400, 1E6*epsilon(), 1E6*epsilon()); // check return value @@ -1218,7 +1218,7 @@ void testNistBoxBOD(void) */ x<< 100., 0.75; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 14000, epsilon(), epsilon()); // check return value @@ -1268,7 +1268,7 @@ void testNistMGH17(void) const int m=33, n=5; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1279,7 +1279,7 @@ void testNistMGH17(void) x<< 50., 150., -100., 1., 2.; // do the computation info = ei_lmder( - x, fvec, nfev, njev, fjac, ipvt, diag, + x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 5000, epsilon(), epsilon()); // check return value @@ -1301,7 +1301,7 @@ void testNistMGH17(void) */ x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -1356,7 +1356,7 @@ void testNistMGH09(void) const int m=11, n=4; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1365,7 +1365,7 @@ void testNistMGH09(void) */ x<< 25., 39, 41.5, 39.; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 5000); // 1, 100., 5000, epsilon(), epsilon()); @@ -1386,7 +1386,7 @@ void testNistMGH09(void) */ x<< 0.25, 0.39, 0.415, 0.39; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -1438,7 +1438,7 @@ void testNistBennett5(void) const int m=154, n=3; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1447,7 +1447,7 @@ void testNistBennett5(void) */ x<< -2000., 50., 0.8; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 5000); // check return value @@ -1465,7 +1465,7 @@ void testNistBennett5(void) */ x<< -1500., 45., 0.85; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -1523,7 +1523,7 @@ void testNistThurber(void) const int m=37, n=7; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1532,7 +1532,7 @@ void testNistThurber(void) */ x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 400, 1.E4*epsilon(), 1.E4*epsilon()); // check return value @@ -1555,7 +1555,7 @@ void testNistThurber(void) */ x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 400, 1.E4*epsilon(), 1.E4*epsilon()); // check return value @@ -1611,7 +1611,7 @@ void testNistRat43(void) const int m=15, n=4; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1620,7 +1620,7 @@ void testNistRat43(void) */ x<< 100., 10., 1., 1.; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 400, 1.E6*epsilon(), 1.E6*epsilon()); // check return value @@ -1640,7 +1640,7 @@ void testNistRat43(void) */ x<< 700., 5., 0.75, 1.3; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag, + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag, 1, 100., 400, 1.E5*epsilon(), 1.E5*epsilon()); // check return value @@ -1694,7 +1694,7 @@ void testNistEckerle4(void) const int m=35, n=3; int info, nfev=0, njev=0; - VectorXd x(n), fvec(m), diag; + VectorXd x(n), fvec(m), diag, qtf; MatrixXd fjac; VectorXi ipvt; @@ -1703,7 +1703,7 @@ void testNistEckerle4(void) */ x<< 1., 10., 500.; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); @@ -1721,7 +1721,7 @@ void testNistEckerle4(void) */ x<< 1.5, 5., 450.; // do the computation - info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, diag); + info = ei_lmder(x, fvec, nfev, njev, fjac, ipvt, qtf, diag); // check return value VERIFY( 1 == info); -- cgit v1.2.3