merge

7 files changed, 201 insertions, 177 deletions

diff --git a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
index 35dc332e0..d75b1407c 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
@@ -194,6 +194,8 @@ template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
 HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(FVectorType  &x)
 {
+    assert(x.size()==n); // check the caller is not cheating us
+
     int j;
     std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
 
@@ -350,6 +352,8 @@ HybridNonLinearSolverSpace::Status
 HybridNonLinearSolver<FunctorType,Scalar>::solve(FVectorType  &x)
 {
     HybridNonLinearSolverSpace::Status status = solveInit(x);
+    if (status==HybridNonLinearSolverSpace::ImproperInputParameters)
+        return status;
     while (status==HybridNonLinearSolverSpace::Running)
         status = solveOneStep(x);
     return status;
@@ -429,6 +433,8 @@ template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
 HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(FVectorType  &x)
 {
+    assert(x.size()==n); // check the caller is not cheating us
+
     int j;
     std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
 
@@ -587,6 +593,8 @@ HybridNonLinearSolverSpace::Status
 HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(FVectorType  &x)
 {
     HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x);
+    if (status==HybridNonLinearSolverSpace::ImproperInputParameters)
+        return status;
     while (status==HybridNonLinearSolverSpace::Running)
         status = solveNumericalDiffOneStep(x);
     return status;
diff --git a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
index 8bae1e131..f99366bbc 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
@@ -161,6 +161,8 @@ LevenbergMarquardtSpace::Status
 LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType  &x)
 {
     LevenbergMarquardtSpace::Status status = minimizeInit(x);
+    if (status==LevenbergMarquardtSpace::ImproperInputParameters)
+        return status;
     do {
         status = minimizeOneStep(x);
     } while (status==LevenbergMarquardtSpace::Running);
@@ -214,7 +216,7 @@ template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
 LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType  &x)
 {
-    int j;
+    assert(x.size()==n); // check the caller is not cheating us
 
     /* calculate the jacobian matrix. */
     int df_ret = functor.df(x, fjac);
@@ -235,7 +237,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType  &x)
     /* to the norms of the columns of the initial jacobian. */
     if (iter == 1) {
         if (!useExternalScaling)
-            for (j = 0; j < n; ++j)
+            for (int j = 0; j < n; ++j)
                 diag[j] = (wa2[j]==0.)? 1. : wa2[j];
 
         /* on the first iteration, calculate the norm of the scaled x */
@@ -255,7 +257,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType  &x)
     /* compute the norm of the scaled gradient. */
     gnorm = 0.;
     if (fnorm != 0.)
-        for (j = 0; j < n; ++j)
+        for (int j = 0; j < n; ++j)
             if (wa2[permutation.indices()[j]] != 0.)
                 gnorm = std::max(gnorm, ei_abs( fjac.col(j).head(j+1).dot(qtf.head(j+1)/fnorm) / wa2[permutation.indices()[j]]));
 
@@ -431,6 +433,8 @@ template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType  &x)
 {
+    assert(x.size()==n); // check the caller is not cheating us
+
     int i, j;
     bool sing;
 
@@ -606,6 +610,8 @@ LevenbergMarquardtSpace::Status
 LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType  &x)
 {
     LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x);
+    if (status==LevenbergMarquardtSpace::ImproperInputParameters)
+        return status;
     do {
         status = minimizeOptimumStorageOneStep(x);
     } while (status==LevenbergMarquardtSpace::Running);
diff --git a/unsupported/doc/examples/MatrixFunction.cpp b/unsupported/doc/examples/MatrixFunction.cpp
index 075fe7361..9b594cf39 100644
--- a/unsupported/doc/examples/MatrixFunction.cpp
+++ b/unsupported/doc/examples/MatrixFunction.cpp
@@ -18,5 +18,5 @@ int main()
 
   std::cout << "The matrix A is:\n" << A << "\n\n";
   std::cout << "The matrix exponential of A is:\n" 
-	    << ei_matrix_function(A, expfn) << "\n\n";
+            << ei_matrix_function(A, expfn) << "\n\n";
 }
diff --git a/unsupported/test/FFTW.cpp b/unsupported/test/FFTW.cpp
index 1cc8b156f..131b4b762 100644
--- a/unsupported/test/FFTW.cpp
+++ b/unsupported/test/FFTW.cpp
@@ -226,29 +226,28 @@ void test_complex2d()
 */
 
 
