merge

9 files changed, 90 insertions, 164 deletions

diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
index ec25106f5..4b7331035 100644
--- a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
+++ b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
@@ -563,6 +563,8 @@ template<typename DerType> struct NumTraits<AutoDiffScalar<DerType> >
     AddCost = 1,
     MulCost = 1
   };
+  inline static Real epsilon() { return std::numeric_limits<Real>::epsilon(); }
+  inline static Real dummy_precision() { return NumTraits<Real>::dummy_precision(); }
 };
 
 }
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
index 39c23cdc5..2c761e648 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
@@ -313,7 +313,7 @@ template<typename Derived> struct MatrixExponentialReturnValue
     inline void evalTo(ResultType& result) const
     {
       const typename ei_eval<Derived>::type srcEvaluated = m_src.eval();
-      MatrixExponential<typename Derived::PlainMatrixType> me(srcEvaluated);
+      MatrixExponential<typename Derived::PlainObject> me(srcEvaluated);
       me.compute(result);
     }
 
@@ -327,7 +327,7 @@ template<typename Derived> struct MatrixExponentialReturnValue
 template<typename Derived>
 struct ei_traits<MatrixExponentialReturnValue<Derived> >
 {
-  typedef typename Derived::PlainMatrixType ReturnMatrixType;
+  typedef typename Derived::PlainObject ReturnType;
 };
 
 /** \ingroup MatrixFunctions_Module
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
index d63bcbce9..e8069154c 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
@@ -178,9 +178,9 @@ class MatrixFunction<MatrixType, 1>
       *
       * This is morally a \c static \c const \c Scalar, but only
       * integers can be static constant class members in C++. The
-      * separation constant is set to 0.01, a value taken from the
+      * separation constant is set to 0.1, a value taken from the
       * paper by Davies and Higham. */
-    static const RealScalar separation() { return static_cast<RealScalar>(0.01); }
+    static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
 };
 
 /** \brief Constructor. 
@@ -492,14 +492,12 @@ typename MatrixFunction<MatrixType,1>::DynMatrixType MatrixFunction<MatrixType,1
 template<typename Derived> class MatrixFunctionReturnValue
 : public ReturnByValue<MatrixFunctionReturnValue<Derived> >
 {
-  private:
+  public:
 
     typedef typename ei_traits<Derived>::Scalar Scalar;
     typedef typename ei_stem_function<Scalar>::type StemFunction;
 
-  public:
-
-    /** \brief Constructor.
+   /** \brief Constructor.
       *
       * \param[in] A  %Matrix (expression) forming the argument of the
       * matrix function.
@@ -516,7 +514,7 @@ template<typename Derived> class MatrixFunctionReturnValue
     inline void evalTo(ResultType& result) const
     {
       const typename ei_eval<Derived>::type Aevaluated = m_A.eval();
-      MatrixFunction<typename Derived::PlainMatrixType> mf(Aevaluated, m_f);
+      MatrixFunction<typename Derived::PlainObject> mf(Aevaluated, m_f);
       mf.compute(result);
     }
 
@@ -531,7 +529,7 @@ template<typename Derived> class MatrixFunctionReturnValue
 template<typename Derived>
 struct ei_traits<MatrixFunctionReturnValue<Derived> >
 {
-  typedef typename Derived::PlainMatrixType ReturnMatrixType;
+  typedef typename Derived::PlainObject ReturnType;
 };
 
 
diff --git a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
index 0754a426b..35dc332e0 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
@@ -57,7 +57,7 @@ class HybridNonLinearSolver
 {
 public:
     HybridNonLinearSolver(FunctorType &_functor)
-        : functor(_functor) { nfev=njev=iter = 0;  fnorm= 0.; }
+        : functor(_functor) { nfev=njev=iter = 0;  fnorm= 0.; useExternalScaling=false;}
 
     struct Parameters {
         Parameters()
@@ -84,36 +84,18 @@ public:
             const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
             );
 
-    HybridNonLinearSolverSpace::Status solveInit(
-            FVectorType  &x,
-            const int mode=1
-            );
-    HybridNonLinearSolverSpace::Status solveOneStep(
-            FVectorType  &x,
-            const int mode=1
-            );
-    HybridNonLinearSolverSpace::Status solve(
-            FVectorType  &x,
-            const int mode=1
-            );
+    HybridNonLinearSolverSpace::Status solveInit(FVectorType  &x);
+    HybridNonLinearSolverSpace::Status solveOneStep(FVectorType  &x);
+    HybridNonLinearSolverSpace::Status solve(FVectorType  &x);
 
     HybridNonLinearSolverSpace::Status hybrd1(
             FVectorType  &x,
             const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
             );
 
-    HybridNonLinearSolverSpace::Status solveNumericalDiffInit(
-            FVectorType  &x,
-            const int mode=1
-            );
-    HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(
-            FVectorType  &x,
-            const int mode=1
-            );
-    HybridNonLinearSolverSpace::Status solveNumericalDiff(
-            FVectorType  &x,
-            const int mode=1
-            );
+    HybridNonLinearSolverSpace::Status solveNumericalDiffInit(FVectorType  &x);
+    HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(FVectorType  &x);
+    HybridNonLinearSolverSpace::Status solveNumericalDiff(FVectorType  &x);
 
     void resetParameters(void) { parameters = Parameters(); }
     Parameters parameters;
@@ -124,6 +106,7 @@ public:
     int njev;
     int iter;
     Scalar fnorm;
+    bool useExternalScaling; 
 private:
     FunctorType &functor;
     int n;
@@ -160,18 +143,13 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrj1(
     parameters.maxfev = 100*(n+1);
     parameters.xtol = tol;
     diag.setConstant(n, 1.);
-    return solve(
-        x,
-        2
-    );
+    useExternalScaling = true;
+    return solve(x);
 }
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveInit(FVectorType  &x)
 {
     n = x.size();
 
@@ -179,9 +157,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
     fvec.resize(n);
     qtf.resize(n);
     fjac.resize(n, n);
-    if (mode != 2)
+    if (!useExternalScaling)
         diag.resize(n);
-    assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
+    assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
 
     /* Function Body */
     nfev = 0;
@@ -190,7 +168,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
     /*     check the input parameters for errors. */
     if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. )
         return HybridNonLinearSolverSpace::ImproperInputParameters;
-    if (mode == 2)
+    if (useExternalScaling)
         for (int j = 0; j < n; ++j)
             if (diag[j] <= 0.)
                 return HybridNonLinearSolverSpace::ImproperInputParameters;
@@ -214,10 +192,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(FVectorType  &x)
 {
     int j;
     std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@@ -231,10 +206,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
 
     wa2 = fjac.colwise().blueNorm();
 
-    /* on the first iteration and if mode is 1, scale according */
+    /* on the first iteration and if external scaling is not used, scale according */
     /* to the norms of the columns of the initial jacobian. */
     if (iter == 1) {
-        if (mode != 2)
+        if (!useExternalScaling)
             for (j = 0; j < n; ++j)
                 diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
 
@@ -260,7 +235,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
     qtf = fjac.transpose() * fvec;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     while (true) {
@@ -372,14 +347,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solve(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solve(FVectorType  &x)
 {
-    HybridNonLinearSolverSpace::Status status = solveInit(x, mode);
+    HybridNonLinearSolverSpace::Status status = solveInit(x);
     while (status==HybridNonLinearSolverSpace::Running)
-        status = solveOneStep(x, mode);
+        status = solveOneStep(x);
     return status;
 }
 
@@ -403,18 +375,13 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrd1(
     parameters.xtol = tol;
 
     diag.setConstant(n, 1.);
-    return solveNumericalDiff(
-        x,
-        2
-    );
+    useExternalScaling = true;
+    return solveNumericalDiff(x);
 }
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(FVectorType  &x)
 {
     n = x.size();
 
@@ -425,10 +392,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
     qtf.resize(n);
     fjac.resize(n, n);
     fvec.resize(n);
-    if (mode != 2)
+    if (!useExternalScaling)
         diag.resize(n);
-    assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
-
+    assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
 
     /* Function Body */
     nfev = 0;
@@ -437,7 +403,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
     /*     check the input parameters for errors. */
     if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. )
         return HybridNonLinearSolverSpace::ImproperInputParameters;
-    if (mode == 2)
+    if (useExternalScaling)
         for (int j = 0; j < n; ++j)
             if (diag[j] <= 0.)
                 return HybridNonLinearSolverSpace::ImproperInputParameters;
@@ -461,10 +427,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(FVectorType  &x)
 {
     int j;
     std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@@ -480,10 +443,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
 
     wa2 = fjac.colwise().blueNorm();
 
-    /* on the first iteration and if mode is 1, scale according */
+    /* on the first iteration and if external scaling is not used, scale according */
     /* to the norms of the columns of the initial jacobian. */
     if (iter == 1) {
-        if (mode != 2)
+        if (!useExternalScaling)
             for (j = 0; j < n; ++j)
                 diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
 
@@ -509,7 +472,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
     qtf = fjac.transpose() * fvec;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     while (true) {
@@ -621,14 +584,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(FVectorType  &x)
 {
-    HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x, mode);
+    HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x);
     while (status==HybridNonLinearSolverSpace::Running)
-        status = solveNumericalDiffOneStep(x, mode);
+        status = solveNumericalDiffOneStep(x);
     return status;
 }
 
diff --git a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
index f21331e3e..8bae1e131 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
@@ -61,7 +61,7 @@ class LevenbergMarquardt
 {
 public:
     LevenbergMarquardt(FunctorType &_functor)
-        : functor(_functor) { nfev = njev = iter = 0;  fnorm=gnorm = 0.; }
+        : functor(_functor) { nfev = njev = iter = 0;  fnorm = gnorm = 0.; useExternalScaling=false; }
 
     struct Parameters {
         Parameters()
@@ -87,18 +87,9 @@ public:
             const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
             );
 
-    LevenbergMarquardtSpace::Status minimize(
-            FVectorType &x,
-            const int mode=1
-            );
-    LevenbergMarquardtSpace::Status minimizeInit(
-            FVectorType &x,
-            const int mode=1
-            );
-    LevenbergMarquardtSpace::Status minimizeOneStep(
-            FVectorType &x,
-            const int mode=1
-            );
+    LevenbergMarquardtSpace::Status minimize(FVectorType &x);
+    LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x);
+    LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x);
 
     static LevenbergMarquardtSpace::Status lmdif1(
             FunctorType &functor,
@@ -112,18 +103,9 @@ public:
             const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
             );
 
-    LevenbergMarquardtSpace::Status minimizeOptimumStorage(
-            FVectorType  &x,
-            const int mode=1
-            );
-    LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(
-            FVectorType  &x,
-            const int mode=1
-            );
-    LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(
-            FVectorType  &x,
-            const int mode=1
-            );
+    LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType  &x);
+    LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType  &x);
+    LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType  &x);
 
     void resetParameters(void) { parameters = Parameters(); }
 
@@ -135,6 +117,7 @@ public:
     int njev;
     int iter;
     Scalar fnorm, gnorm;
+    bool useExternalScaling; 
 
     Scalar lm_param(void) { return par; }
 private:
@@ -175,24 +158,18 @@ LevenbergMarquardt<FunctorType,Scalar>::lmder1(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimize(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType  &x)
 {
-    LevenbergMarquardtSpace::Status status = minimizeInit(x, mode);
+    LevenbergMarquardtSpace::Status status = minimizeInit(x);
     do {
-        status = minimizeOneStep(x, mode);
+        status = minimizeOneStep(x);
     } while (status==LevenbergMarquardtSpace::Running);
     return status;
 }
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType  &x)
 {
     n = x.size();
     m = functor.values();
@@ -201,9 +178,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
     wa4.resize(m);
     fvec.resize(m);
     fjac.resize(m, n);
-    if (mode != 2)
+    if (!useExternalScaling)
         diag.resize(n);
-    assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
+    assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
     qtf.resize(n);
 
     /* Function Body */
@@ -214,7 +191,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
     if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
         return LevenbergMarquardtSpace::ImproperInputParameters;
 
-    if (mode == 2)
+    if (useExternalScaling)
         for (int j = 0; j < n; ++j)
             if (diag[j] <= 0.)
                 return LevenbergMarquardtSpace::ImproperInputParameters;
@@ -235,10 +212,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType  &x)
 {
     int j;
 
@@ -257,10 +231,10 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
     fjac = qrfac.matrixQR();
     permutation = qrfac.colsPermutation();
 
-    /* on the first iteration and if mode is 1, scale according */
+    /* on the first iteration and if external scaling is not used, scale according */
     /* to the norms of the columns of the initial jacobian. */
     if (iter == 1) {
-        if (mode != 2)
+        if (!useExternalScaling)
             for (j = 0; j < n; ++j)
                 diag[j] = (wa2[j]==0.)? 1. : wa2[j];
 
@@ -290,7 +264,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
         return LevenbergMarquardtSpace::CosinusTooSmall;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     do {
@@ -406,10 +380,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType  &x)
 {
     n = x.size();
     m = functor.values();
@@ -423,9 +394,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
     // The purpose it to only use a nxn matrix, instead of mxn here, so
     // that we can handle cases where m>>n :
     fjac.resize(n, n);
-    if (mode != 2)
+    if (!useExternalScaling)
         diag.resize(n);
-    assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
+    assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
     qtf.resize(n);
 
     /* Function Body */
@@ -436,7 +407,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
     if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
         return LevenbergMarquardtSpace::ImproperInputParameters;
 
-    if (mode == 2)
+    if (useExternalScaling)
         for (int j = 0; j < n; ++j)
             if (diag[j] <= 0.)
                 return LevenbergMarquardtSpace::ImproperInputParameters;
@@ -458,10 +429,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType  &x)
 {
     int i, j;
     bool sing;
@@ -514,10 +482,10 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
         }
     }
 
-    /* on the first iteration and if mode is 1, scale according */
+    /* on the first iteration and if external scaling is not used, scale according */
     /* to the norms of the columns of the initial jacobian. */
     if (iter == 1) {
-        if (mode != 2)
+        if (!useExternalScaling)
             for (j = 0; j < n; ++j)
                 diag[j] = (wa2[j]==0.)? 1. : wa2[j];
 
@@ -541,7 +509,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
         return LevenbergMarquardtSpace::CosinusTooSmall;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     do {
@@ -635,14 +603,11 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType  &x)
 {
-    LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x, mode);
+    LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x);
     do {
-        status = minimizeOptimumStorageOneStep(x, mode);
+        status = minimizeOptimumStorageOneStep(x);
     } while (status==LevenbergMarquardtSpace::Running);
     return status;
 }
diff --git a/unsupported/Eigen/src/Skyline/SkylineMatrixBase.h b/unsupported/Eigen/src/Skyline/SkylineMatrixBase.h
index b90a6f9e9..ff20b830c 100644
--- a/unsupported/Eigen/src/Skyline/SkylineMatrixBase.h
+++ b/unsupported/Eigen/src/Skyline/SkylineMatrixBase.h
@@ -36,7 +36,7 @@
  * \param Derived
  *
  */
-template<typename Derived> class SkylineMatrixBase : public AnyMatrixBase<Derived> {
+template<typename Derived> class SkylineMatrixBase : public EigenBase<Derived> {
 public:
 
     typedef typename ei_traits<Derived>::Scalar Scalar;
diff --git a/unsupported/test/NonLinearOptimization.cpp b/unsupported/test/NonLinearOptimization.cpp
index 1313726a1..7aea7b361 100644
--- a/unsupported/test/NonLinearOptimization.cpp
+++ b/unsupported/test/NonLinearOptimization.cpp
@@ -317,7 +317,8 @@ void testHybrj()
   hybrj_functor functor;
   HybridNonLinearSolver<hybrj_functor> solver(functor);
   solver.diag.setConstant(n, 1.);
-  info = solver.solve(x, 2);
+  solver.useExternalScaling = true;
+  info = solver.solve(x);
 
   // check return value
   VERIFY( 1 == info);
@@ -401,7 +402,8 @@ void testHybrd()
   solver.parameters.nb_of_subdiagonals = 1;
   solver.parameters.nb_of_superdiagonals = 1;
   solver.diag.setConstant(n, 1.);
-  info = solver.solveNumericalDiff(x, 2);
+  solver.useExternalScaling = true;
+  info = solver.solveNumericalDiff(x);
 
   // check return value
   VERIFY( 1 == info);
diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp
index 6150439c5..86e942edb 100644
--- a/unsupported/test/matrix_exponential.cpp
+++ b/unsupported/test/matrix_exponential.cpp
@@ -57,7 +57,7 @@ void test2dRotation(double tol)
     angle = static_cast<T>(pow(10, i / 5. - 2));
     B << cos(angle), sin(angle), -sin(angle), cos(angle);
 
-    ei_matrix_function(angle*A, expfn, &C);
+    C = ei_matrix_function(angle*A, expfn);
     std::cout << "test2dRotation: i = " << i << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
@@ -82,7 +82,7 @@ void test2dHyperbolicRotation(double tol)
     A << 0, angle*imagUnit, -angle*imagUnit, 0;
     B << ch, sh*imagUnit, -sh*imagUnit, ch;
 
-    ei_matrix_function(A, expfn, &C);
+    C = ei_matrix_function(A, expfn);
     std::cout << "test2dHyperbolicRotation: i = " << i << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
@@ -106,7 +106,7 @@ void testPascal(double tol)
       for (int j=0; j<=i; j++)
     B(i,j) = static_cast<T>(binom(i,j));
 
-    ei_matrix_function(A, expfn, &C);
+    C = ei_matrix_function(A, expfn);
     std::cout << "testPascal: size = " << size << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
@@ -132,10 +132,9 @@ void randomTest(const MatrixType& m, double tol)
   for(int i = 0; i < g_repeat; i++) {
     m1 = MatrixType::Random(rows, cols);
 
-    ei_matrix_function(m1, expfn, &m2);
-    ei_matrix_function(-m1, expfn, &m3);
+    m2 = ei_matrix_function(m1, expfn) * ei_matrix_function(-m1, expfn);
     std::cout << "randomTest: error funm = " << relerr(identity, m2 * m3);
-    VERIFY(identity.isApprox(m2 * m3, static_cast<RealScalar>(tol)));
+    VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
 
     m2 = ei_matrix_exponential(m1) * ei_matrix_exponential(-m1);
     std::cout << "   error expm = " << relerr(identity, m2) << "\n";
diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp
index 4ff6d7f1e..7a1501da2 100644
--- a/unsupported/test/matrix_function.cpp
+++ b/unsupported/test/matrix_function.cpp
@@ -109,11 +109,10 @@ template<typename MatrixType>
 void testHyperbolicFunctions(const MatrixType& A)
 {
   for (int i = 0; i < g_repeat; i++) {
-    MatrixType sinhA = ei_matrix_sinh(A);
-    MatrixType coshA = ei_matrix_cosh(A);
     MatrixType expA = ei_matrix_exponential(A);
-    VERIFY_IS_APPROX(sinhA, (expA - expA.inverse())/2);
-    VERIFY_IS_APPROX(coshA, (expA + expA.inverse())/2);
+    MatrixType expmA = ei_matrix_exponential(-A);
+    VERIFY_IS_APPROX(ei_matrix_sinh(A), (expA - expmA) / 2);
+    VERIFY_IS_APPROX(ei_matrix_cosh(A), (expA + expmA) / 2);
   }
 }
 
@@ -134,14 +133,15 @@ void testGonioFunctions(const MatrixType& A)
     ComplexMatrix Ac = A.template cast<ComplexScalar>();
 
     ComplexMatrix exp_iA = ei_matrix_exponential(imagUnit * Ac);
+    ComplexMatrix exp_miA = ei_matrix_exponential(-imagUnit * Ac);
 
     MatrixType sinA = ei_matrix_sin(A);
     ComplexMatrix sinAc = sinA.template cast<ComplexScalar>();
-    VERIFY_IS_APPROX(sinAc, (exp_iA - exp_iA.inverse()) / (two*imagUnit));
+    VERIFY_IS_APPROX(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
 
     MatrixType cosA = ei_matrix_cos(A);
     ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
-    VERIFY_IS_APPROX(cosAc, (exp_iA + exp_iA.inverse()) / 2);
+    VERIFY_IS_APPROX(cosAc, (exp_iA + exp_miA) / 2);
   }
 }