merge

5 files changed, 147 insertions, 328 deletions

diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp
index 02cd5a48f..45c87f5a7 100644
--- a/unsupported/test/FFT.cpp
+++ b/unsupported/test/FFT.cpp
@@ -1,245 +1,2 @@
-#if 0
-
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra. Eigen itself is part of the KDE project.
-//
-// Copyright (C) 2009 Mark Borgerding mark a borgerding net
-//
-// Eigen is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 3 of the License, or (at your option) any later version.
-//
-// Alternatively, you can redistribute it and/or
-// modify it under the terms of the GNU General Public License as
-// published by the Free Software Foundation; either version 2 of
-// the License, or (at your option) any later version.
-//
-// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
-// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License and a copy of the GNU General Public License along with
-// Eigen. If not, see <http://www.gnu.org/licenses/>.
-
-#include "main.h"
-#include <unsupported/Eigen/FFT>
-
-template <typename T>
-std::complex<T> RandomCpx() { return std::complex<T>( (T)(rand()/(T)RAND_MAX - .5), (T)(rand()/(T)RAND_MAX - .5) ); }
-
-using namespace std;
-using namespace Eigen;
-
-float norm(float x) {return x*x;}
-double norm(double x) {return x*x;}
-long double norm(long double x) {return x*x;}
-
-template < typename T>
-complex<long double>  promote(complex<T> x) { return complex<long double>(x.real(),x.imag()); }
-
-complex<long double>  promote(float x) { return complex<long double>( x); }
-complex<long double>  promote(double x) { return complex<long double>( x); }
-complex<long double>  promote(long double x) { return complex<long double>( x); }
-    
-
-    template <typename T1,typename T2>
-    long double fft_rmse( const vector<T1> & fftbuf,const vector<T2> & timebuf)
-    {
-        long double totalpower=0;
-        long double difpower=0;
-        long double pi = acos((long double)-1 );
-        for (size_t k0=0;k0<fftbuf.size();++k0) {
-            complex<long double> acc = 0;
-            long double phinc = -2.*k0* pi / timebuf.size();
-            for (size_t k1=0;k1<timebuf.size();++k1) {
-                acc +=  promote( timebuf[k1] ) * exp( complex<long double>(0,k1*phinc) );
-            }
-            totalpower += norm(acc);
-            complex<long double> x = promote(fftbuf[k0]); 
-            complex<long double> dif = acc - x;
-            difpower += norm(dif);
-            cerr << k0 << "\t" << acc << "\t" <<  x << "\t" << sqrt(norm(dif)) << endl;
-        }
-        cerr << "rmse:" << sqrt(difpower/totalpower) << endl;
-        return sqrt(difpower/totalpower);
-    }
-
-    template <typename T1,typename T2>
-    long double dif_rmse( const vector<T1> buf1,const vector<T2> buf2)
-    {
-        long double totalpower=0;
-        long double difpower=0;
-        size_t n = min( buf1.size(),buf2.size() );
-        for (size_t k=0;k<n;++k) {
-            totalpower += (norm( buf1[k] ) + norm(buf2[k]) )/2.;
-            difpower += norm(buf1[k] - buf2[k]);
-        }
-        return sqrt(difpower/totalpower);
-    }
-
-enum { StdVectorContainer, EigenVectorContainer };
-
-template<int Container, typename Scalar> struct VectorType;
-
-template<typename Scalar> struct VectorType<StdVectorContainer,Scalar>
-{
-  typedef vector<Scalar> type;
-};
-
-template<typename Scalar> struct VectorType<EigenVectorContainer,Scalar>
-{
-  typedef Matrix<Scalar,Dynamic,1> type;
-};
-
-template <int Container, typename T>
-void test_scalar_generic(int nfft)
-{
-    typedef typename FFT<T>::Complex Complex;
-    typedef typename FFT<T>::Scalar Scalar;
-    typedef typename VectorType<Container,Scalar>::type ScalarVector;
-    typedef typename VectorType<Container,Complex>::type ComplexVector;
-
-    FFT<T> fft;
-    ScalarVector inbuf(nfft);
-    ComplexVector outbuf;
-    for (int k=0;k<nfft;++k)
-        inbuf[k]= (T)(rand()/(double)RAND_MAX - .5);
-
-    // make sure it DOESN'T give the right full spectrum answer
-    // if we've asked for half-spectrum
-    fft.SetFlag(fft.HalfSpectrum );
-    fft.fwd( outbuf,inbuf);
-    VERIFY(outbuf.size() == (size_t)( (nfft>>1)+1) );
-    VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>()  );// gross check
-
-    fft.ClearFlag(fft.HalfSpectrum );
-    fft.fwd( outbuf,inbuf);
-    VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>()  );// gross check
-
-    ScalarVector buf3;
-    fft.inv( buf3 , outbuf);
-    VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check
-
-    // verify that the Unscaled flag takes effect
-    ComplexVector buf4;
-    fft.SetFlag(fft.Unscaled);
-    fft.inv( buf4 , outbuf);
-    for (int k=0;k<nfft;++k)
-        buf4[k] *= T(1./nfft);
-    VERIFY( dif_rmse(inbuf,buf4) < test_precision<T>()  );// gross check
-
-    // verify that ClearFlag works
-    fft.ClearFlag(fft.Unscaled);
-    fft.inv( buf3 , outbuf);
-    VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check
-}
-
-template <typename T>
-void test_scalar(int nfft)
-{
-  test_scalar_generic<StdVectorContainer,T>(nfft);
-  test_scalar_generic<EigenVectorContainer,T>(nfft);
-}
-
-template <int Container, typename T>
-void test_complex_generic(int nfft)
-{
-    typedef typename FFT<T>::Complex Complex;
-    typedef typename VectorType<Container,Complex>::type ComplexVector;
-
-    FFT<T> fft;
-
-    ComplexVector inbuf(nfft);
-    ComplexVector outbuf;
-    ComplexVector buf3;
-    for (int k=0;k<nfft;++k)
-        inbuf[k]= Complex( (T)(rand()/(double)RAND_MAX - .5), (T)(rand()/(double)RAND_MAX - .5) );
-    fft.fwd( outbuf , inbuf);
-
-    VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>()  );// gross check
-
-    fft.inv( buf3 , outbuf);
-
-    VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check
-
-    // verify that the Unscaled flag takes effect
-    ComplexVector buf4;
-    fft.SetFlag(fft.Unscaled);
-    fft.inv( buf4 , outbuf);
-    for (int k=0;k<nfft;++k)
-        buf4[k] *= T(1./nfft);
-    VERIFY( dif_rmse(inbuf,buf4) < test_precision<T>()  );// gross check
-
-    // verify that ClearFlag works
-    fft.ClearFlag(fft.Unscaled);
-    fft.inv( buf3 , outbuf);
-    VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check
-}
-
-template <typename T>
-void test_complex(int nfft)
-{
-  test_complex_generic<StdVectorContainer,T>(nfft);
-  test_complex_generic<EigenVectorContainer,T>(nfft);
-}
-
-void test_FFT()
-{
-
-  CALL_SUBTEST( test_complex<float>(32) );
-  CALL_SUBTEST( test_complex<double>(32) );
-  CALL_SUBTEST( test_complex<long double>(32) );
-  
-  CALL_SUBTEST( test_complex<float>(256) );
-  CALL_SUBTEST( test_complex<double>(256) );
-  CALL_SUBTEST( test_complex<long double>(256) );
-  
-  CALL_SUBTEST( test_complex<float>(3*8) );
-  CALL_SUBTEST( test_complex<double>(3*8) );
-  CALL_SUBTEST( test_complex<long double>(3*8) );
-  
-  CALL_SUBTEST( test_complex<float>(5*32) );
-  CALL_SUBTEST( test_complex<double>(5*32) );
-  CALL_SUBTEST( test_complex<long double>(5*32) );
-  
-  CALL_SUBTEST( test_complex<float>(2*3*4) );
-  CALL_SUBTEST( test_complex<double>(2*3*4) );
-  CALL_SUBTEST( test_complex<long double>(2*3*4) );
-  
-  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) );
-
-
-
-  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) );
-  
-  CALL_SUBTEST( test_scalar<float>(256) );
-  CALL_SUBTEST( test_scalar<double>(256) );
-  CALL_SUBTEST( test_scalar<long double>(256) );
-  
-  CALL_SUBTEST( test_scalar<float>(2*3*4*5*7) );
-  CALL_SUBTEST( test_scalar<double>(2*3*4*5*7) );
-  CALL_SUBTEST( test_scalar<long double>(2*3*4*5*7) );
-}
-#else
 #define test_FFTW test_FFT
 #include "FFTW.cpp"
