From de7d14b2b3eb3c9117d85ad940d1af2fbf41b54b Mon Sep 17 00:00:00 2001
From: Thomas Capricelli <orzel@freehackers.org>
Date: Thu, 20 Aug 2009 22:16:30 +0200
Subject: [PATCH] porting lmstr1 to eigen

---
 unsupported/Eigen/NonLinear                   |  2 +-
 .../Eigen/src/NonLinear/MathFunctions.h       | 24 ------
 unsupported/Eigen/src/NonLinear/lmder1.h      |  5 +-
 unsupported/Eigen/src/NonLinear/lmstr1.h      | 82 ++++++-------------
 4 files changed, 30 insertions(+), 83 deletions(-)

diff --git a/unsupported/Eigen/NonLinear b/unsupported/Eigen/NonLinear
index ac0f787e1..dbb9e8703 100644
--- a/unsupported/Eigen/NonLinear
+++ b/unsupported/Eigen/NonLinear
@@ -130,7 +130,6 @@ Scalar ei_enorm ( int n, const Scalar *x ){
 #include "src/NonLinear/lmder.h"
 #include "src/NonLinear/hybrd1.h"
 #include "src/NonLinear/hybrd.h"
-#include "src/NonLinear/lmstr1.h"
 #include "src/NonLinear/lmstr.h"
 #include "src/NonLinear/lmdif1.h"
 #include "src/NonLinear/lmdif.h"
@@ -139,6 +138,7 @@ Scalar ei_enorm ( int n, const Scalar *x ){
 #include "src/NonLinear/chkder.h"
 #include "src/NonLinear/MathFunctions.h"
 #include "src/NonLinear/lmder1.h"
+#include "src/NonLinear/lmstr1.h"
 
 //@}
 
diff --git a/unsupported/Eigen/src/NonLinear/MathFunctions.h b/unsupported/Eigen/src/NonLinear/MathFunctions.h
index ad16d77e8..d7e221ba1 100644
--- a/unsupported/Eigen/src/NonLinear/MathFunctions.h
+++ b/unsupported/Eigen/src/NonLinear/MathFunctions.h
@@ -149,30 +149,6 @@ int ei_hybrj(
     );
 }
 
-template<typename Functor, typename Scalar>
-int ei_lmstr1(
-        Matrix< Scalar, Dynamic, 1 >  &x,
-        Matrix< Scalar, Dynamic, 1 >  &fvec,
-        VectorXi &ipvt,
-        Scalar tol = ei_sqrt(epsilon<Scalar>())
-        )
-{
-    int lwa = 5*x.size()+fvec.size();
-    int ldfjac = fvec.size();
-    Matrix< Scalar, Dynamic, 1 > wa(lwa);
-    Matrix< Scalar, Dynamic, Dynamic > fjac(ldfjac, x.size());
-
-    ipvt.resize(x.size());
-    return lmstr1_template<Scalar>(
-            Functor::f, 0,
-            fvec.size(), x.size(), x.data(), fvec.data(),
-            fjac.data() , ldfjac,
-            tol,
-            ipvt.data(),
-            wa.data(), lwa
-    );
-}
-
 template<typename Functor, typename Scalar>
 int ei_lmstr(
         Matrix< Scalar, Dynamic, 1 >  &x,
diff --git a/unsupported/Eigen/src/NonLinear/lmder1.h b/unsupported/Eigen/src/NonLinear/lmder1.h
index f1322261d..b12f7595a 100644
--- a/unsupported/Eigen/src/NonLinear/lmder1.h
+++ b/unsupported/Eigen/src/NonLinear/lmder1.h
@@ -11,15 +11,14 @@ int ei_lmder1(
     int info, nfev, njev;
     Matrix< Scalar, Dynamic, Dynamic > fjac(m, n);
     Matrix< Scalar, Dynamic, 1> diag;
-    ipvt.resize(n);
 
-
-/*     check the input parameters for errors. */
+    /* check the input parameters for errors. */
     if (n <= 0 || m < n || tol < 0.) {
         printf("ei_lmder1 bad args : m,n,tol,...");
         return 0;
     }
 
+    ipvt.resize(n);
     info = ei_lmder<Functor,Scalar>(
         x, fvec,
         nfev, njev,
diff --git a/unsupported/Eigen/src/NonLinear/lmstr1.h b/unsupported/Eigen/src/NonLinear/lmstr1.h
index 58c248f11..b1e1827b3 100644
--- a/unsupported/Eigen/src/NonLinear/lmstr1.h
+++ b/unsupported/Eigen/src/NonLinear/lmstr1.h
@@ -1,61 +1,33 @@
 
-template<typename Scalar>
-int lmstr1_template(minpack_funcderstr_mn fcn, void *p, int m, int n, Scalar *x, 
-	Scalar *fvec, Scalar *fjac, int ldfjac, Scalar tol, 
-	int *ipvt, Scalar *wa, int lwa)
+template<typename Functor, typename Scalar>
+int ei_lmstr1(
+        Matrix< Scalar, Dynamic, 1 >  &x,
+        Matrix< Scalar, Dynamic, 1 >  &fvec,
+        VectorXi &ipvt,
+        Scalar tol = ei_sqrt(epsilon<Scalar>())
+        )
 {
-    /* Initialized data */
+    const int n = x.size(), m=fvec.size();
+    int info, nfev, njev;
+    Matrix< Scalar, Dynamic, Dynamic > fjac(m, n);
+    Matrix< Scalar, Dynamic, 1> diag;
 
-    const Scalar factor = 100.;
-
-    /* System generated locals */
-    int fjac_dim1, fjac_offset;
-
-    /* Local variables */
-    int mode, nfev, njev;
-    Scalar ftol, gtol, xtol;
-    int maxfev, nprint;
-    int info;
-
-    /* Parameter adjustments */
-    --fvec;
-    --ipvt;
-    --x;
-    fjac_dim1 = ldfjac;
-    fjac_offset = 1 + fjac_dim1 * 1;
-    fjac -= fjac_offset;
-    --wa;
-
-    /* Function Body */
-    info = 0;
-
-/*     check the input parameters for errors. */
-
-    if (n <= 0 || m < n || ldfjac < n || tol < 0. || lwa < n * 5 +
-	    m) {
-	/* goto L10; */
-        return info;
+    /* check the input parameters for errors. */
+    if (n <= 0 || m < n || tol < 0.) {
+        printf("ei_lmstr1 bad args : m,n,tol,...");
+        return 0;
     }
 
-/*     call lmstr. */
-
-    maxfev = (n + 1) * 100;
-    ftol = tol;
-    xtol = tol;
-    gtol = 0.;
-    mode = 1;
-    nprint = 0;
-    info = lmstr(fcn, p, m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac,
-	    ftol, xtol, gtol, maxfev, &wa[1], mode, factor, nprint, 
-	    &nfev, &njev, &ipvt[1], &wa[n + 1], &wa[(n << 1) + 1], &
-	    wa[n * 3 + 1], &wa[(n << 2) + 1], &wa[n * 5 + 1]);
-    if (info == 8) {
-	info = 4;
-    }
-/* L10: */
-    return info;
-
-/*     last card of subroutine lmstr1. */
-
-} /* lmstr1_ */
+    ipvt.resize(n);
+    info = ei_lmstr<Functor,Scalar>(
+        x, fvec,
+        nfev, njev,
+        fjac, ipvt, diag,
+        1,
+        100.,
+        (n+1)*100,
+        tol, tol, Scalar(0.)
+    );
+    return (info==8)?4:info;
+}