From f2fcbb0207dc1cd3d76c8cffacf92af36ab858ba Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Sun, 23 Aug 2009 03:54:40 +0200 Subject: [PATCH] use eigen objects for ei_fdjac*(), this is a prerequisite before porting hybrd/lmdif.. --- unsupported/Eigen/src/NonLinear/fdjac1.h | 74 ++++++++++-------------- unsupported/Eigen/src/NonLinear/fdjac2.h | 40 ++++++------- unsupported/Eigen/src/NonLinear/hybrd.h | 4 +- unsupported/Eigen/src/NonLinear/lmdif.h | 3 +- 4 files changed, 51 insertions(+), 70 deletions(-) diff --git a/unsupported/Eigen/src/NonLinear/fdjac1.h b/unsupported/Eigen/src/NonLinear/fdjac1.h index 9d5b68ce2..9388bd347 100644 --- a/unsupported/Eigen/src/NonLinear/fdjac1.h +++ b/unsupported/Eigen/src/NonLinear/fdjac1.h @@ -1,30 +1,24 @@ template -int ei_fdjac1(minpack_func_nn fcn, int n, Scalar *x, const Scalar * - fvec, Scalar *fjac, int ldfjac, int ml, - int mu, Scalar epsfcn, Scalar *wa1, Scalar *wa2) +int ei_fdjac1(minpack_func_nn fcn, + Matrix< Scalar, Dynamic, 1 > &x, + Matrix< Scalar, Dynamic, 1 > &fvec, + Matrix< Scalar, Dynamic, Dynamic > &fjac, + int ml, int mu, + Scalar epsfcn, + Matrix< Scalar, Dynamic, 1 > &wa1, + Matrix< Scalar, Dynamic, 1 > &wa2) { - /* System generated locals */ - int fjac_dim1, fjac_offset; - /* Local variables */ - Scalar h__; + Scalar h; int i, j, k; Scalar eps, temp; int msum; int iflag = 0; - /* Parameter adjustments */ - --wa2; - --wa1; - --fvec; - --x; - fjac_dim1 = ldfjac; - fjac_offset = 1 + fjac_dim1 * 1; - fjac -= fjac_offset; - /* Function Body */ const Scalar epsmch = epsilon(); + const int n = x.size(); eps = ei_sqrt((std::max(epsfcn,epsmch))); msum = ml + mu + 1; @@ -34,20 +28,18 @@ int ei_fdjac1(minpack_func_nn fcn, int n, Scalar *x, const Scalar * /* computation of dense approximate jacobian. */ - for (j = 1; j <= n; ++j) { + for (j = 0; j < n; ++j) { temp = x[j]; - h__ = eps * ei_abs(temp); - if (h__ == 0.) { - h__ = eps; - } - x[j] = temp + h__; - iflag = (*fcn)(n, &x[1], &wa1[1], 1); - if (iflag < 0) { + h = eps * ei_abs(temp); + if (h == 0.) + h = eps; + x[j] = temp + h; + iflag = (*fcn)(n, x.data(), wa1.data(), 1); + if (iflag < 0) goto L30; - } x[j] = temp; - for (i = 1; i <= n; ++i) { - fjac[i + j * fjac_dim1] = (wa1[i] - fvec[i]) / h__; + for (i = 0; i < n; ++i) { + fjac(i,j) = (wa1[i] - fvec[i]) / h; /* L10: */ } /* L20: */ @@ -59,31 +51,27 @@ L40: /* computation of banded approximate jacobian. */ - for (k = 1; k <= msum; ++k) { - for (j = k; msum< 0 ? j >= n: j <= n; j += msum) { + for (k = 0; k < msum; ++k) { + for (j = k; msum< 0 ? j > n: j < n; j += msum) { wa2[j] = x[j]; - h__ = eps * ei_abs(wa2[j]); - if (h__ == 0.) { - h__ = eps; - } - x[j] = wa2[j] + h__; + h = eps * ei_abs(wa2[j]); + if (h == 0.) h = eps; + x[j] = wa2[j] + h; /* L60: */ } - iflag = (*fcn)(n, &x[1], &wa1[1], 1); + iflag = (*fcn)(n, x.data(), wa1.data(), 1); if (iflag < 0) { /* goto L100; */ return iflag; } - for (j = k; msum< 0 ? j >= n: j <= n; j += msum) { + for (j = k; msum< 0 ? j > n: j < n; j += msum) { x[j] = wa2[j]; - h__ = eps * ei_abs(wa2[j]); - if (h__ == 0.) { - h__ = eps; - } - for (i = 1; i <= n; ++i) { - fjac[i + j * fjac_dim1] = 0.; + h = eps * ei_abs(wa2[j]); + if (h == 0.) h = eps; + for (i = 0; i < n; ++i) { + fjac(i,j) = 0.; if (i >= j - mu && i <= j + ml) { - fjac[i + j * fjac_dim1] = (wa1[i] - fvec[i]) / h__; + fjac(i,j) = (wa1[i] - fvec[i]) / h; } /* L70: */ } diff --git a/unsupported/Eigen/src/NonLinear/fdjac2.h b/unsupported/Eigen/src/NonLinear/fdjac2.h index 730bc1dfc..c8cd74b20 100644 --- a/unsupported/Eigen/src/NonLinear/fdjac2.h +++ b/unsupported/Eigen/src/NonLinear/fdjac2.h @@ -1,45 +1,39 @@ template -int ei_fdjac2(minpack_func_mn fcn, int m, int n, Scalar *x, - const Scalar *fvec, Scalar *fjac, int ldfjac, - Scalar epsfcn, Scalar *wa) +int ei_fdjac2(minpack_func_mn fcn, + Matrix< Scalar, Dynamic, 1 > &x, + Matrix< Scalar, Dynamic, 1 > &fvec, + Matrix< Scalar, Dynamic, Dynamic > &fjac, + Scalar epsfcn, + Matrix< Scalar, Dynamic, 1 > &wa) { - /* System generated locals */ - int fjac_dim1, fjac_offset; - /* Local variables */ - Scalar h__; + Scalar h; int i, j; Scalar eps, temp; int iflag; - /* Parameter adjustments */ - --wa; - --fvec; - --x; - fjac_dim1 = ldfjac; - fjac_offset = 1 + fjac_dim1 * 1; - fjac -= fjac_offset; - /* Function Body */ const Scalar epsmch = epsilon(); + const int n = x.size(); + const int m = fvec.size(); eps = ei_sqrt((std::max(epsfcn,epsmch))); - for (j = 1; j <= n; ++j) { + for (j = 0; j < n; ++j) { temp = x[j]; - h__ = eps * ei_abs(temp); - if (h__ == 0.) { - h__ = eps; + h = eps * ei_abs(temp); + if (h == 0.) { + h = eps; } - x[j] = temp + h__; - iflag = (*fcn)(m, n, &x[1], &wa[1], 1); + x[j] = temp + h; + iflag = (*fcn)(m, n, x.data(), wa.data(), 1); if (iflag < 0) { /* goto L30; */ return iflag; } x[j] = temp; - for (i = 1; i <= m; ++i) { - fjac[i + j * fjac_dim1] = (wa[i] - fvec[i]) / h__; + for (i = 0; i < m; ++i) { + fjac(i,j) = (wa[i] - fvec[i]) / h; /* L10: */ } /* L20: */ diff --git a/unsupported/Eigen/src/NonLinear/hybrd.h b/unsupported/Eigen/src/NonLinear/hybrd.h index 82554cdcf..b3cf51727 100644 --- a/unsupported/Eigen/src/NonLinear/hybrd.h +++ b/unsupported/Eigen/src/NonLinear/hybrd.h @@ -103,8 +103,8 @@ L30: /* calculate the jacobian matrix. */ - iflag = ei_fdjac1(Functor::f, n, x.data(), fvec.data(), fjac.data(), ldfjac, - nb_of_subdiagonals, nb_of_superdiagonals, epsfcn, wa1.data(), wa2.data()); + iflag = ei_fdjac1(Functor::f, x, fvec, fjac, + nb_of_subdiagonals, nb_of_superdiagonals, epsfcn, wa1, wa2); nfev += msum; if (iflag < 0) { goto L300; diff --git a/unsupported/Eigen/src/NonLinear/lmdif.h b/unsupported/Eigen/src/NonLinear/lmdif.h index 3df6effd9..ae50edd58 100644 --- a/unsupported/Eigen/src/NonLinear/lmdif.h +++ b/unsupported/Eigen/src/NonLinear/lmdif.h @@ -84,8 +84,7 @@ L30: /* calculate the jacobian matrix. */ - iflag = ei_fdjac2(Functor::f, m, n, x.data(), fvec.data(), fjac.data(), ldfjac, - epsfcn, wa4.data()); + iflag = ei_fdjac2(Functor::f, x, fvec, fjac, epsfcn, wa4); nfev += n; if (iflag < 0) { goto L300;