-#endif
diff --git a/unsupported/test/FFTW.cpp b/unsupported/test/FFTW.cpp
index 94de53e36..838ddb238 100644
--- a/unsupported/test/FFTW.cpp
+++ b/unsupported/test/FFTW.cpp
@@ -23,7 +23,6 @@
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
 #include "main.h"
-#include <iostream>
 #include <unsupported/Eigen/FFT>
 
 template <typename T> 
@@ -107,8 +106,6 @@ void test_scalar_generic(int nfft)
     for (int k=0;k<nfft;++k)
         tbuf[k]= (T)( rand()/(double)RAND_MAX - .5);
 
-    cout << "tbuf=["; for (size_t i=0;i<(size_t) tbuf.size();++i) {cout << tbuf[i] << " ";} cout << "];\n";
-
     // make sure it DOESN'T give the right full spectrum answer
     // if we've asked for half-spectrum
     fft.SetFlag(fft.HalfSpectrum );
@@ -125,9 +122,7 @@ void test_scalar_generic(int nfft)
         return; // odd FFTs get the wrong size inverse FFT
 
     ScalarVector tbuf2;
-    cout << "freqBuf=["; for (size_t i=0;i<(size_t) freqBuf.size();++i) {cout << freqBuf[i] << " ";} cout << "];\n";
     fft.inv( tbuf2 , freqBuf);
-    cout << "tbuf2=["; for (size_t i=0;i<(size_t) tbuf2.size();++i) {cout << tbuf2[i] << " ";} cout << "];\n";
     VERIFY( dif_rmse(tbuf,tbuf2) < test_precision<T>()  );// gross check
 
 
@@ -135,9 +130,7 @@ void test_scalar_generic(int nfft)
     ScalarVector tbuf3;
     fft.SetFlag(fft.Unscaled);
 
-    cout << "freqBuf=["; for (size_t i=0;i<(size_t) freqBuf.size();++i) {cout << freqBuf[i] << " ";} cout << "];\n";
     fft.inv( tbuf3 , freqBuf);
-    cout << "tbuf3=["; for (size_t i=0;i<(size_t) tbuf3.size();++i) {cout << tbuf3[i] << " ";} cout << "];\n";
 
     for (int k=0;k<nfft;++k)
         tbuf3[k] *= T(1./nfft);
@@ -146,8 +139,6 @@ void test_scalar_generic(int nfft)
     //for (size_t i=0;i<(size_t) tbuf.size();++i)
     //    cout << "freqBuf=" << freqBuf[i] << " in2=" << tbuf3[i] << " -  in=" << tbuf[i] << " => " << (tbuf3[i] - tbuf[i] ) <<  endl;
 
-    cout << "dif_rmse = " <<  dif_rmse(tbuf,tbuf3)  << endl;
-    cout << "test_precision = " << test_precision<T>()  << endl;
     VERIFY( dif_rmse(tbuf,tbuf3) < test_precision<T>()  );// gross check
 
     // verify that ClearFlag works
diff --git a/unsupported/test/NonLinearOptimization.cpp b/unsupported/test/NonLinearOptimization.cpp
index ae587f016..1313726a1 100644
--- a/unsupported/test/NonLinearOptimization.cpp
+++ b/unsupported/test/NonLinearOptimization.cpp
@@ -216,7 +216,7 @@ void testLmder()
 
   // check covariance
   covfac = fnorm*fnorm/(m-n);
-  ei_covar(lm.fjac, lm.ipvt); // TODO : move this as a function of lm
+  ei_covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm
 
   MatrixXd cov_ref(n,n);
   cov_ref <<
