From bb6ffafdb9eaed3ed4d398a510c52e026f19e434 Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Sat, 22 Aug 2009 07:27:17 +0200 Subject: [PATCH] keep on cleaning f2c mess --- unsupported/Eigen/src/NonLinear/dogleg.h | 60 ++++++-------------- unsupported/Eigen/src/NonLinear/lmpar.h | 70 +++++++++--------------- 2 files changed, 43 insertions(+), 87 deletions(-) diff --git a/unsupported/Eigen/src/NonLinear/dogleg.h b/unsupported/Eigen/src/NonLinear/dogleg.h index 0909b9621..3ca406c6a 100644 --- a/unsupported/Eigen/src/NonLinear/dogleg.h +++ b/unsupported/Eigen/src/NonLinear/dogleg.h @@ -4,12 +4,8 @@ void ei_dogleg(int n, const Scalar *r__, int /* lr*/ , const Scalar *diag, const Scalar *qtb, Scalar delta, Scalar *x, Scalar *wa1, Scalar *wa2) { - /* System generated locals */ - int i__1, i__2; - Scalar d__1, d__2, d__3, d__4; - /* Local variables */ - int i__, j, k, l, jj, jp1; + int i, j, k, l, jj, jp1; Scalar sum, temp, alpha, bnorm; Scalar gnorm, qnorm, epsmch; Scalar sgnorm; @@ -31,8 +27,7 @@ void ei_dogleg(int n, const Scalar *r__, int /* lr*/ , /* first, calculate the gauss-newton direction. */ jj = n * (n + 1) / 2 + 1; - i__1 = n; - for (k = 1; k <= i__1; ++k) { + for (k = 1; k <= n; ++k) { j = n - k + 1; jp1 = j + 1; jj -= k; @@ -41,9 +36,8 @@ void ei_dogleg(int n, const Scalar *r__, int /* lr*/ , if (n < jp1) { goto L20; } - i__2 = n; - for (i__ = jp1; i__ <= i__2; ++i__) { - sum += r__[l] * x[i__]; + for (i = jp1; i <= n; ++i) { + sum += r__[l] * x[i]; ++l; /* L10: */ } @@ -53,12 +47,10 @@ L20: goto L40; } l = j; - i__2 = j; - for (i__ = 1; i__ <= i__2; ++i__) { + for (i = 1; i <= j; ++i) { /* Computing MAX */ - d__2 = temp, d__3 = fabs(r__[l]); - temp = std::max(d__2,d__3); - l = l + n - i__; + temp = std::max(temp,ei_abs(r__[l])); + l = l + n - i; /* L30: */ } temp = epsmch * temp; @@ -72,8 +64,7 @@ L40: /* test whether the gauss-newton direction is acceptable. */ - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { wa1[j] = 0.; wa2[j] = diag[j] * x[j]; /* L60: */ @@ -88,12 +79,10 @@ L40: /* next, calculate the scaled gradient direction. */ l = 1; - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { temp = qtb[j]; - i__2 = n; - for (i__ = j; i__ <= i__2; ++i__) { - wa1[i__] += r__[l] * temp; + for (i = j; i <= n; ++i) { + wa1[i] += r__[l] * temp; ++l; /* L70: */ } @@ -114,18 +103,15 @@ L40: /* calculate the point along the scaled gradient */ /* at which the quadratic is minimized. */ - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { wa1[j] = wa1[j] / gnorm / diag[j]; /* L90: */ } l = 1; - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { sum = 0.; - i__2 = n; - for (i__ = j; i__ <= i__2; ++i__) { - sum += r__[l] * wa1[i__]; + for (i = j; i <= n; ++i) { + sum += r__[l] * wa1[i]; ++l; /* L100: */ } @@ -149,26 +135,16 @@ L40: bnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&qtb[1],n).stableNorm(); temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta); /* Computing 2nd power */ - d__1 = sgnorm / delta; + temp = temp - delta / qnorm * ei_abs2(sgnorm / delta) + ei_sqrt(ei_abs2(temp - delta / qnorm) + (1.-ei_abs2(delta / qnorm)) * (1.-ei_abs2(sgnorm / delta))); /* Computing 2nd power */ - d__2 = temp - delta / qnorm; -/* Computing 2nd power */ - d__3 = delta / qnorm; -/* Computing 2nd power */ - d__4 = sgnorm / delta; - temp = temp - delta / qnorm * (d__1 * d__1) + sqrt(d__2 * d__2 + (1. - - d__3 * d__3) * (1. - d__4 * d__4)); -/* Computing 2nd power */ - d__1 = sgnorm / delta; - alpha = delta / qnorm * (1. - d__1 * d__1) / temp; + alpha = delta / qnorm * (1. - ei_abs2(sgnorm / delta)) / temp; L120: /* form appropriate convex combination of the gauss-newton */ /* direction and the scaled gradient direction. */ temp = (1. - alpha) * std::min(sgnorm,delta); - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { x[j] = temp * wa1[j] + alpha * x[j]; /* L130: */ } diff --git a/unsupported/Eigen/src/NonLinear/lmpar.h b/unsupported/Eigen/src/NonLinear/lmpar.h index 2d75b3b15..8fe5dcff6 100644 --- a/unsupported/Eigen/src/NonLinear/lmpar.h +++ b/unsupported/Eigen/src/NonLinear/lmpar.h @@ -5,17 +5,12 @@ void ei_lmpar(int n, Scalar *r__, int ldr, Scalar *par, Scalar *x, Scalar *sdiag, Scalar *wa1, Scalar *wa2) { - /* Initialized data */ - -#define p1 .1 -#define p001 .001 - /* System generated locals */ - int r_dim1, r_offset, i__1, i__2; + int r_dim1, r_offset; Scalar d__1, d__2; /* Local variables */ - int i__, j, k, l; + int i, j, k, l; Scalar fp; int jm1, jp1; Scalar sum, parc, parl; @@ -47,8 +42,7 @@ void ei_lmpar(int n, Scalar *r__, int ldr, /* jacobian is rank-deficient, obtain a least squares solution. */ nsing = n; - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { wa1[j] = qtb[j]; if (r__[j + j * r_dim1] == 0. && nsing == n) { nsing = j - 1; @@ -61,8 +55,7 @@ void ei_lmpar(int n, Scalar *r__, int ldr, if (nsing < 1) { goto L50; } - i__1 = nsing; - for (k = 1; k <= i__1; ++k) { + for (k = 1; k <= nsing; ++k) { j = nsing - k + 1; wa1[j] /= r__[j + j * r_dim1]; temp = wa1[j]; @@ -70,9 +63,8 @@ void ei_lmpar(int n, Scalar *r__, int ldr, if (jm1 < 1) { goto L30; } - i__2 = jm1; - for (i__ = 1; i__ <= i__2; ++i__) { - wa1[i__] -= r__[i__ + j * r_dim1] * temp; + for (i = 1; i <= jm1; ++i) { + wa1[i] -= r__[i + j * r_dim1] * temp; /* L20: */ } L30: @@ -80,8 +72,7 @@ L30: ; } L50: - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { l = ipvt[j]; x[l] = wa1[j]; /* L60: */ @@ -92,14 +83,13 @@ L50: /* for acceptance of the gauss-newton direction. */ iter = 0; - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { wa2[j] = diag[j] * x[j]; /* L70: */ } dxnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&wa2[1],n).blueNorm(); fp = dxnorm - delta; - if (fp <= p1 * delta) { + if (fp <= Scalar(0.1) * delta) { goto L220; } @@ -111,22 +101,19 @@ L50: if (nsing < n) { goto L120; } - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { l = ipvt[j]; wa1[j] = diag[l] * (wa2[l] / dxnorm); /* L80: */ } - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { sum = 0.; jm1 = j - 1; if (jm1 < 1) { goto L100; } - i__2 = jm1; - for (i__ = 1; i__ <= i__2; ++i__) { - sum += r__[i__ + j * r_dim1] * wa1[i__]; + for (i = 1; i <= jm1; ++i) { + sum += r__[i + j * r_dim1] * wa1[i]; /* L90: */ } L100: @@ -139,12 +126,10 @@ L120: /* calculate an upper bound, paru, for the zero of the function. */ - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { sum = 0.; - i__2 = j; - for (i__ = 1; i__ <= i__2; ++i__) { - sum += r__[i__ + j * r_dim1] * qtb[i__]; + for (i = 1; i <= j; ++i) { + sum += r__[i + j * r_dim1] * qtb[i]; /* L130: */ } l = ipvt[j]; @@ -154,7 +139,7 @@ L120: gnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&wa1[1],n).stableNorm(); paru = gnorm / delta; if (paru == 0.) { - paru = dwarf / std::min(delta,p1); + paru = dwarf / std::min(delta,Scalar(0.1)); } /* if the input par lies outside of the interval (parl,paru), */ @@ -175,18 +160,16 @@ L150: if (*par == 0.) { /* Computing MAX */ - d__1 = dwarf, d__2 = p001 * paru; + d__1 = dwarf, d__2 = Scalar(.001) * paru; *par = std::max(d__1,d__2); } temp = ei_sqrt(*par); - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { wa1[j] = temp * diag[j]; /* L160: */ } ei_qrsolv(n, &r__[r_offset], ldr, &ipvt[1], &wa1[1], &qtb[1], &x[1], &sdiag[1], &wa2[1]); - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { wa2[j] = diag[j] * x[j]; /* L170: */ } @@ -198,30 +181,27 @@ L150: /* of par. also test for the exceptional cases where parl */ /* is zero or the number of iterations has reached 10. */ - if (ei_abs(fp) <= p1 * delta || (parl == 0. && fp <= temp && temp < 0.) || + if (ei_abs(fp) <= Scalar(0.1) * delta || (parl == 0. && fp <= temp && temp < 0.) || iter == 10) { goto L220; } /* compute the newton correction. */ - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { l = ipvt[j]; wa1[j] = diag[l] * (wa2[l] / dxnorm); /* L180: */ } - i__1 = n; - for (j = 1; j <= i__1; ++j) { + for (j = 1; j <= n; ++j) { wa1[j] /= sdiag[j]; temp = wa1[j]; jp1 = j + 1; if (n < jp1) { goto L200; } - i__2 = n; - for (i__ = jp1; i__ <= i__2; ++i__) { - wa1[i__] -= r__[i__ + j * r_dim1] * temp; + for (i = jp1; i <= n; ++i) { + wa1[i] -= r__[i + j * r_dim1] * temp; /* L190: */ } L200: