diff --git a/unsupported/doc/examples/BVH_Example.cpp b/unsupported/doc/examples/BVH_Example.cpp
new file mode 100644
index 0000000..6b6fac0
--- /dev/null
+++ b/unsupported/doc/examples/BVH_Example.cpp
@@ -0,0 +1,52 @@
+#include <Eigen/StdVector>
+#include <unsupported/Eigen/BVH>
+#include <iostream>
+
+using namespace Eigen;
+typedef AlignedBox<double, 2> Box2d;
+
+namespace Eigen {
+    namespace internal {
+        Box2d bounding_box(const Vector2d &v) { return Box2d(v, v); } //compute the bounding box of a single point
+    }
+}
+
+struct PointPointMinimizer //how to compute squared distances between points and rectangles
+{
+  PointPointMinimizer() : calls(0) {}
+  typedef double Scalar;
+
+  double minimumOnVolumeVolume(const Box2d &r1, const Box2d &r2) { ++calls; return r1.squaredExteriorDistance(r2); }
+  double minimumOnVolumeObject(const Box2d &r, const Vector2d &v) { ++calls; return r.squaredExteriorDistance(v); }
+  double minimumOnObjectVolume(const Vector2d &v, const Box2d &r) { ++calls; return r.squaredExteriorDistance(v); }
+  double minimumOnObjectObject(const Vector2d &v1, const Vector2d &v2) { ++calls; return (v1 - v2).squaredNorm(); }
+
+  int calls;
+};
+
+int main()
+{
+  typedef std::vector<Vector2d, aligned_allocator<Vector2d> > StdVectorOfVector2d;
+  StdVectorOfVector2d redPoints, bluePoints;
+  for(int i = 0; i < 100; ++i) { //initialize random set of red points and blue points
+    redPoints.push_back(Vector2d::Random());
+    bluePoints.push_back(Vector2d::Random());
+  }
+
+  PointPointMinimizer minimizer;
+  double minDistSq = std::numeric_limits<double>::max();
+
+  //brute force to find closest red-blue pair
+  for(int i = 0; i < (int)redPoints.size(); ++i)
+    for(int j = 0; j < (int)bluePoints.size(); ++j)
+      minDistSq = std::min(minDistSq, minimizer.minimumOnObjectObject(redPoints[i], bluePoints[j]));
+  std::cout << "Brute force distance = " << sqrt(minDistSq) << ", calls = " << minimizer.calls << std::endl;
+
+  //using BVH to find closest red-blue pair
+  minimizer.calls = 0;
+  KdBVH<double, 2, Vector2d> redTree(redPoints.begin(), redPoints.end()), blueTree(bluePoints.begin(), bluePoints.end()); //construct the trees
+  minDistSq = BVMinimize(redTree, blueTree, minimizer); //actual BVH minimization call
+  std::cout << "BVH distance         = " << sqrt(minDistSq) << ", calls = " << minimizer.calls << std::endl;
+
+  return 0;
+}
diff --git a/unsupported/doc/examples/CMakeLists.txt b/unsupported/doc/examples/CMakeLists.txt
new file mode 100644
index 0000000..c47646d
--- /dev/null
+++ b/unsupported/doc/examples/CMakeLists.txt
@@ -0,0 +1,20 @@
+FILE(GLOB examples_SRCS "*.cpp")
+
+ADD_CUSTOM_TARGET(unsupported_examples)
+
+INCLUDE_DIRECTORIES(../../../unsupported ../../../unsupported/test)
+
+FOREACH(example_src ${examples_SRCS})
+  GET_FILENAME_COMPONENT(example ${example_src} NAME_WE)
+  ADD_EXECUTABLE(example_${example} ${example_src})
+  if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
+    target_link_libraries(example_${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
+  endif()
+  ADD_CUSTOM_COMMAND(
+    TARGET example_${example}
+    POST_BUILD
+    COMMAND example_${example}
+    ARGS >${CMAKE_CURRENT_BINARY_DIR}/${example}.out
+  )
+  ADD_DEPENDENCIES(unsupported_examples example_${example})
+ENDFOREACH(example_src)
diff --git a/unsupported/doc/examples/FFT.cpp b/unsupported/doc/examples/FFT.cpp
new file mode 100644
index 0000000..fcbf812
--- /dev/null
+++ b/unsupported/doc/examples/FFT.cpp
@@ -0,0 +1,118 @@
+//  To use the simple FFT implementation
+//  g++ -o demofft -I.. -Wall -O3 FFT.cpp 
+
+//  To use the FFTW implementation
+//  g++ -o demofft -I.. -DUSE_FFTW -Wall -O3 FFT.cpp -lfftw3 -lfftw3f -lfftw3l
+
+#ifdef USE_FFTW
+#include <fftw3.h>
+#endif
+
+#include <vector>
+#include <complex>
+#include <algorithm>
+#include <iterator>
+#include <iostream>
+#include <Eigen/Core>
+#include <unsupported/Eigen/FFT>
+
+using namespace std;
+using namespace Eigen;
+
+template <typename T>
+T mag2(T a)
+{
+    return a*a;
+}
+template <typename T>
+T mag2(std::complex<T> a)
+{
+    return norm(a);
+}
+
+template <typename T>
+T mag2(const std::vector<T> & vec)
+{
+    T out=0;
+    for (size_t k=0;k<vec.size();++k)
+        out += mag2(vec[k]);
+    return out;
+}
+
+template <typename T>
+T mag2(const std::vector<std::complex<T> > & vec)
+{
+    T out=0;
+    for (size_t k=0;k<vec.size();++k)
+        out += mag2(vec[k]);
+    return out;
+}
+
+template <typename T>
+vector<T> operator-(const vector<T> & a,const vector<T> & b )
+{
+    vector<T> c(a);
+    for (size_t k=0;k<b.size();++k) 
+        c[k] -= b[k];
+    return c;
+}
+
+template <typename T>
+void RandomFill(std::vector<T> & vec)
+{
+    for (size_t k=0;k<vec.size();++k)
+        vec[k] = T( rand() )/T(RAND_MAX) - .5;
+}
+
+template <typename T>
+void RandomFill(std::vector<std::complex<T> > & vec)
+{
+    for (size_t k=0;k<vec.size();++k)
+        vec[k] = std::complex<T> ( T( rand() )/T(RAND_MAX) - .5, T( rand() )/T(RAND_MAX) - .5);
+}
+
+template <typename T_time,typename T_freq>
+void fwd_inv(size_t nfft)
+{
+    typedef typename NumTraits<T_freq>::Real Scalar;
+    vector<T_time> timebuf(nfft);
+    RandomFill(timebuf);
+
+    vector<T_freq> freqbuf;
+    static FFT<Scalar> fft;
+    fft.fwd(freqbuf,timebuf);
+
+    vector<T_time> timebuf2;
+    fft.inv(timebuf2,freqbuf);
+
+    long double rmse = mag2(timebuf - timebuf2) / mag2(timebuf);
+    cout << "roundtrip rmse: " << rmse << endl;
+}
+
+template <typename T_scalar>
+void two_demos(int nfft)
+{
+    cout << "     scalar ";
+    fwd_inv<T_scalar,std::complex<T_scalar> >(nfft);
+    cout << "    complex ";
+    fwd_inv<std::complex<T_scalar>,std::complex<T_scalar> >(nfft);
+}
+
+void demo_all_types(int nfft)
+{
+    cout << "nfft=" << nfft << endl;
+    cout << "   float" << endl;
+    two_demos<float>(nfft);
+    cout << "   double" << endl;
+    two_demos<double>(nfft);
+    cout << "   long double" << endl;
+    two_demos<long double>(nfft);
+}
+
+int main()
+{
+    demo_all_types( 2*3*4*5*7 );
+    demo_all_types( 2*9*16*25 );
+    demo_all_types( 1024 );
+    return 0;
+}
diff --git a/unsupported/doc/examples/MatrixExponential.cpp b/unsupported/doc/examples/MatrixExponential.cpp
new file mode 100644
index 0000000..ebd3b96
--- /dev/null
+++ b/unsupported/doc/examples/MatrixExponential.cpp
@@ -0,0 +1,16 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  const double pi = std::acos(-1.0);
+
+  MatrixXd A(3,3);
+  A << 0,    -pi/4, 0,
+       pi/4, 0,     0,
+       0,    0,     0;
+  std::cout << "The matrix A is:\n" << A << "\n\n";
+  std::cout << "The matrix exponential of A is:\n" << A.exp() << "\n\n";
+}
diff --git a/unsupported/doc/examples/MatrixFunction.cpp b/unsupported/doc/examples/MatrixFunction.cpp
new file mode 100644
index 0000000..a4172e4
--- /dev/null
+++ b/unsupported/doc/examples/MatrixFunction.cpp
@@ -0,0 +1,23 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+std::complex<double> expfn(std::complex<double> x, int)
+{
+  return std::exp(x);
+}
+
+int main()
+{
+  const double pi = std::acos(-1.0);
+
+  MatrixXd A(3,3);
+  A << 0,    -pi/4, 0,
+       pi/4, 0,     0,
+       0,    0,     0;
+
+  std::cout << "The matrix A is:\n" << A << "\n\n";
+  std::cout << "The matrix exponential of A is:\n" 
+            << A.matrixFunction(expfn) << "\n\n";
+}
diff --git a/unsupported/doc/examples/MatrixLogarithm.cpp b/unsupported/doc/examples/MatrixLogarithm.cpp
new file mode 100644
index 0000000..8c5d970
--- /dev/null
+++ b/unsupported/doc/examples/MatrixLogarithm.cpp
@@ -0,0 +1,15 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  using std::sqrt;
+  MatrixXd A(3,3);
+  A << 0.5*sqrt(2), -0.5*sqrt(2), 0,
+       0.5*sqrt(2),  0.5*sqrt(2), 0,
+       0,            0,           1;
+  std::cout << "The matrix A is:\n" << A << "\n\n";
+  std::cout << "The matrix logarithm of A is:\n" << A.log() << "\n";
+}
diff --git a/unsupported/doc/examples/MatrixPower.cpp b/unsupported/doc/examples/MatrixPower.cpp
new file mode 100644
index 0000000..2224524
--- /dev/null
+++ b/unsupported/doc/examples/MatrixPower.cpp
@@ -0,0 +1,16 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  const double pi = std::acos(-1.0);
+  Matrix3d A;
+  A << cos(1), -sin(1), 0,
+       sin(1),  cos(1), 0,
+	   0 ,      0 , 1;
+  std::cout << "The matrix A is:\n" << A << "\n\n"
+	       "The matrix power A^(pi/4) is:\n" << A.pow(pi/4) << std::endl;
+  return 0;
+}
diff --git a/unsupported/doc/examples/MatrixPower_optimal.cpp b/unsupported/doc/examples/MatrixPower_optimal.cpp
new file mode 100644
index 0000000..86470ba
--- /dev/null
+++ b/unsupported/doc/examples/MatrixPower_optimal.cpp
@@ -0,0 +1,17 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  Matrix4cd A = Matrix4cd::Random();
+  MatrixPower<Matrix4cd> Apow(A);
+
+  std::cout << "The matrix A is:\n" << A << "\n\n"
+	       "A^3.1 is:\n" << Apow(3.1) << "\n\n"
+	       "A^3.3 is:\n" << Apow(3.3) << "\n\n"
+	       "A^3.7 is:\n" << Apow(3.7) << "\n\n"
+	       "A^3.9 is:\n" << Apow(3.9) << std::endl;
+  return 0;
+}
diff --git a/unsupported/doc/examples/MatrixSine.cpp b/unsupported/doc/examples/MatrixSine.cpp
new file mode 100644
index 0000000..9eea9a0
--- /dev/null
+++ b/unsupported/doc/examples/MatrixSine.cpp
@@ -0,0 +1,20 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  MatrixXd A = MatrixXd::Random(3,3);
+  std::cout << "A = \n" << A << "\n\n";
+
+  MatrixXd sinA = A.sin();
+  std::cout << "sin(A) = \n" << sinA << "\n\n";
+
+  MatrixXd cosA = A.cos();
+  std::cout << "cos(A) = \n" << cosA << "\n\n";
+  
+  // The matrix functions satisfy sin^2(A) + cos^2(A) = I, 
+  // like the scalar functions.
+  std::cout << "sin^2(A) + cos^2(A) = \n" << sinA*sinA + cosA*cosA << "\n\n";
+}
diff --git a/unsupported/doc/examples/MatrixSinh.cpp b/unsupported/doc/examples/MatrixSinh.cpp
new file mode 100644
index 0000000..f771867
--- /dev/null
+++ b/unsupported/doc/examples/MatrixSinh.cpp
@@ -0,0 +1,20 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  MatrixXf A = MatrixXf::Random(3,3);
+  std::cout << "A = \n" << A << "\n\n";
+
+  MatrixXf sinhA = A.sinh();
+  std::cout << "sinh(A) = \n" << sinhA << "\n\n";
+
+  MatrixXf coshA = A.cosh();
+  std::cout << "cosh(A) = \n" << coshA << "\n\n";
+  
+  // The matrix functions satisfy cosh^2(A) - sinh^2(A) = I, 
+  // like the scalar functions.
+  std::cout << "cosh^2(A) - sinh^2(A) = \n" << coshA*coshA - sinhA*sinhA << "\n\n";
+}
diff --git a/unsupported/doc/examples/MatrixSquareRoot.cpp b/unsupported/doc/examples/MatrixSquareRoot.cpp
new file mode 100644
index 0000000..88e7557
--- /dev/null
+++ b/unsupported/doc/examples/MatrixSquareRoot.cpp
@@ -0,0 +1,16 @@
+#include <unsupported/Eigen/MatrixFunctions>
+#include <iostream>
+
+using namespace Eigen;
+
+int main()
+{
+  const double pi = std::acos(-1.0);
+
+  MatrixXd A(2,2);
+  A << cos(pi/3), -sin(pi/3), 
+       sin(pi/3),  cos(pi/3);
+  std::cout << "The matrix A is:\n" << A << "\n\n";
+  std::cout << "The matrix square root of A is:\n" << A.sqrt() << "\n\n";
+  std::cout << "The square of the last matrix is:\n" << A.sqrt() * A.sqrt() << "\n";
+}
diff --git a/unsupported/doc/examples/PolynomialSolver1.cpp b/unsupported/doc/examples/PolynomialSolver1.cpp
new file mode 100644
index 0000000..cd777a4
--- /dev/null
+++ b/unsupported/doc/examples/PolynomialSolver1.cpp
@@ -0,0 +1,53 @@
+#include <unsupported/Eigen/Polynomials>
+#include <vector>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+  typedef Matrix<double,5,1> Vector5d;
+
+  Vector5d roots = Vector5d::Random();
+  cout << "Roots: " << roots.transpose() << endl;
+  Eigen::Matrix<double,6,1> polynomial;
+  roots_to_monicPolynomial( roots, polynomial );
+
+  PolynomialSolver<double,5> psolve( polynomial );
+  cout << "Complex roots: " << psolve.roots().transpose() << endl;
+
+  std::vector<double> realRoots;
+  psolve.realRoots( realRoots );
+  Map<Vector5d> mapRR( &realRoots[0] );
+  cout << "Real roots: " << mapRR.transpose() << endl;
+
+  cout << endl;
+  cout << "Illustration of the convergence problem with the QR algorithm: " << endl;
+  cout << "---------------------------------------------------------------" << endl;
+  Eigen::Matrix<float,7,1> hardCase_polynomial;
+  hardCase_polynomial <<
+  -0.957, 0.9219, 0.3516, 0.9453, -0.4023, -0.5508, -0.03125;
+  cout << "Hard case polynomial defined by floats: " << hardCase_polynomial.transpose() << endl;
+  PolynomialSolver<float,6> psolvef( hardCase_polynomial );
+  cout << "Complex roots: " << psolvef.roots().transpose() << endl;
+  Eigen::Matrix<float,6,1> evals;
+  for( int i=0; i<6; ++i ){ evals[i] = std::abs( poly_eval( hardCase_polynomial, psolvef.roots()[i] ) ); }
+  cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
+
+  cout << "Using double's almost always solves the problem for small degrees: " << endl;
+  cout << "-------------------------------------------------------------------" << endl;
+  PolynomialSolver<double,6> psolve6d( hardCase_polynomial.cast<double>() );
+  cout << "Complex roots: " << psolve6d.roots().transpose() << endl;
+  for( int i=0; i<6; ++i )
+  {
+    std::complex<float> castedRoot( psolve6d.roots()[i].real(), psolve6d.roots()[i].imag() );
+    evals[i] = std::abs( poly_eval( hardCase_polynomial, castedRoot ) );
+  }
+  cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
+
+  cout.precision(10);
+  cout << "The last root in float then in double: " << psolvef.roots()[5] << "\t" << psolve6d.roots()[5] << endl;
+  std::complex<float> castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() );
+  cout << "Norm of the difference: " << std::abs( psolvef.roots()[5] - castedRoot ) << endl;
+}
diff --git a/unsupported/doc/examples/PolynomialUtils1.cpp b/unsupported/doc/examples/PolynomialUtils1.cpp
new file mode 100644
index 0000000..dbfe520
--- /dev/null
+++ b/unsupported/doc/examples/PolynomialUtils1.cpp
@@ -0,0 +1,20 @@
+#include <unsupported/Eigen/Polynomials>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+  Vector4d roots = Vector4d::Random();
+  cout << "Roots: " << roots.transpose() << endl;
+  Eigen::Matrix<double,5,1> polynomial;
+  roots_to_monicPolynomial( roots, polynomial );
+  cout << "Polynomial: ";
+  for( int i=0; i<4; ++i ){ cout << polynomial[i] << ".x^" << i << "+ "; }
+  cout << polynomial[4] << ".x^4" << endl;
+  Vector4d evaluation;
+  for( int i=0; i<4; ++i ){
+    evaluation[i] = poly_eval( polynomial, roots[i] ); }
+  cout << "Evaluation of the polynomial at the roots: " << evaluation.transpose();
+}