@@ -605,7 +605,7 @@ void testLmdif()
 
   // check covariance
   covfac = fnorm*fnorm/(m-n);
-  ei_covar(lm.fjac, lm.ipvt);
+  ei_covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm
 
   MatrixXd cov_ref(n,n);
   cov_ref <<
@@ -692,8 +692,8 @@ void testNistChwirut2(void)
   x<< 0.15, 0.008, 0.010;
   // do the computation
   lm.resetParameters();
-  lm.parameters.ftol = 1.E6*epsilon<double>();
-  lm.parameters.xtol = 1.E6*epsilon<double>();
+  lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
+  lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
   info = lm.minimize(x);
 
   // check return value
@@ -1010,7 +1010,7 @@ void testNistLanczos1(void)
   VERIFY( 79 == lm.nfev);
   VERIFY( 72 == lm.njev);
   // check norm^2
-  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.429604433690E-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
+  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 );
@@ -1031,7 +1031,7 @@ void testNistLanczos1(void)
   VERIFY( 9 == lm.nfev);
   VERIFY( 8 == lm.njev);
   // check norm^2
-  VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.43049947737308E-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
+  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 );
@@ -1170,9 +1170,9 @@ void testNistMGH10(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 2 == info);
-  VERIFY( 285 == lm.nfev);
-  VERIFY( 250 == lm.njev);
+  VERIFY( 2 == info); 
+  VERIFY( 284 == lm.nfev); 
+  VERIFY( 249 == lm.njev); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01);
   // check x
@@ -1188,7 +1188,7 @@ void testNistMGH10(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 2 == info);
+  VERIFY( 3 == info);
   VERIFY( 126 == lm.nfev);
   VERIFY( 116 == lm.njev);
   // check norm^2
@@ -1243,8 +1243,8 @@ void testNistBoxBOD(void)
   // do the computation
   BoxBOD_functor functor;
   LevenbergMarquardt<BoxBOD_functor> lm(functor);
-  lm.parameters.ftol = 1.E6*epsilon<double>();
-  lm.parameters.xtol = 1.E6*epsilon<double>();
+  lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
+  lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
   lm.parameters.factor = 10.;
   info = lm.minimize(x);
 
@@ -1264,14 +1264,14 @@ void testNistBoxBOD(void)
   x<< 100., 0.75;
   // do the computation
   lm.resetParameters();
-  lm.parameters.ftol = epsilon<double>();
-  lm.parameters.xtol = epsilon<double>();
+  lm.parameters.ftol = NumTraits<double>::epsilon();
+  lm.parameters.xtol = NumTraits<double>::epsilon();
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 15 == lm.nfev);
-  VERIFY( 14 == lm.njev);
+  VERIFY( 1 == info); 
+  VERIFY( 15 == lm.nfev); 
+  VERIFY( 14 == lm.njev); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03);
   // check x
@@ -1325,15 +1325,15 @@ void testNistMGH17(void)
   // do the computation
   MGH17_functor functor;
   LevenbergMarquardt<MGH17_functor> lm(functor);
-  lm.parameters.ftol = epsilon<double>();
-  lm.parameters.xtol = epsilon<double>();
+  lm.parameters.ftol = NumTraits<double>::epsilon();
+  lm.parameters.xtol = NumTraits<double>::epsilon();
   lm.parameters.maxfev = 1000;
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 599 == lm.nfev);
-  VERIFY( 544 == lm.njev);
+  VERIFY( 2 == info); 
+  VERIFY( 602 == lm.nfev); 
+  VERIFY( 545 == lm.njev); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05);
   // check x
@@ -1418,16 +1418,16 @@ void testNistMGH09(void)
   info = lm.minimize(x);
 
   // check return value
-  VERIFY( 1 == info);
-  VERIFY( 503== lm.nfev);
-  VERIFY( 385 == lm.njev);
+  VERIFY( 1 == info); 
+  VERIFY( 490 == lm.nfev); 
+  VERIFY( 376 == lm.njev); 
   // check norm^2
   VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04);
   // check x