-template <typename T,int nrows,int ncols>
-void test_return_by_value()
+void test_return_by_value(int len)
 {
-    Matrix<complex<T>,nrows,ncols> in;
-    Matrix<complex<T>,nrows,ncols> in1;
-    in.Random();
-    Matrix<complex<T>,nrows,ncols> out1;
-    Matrix<complex<T>,nrows,ncols> out2;
-    FFT<T> fft;
+    VectorXf in;
+    VectorXf in1;
+    in = in.Random( len );
+    VectorXcf out1,out2;
+    FFT<float> fft;
 
     fft.SetFlag(fft.HalfSpectrum );
 
     fft.fwd(out1,in);
     out2 = fft.fwd(in);
-    VERIFY( (out1-out2).norm() < test_precision<T>() );
+    VERIFY( (out1-out2).norm() < test_precision<float>() );
     in1 = fft.inv(out1);
-    VERIFY( (in1-in).norm() < test_precision<T>() );
+    VERIFY( (in1-in).norm() < test_precision<float>() );
 }
 
 void test_FFTW()
 {
-    test_return_by_value<float,1,32>();
-    test_return_by_value<double,1,32>();
+    cout << "testing return-by-value\n";
+    CALL_SUBTEST( test_return_by_value(32) );
+    cout << "testing complex\n";
   //CALL_SUBTEST( ( test_complex2d<float,4,8> () ) ); CALL_SUBTEST( ( test_complex2d<double,4,8> () ) );
   //CALL_SUBTEST( ( test_complex2d<long double,4,8> () ) );
   CALL_SUBTEST( test_complex<float>(32) ); CALL_SUBTEST( test_complex<double>(32) ); CALL_SUBTEST( test_complex<long double>(32) );
@@ -259,6 +258,7 @@ void test_FFTW()
   CALL_SUBTEST( test_complex<float>(2*3*4*5) ); CALL_SUBTEST( test_complex<double>(2*3*4*5) ); CALL_SUBTEST( test_complex<long double>(2*3*4*5) );
   CALL_SUBTEST( test_complex<float>(2*3*4*5*7) ); CALL_SUBTEST( test_complex<double>(2*3*4*5*7) ); CALL_SUBTEST( test_complex<long double>(2*3*4*5*7) );
 
+    cout << "testing scalar\n";
   CALL_SUBTEST( test_scalar<float>(32) ); CALL_SUBTEST( test_scalar<double>(32) ); CALL_SUBTEST( test_scalar<long double>(32) );
   CALL_SUBTEST( test_scalar<float>(45) ); CALL_SUBTEST( test_scalar<double>(45) ); CALL_SUBTEST( test_scalar<long double>(45) );
   CALL_SUBTEST( test_scalar<float>(50) ); CALL_SUBTEST( test_scalar<double>(50) ); CALL_SUBTEST( test_scalar<long double>(50) );
diff --git a/unsupported/test/NonLinearOptimization.cpp b/unsupported/test/NonLinearOptimization.cpp
index 7aea7b361..e68745ad1 100644
--- a/unsupported/test/NonLinearOptimization.cpp
+++ b/unsupported/test/NonLinearOptimization.cpp
@@ -172,9 +172,9 @@ void testLmder1()
   info = lm.lmder1(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY(lm.nfev==6);
-  VERIFY(lm.njev==5);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 6);
+  VERIFY_IS_EQUAL(lm.njev, 5);
 
   // check norm
   VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596);
@@ -201,9 +201,9 @@ void testLmder()
   info = lm.minimize(x);
 
   // check return values
-  VERIFY( 1 == info);
-  VERIFY(lm.nfev==6);
-  VERIFY(lm.njev==5);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 6);
+  VERIFY_IS_EQUAL(lm.njev, 5);
 
   // check norm
   fnorm = lm.fvec.blueNorm();
@@ -286,9 +286,9 @@ void testHybrj1()
   info = solver.hybrj1(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY(solver.nfev==11);
-  VERIFY(solver.njev==1);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(solver.nfev, 11);
+  VERIFY_IS_EQUAL(solver.njev, 1);
 
   // check norm
   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
@@ -321,9 +321,9 @@ void testHybrj()
   info = solver.solve(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY(solver.nfev==11);
-  VERIFY(solver.njev==1);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(solver.nfev, 11);
+  VERIFY_IS_EQUAL(solver.njev, 1);
 
   // check norm
   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
@@ -375,8 +375,8 @@ void testHybrd1()
   info = solver.hybrd1(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY(solver.nfev==20);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(solver.nfev, 20);
 
   // check norm
   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
@@ -406,8 +406,8 @@ void testHybrd()
   info = solver.solveNumericalDiff(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY(solver.nfev==14);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(solver.nfev, 14);
 
   // check norm
   VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08);
@@ -477,9 +477,9 @@ void testLmstr1()
   info = lm.lmstr1(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY(lm.nfev==6);
-  VERIFY(lm.njev==5);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 6);
+  VERIFY_IS_EQUAL(lm.njev, 5);
 
   // check norm
   VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596);
@@ -506,9 +506,9 @@ void testLmstr()
   info = lm.minimizeOptimumStorage(x);
 
   // check return values
-  VERIFY( 1 == info);
-  VERIFY(lm.nfev==6);
-  VERIFY(lm.njev==5);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 6);
+  VERIFY_IS_EQUAL(lm.njev, 5);
 
   // check norm
   fnorm = lm.fvec.blueNorm();
@@ -562,8 +562,8 @@ void testLmdif1()
   info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY(nfev==26);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(nfev, 26);
 
   // check norm
   functor(x, fvec);
@@ -593,8 +593,8 @@ void testLmdif()
   info = lm.minimize(x);
 
   // check return values
-  VERIFY( 1 == info);
-  VERIFY(lm.nfev==26);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 26);
 
   // check norm
   fnorm = lm.fvec.blueNorm();
@@ -678,9 +678,9 @@ void testNistChwirut2(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 10 == lm.nfev);
-  VERIFY( 8 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 10);
+  VERIFY_IS_EQUAL(lm.njev, 8);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02);
   // check x
@@ -699,9 +699,9 @@ void testNistChwirut2(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 7 == lm.nfev);
-  VERIFY( 6 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 7);
+  VERIFY_IS_EQUAL(lm.njev, 6);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02);
   // check x
@@ -758,9 +758,9 @@ void testNistMisra1a(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 19 == lm.nfev);
-  VERIFY( 15 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 19);
+  VERIFY_IS_EQUAL(lm.njev, 15);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01);
   // check x
@@ -775,9 +775,9 @@ void testNistMisra1a(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 5 == lm.nfev);
-  VERIFY( 4 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 5);
+  VERIFY_IS_EQUAL(lm.njev, 4);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01);
   // check x
@@ -844,19 +844,19 @@ void testNistHahn1(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 11== lm.nfev);
-  VERIFY( 10== lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 11);
+  VERIFY_IS_EQUAL(lm.njev, 10);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00);
   // check x
-  VERIFY_IS_APPROX(x[0], 1.0776351733E+00  );
-  VERIFY_IS_APPROX(x[1],-1.2269296921E-01  );
-  VERIFY_IS_APPROX(x[2], 4.0863750610E-03  );
+  VERIFY_IS_APPROX(x[0], 1.0776351733E+00);
+  VERIFY_IS_APPROX(x[1],-1.2269296921E-01);
+  VERIFY_IS_APPROX(x[2], 4.0863750610E-03);
   VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06