-  VERIFY_IS_APPROX(x[0], 0.19280624); // should be 1.9280693458E-01
-  VERIFY_IS_APPROX(x[1], 0.19129774); // should be 1.9128232873E-01
-  VERIFY_IS_APPROX(x[2], 0.12305940); // should be 1.2305650693E-01
-  VERIFY_IS_APPROX(x[3], 0.13606946); // should be 1.3606233068E-01
+  VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01
+  VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01
+  VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01
+  VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01
 
   /*
    * Second try
@@ -1584,8 +1584,8 @@ void testNistThurber(void)
   // do the computation
   thurber_functor functor;
   LevenbergMarquardt<thurber_functor> lm(functor);
-  lm.parameters.ftol = 1.E4*epsilon<double>();
-  lm.parameters.xtol = 1.E4*epsilon<double>();
+  lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon();
+  lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon();
   info = lm.minimize(x);
 
   // check return value
@@ -1609,8 +1609,8 @@ void testNistThurber(void)
   x<< 1300 ,1500 ,500  ,75   ,1    ,0.4  ,0.05  ;
   // do the computation
   lm.resetParameters();
-  lm.parameters.ftol = 1.E4*epsilon<double>();
-  lm.parameters.xtol = 1.E4*epsilon<double>();
+  lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon();
+  lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon();
   info = lm.minimize(x);
 
   // check return value
@@ -1676,8 +1676,8 @@ void testNistRat43(void)
   // do the computation
   rat43_functor functor;
   LevenbergMarquardt<rat43_functor> lm(functor);
-  lm.parameters.ftol = 1.E6*epsilon<double>();
-  lm.parameters.xtol = 1.E6*epsilon<double>();
+  lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon();
+  lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon();
   info = lm.minimize(x);
 
   // check return value
@@ -1698,8 +1698,8 @@ void testNistRat43(void)
   x<< 700., 5., 0.75, 1.3;
   // do the computation
   lm.resetParameters();
-  lm.parameters.ftol = 1.E5*epsilon<double>();
-  lm.parameters.xtol = 1.E5*epsilon<double>();
+  lm.parameters.ftol = 1.E5*NumTraits<double>::epsilon();
+  lm.parameters.xtol = 1.E5*NumTraits<double>::epsilon();
   info = lm.minimize(x);
 
   // check return value
@@ -1833,7 +1833,6 @@ void test_NonLinearOptimization()
 
 /*
  * Can be useful for debugging...
-  printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev);
   printf("info, nfev : %d, %d\n", info, lm.nfev);
   printf("info, nfev, njev : %d, %d, %d\n", info, solver.nfev, solver.njev);
   printf("info, nfev : %d, %d\n", info, solver.nfev);
@@ -1843,5 +1842,14 @@ void test_NonLinearOptimization()
   printf("x[3] : %.32g\n", x[3]);
   printf("fvec.blueNorm() : %.32g\n", solver.fvec.blueNorm());
   printf("fvec.blueNorm() : %.32g\n", lm.fvec.blueNorm());
+
+  printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev);
+  printf("fvec.squaredNorm() : %.13g\n", lm.fvec.squaredNorm());
+  std::cout << x << std::endl;
+  std::cout.precision(9);
+  std::cout << x[0] << std::endl;
+  std::cout << x[1] << std::endl;
+  std::cout << x[2] << std::endl;
+  std::cout << x[3] << std::endl;
 */
 
diff --git a/unsupported/test/matrix_exponential.cpp b/unsupported/test/matrix_exponential.cpp
index a5b40adde..6150439c5 100644
--- a/unsupported/test/matrix_exponential.cpp
+++ b/unsupported/test/matrix_exponential.cpp
@@ -61,7 +61,7 @@ void test2dRotation(double tol)
     std::cout << "test2dRotation: i = " << i << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
-    ei_matrix_exponential(angle*A, &C);
+    C = ei_matrix_exponential(angle*A);
     std::cout << "   error expm = " << relerr(C, B) << "\n";
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
   }
@@ -86,7 +86,7 @@ void test2dHyperbolicRotation(double tol)
     std::cout << "test2dHyperbolicRotation: i = " << i << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
-    ei_matrix_exponential(A, &C);
+    C = ei_matrix_exponential(A);
     std::cout << "   error expm = " << relerr(C, B) << "\n";
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
   }
@@ -110,7 +110,7 @@ void testPascal(double tol)
     std::cout << "testPascal: size = " << size << "   error funm = " << relerr(C, B);
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
 
-    ei_matrix_exponential(A, &C);
+    C = ei_matrix_exponential(A);
     std::cout << "   error expm = " << relerr(C, B) << "\n";
     VERIFY(C.isApprox(B, static_cast<T>(tol)));
   }
@@ -137,10 +137,9 @@ void randomTest(const MatrixType& m, double tol)
     std::cout << "randomTest: error funm = " << relerr(identity, m2 * m3);
     VERIFY(identity.isApprox(m2 * m3, static_cast<RealScalar>(tol)));
 
-    ei_matrix_exponential(m1, &m2);
-    ei_matrix_exponential(-m1, &m3);
-    std::cout << "   error expm = " << relerr(identity, m2 * m3) << "\n";
-    VERIFY(identity.isApprox(m2 * m3, static_cast<RealScalar>(tol)));
+    m2 = ei_matrix_exponential(m1) * ei_matrix_exponential(-m1);
+    std::cout << "   error expm = " << relerr(identity, m2) << "\n";
+    VERIFY(identity.isApprox(m2, static_cast<RealScalar>(tol)));
   }
 }
 
diff --git a/unsupported/test/matrix_function.cpp b/unsupported/test/matrix_function.cpp
index 0eb30eecb..4ff6d7f1e 100644
--- a/unsupported/test/matrix_function.cpp
+++ b/unsupported/test/matrix_function.cpp
@@ -25,44 +25,100 @@
 #include "main.h"
 #include <unsupported/Eigen/MatrixFunctions>
 
+// Returns a matrix with eigenvalues clustered around 0, 1 and 2.
 template<typename MatrixType>
-void testMatrixExponential(const MatrixType& m)
+MatrixType randomMatrixWithRealEivals(const int size)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  MatrixType diag = MatrixType::Zero(size, size);
+  for (int i = 0; i < size; ++i) {
+    diag(i, i) = Scalar(RealScalar(ei_random<int>(0,2)))
+      + ei_random<Scalar>() * Scalar(RealScalar(0.01));
+  }
+  MatrixType A = MatrixType::Random(size, size);
+  return A.inverse() * diag * A;
+}
+
+template <typename MatrixType, int IsComplex = NumTraits<typename ei_traits<MatrixType>::Scalar>::IsComplex>
+struct randomMatrixWithImagEivals
+{
+  // Returns a matrix with eigenvalues clustered around 0 and +/- i.
+  static MatrixType run(const int size);
+};
+
+// Partial specialization for real matrices
+template<typename MatrixType>
+struct randomMatrixWithImagEivals<MatrixType, 0>
+{
+  static MatrixType run(const int size)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    MatrixType diag = MatrixType::Zero(size, size);
+    int i = 0;
+    while (i < size) {
+      int randomInt = ei_random<int>(-1, 1);
+      if (randomInt == 0 || i == size-1) {
+        diag(i, i) = ei_random<Scalar>() * Scalar(0.01);
+        ++i;
+      } else {
+        Scalar alpha = Scalar(randomInt) + ei_random<Scalar>() * Scalar(0.01);
+        diag(i, i+1) = alpha;
+        diag(i+1, i) = -alpha;
+        i += 2;
+      }
+    }
+    MatrixType A = MatrixType::Random(size, size);
+    return A.inverse() * diag * A;
+  }
+};
+
+// Partial specialization for complex matrices
+template<typename MatrixType>
+struct randomMatrixWithImagEivals<MatrixType, 1>
+{
+  static MatrixType run(const int size)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    const Scalar imagUnit(0, 1);
+    MatrixType diag = MatrixType::Zero(size, size);
+    for (int i = 0; i < size; ++i) {
+      diag(i, i) = Scalar(RealScalar(ei_random<int>(-1, 1))) * imagUnit
+        + ei_random<Scalar>() * Scalar(RealScalar(0.01));
+    }
+    MatrixType A = MatrixType::Random(size, size);
+    return A.inverse() * diag * A;
+  }
+};
+
+template<typename MatrixType>
+void testMatrixExponential(const MatrixType& A)
 {
   typedef typename ei_traits<MatrixType>::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef std::complex<RealScalar> ComplexScalar;
 
-  const int rows = m.rows();
-  const int cols = m.cols();
-
   for (int i = 0; i < g_repeat; i++) {
-    MatrixType A = MatrixType::Random(rows, cols);
-    MatrixType expA1, expA2;
-    ei_matrix_exponential(A, &expA1);
-    ei_matrix_function(A, StdStemFunctions<ComplexScalar>::exp, &expA2);
-    VERIFY_IS_APPROX(expA1, expA2);
+    VERIFY_IS_APPROX(ei_matrix_exponential(A), 
+		     ei_matrix_function(A, StdStemFunctions<ComplexScalar>::exp));
   }
 }
 
 template<typename MatrixType>