-  VERIFY_IS_APPROX(x[4],-5.7609940901E-03  );
-  VERIFY_IS_APPROX(x[5], 2.4053735503E-04  );
-  VERIFY_IS_APPROX(x[6],-1.2314450199E-07  );
+  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
+  VERIFY_IS_APPROX(x[5], 2.4053735503E-04);
+  VERIFY_IS_APPROX(x[6],-1.2314450199E-07);
 
   /*
    * Second try
@@ -866,9 +866,9 @@ void testNistHahn1(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 11 == lm.nfev);
-  VERIFY( 10 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 11);
+  VERIFY_IS_EQUAL(lm.njev, 10);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00);
   // check x
@@ -876,7 +876,7 @@ void testNistHahn1(void)
   VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01
   VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03
   VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06
-  VERIFY_IS_APPROX(x[4],-5.7609940901E-03  );
+  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
   VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04
   VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07
 
@@ -930,9 +930,9 @@ void testNistMisra1d(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 3 == info);
-  VERIFY( 9 == lm.nfev);
-  VERIFY( 7 == lm.njev);
+  VERIFY_IS_EQUAL(info, 3);
+  VERIFY_IS_EQUAL(lm.nfev, 9);
+  VERIFY_IS_EQUAL(lm.njev, 7);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02);
   // check x
@@ -947,9 +947,9 @@ void testNistMisra1d(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 4 == lm.nfev);
-  VERIFY( 3 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 4);
+  VERIFY_IS_EQUAL(lm.njev, 3);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02);
   // check x
@@ -1008,18 +1008,18 @@ void testNistLanczos1(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 2 == info);
-  VERIFY( 79 == lm.nfev);
-  VERIFY( 72 == lm.njev);
+  VERIFY_IS_EQUAL(info, 2);
+  VERIFY_IS_EQUAL(lm.nfev, 79);
+  VERIFY_IS_EQUAL(lm.njev, 72);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430899764097e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
   // check x
-  VERIFY_IS_APPROX(x[0], 9.5100000027E-02 );
-  VERIFY_IS_APPROX(x[1], 1.0000000001E+00 );
-  VERIFY_IS_APPROX(x[2], 8.6070000013E-01 );
-  VERIFY_IS_APPROX(x[3], 3.0000000002E+00 );
-  VERIFY_IS_APPROX(x[4], 1.5575999998E+00 );
-  VERIFY_IS_APPROX(x[5], 5.0000000001E+00 );
+  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
+  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
+  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
+  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
+  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
+  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
 
   /*
    * Second try
@@ -1029,18 +1029,18 @@ void testNistLanczos1(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 2 == info);
-  VERIFY( 9 == lm.nfev);
-  VERIFY( 8 == lm.njev);
+  VERIFY_IS_EQUAL(info, 2);
+  VERIFY_IS_EQUAL(lm.nfev, 9);
+  VERIFY_IS_EQUAL(lm.njev, 8);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.428595533845e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
   // check x
-  VERIFY_IS_APPROX(x[0], 9.5100000027E-02 );
-  VERIFY_IS_APPROX(x[1], 1.0000000001E+00 );
-  VERIFY_IS_APPROX(x[2], 8.6070000013E-01 );
-  VERIFY_IS_APPROX(x[3], 3.0000000002E+00 );
-  VERIFY_IS_APPROX(x[4], 1.5575999998E+00 );
-  VERIFY_IS_APPROX(x[5], 5.0000000001E+00 );
+  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
+  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
+  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
+  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
+  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
+  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
 
 }
 
@@ -1094,9 +1094,9 @@ void testNistRat42(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 10 == lm.nfev);
-  VERIFY( 8 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 10);
+  VERIFY_IS_EQUAL(lm.njev, 8);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00);
   // check x
@@ -1112,9 +1112,9 @@ void testNistRat42(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 6 == lm.nfev);
-  VERIFY( 5 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 6);
+  VERIFY_IS_EQUAL(lm.njev, 5);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00);
   // check x
@@ -1172,9 +1172,9 @@ void testNistMGH10(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 2 == info); 
-  VERIFY( 284 == lm.nfev); 
-  VERIFY( 249 == lm.njev); 
+  VERIFY_IS_EQUAL(info, 2); 
+  VERIFY_IS_EQUAL(lm.nfev, 284 ); 
+  VERIFY_IS_EQUAL(lm.njev, 249 ); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01);
   // check x
@@ -1190,9 +1190,9 @@ void testNistMGH10(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 3 == info);
-  VERIFY( 126 == lm.nfev);
-  VERIFY( 116 == lm.njev);
+  VERIFY_IS_EQUAL(info, 3);
+  VERIFY_IS_EQUAL(lm.nfev, 126);
+  VERIFY_IS_EQUAL(lm.njev, 116);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01);
   // check x
@@ -1251,9 +1251,9 @@ void testNistBoxBOD(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 31 == lm.nfev);
-  VERIFY( 25  == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 31);
+  VERIFY_IS_EQUAL(lm.njev, 25);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
   // check x
@@ -1271,9 +1271,9 @@ void testNistBoxBOD(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info); 
-  VERIFY( 15 == lm.nfev); 
-  VERIFY( 14 == lm.njev); 
+  VERIFY_IS_EQUAL(info, 1); 
+  VERIFY_IS_EQUAL(lm.nfev, 15 ); 
+  VERIFY_IS_EQUAL(lm.njev, 14 ); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
   // check x
@@ -1333,9 +1333,9 @@ void testNistMGH17(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 2 == info); 
-  VERIFY( 602 == lm.nfev); 
-  VERIFY( 545 == lm.njev); 
+  VERIFY_IS_EQUAL(info, 2); 
+  VERIFY_IS_EQUAL(lm.nfev, 602 ); 
+  VERIFY_IS_EQUAL(lm.njev, 545 ); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
   // check x
@@ -1354,9 +1354,9 @@ void testNistMGH17(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 18 == lm.nfev);
-  VERIFY( 15 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 18);
+  VERIFY_IS_EQUAL(lm.njev, 15);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
   // check x
@@ -1420,9 +1420,9 @@ void testNistMGH09(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info); 
-  VERIFY( 490 == lm.nfev); 
-  VERIFY( 376 == lm.njev); 
+  VERIFY_IS_EQUAL(info, 1); 
+  VERIFY_IS_EQUAL(lm.nfev, 490 ); 
+  VERIFY_IS_EQUAL(lm.njev, 376 ); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04);
   // check x
@@ -1440,9 +1440,9 @@ void testNistMGH09(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 18 == lm.nfev);
-  VERIFY( 16 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 18);
+  VERIFY_IS_EQUAL(lm.njev, 16);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04);
   // check x
@@ -1503,9 +1503,9 @@ void testNistBennett5(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 758 == lm.nfev);
-  VERIFY( 744 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 758);
+  VERIFY_IS_EQUAL(lm.njev, 744);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04);
   // check x
@@ -1521,9 +1521,9 @@ void testNistBennett5(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 203 == lm.nfev);
-  VERIFY( 192 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 203);
+  VERIFY_IS_EQUAL(lm.njev, 192);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04);
   // check x
@@ -1591,9 +1591,9 @@ void testNistThurber(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 39 == lm.nfev);
-  VERIFY( 36== lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 39);
+  VERIFY_IS_EQUAL(lm.njev, 36);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03);
   // check x
@@ -1616,9 +1616,9 @@ void testNistThurber(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 29 == lm.nfev);
-  VERIFY( 28 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 29);
+  VERIFY_IS_EQUAL(lm.njev, 28);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03);
   // check x
@@ -1683,9 +1683,9 @@ void testNistRat43(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 27 == lm.nfev);
-  VERIFY( 20 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 27);
+  VERIFY_IS_EQUAL(lm.njev, 20);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03);
   // check x
@@ -1705,9 +1705,9 @@ void testNistRat43(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 9 == lm.nfev);
-  VERIFY( 8 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 9);
+  VERIFY_IS_EQUAL(lm.njev, 8);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03);
   // check x
@@ -1768,9 +1768,9 @@ void testNistEckerle4(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 18 == lm.nfev);
-  VERIFY( 15 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 18);
+  VERIFY_IS_EQUAL(lm.njev, 15);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03);
   // check x
@@ -1786,9 +1786,9 @@ void testNistEckerle4(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 7 == lm.nfev);
-  VERIFY( 6 == lm.njev);
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev, 7);
+  VERIFY_IS_EQUAL(lm.njev, 6);
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03);
   // check x
diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp
index 86e942edb..61f30334d 100644
--- a/unsupported/test/matrix_exponential.cpp
+++ b/unsupported/test/matrix_exponential.cpp
@@ -133,7 +133,7 @@ void randomTest(const MatrixType& m, double tol)
     m1 = MatrixType::Random(rows, cols);
 
     m2 = ei_matrix_function(m1, expfn) * ei_matrix_function(-m1, expfn);
-    std::cout << "randomTest: error funm = " << relerr(identity, m2 * m3);
+    std::cout << "randomTest: error funm = " << relerr(identity, m2);
     VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
 
     m2 = ei_matrix_exponential(m1) * ei_matrix_exponential(-m1);
diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp
index 7a1501da2..e40af7b4e 100644
--- a/unsupported/test/matrix_function.cpp
+++ b/unsupported/test/matrix_function.cpp
@@ -25,6 +25,17 @@
 #include "main.h"
 #include <unsupported/Eigen/MatrixFunctions>
 
+// Variant of VERIFY_IS_APPROX which uses absolute error instead of
+// relative error.
+#define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b))
+
+template<typename Type1, typename Type2>
+inline bool test_isApprox_abs(const Type1& a, const Type2& b)
+{
+  return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all();
+}
+
+
 // Returns a matrix with eigenvalues clustered around 0, 1 and 2.
 template<typename MatrixType>
 MatrixType randomMatrixWithRealEivals(const int size)
@@ -37,7 +48,8 @@ MatrixType randomMatrixWithRealEivals(const int size)
       + ei_random<Scalar>() * Scalar(RealScalar(0.01));
   }
   MatrixType A = MatrixType::Random(size, size);
-  return A.inverse() * diag * A;
+  HouseholderQR<MatrixType> QRofA(A);
+  return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
 }
 
 template <typename MatrixType, int IsComplex = NumTraits<typename ei_traits<MatrixType>::Scalar>::IsComplex>
@@ -69,7 +81,8 @@ struct randomMatrixWithImagEivals<MatrixType, 0>
       }
     }
     MatrixType A = MatrixType::Random(size, size);
-    return A.inverse() * diag * A;
+    HouseholderQR<MatrixType> QRofA(A);
+    return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
   }
 };
 
@@ -88,10 +101,12 @@ struct randomMatrixWithImagEivals<MatrixType, 1>
         + ei_random<Scalar>() * Scalar(RealScalar(0.01));
     }
     MatrixType A = MatrixType::Random(size, size);
-    return A.inverse() * diag * A;
+    HouseholderQR<MatrixType> QRofA(A);
+    return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
   }
 };
 
+
 template<typename MatrixType>
 void testMatrixExponential(const MatrixType& A)
 {
@@ -99,50 +114,45 @@ void testMatrixExponential(const MatrixType& A)
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef std::complex<RealScalar> ComplexScalar;
 
-  for (int i = 0; i < g_repeat; i++) {
-    VERIFY_IS_APPROX(ei_matrix_exponential(A), 
-		     ei_matrix_function(A, StdStemFunctions<ComplexScalar>::exp));
-  }
+  VERIFY_IS_APPROX(ei_matrix_exponential(A), 
+		   ei_matrix_function(A, StdStemFunctions<ComplexScalar>::exp));
 }
 
 template<typename MatrixType>
 void testHyperbolicFunctions(const MatrixType& A)
 {
-  for (int i = 0; i < g_repeat; i++) {
-    MatrixType expA = ei_matrix_exponential(A);
-    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);
-  }
+  // Need to use absolute error because of possible cancellation when
+  // adding/subtracting expA and expmA.
+  MatrixType expA = ei_matrix_exponential(A);
+  MatrixType expmA = ei_matrix_exponential(-A);
+  VERIFY_IS_APPROX_ABS(ei_matrix_sinh(A), (expA - expmA) / 2);
+  VERIFY_IS_APPROX_ABS(ei_matrix_cosh(A), (expA + expmA) / 2);
 }
 
 template<typename MatrixType>
 void testGonioFunctions(const MatrixType& A)
 {
-  typedef ei_traits<MatrixType> Traits;
-  typedef typename Traits::Scalar Scalar;
+  typedef typename MatrixType::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef std::complex<RealScalar> ComplexScalar;
-  typedef Matrix<ComplexScalar, Traits::RowsAtCompileTime, 
-                 Traits::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
+  typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime, 
+                 MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
 
   ComplexScalar imagUnit(0,1);
   ComplexScalar two(2,0);
 
-  for (int i = 0; i < g_repeat; i++) {
-    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_miA) / (two*imagUnit));
-
-    MatrixType cosA = ei_matrix_cos(A);
-    ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
-    VERIFY_IS_APPROX(cosAc, (exp_iA + exp_miA) / 2);
-  }
+  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_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
+  
+  MatrixType cosA = ei_matrix_cos(A);
+  ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
+  VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
 }
 
 template<typename MatrixType>