-void testHyperbolicFunctions(const MatrixType& m)
+void testHyperbolicFunctions(const MatrixType& A)
 {
-  const int rows = m.rows();
-  const int cols = m.cols();
-
   for (int i = 0; i < g_repeat; i++) {
-    MatrixType A = MatrixType::Random(rows, cols);
-    MatrixType sinhA, coshA, expA;
-    ei_matrix_sinh(A, &sinhA);
-    ei_matrix_cosh(A, &coshA);
-    ei_matrix_exponential(A, &expA);
+    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);
   }
 }
 
 template<typename MatrixType>
-void testGonioFunctions(const MatrixType& m)
+void testGonioFunctions(const MatrixType& A)
 {
   typedef ei_traits<MatrixType> Traits;
   typedef typename Traits::Scalar Scalar;
@@ -71,36 +127,44 @@ void testGonioFunctions(const MatrixType& m)
   typedef Matrix<ComplexScalar, Traits::RowsAtCompileTime, 
                  Traits::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;
 
-  const int rows = m.rows();
-  const int cols = m.cols();
   ComplexScalar imagUnit(0,1);
   ComplexScalar two(2,0);
 
   for (int i = 0; i < g_repeat; i++) {
-    MatrixType A = MatrixType::Random(rows, cols);
     ComplexMatrix Ac = A.template cast<ComplexScalar>();
 
-    ComplexMatrix exp_iA;
-    ei_matrix_exponential(imagUnit * Ac, &exp_iA);
+    ComplexMatrix exp_iA = ei_matrix_exponential(imagUnit * Ac);
 
-    MatrixType sinA;
-    ei_matrix_sin(A, &sinA);
+    MatrixType sinA = ei_matrix_sin(A);
     ComplexMatrix sinAc = sinA.template cast<ComplexScalar>();
     VERIFY_IS_APPROX(sinAc, (exp_iA - exp_iA.inverse()) / (two*imagUnit));
 
-    MatrixType cosA;
-    ei_matrix_cos(A, &cosA);
+    MatrixType cosA = ei_matrix_cos(A);
     ComplexMatrix cosAc = cosA.template cast<ComplexScalar>();
     VERIFY_IS_APPROX(cosAc, (exp_iA + exp_iA.inverse()) / 2);
   }
 }
 
 template<typename MatrixType>
+void testMatrix(const MatrixType& A)
+{
+  testMatrixExponential(A);
+  testHyperbolicFunctions(A);
+  testGonioFunctions(A);
+}
+
+template<typename MatrixType>
 void testMatrixType(const MatrixType& m)
 {
-  testMatrixExponential(m);
-  testHyperbolicFunctions(m);
-  testGonioFunctions(m);
+  // Matrices with clustered eigenvalue lead to different code paths
+  // in MatrixFunction.h and are thus useful for testing.
+
+  const int size = m.rows();
+  for (int i = 0; i < g_repeat; i++) {
+    testMatrix(MatrixType::Random(size, size).eval());
+    testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
+    testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
+  }
 }
 
 void test_matrix_function()