diff --git a/doc/snippets/.krazy b/doc/snippets/.krazy
new file mode 100644
index 0000000..00b9940
--- /dev/null
+++ b/doc/snippets/.krazy
@@ -0,0 +1,2 @@
+EXCLUDE copyright
+EXCLUDE license
diff --git a/doc/snippets/AngleAxis_mimic_euler.cpp b/doc/snippets/AngleAxis_mimic_euler.cpp
new file mode 100644
index 0000000..456de7f
--- /dev/null
+++ b/doc/snippets/AngleAxis_mimic_euler.cpp
@@ -0,0 +1,5 @@
+Matrix3f m;
+m = AngleAxisf(0.25*M_PI, Vector3f::UnitX())
+  * AngleAxisf(0.5*M_PI,  Vector3f::UnitY())
+  * AngleAxisf(0.33*M_PI, Vector3f::UnitZ());
+cout << m << endl << "is unitary: " << m.isUnitary() << endl;
diff --git a/doc/snippets/CMakeLists.txt b/doc/snippets/CMakeLists.txt
new file mode 100644
index 0000000..1135900
--- /dev/null
+++ b/doc/snippets/CMakeLists.txt
@@ -0,0 +1,28 @@
+file(GLOB snippets_SRCS "*.cpp")
+
+add_custom_target(all_snippets)
+
+foreach(snippet_src ${snippets_SRCS})
+  get_filename_component(snippet ${snippet_src} NAME_WE)
+  set(compile_snippet_target compile_${snippet})
+  set(compile_snippet_src ${compile_snippet_target}.cpp)
+  file(READ ${snippet_src} snippet_source_code)
+  configure_file(${CMAKE_CURRENT_SOURCE_DIR}/compile_snippet.cpp.in
+                 ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src})
+  add_executable(${compile_snippet_target}
+                 ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src})
+  if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
+    target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
+  endif()
+  add_custom_command(
+    TARGET ${compile_snippet_target}
+    POST_BUILD
+    COMMAND ${compile_snippet_target}
+    ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out
+  )
+  add_dependencies(all_snippets ${compile_snippet_target})
+  set_source_files_properties(${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src}
+                              PROPERTIES OBJECT_DEPENDS ${snippet_src})
+endforeach(snippet_src)
+
+ei_add_target_property(compile_tut_arithmetic_transpose_aliasing COMPILE_FLAGS -DEIGEN_NO_DEBUG)
diff --git a/doc/snippets/ColPivHouseholderQR_solve.cpp b/doc/snippets/ColPivHouseholderQR_solve.cpp
new file mode 100644
index 0000000..b7b204a
--- /dev/null
+++ b/doc/snippets/ColPivHouseholderQR_solve.cpp
@@ -0,0 +1,8 @@
+Matrix3f m = Matrix3f::Random();
+Matrix3f y = Matrix3f::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the matrix y:" << endl << y << endl;
+Matrix3f x;
+x = m.colPivHouseholderQr().solve(y);
+assert(y.isApprox(m*x));
+cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
diff --git a/doc/snippets/ComplexEigenSolver_compute.cpp b/doc/snippets/ComplexEigenSolver_compute.cpp
new file mode 100644
index 0000000..11d6bd3
--- /dev/null
+++ b/doc/snippets/ComplexEigenSolver_compute.cpp
@@ -0,0 +1,16 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+cout << "Here is a random 4x4 matrix, A:" << endl << A << endl << endl;
+
+ComplexEigenSolver<MatrixXcf> ces;
+ces.compute(A);
+cout << "The eigenvalues of A are:" << endl << ces.eigenvalues() << endl;
+cout << "The matrix of eigenvectors, V, is:" << endl << ces.eigenvectors() << endl << endl;
+
+complex<float> lambda = ces.eigenvalues()[0];
+cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
+VectorXcf v = ces.eigenvectors().col(0);
+cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
+cout << "... and A * v = " << endl << A * v << endl << endl;
+
+cout << "Finally, V * D * V^(-1) = " << endl
+     << ces.eigenvectors() * ces.eigenvalues().asDiagonal() * ces.eigenvectors().inverse() << endl;
diff --git a/doc/snippets/ComplexEigenSolver_eigenvalues.cpp b/doc/snippets/ComplexEigenSolver_eigenvalues.cpp
new file mode 100644
index 0000000..5509bd8
--- /dev/null
+++ b/doc/snippets/ComplexEigenSolver_eigenvalues.cpp
@@ -0,0 +1,4 @@
+MatrixXcf ones = MatrixXcf::Ones(3,3);
+ComplexEigenSolver<MatrixXcf> ces(ones, /* computeEigenvectors = */ false);
+cout << "The eigenvalues of the 3x3 matrix of ones are:" 
+     << endl << ces.eigenvalues() << endl;
diff --git a/doc/snippets/ComplexEigenSolver_eigenvectors.cpp b/doc/snippets/ComplexEigenSolver_eigenvectors.cpp
new file mode 100644
index 0000000..bb1c2cc
--- /dev/null
+++ b/doc/snippets/ComplexEigenSolver_eigenvectors.cpp
@@ -0,0 +1,4 @@
+MatrixXcf ones = MatrixXcf::Ones(3,3);
+ComplexEigenSolver<MatrixXcf> ces(ones);
+cout << "The first eigenvector of the 3x3 matrix of ones is:" 
+     << endl << ces.eigenvectors().col(1) << endl;
diff --git a/doc/snippets/ComplexSchur_compute.cpp b/doc/snippets/ComplexSchur_compute.cpp
new file mode 100644
index 0000000..3a51701
--- /dev/null
+++ b/doc/snippets/ComplexSchur_compute.cpp
@@ -0,0 +1,6 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+ComplexSchur<MatrixXcf> schur(4);
+schur.compute(A);
+cout << "The matrix T in the decomposition of A is:" << endl << schur.matrixT() << endl;
+schur.compute(A.inverse());
+cout << "The matrix T in the decomposition of A^(-1) is:" << endl << schur.matrixT() << endl;
diff --git a/doc/snippets/ComplexSchur_matrixT.cpp b/doc/snippets/ComplexSchur_matrixT.cpp
new file mode 100644
index 0000000..8380571
--- /dev/null
+++ b/doc/snippets/ComplexSchur_matrixT.cpp
@@ -0,0 +1,4 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+cout << "Here is a random 4x4 matrix, A:" << endl << A << endl << endl;
+ComplexSchur<MatrixXcf> schurOfA(A, false); // false means do not compute U
+cout << "The triangular matrix T is:" << endl << schurOfA.matrixT() << endl;
diff --git a/doc/snippets/ComplexSchur_matrixU.cpp b/doc/snippets/ComplexSchur_matrixU.cpp
new file mode 100644
index 0000000..ba3d9c2
--- /dev/null
+++ b/doc/snippets/ComplexSchur_matrixU.cpp
@@ -0,0 +1,4 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+cout << "Here is a random 4x4 matrix, A:" << endl << A << endl << endl;
+ComplexSchur<MatrixXcf> schurOfA(A);
+cout << "The unitary matrix U is:" << endl << schurOfA.matrixU() << endl;
diff --git a/doc/snippets/Cwise_abs.cpp b/doc/snippets/Cwise_abs.cpp
new file mode 100644
index 0000000..0aeec3a
--- /dev/null
+++ b/doc/snippets/Cwise_abs.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,-2,-3);
+cout << v.abs() << endl;
diff --git a/doc/snippets/Cwise_abs2.cpp b/doc/snippets/Cwise_abs2.cpp
new file mode 100644
index 0000000..2c4f9b3
--- /dev/null
+++ b/doc/snippets/Cwise_abs2.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,-2,-3);
+cout << v.abs2() << endl;
diff --git a/doc/snippets/Cwise_acos.cpp b/doc/snippets/Cwise_acos.cpp
new file mode 100644
index 0000000..34432cb
--- /dev/null
+++ b/doc/snippets/Cwise_acos.cpp
@@ -0,0 +1,2 @@
+Array3d v(0, sqrt(2.)/2, 1);
+cout << v.acos() << endl;
diff --git a/doc/snippets/Cwise_asin.cpp b/doc/snippets/Cwise_asin.cpp
new file mode 100644
index 0000000..8dad838
--- /dev/null
+++ b/doc/snippets/Cwise_asin.cpp
@@ -0,0 +1,2 @@
+Array3d v(0, sqrt(2.)/2, 1);
+cout << v.asin() << endl;
diff --git a/doc/snippets/Cwise_boolean_and.cpp b/doc/snippets/Cwise_boolean_and.cpp
new file mode 100644
index 0000000..df6b60d
--- /dev/null
+++ b/doc/snippets/Cwise_boolean_and.cpp
@@ -0,0 +1,2 @@
+Array3d v(-1,2,1), w(-3,2,3);
+cout << ((v<w) && (v<0)) << endl;
diff --git a/doc/snippets/Cwise_boolean_or.cpp b/doc/snippets/Cwise_boolean_or.cpp
new file mode 100644
index 0000000..83eb006
--- /dev/null
+++ b/doc/snippets/Cwise_boolean_or.cpp
@@ -0,0 +1,2 @@
+Array3d v(-1,2,1), w(-3,2,3);
+cout << ((v<w) || (v<0)) << endl;
diff --git a/doc/snippets/Cwise_cos.cpp b/doc/snippets/Cwise_cos.cpp
new file mode 100644
index 0000000..f589f07
--- /dev/null
+++ b/doc/snippets/Cwise_cos.cpp
@@ -0,0 +1,2 @@
+Array3d v(M_PI, M_PI/2, M_PI/3);
+cout << v.cos() << endl;
diff --git a/doc/snippets/Cwise_cube.cpp b/doc/snippets/Cwise_cube.cpp
new file mode 100644
index 0000000..85e41dc
--- /dev/null
+++ b/doc/snippets/Cwise_cube.cpp
@@ -0,0 +1,2 @@
+Array3d v(2,3,4);
+cout << v.cube() << endl;
diff --git a/doc/snippets/Cwise_equal_equal.cpp b/doc/snippets/Cwise_equal_equal.cpp
new file mode 100644
index 0000000..0ba96f6
--- /dev/null
+++ b/doc/snippets/Cwise_equal_equal.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3), w(3,2,1);
+cout << (v==w) << endl;
diff --git a/doc/snippets/Cwise_exp.cpp b/doc/snippets/Cwise_exp.cpp
new file mode 100644
index 0000000..db23618
--- /dev/null
+++ b/doc/snippets/Cwise_exp.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3);
+cout << v.exp() << endl;
diff --git a/doc/snippets/Cwise_greater.cpp b/doc/snippets/Cwise_greater.cpp
new file mode 100644
index 0000000..40ad029
--- /dev/null
+++ b/doc/snippets/Cwise_greater.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3), w(3,2,1);
+cout << (v>w) << endl;
diff --git a/doc/snippets/Cwise_greater_equal.cpp b/doc/snippets/Cwise_greater_equal.cpp
new file mode 100644
index 0000000..6a08f89
--- /dev/null
+++ b/doc/snippets/Cwise_greater_equal.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3), w(3,2,1);
+cout << (v>=w) << endl;
diff --git a/doc/snippets/Cwise_inverse.cpp b/doc/snippets/Cwise_inverse.cpp
new file mode 100644
index 0000000..3967a7e
--- /dev/null
+++ b/doc/snippets/Cwise_inverse.cpp
@@ -0,0 +1,2 @@
+Array3d v(2,3,4);
+cout << v.inverse() << endl;
diff --git a/doc/snippets/Cwise_less.cpp b/doc/snippets/Cwise_less.cpp
new file mode 100644
index 0000000..cafd3b6
--- /dev/null
+++ b/doc/snippets/Cwise_less.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3), w(3,2,1);
+cout << (v<w) << endl;
diff --git a/doc/snippets/Cwise_less_equal.cpp b/doc/snippets/Cwise_less_equal.cpp
new file mode 100644
index 0000000..1600e39
--- /dev/null
+++ b/doc/snippets/Cwise_less_equal.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3), w(3,2,1);
+cout << (v<=w) << endl;
diff --git a/doc/snippets/Cwise_log.cpp b/doc/snippets/Cwise_log.cpp
new file mode 100644
index 0000000..f7aca72
--- /dev/null
+++ b/doc/snippets/Cwise_log.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3);
+cout << v.log() << endl;
diff --git a/doc/snippets/Cwise_max.cpp b/doc/snippets/Cwise_max.cpp
new file mode 100644
index 0000000..6602881
--- /dev/null
+++ b/doc/snippets/Cwise_max.cpp
@@ -0,0 +1,2 @@
+Array3d v(2,3,4), w(4,2,3);
+cout << v.max(w) << endl;
diff --git a/doc/snippets/Cwise_min.cpp b/doc/snippets/Cwise_min.cpp
new file mode 100644
index 0000000..1c01c76
--- /dev/null
+++ b/doc/snippets/Cwise_min.cpp
@@ -0,0 +1,2 @@
+Array3d v(2,3,4), w(4,2,3);
+cout << v.min(w) << endl;
diff --git a/doc/snippets/Cwise_minus.cpp b/doc/snippets/Cwise_minus.cpp
new file mode 100644
index 0000000..b89b9fb
--- /dev/null
+++ b/doc/snippets/Cwise_minus.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3);
+cout << v-5 << endl;
diff --git a/doc/snippets/Cwise_minus_equal.cpp b/doc/snippets/Cwise_minus_equal.cpp
new file mode 100644
index 0000000..dfde49d
--- /dev/null
+++ b/doc/snippets/Cwise_minus_equal.cpp
@@ -0,0 +1,3 @@
+Array3d v(1,2,3);
+v -= 5;
+cout << v << endl;
diff --git a/doc/snippets/Cwise_not_equal.cpp b/doc/snippets/Cwise_not_equal.cpp
new file mode 100644
index 0000000..57a407a
--- /dev/null
+++ b/doc/snippets/Cwise_not_equal.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3), w(3,2,1);
+cout << (v!=w) << endl;
diff --git a/doc/snippets/Cwise_plus.cpp b/doc/snippets/Cwise_plus.cpp
new file mode 100644
index 0000000..9d47327
--- /dev/null
+++ b/doc/snippets/Cwise_plus.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,3);
+cout << v+5 << endl;
diff --git a/doc/snippets/Cwise_plus_equal.cpp b/doc/snippets/Cwise_plus_equal.cpp
new file mode 100644
index 0000000..d744b1e
--- /dev/null
+++ b/doc/snippets/Cwise_plus_equal.cpp
@@ -0,0 +1,3 @@
+Array3d v(1,2,3);
+v += 5;
+cout << v << endl;
diff --git a/doc/snippets/Cwise_pow.cpp b/doc/snippets/Cwise_pow.cpp
new file mode 100644
index 0000000..a723ed8
--- /dev/null
+++ b/doc/snippets/Cwise_pow.cpp
@@ -0,0 +1,2 @@
+Array3d v(8,27,64);
+cout << v.pow(0.333333) << endl;
diff --git a/doc/snippets/Cwise_product.cpp b/doc/snippets/Cwise_product.cpp
new file mode 100644
index 0000000..714d66d
--- /dev/null
+++ b/doc/snippets/Cwise_product.cpp
@@ -0,0 +1,4 @@
+Array33i a = Array33i::Random(), b = Array33i::Random();
+Array33i c = a * b;
+cout << "a:\n" << a << "\nb:\n" << b << "\nc:\n" << c << endl;
+
diff --git a/doc/snippets/Cwise_quotient.cpp b/doc/snippets/Cwise_quotient.cpp
new file mode 100644
index 0000000..7cb9f7f
--- /dev/null
+++ b/doc/snippets/Cwise_quotient.cpp
@@ -0,0 +1,2 @@
+Array3d v(2,3,4), w(4,2,3);
+cout << v/w << endl;
diff --git a/doc/snippets/Cwise_sin.cpp b/doc/snippets/Cwise_sin.cpp
new file mode 100644
index 0000000..46fa908
--- /dev/null
+++ b/doc/snippets/Cwise_sin.cpp
@@ -0,0 +1,2 @@
+Array3d v(M_PI, M_PI/2, M_PI/3);
+cout << v.sin() << endl;
diff --git a/doc/snippets/Cwise_slash_equal.cpp b/doc/snippets/Cwise_slash_equal.cpp
new file mode 100644
index 0000000..2efd32d
--- /dev/null
+++ b/doc/snippets/Cwise_slash_equal.cpp
@@ -0,0 +1,3 @@
+Array3d v(3,2,4), w(5,4,2);
+v /= w;
+cout << v << endl;
diff --git a/doc/snippets/Cwise_sqrt.cpp b/doc/snippets/Cwise_sqrt.cpp
new file mode 100644
index 0000000..97bafe8
--- /dev/null
+++ b/doc/snippets/Cwise_sqrt.cpp
@@ -0,0 +1,2 @@
+Array3d v(1,2,4);
+cout << v.sqrt() << endl;
diff --git a/doc/snippets/Cwise_square.cpp b/doc/snippets/Cwise_square.cpp
new file mode 100644
index 0000000..f704c5e
--- /dev/null
+++ b/doc/snippets/Cwise_square.cpp
@@ -0,0 +1,2 @@
+Array3d v(2,3,4);
+cout << v.square() << endl;
diff --git a/doc/snippets/Cwise_tan.cpp b/doc/snippets/Cwise_tan.cpp
new file mode 100644
index 0000000..b758ef0
--- /dev/null
+++ b/doc/snippets/Cwise_tan.cpp
@@ -0,0 +1,2 @@
+Array3d v(M_PI, M_PI/2, M_PI/3);
+cout << v.tan() << endl;
diff --git a/doc/snippets/Cwise_times_equal.cpp b/doc/snippets/Cwise_times_equal.cpp
new file mode 100644
index 0000000..147556c
--- /dev/null
+++ b/doc/snippets/Cwise_times_equal.cpp
@@ -0,0 +1,3 @@
+Array3d v(1,2,3), w(2,3,0);
+v *= w;
+cout << v << endl;
diff --git a/doc/snippets/DenseBase_LinSpaced.cpp b/doc/snippets/DenseBase_LinSpaced.cpp
new file mode 100644
index 0000000..8e54b17
--- /dev/null
+++ b/doc/snippets/DenseBase_LinSpaced.cpp
@@ -0,0 +1,2 @@
+cout << VectorXi::LinSpaced(4,7,10).transpose() << endl;
+cout << VectorXd::LinSpaced(5,0.0,1.0).transpose() << endl;
diff --git a/doc/snippets/DenseBase_LinSpaced_seq.cpp b/doc/snippets/DenseBase_LinSpaced_seq.cpp
new file mode 100644
index 0000000..f55c508
--- /dev/null
+++ b/doc/snippets/DenseBase_LinSpaced_seq.cpp
@@ -0,0 +1,2 @@
+cout << VectorXi::LinSpaced(Sequential,4,7,10).transpose() << endl;
+cout << VectorXd::LinSpaced(Sequential,5,0.0,1.0).transpose() << endl;
diff --git a/doc/snippets/DenseBase_setLinSpaced.cpp b/doc/snippets/DenseBase_setLinSpaced.cpp
new file mode 100644
index 0000000..46054f2
--- /dev/null
+++ b/doc/snippets/DenseBase_setLinSpaced.cpp
@@ -0,0 +1,3 @@
+VectorXf v;
+v.setLinSpaced(5,0.5f,1.5f);
+cout << v << endl;
diff --git a/doc/snippets/DirectionWise_replicate.cpp b/doc/snippets/DirectionWise_replicate.cpp
new file mode 100644
index 0000000..d92d4a3
--- /dev/null
+++ b/doc/snippets/DirectionWise_replicate.cpp
@@ -0,0 +1,4 @@
+MatrixXi m = MatrixXi::Random(2,3);
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "m.colwise().replicate<3>() = ..." << endl;
+cout << m.colwise().replicate<3>() << endl;
diff --git a/doc/snippets/DirectionWise_replicate_int.cpp b/doc/snippets/DirectionWise_replicate_int.cpp
new file mode 100644
index 0000000..f9b1b53
--- /dev/null
+++ b/doc/snippets/DirectionWise_replicate_int.cpp
@@ -0,0 +1,4 @@
+Vector3i v = Vector3i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "v.rowwise().replicate(5) = ..." << endl;
+cout << v.rowwise().replicate(5) << endl;
diff --git a/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp b/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp
new file mode 100644
index 0000000..c1d9fa8
--- /dev/null
+++ b/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp
@@ -0,0 +1,16 @@
+MatrixXd A = MatrixXd::Random(6,6);
+cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
+
+EigenSolver<MatrixXd> es(A);
+cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl;
+cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;
+
+complex<double> lambda = es.eigenvalues()[0];
+cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
+VectorXcd v = es.eigenvectors().col(0);
+cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
+cout << "... and A * v = " << endl << A.cast<complex<double> >() * v << endl << endl;
+
+MatrixXcd D = es.eigenvalues().asDiagonal();
+MatrixXcd V = es.eigenvectors();
+cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
diff --git a/doc/snippets/EigenSolver_compute.cpp b/doc/snippets/EigenSolver_compute.cpp
new file mode 100644
index 0000000..a5c96e9
--- /dev/null
+++ b/doc/snippets/EigenSolver_compute.cpp
@@ -0,0 +1,6 @@
+EigenSolver<MatrixXf> es;
+MatrixXf A = MatrixXf::Random(4,4);
+es.compute(A, /* computeEigenvectors = */ false);
+cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl;
+es.compute(A + MatrixXf::Identity(4,4), false); // re-use es to compute eigenvalues of A+I
+cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl;
diff --git a/doc/snippets/EigenSolver_eigenvalues.cpp b/doc/snippets/EigenSolver_eigenvalues.cpp
new file mode 100644
index 0000000..ed28869
--- /dev/null
+++ b/doc/snippets/EigenSolver_eigenvalues.cpp
@@ -0,0 +1,4 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+EigenSolver<MatrixXd> es(ones, false);
+cout << "The eigenvalues of the 3x3 matrix of ones are:" 
+     << endl << es.eigenvalues() << endl;
diff --git a/doc/snippets/EigenSolver_eigenvectors.cpp b/doc/snippets/EigenSolver_eigenvectors.cpp
new file mode 100644
index 0000000..0fad4da
--- /dev/null
+++ b/doc/snippets/EigenSolver_eigenvectors.cpp
@@ -0,0 +1,4 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+EigenSolver<MatrixXd> es(ones);
+cout << "The first eigenvector of the 3x3 matrix of ones is:" 
+     << endl << es.eigenvectors().col(1) << endl;
diff --git a/doc/snippets/EigenSolver_pseudoEigenvectors.cpp b/doc/snippets/EigenSolver_pseudoEigenvectors.cpp
new file mode 100644
index 0000000..85e2569
--- /dev/null
+++ b/doc/snippets/EigenSolver_pseudoEigenvectors.cpp
@@ -0,0 +1,9 @@
+MatrixXd A = MatrixXd::Random(6,6);
+cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
+
+EigenSolver<MatrixXd> es(A);
+MatrixXd D = es.pseudoEigenvalueMatrix();
+MatrixXd V = es.pseudoEigenvectors();
+cout << "The pseudo-eigenvalue matrix D is:" << endl << D << endl;
+cout << "The pseudo-eigenvector matrix V is:" << endl << V << endl;
+cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
diff --git a/doc/snippets/FullPivHouseholderQR_solve.cpp b/doc/snippets/FullPivHouseholderQR_solve.cpp
new file mode 100644
index 0000000..23bc074
--- /dev/null
+++ b/doc/snippets/FullPivHouseholderQR_solve.cpp
@@ -0,0 +1,8 @@
+Matrix3f m = Matrix3f::Random();
+Matrix3f y = Matrix3f::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the matrix y:" << endl << y << endl;
+Matrix3f x;
+x = m.fullPivHouseholderQr().solve(y);
+assert(y.isApprox(m*x));
+cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
diff --git a/doc/snippets/FullPivLU_image.cpp b/doc/snippets/FullPivLU_image.cpp
new file mode 100644
index 0000000..817bc1e
--- /dev/null
+++ b/doc/snippets/FullPivLU_image.cpp
@@ -0,0 +1,9 @@
+Matrix3d m;
+m << 1,1,0,
+     1,3,2,
+     0,1,1;
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Notice that the middle column is the sum of the two others, so the "
+     << "columns are linearly dependent." << endl;
+cout << "Here is a matrix whose columns have the same span but are linearly independent:"
+     << endl << m.fullPivLu().image(m) << endl;
diff --git a/doc/snippets/FullPivLU_kernel.cpp b/doc/snippets/FullPivLU_kernel.cpp
new file mode 100644
index 0000000..7086e01
--- /dev/null
+++ b/doc/snippets/FullPivLU_kernel.cpp
@@ -0,0 +1,7 @@
+MatrixXf m = MatrixXf::Random(3,5);
+cout << "Here is the matrix m:" << endl << m << endl;
+MatrixXf ker = m.fullPivLu().kernel();
+cout << "Here is a matrix whose columns form a basis of the kernel of m:"
+     << endl << ker << endl;
+cout << "By definition of the kernel, m*ker is zero:"
+     << endl << m*ker << endl;
diff --git a/doc/snippets/FullPivLU_solve.cpp b/doc/snippets/FullPivLU_solve.cpp
new file mode 100644
index 0000000..c1f8823
--- /dev/null
+++ b/doc/snippets/FullPivLU_solve.cpp
@@ -0,0 +1,11 @@
+Matrix<float,2,3> m = Matrix<float,2,3>::Random();
+Matrix2f y = Matrix2f::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the matrix y:" << endl << y << endl;
+Matrix<float,3,2> x = m.fullPivLu().solve(y);
+if((m*x).isApprox(y))
+{
+  cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
+}
+else
+  cout << "The equation mx=y does not have any solution." << endl;
diff --git a/doc/snippets/GeneralizedEigenSolver.cpp b/doc/snippets/GeneralizedEigenSolver.cpp
new file mode 100644
index 0000000..2acda45
--- /dev/null
+++ b/doc/snippets/GeneralizedEigenSolver.cpp
@@ -0,0 +1,7 @@
+GeneralizedEigenSolver<MatrixXf> ges;
+MatrixXf A = MatrixXf::Random(4,4);
+MatrixXf B = MatrixXf::Random(4,4);
+ges.compute(A, B);
+cout << "The (complex) numerators of the generalzied eigenvalues are: " << ges.alphas().transpose() << endl;
+cout << "The (real) denominatore of the generalzied eigenvalues are: " << ges.betas().transpose() << endl;
+cout << "The (complex) generalzied eigenvalues are (alphas./beta): " << ges.eigenvalues().transpose() << endl;
diff --git a/doc/snippets/HessenbergDecomposition_compute.cpp b/doc/snippets/HessenbergDecomposition_compute.cpp
new file mode 100644
index 0000000..50e3783
--- /dev/null
+++ b/doc/snippets/HessenbergDecomposition_compute.cpp
@@ -0,0 +1,6 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+HessenbergDecomposition<MatrixXcf> hd(4);
+hd.compute(A);
+cout << "The matrix H in the decomposition of A is:" << endl << hd.matrixH() << endl;
+hd.compute(2*A); // re-use hd to compute and store decomposition of 2A
+cout << "The matrix H in the decomposition of 2A is:" << endl << hd.matrixH() << endl;
diff --git a/doc/snippets/HessenbergDecomposition_matrixH.cpp b/doc/snippets/HessenbergDecomposition_matrixH.cpp
new file mode 100644
index 0000000..af01366
--- /dev/null
+++ b/doc/snippets/HessenbergDecomposition_matrixH.cpp
@@ -0,0 +1,8 @@
+Matrix4f A = MatrixXf::Random(4,4);
+cout << "Here is a random 4x4 matrix:" << endl << A << endl;
+HessenbergDecomposition<MatrixXf> hessOfA(A);
+MatrixXf H = hessOfA.matrixH();
+cout << "The Hessenberg matrix H is:" << endl << H << endl;
+MatrixXf Q = hessOfA.matrixQ();
+cout << "The orthogonal matrix Q is:" << endl << Q << endl;
+cout << "Q H Q^T is:" << endl << Q * H * Q.transpose() << endl;
diff --git a/doc/snippets/HessenbergDecomposition_packedMatrix.cpp b/doc/snippets/HessenbergDecomposition_packedMatrix.cpp
new file mode 100644
index 0000000..4fa5957
--- /dev/null
+++ b/doc/snippets/HessenbergDecomposition_packedMatrix.cpp
@@ -0,0 +1,9 @@
+Matrix4d A = Matrix4d::Random(4,4);
+cout << "Here is a random 4x4 matrix:" << endl << A << endl;
+HessenbergDecomposition<Matrix4d> hessOfA(A);
+Matrix4d pm = hessOfA.packedMatrix();
+cout << "The packed matrix M is:" << endl << pm << endl;
+cout << "The upper Hessenberg part corresponds to the matrix H, which is:" 
+     << endl << hessOfA.matrixH() << endl;
+Vector3d hc = hessOfA.householderCoefficients();
+cout << "The vector of Householder coefficients is:" << endl << hc << endl;
diff --git a/doc/snippets/HouseholderQR_householderQ.cpp b/doc/snippets/HouseholderQR_householderQ.cpp
new file mode 100644
index 0000000..e859ce5
--- /dev/null
+++ b/doc/snippets/HouseholderQR_householderQ.cpp
@@ -0,0 +1,7 @@
+MatrixXf A(MatrixXf::Random(5,3)), thinQ(MatrixXf::Identity(5,3)), Q;
+A.setRandom();
+HouseholderQR<MatrixXf> qr(A);
+Q = qr.householderQ();
+thinQ = qr.householderQ() * thinQ;
+std::cout << "The complete unitary matrix Q is:\n" << Q << "\n\n";
+std::cout << "The thin matrix Q is:\n" << thinQ << "\n\n";
diff --git a/doc/snippets/HouseholderQR_solve.cpp b/doc/snippets/HouseholderQR_solve.cpp
new file mode 100644
index 0000000..8cce6ce
--- /dev/null
+++ b/doc/snippets/HouseholderQR_solve.cpp
@@ -0,0 +1,9 @@
+typedef Matrix<float,3,3> Matrix3x3;
+Matrix3x3 m = Matrix3x3::Random();
+Matrix3f y = Matrix3f::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the matrix y:" << endl << y << endl;
+Matrix3f x;
+x = m.householderQr().solve(y);
+assert(y.isApprox(m*x));
+cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
diff --git a/doc/snippets/HouseholderSequence_HouseholderSequence.cpp b/doc/snippets/HouseholderSequence_HouseholderSequence.cpp
new file mode 100644
index 0000000..2632b83
--- /dev/null
+++ b/doc/snippets/HouseholderSequence_HouseholderSequence.cpp
@@ -0,0 +1,31 @@
+Matrix3d v = Matrix3d::Random();
+cout << "The matrix v is:" << endl;
+cout << v << endl;
+
+Vector3d v0(1, v(1,0), v(2,0));
+cout << "The first Householder vector is: v_0 = " << v0.transpose() << endl;
+Vector3d v1(0, 1, v(2,1));
+cout << "The second Householder vector is: v_1 = " << v1.transpose()  << endl;
+Vector3d v2(0, 0, 1);
+cout << "The third Householder vector is: v_2 = " << v2.transpose() << endl;
+
+Vector3d h = Vector3d::Random();
+cout << "The Householder coefficients are: h = " << h.transpose() << endl;
+
+Matrix3d H0 = Matrix3d::Identity() - h(0) * v0 * v0.adjoint();
+cout << "The first Householder reflection is represented by H_0 = " << endl;
+cout << H0 << endl;
+Matrix3d H1 = Matrix3d::Identity() - h(1) * v1 * v1.adjoint();
+cout << "The second Householder reflection is represented by H_1 = " << endl;
+cout << H1 << endl;
+Matrix3d H2 = Matrix3d::Identity() - h(2) * v2 * v2.adjoint();
+cout << "The third Householder reflection is represented by H_2 = " << endl;
+cout << H2 << endl;
+cout << "Their product is H_0 H_1 H_2 = " << endl;
+cout << H0 * H1 * H2 << endl;
+
+HouseholderSequence<Matrix3d, Vector3d> hhSeq(v, h);
+Matrix3d hhSeqAsMatrix(hhSeq);
+cout << "If we construct a HouseholderSequence from v and h" << endl;
+cout << "and convert it to a matrix, we get:" << endl;
+cout << hhSeqAsMatrix << endl;
diff --git a/doc/snippets/IOFormat.cpp b/doc/snippets/IOFormat.cpp
new file mode 100644
index 0000000..735f5dd
--- /dev/null
+++ b/doc/snippets/IOFormat.cpp
@@ -0,0 +1,14 @@
+std::string sep = "\n----------------------------------------\n";
+Matrix3d m1;
+m1 << 1.111111, 2, 3.33333, 4, 5, 6, 7, 8.888888, 9;
+
+IOFormat CommaInitFmt(StreamPrecision, DontAlignCols, ", ", ", ", "", "", " << ", ";");
+IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
+IOFormat OctaveFmt(StreamPrecision, 0, ", ", ";\n", "", "", "[", "]");
+IOFormat HeavyFmt(FullPrecision, 0, ", ", ";\n", "[", "]", "[", "]");
+
+std::cout << m1 << sep;
+std::cout << m1.format(CommaInitFmt) << sep;
+std::cout << m1.format(CleanFmt) << sep;
+std::cout << m1.format(OctaveFmt) << sep;
+std::cout << m1.format(HeavyFmt) << sep;
diff --git a/doc/snippets/JacobiSVD_basic.cpp b/doc/snippets/JacobiSVD_basic.cpp
new file mode 100644
index 0000000..ab24b9b
--- /dev/null
+++ b/doc/snippets/JacobiSVD_basic.cpp
@@ -0,0 +1,9 @@
+MatrixXf m = MatrixXf::Random(3,2);
+cout << "Here is the matrix m:" << endl << m << endl;
+JacobiSVD<MatrixXf> svd(m, ComputeThinU | ComputeThinV);
+cout << "Its singular values are:" << endl << svd.singularValues() << endl;
+cout << "Its left singular vectors are the columns of the thin U matrix:" << endl << svd.matrixU() << endl;
+cout << "Its right singular vectors are the columns of the thin V matrix:" << endl << svd.matrixV() << endl;
+Vector3f rhs(1, 0, 0);
+cout << "Now consider this rhs vector:" << endl << rhs << endl;
+cout << "A least-squares solution of m*x = rhs is:" << endl << svd.solve(rhs) << endl;
diff --git a/doc/snippets/Jacobi_makeGivens.cpp b/doc/snippets/Jacobi_makeGivens.cpp
new file mode 100644
index 0000000..4b733c3
--- /dev/null
+++ b/doc/snippets/Jacobi_makeGivens.cpp
@@ -0,0 +1,6 @@
+Vector2f v = Vector2f::Random();
+JacobiRotation<float> G;
+G.makeGivens(v.x(), v.y());
+cout << "Here is the vector v:" << endl << v << endl;
+v.applyOnTheLeft(0, 1, G.adjoint());
+cout << "Here is the vector J' * v:" << endl << v << endl;
\ No newline at end of file
diff --git a/doc/snippets/Jacobi_makeJacobi.cpp b/doc/snippets/Jacobi_makeJacobi.cpp
new file mode 100644
index 0000000..0cc331d
--- /dev/null
+++ b/doc/snippets/Jacobi_makeJacobi.cpp
@@ -0,0 +1,8 @@
+Matrix2f m = Matrix2f::Random();
+m = (m + m.adjoint()).eval();
+JacobiRotation<float> J;
+J.makeJacobi(m, 0, 1);
+cout << "Here is the matrix m:" << endl << m << endl;
+m.applyOnTheLeft(0, 1, J.adjoint());
+m.applyOnTheRight(0, 1, J);
+cout << "Here is the matrix J' * m * J:" << endl << m << endl;
\ No newline at end of file
diff --git a/doc/snippets/LLT_example.cpp b/doc/snippets/LLT_example.cpp
new file mode 100644
index 0000000..46fb407
--- /dev/null
+++ b/doc/snippets/LLT_example.cpp
@@ -0,0 +1,12 @@
+MatrixXd A(3,3);
+A << 4,-1,2, -1,6,0, 2,0,5;
+cout << "The matrix A is" << endl << A << endl;
+
+LLT<MatrixXd> lltOfA(A); // compute the Cholesky decomposition of A
+MatrixXd L = lltOfA.matrixL(); // retrieve factor L  in the decomposition
+// The previous two lines can also be written as "L = A.llt().matrixL()"
+
+cout << "The Cholesky factor L is" << endl << L << endl;
+cout << "To check this, let us compute L * L.transpose()" << endl;
+cout << L * L.transpose() << endl;
+cout << "This should equal the matrix A" << endl;
diff --git a/doc/snippets/LLT_solve.cpp b/doc/snippets/LLT_solve.cpp
new file mode 100644
index 0000000..7095d2c
--- /dev/null
+++ b/doc/snippets/LLT_solve.cpp
@@ -0,0 +1,8 @@
+typedef Matrix<float,Dynamic,2> DataMatrix;
+// let's generate some samples on the 3D plane of equation z = 2x+3y (with some noise)
+DataMatrix samples = DataMatrix::Random(12,2);
+VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1;
+// and let's solve samples * [x y]^T = elevations in least square sense:
+Matrix<float,2,1> xy
+ = (samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations));
+cout << xy << endl;
diff --git a/doc/snippets/Map_general_stride.cpp b/doc/snippets/Map_general_stride.cpp
new file mode 100644
index 0000000..0657e7f
--- /dev/null
+++ b/doc/snippets/Map_general_stride.cpp
@@ -0,0 +1,5 @@
+int array[24];
+for(int i = 0; i < 24; ++i) array[i] = i;
+cout << Map<MatrixXi, 0, Stride<Dynamic,2> >
+         (array, 3, 3, Stride<Dynamic,2>(8, 2))
+     << endl;
diff --git a/doc/snippets/Map_inner_stride.cpp b/doc/snippets/Map_inner_stride.cpp
new file mode 100644
index 0000000..d95ae9b
--- /dev/null
+++ b/doc/snippets/Map_inner_stride.cpp
@@ -0,0 +1,5 @@
+int array[12];
+for(int i = 0; i < 12; ++i) array[i] = i;
+cout << Map<VectorXi, 0, InnerStride<2> >
+         (array, 6) // the inner stride has already been passed as template parameter
+     << endl;
diff --git a/doc/snippets/Map_outer_stride.cpp b/doc/snippets/Map_outer_stride.cpp
new file mode 100644
index 0000000..2f6f052
--- /dev/null
+++ b/doc/snippets/Map_outer_stride.cpp
@@ -0,0 +1,3 @@
+int array[12];
+for(int i = 0; i < 12; ++i) array[i] = i;
+cout << Map<MatrixXi, 0, OuterStride<> >(array, 3, 3, OuterStride<>(4)) << endl;
diff --git a/doc/snippets/Map_placement_new.cpp b/doc/snippets/Map_placement_new.cpp
new file mode 100644
index 0000000..2e40eca
--- /dev/null
+++ b/doc/snippets/Map_placement_new.cpp
@@ -0,0 +1,5 @@
+int data[] = {1,2,3,4,5,6,7,8,9};
+Map<RowVectorXi> v(data,4);
+cout << "The mapped vector v is: " << v << "\n";
+new (&v) Map<RowVectorXi>(data+4,5);
+cout << "Now v is: " << v << "\n";
\ No newline at end of file
diff --git a/doc/snippets/Map_simple.cpp b/doc/snippets/Map_simple.cpp
new file mode 100644
index 0000000..423bb52
--- /dev/null
+++ b/doc/snippets/Map_simple.cpp
@@ -0,0 +1,3 @@
+int array[9];
+for(int i = 0; i < 9; ++i) array[i] = i;
+cout << Map<Matrix3i>(array) << endl;
diff --git a/doc/snippets/MatrixBase_adjoint.cpp b/doc/snippets/MatrixBase_adjoint.cpp
new file mode 100644
index 0000000..4680d59
--- /dev/null
+++ b/doc/snippets/MatrixBase_adjoint.cpp
@@ -0,0 +1,3 @@
+Matrix2cf m = Matrix2cf::Random();
+cout << "Here is the 2x2 complex matrix m:" << endl << m << endl;
+cout << "Here is the adjoint of m:" << endl << m.adjoint() << endl;
diff --git a/doc/snippets/MatrixBase_all.cpp b/doc/snippets/MatrixBase_all.cpp
new file mode 100644
index 0000000..46f26f1
--- /dev/null
+++ b/doc/snippets/MatrixBase_all.cpp
@@ -0,0 +1,7 @@
+Vector3f boxMin(Vector3f::Zero()), boxMax(Vector3f::Ones());
+Vector3f p0 = Vector3f::Random(), p1 = Vector3f::Random().cwiseAbs();
+// let's check if p0 and p1 are inside the axis aligned box defined by the corners boxMin,boxMax:
+cout << "Is (" << p0.transpose() << ") inside the box: "
+     << ((boxMin.array()<p0.array()).all() && (boxMax.array()>p0.array()).all()) << endl;
+cout << "Is (" << p1.transpose() << ") inside the box: "
+     << ((boxMin.array()<p1.array()).all() && (boxMax.array()>p1.array()).all()) << endl;
diff --git a/doc/snippets/MatrixBase_applyOnTheLeft.cpp b/doc/snippets/MatrixBase_applyOnTheLeft.cpp
new file mode 100644
index 0000000..6398c87
--- /dev/null
+++ b/doc/snippets/MatrixBase_applyOnTheLeft.cpp
@@ -0,0 +1,7 @@
+Matrix3f A = Matrix3f::Random(3,3), B;
+B << 0,1,0,  
+     0,0,1,  
+     1,0,0;
+cout << "At start, A = " << endl << A << endl;
+A.applyOnTheLeft(B); 
+cout << "After applyOnTheLeft, A = " << endl << A << endl;
diff --git a/doc/snippets/MatrixBase_applyOnTheRight.cpp b/doc/snippets/MatrixBase_applyOnTheRight.cpp
new file mode 100644
index 0000000..e4b71b2
--- /dev/null
+++ b/doc/snippets/MatrixBase_applyOnTheRight.cpp
@@ -0,0 +1,9 @@
+Matrix3f A = Matrix3f::Random(3,3), B;
+B << 0,1,0,  
+     0,0,1,  
+     1,0,0;
+cout << "At start, A = " << endl << A << endl;
+A *= B;
+cout << "After A *= B, A = " << endl << A << endl;
+A.applyOnTheRight(B);  // equivalent to A *= B
+cout << "After applyOnTheRight, A = " << endl << A << endl;
diff --git a/doc/snippets/MatrixBase_array.cpp b/doc/snippets/MatrixBase_array.cpp
new file mode 100644
index 0000000..f215086
--- /dev/null
+++ b/doc/snippets/MatrixBase_array.cpp
@@ -0,0 +1,4 @@
+Vector3d v(1,2,3);
+v.array() += 3;
+v.array() -= 2;
+cout << v << endl;
diff --git a/doc/snippets/MatrixBase_array_const.cpp b/doc/snippets/MatrixBase_array_const.cpp
new file mode 100644
index 0000000..cd3b26a
--- /dev/null
+++ b/doc/snippets/MatrixBase_array_const.cpp
@@ -0,0 +1,4 @@
+Vector3d v(-1,2,-3);
+cout << "the absolute values:" << endl << v.array().abs() << endl;
+cout << "the absolute values plus one:" << endl << v.array().abs()+1 << endl;
+cout << "sum of the squares: " << v.array().square().sum() << endl;
diff --git a/doc/snippets/MatrixBase_asDiagonal.cpp b/doc/snippets/MatrixBase_asDiagonal.cpp
new file mode 100644
index 0000000..b01082d
--- /dev/null
+++ b/doc/snippets/MatrixBase_asDiagonal.cpp
@@ -0,0 +1 @@
+cout << Matrix3i(Vector3i(2,5,6).asDiagonal()) << endl;
diff --git a/doc/snippets/MatrixBase_block_int_int.cpp b/doc/snippets/MatrixBase_block_int_int.cpp
new file mode 100644
index 0000000..f99b6d4
--- /dev/null
+++ b/doc/snippets/MatrixBase_block_int_int.cpp
@@ -0,0 +1,5 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.block<2,2>(1,1):" << endl << m.block<2,2>(1,1) << endl;
+m.block<2,2>(1,1).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_block_int_int_int_int.cpp b/doc/snippets/MatrixBase_block_int_int_int_int.cpp
new file mode 100644
index 0000000..7238cbb
--- /dev/null
+++ b/doc/snippets/MatrixBase_block_int_int_int_int.cpp
@@ -0,0 +1,5 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.block(1, 1, 2, 2):" << endl << m.block(1, 1, 2, 2) << endl;
+m.block(1, 1, 2, 2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_bottomLeftCorner_int_int.cpp b/doc/snippets/MatrixBase_bottomLeftCorner_int_int.cpp
new file mode 100644
index 0000000..ebae95e
--- /dev/null
+++ b/doc/snippets/MatrixBase_bottomLeftCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.bottomLeftCorner(2, 2):" << endl;
+cout << m.bottomLeftCorner(2, 2) << endl;
+m.bottomLeftCorner(2, 2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_bottomRightCorner_int_int.cpp b/doc/snippets/MatrixBase_bottomRightCorner_int_int.cpp
new file mode 100644
index 0000000..bf05093
--- /dev/null
+++ b/doc/snippets/MatrixBase_bottomRightCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.bottomRightCorner(2, 2):" << endl;
+cout << m.bottomRightCorner(2, 2) << endl;
+m.bottomRightCorner(2, 2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_bottomRows_int.cpp b/doc/snippets/MatrixBase_bottomRows_int.cpp
new file mode 100644
index 0000000..47ca92e
--- /dev/null
+++ b/doc/snippets/MatrixBase_bottomRows_int.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.bottomRows(2):" << endl;
+cout << a.bottomRows(2) << endl;
+a.bottomRows(2).setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_cast.cpp b/doc/snippets/MatrixBase_cast.cpp
new file mode 100644
index 0000000..016880b
--- /dev/null
+++ b/doc/snippets/MatrixBase_cast.cpp
@@ -0,0 +1,3 @@
+Matrix2d md = Matrix2d::Identity() * 0.45;
+Matrix2f mf = Matrix2f::Identity();
+cout << md + mf.cast<double>() << endl;
diff --git a/doc/snippets/MatrixBase_col.cpp b/doc/snippets/MatrixBase_col.cpp
new file mode 100644
index 0000000..87c91b1
--- /dev/null
+++ b/doc/snippets/MatrixBase_col.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Identity();
+m.col(1) = Vector3d(4,5,6);
+cout << m << endl;
diff --git a/doc/snippets/MatrixBase_colwise.cpp b/doc/snippets/MatrixBase_colwise.cpp
new file mode 100644
index 0000000..a048bef
--- /dev/null
+++ b/doc/snippets/MatrixBase_colwise.cpp
@@ -0,0 +1,5 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the sum of each column:" << endl << m.colwise().sum() << endl;
+cout << "Here is the maximum absolute value of each column:"
+     << endl << m.cwiseAbs().colwise().maxCoeff() << endl;
diff --git a/doc/snippets/MatrixBase_computeInverseAndDetWithCheck.cpp b/doc/snippets/MatrixBase_computeInverseAndDetWithCheck.cpp
new file mode 100644
index 0000000..a7b084f
--- /dev/null
+++ b/doc/snippets/MatrixBase_computeInverseAndDetWithCheck.cpp
@@ -0,0 +1,13 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+Matrix3d inverse;
+bool invertible;
+double determinant;
+m.computeInverseAndDetWithCheck(inverse,determinant,invertible);
+cout << "Its determinant is " << determinant << endl;
+if(invertible) {
+  cout << "It is invertible, and its inverse is:" << endl << inverse << endl;
+}
+else {
+  cout << "It is not invertible." << endl;
+}
diff --git a/doc/snippets/MatrixBase_computeInverseWithCheck.cpp b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp
new file mode 100644
index 0000000..873a9f8
--- /dev/null
+++ b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp
@@ -0,0 +1,11 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+Matrix3d inverse;
+bool invertible;
+m.computeInverseWithCheck(inverse,invertible);
+if(invertible) {
+  cout << "It is invertible, and its inverse is:" << endl << inverse << endl;
+}
+else {
+  cout << "It is not invertible." << endl;
+}
diff --git a/doc/snippets/MatrixBase_cwiseAbs.cpp b/doc/snippets/MatrixBase_cwiseAbs.cpp
new file mode 100644
index 0000000..28a3160
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseAbs.cpp
@@ -0,0 +1,4 @@
+MatrixXd m(2,3);
+m << 2, -4, 6,   
+     -5, 1, 0;
+cout << m.cwiseAbs() << endl;
diff --git a/doc/snippets/MatrixBase_cwiseAbs2.cpp b/doc/snippets/MatrixBase_cwiseAbs2.cpp
new file mode 100644
index 0000000..889a2e2
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseAbs2.cpp
@@ -0,0 +1,4 @@
+MatrixXd m(2,3);
+m << 2, -4, 6,   
+     -5, 1, 0;
+cout << m.cwiseAbs2() << endl;
diff --git a/doc/snippets/MatrixBase_cwiseEqual.cpp b/doc/snippets/MatrixBase_cwiseEqual.cpp
new file mode 100644
index 0000000..eb3656f
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseEqual.cpp
@@ -0,0 +1,7 @@
+MatrixXi m(2,2);
+m << 1, 0,
+     1, 1;
+cout << "Comparing m with identity matrix:" << endl;
+cout << m.cwiseEqual(MatrixXi::Identity(2,2)) << endl;
+int count = m.cwiseEqual(MatrixXi::Identity(2,2)).count();
+cout << "Number of coefficients that are equal: " << count << endl;
diff --git a/doc/snippets/MatrixBase_cwiseInverse.cpp b/doc/snippets/MatrixBase_cwiseInverse.cpp
new file mode 100644
index 0000000..23e08f7
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseInverse.cpp
@@ -0,0 +1,4 @@
+MatrixXd m(2,3);
+m << 2, 0.5, 1,   
+     3, 0.25, 1;
+cout << m.cwiseInverse() << endl;
diff --git a/doc/snippets/MatrixBase_cwiseMax.cpp b/doc/snippets/MatrixBase_cwiseMax.cpp
new file mode 100644
index 0000000..3c95681
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseMax.cpp
@@ -0,0 +1,2 @@
+Vector3d v(2,3,4), w(4,2,3);
+cout << v.cwiseMax(w) << endl;
diff --git a/doc/snippets/MatrixBase_cwiseMin.cpp b/doc/snippets/MatrixBase_cwiseMin.cpp
new file mode 100644
index 0000000..82fc761
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseMin.cpp
@@ -0,0 +1,2 @@
+Vector3d v(2,3,4), w(4,2,3);
+cout << v.cwiseMin(w) << endl;
diff --git a/doc/snippets/MatrixBase_cwiseNotEqual.cpp b/doc/snippets/MatrixBase_cwiseNotEqual.cpp
new file mode 100644
index 0000000..6a2e4fb
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseNotEqual.cpp
@@ -0,0 +1,7 @@
+MatrixXi m(2,2);
+m << 1, 0,
+     1, 1;
+cout << "Comparing m with identity matrix:" << endl;
+cout << m.cwiseNotEqual(MatrixXi::Identity(2,2)) << endl;
+int count = m.cwiseNotEqual(MatrixXi::Identity(2,2)).count();
+cout << "Number of coefficients that are not equal: " << count << endl;
diff --git a/doc/snippets/MatrixBase_cwiseProduct.cpp b/doc/snippets/MatrixBase_cwiseProduct.cpp
new file mode 100644
index 0000000..1db3a11
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseProduct.cpp
@@ -0,0 +1,4 @@
+Matrix3i a = Matrix3i::Random(), b = Matrix3i::Random();
+Matrix3i c = a.cwiseProduct(b);
+cout << "a:\n" << a << "\nb:\n" << b << "\nc:\n" << c << endl;
+
diff --git a/doc/snippets/MatrixBase_cwiseQuotient.cpp b/doc/snippets/MatrixBase_cwiseQuotient.cpp
new file mode 100644
index 0000000..9691212
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseQuotient.cpp
@@ -0,0 +1,2 @@
+Vector3d v(2,3,4), w(4,2,3);
+cout << v.cwiseQuotient(w) << endl;
diff --git a/doc/snippets/MatrixBase_cwiseSqrt.cpp b/doc/snippets/MatrixBase_cwiseSqrt.cpp
new file mode 100644
index 0000000..4bfd75d
--- /dev/null
+++ b/doc/snippets/MatrixBase_cwiseSqrt.cpp
@@ -0,0 +1,2 @@
+Vector3d v(1,2,4);
+cout << v.cwiseSqrt() << endl;
diff --git a/doc/snippets/MatrixBase_diagonal.cpp b/doc/snippets/MatrixBase_diagonal.cpp
new file mode 100644
index 0000000..cd63413
--- /dev/null
+++ b/doc/snippets/MatrixBase_diagonal.cpp
@@ -0,0 +1,4 @@
+Matrix3i m = Matrix3i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here are the coefficients on the main diagonal of m:" << endl
+     << m.diagonal() << endl;
diff --git a/doc/snippets/MatrixBase_diagonal_int.cpp b/doc/snippets/MatrixBase_diagonal_int.cpp
new file mode 100644
index 0000000..7b66abf
--- /dev/null
+++ b/doc/snippets/MatrixBase_diagonal_int.cpp
@@ -0,0 +1,5 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here are the coefficients on the 1st super-diagonal and 2nd sub-diagonal of m:" << endl
+     << m.diagonal(1).transpose() << endl
+     << m.diagonal(-2).transpose() << endl;
diff --git a/doc/snippets/MatrixBase_diagonal_template_int.cpp b/doc/snippets/MatrixBase_diagonal_template_int.cpp
new file mode 100644
index 0000000..0e73d1c
--- /dev/null
+++ b/doc/snippets/MatrixBase_diagonal_template_int.cpp
@@ -0,0 +1,5 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here are the coefficients on the 1st super-diagonal and 2nd sub-diagonal of m:" << endl
+     << m.diagonal<1>().transpose() << endl
+     << m.diagonal<-2>().transpose() << endl;
diff --git a/doc/snippets/MatrixBase_eigenvalues.cpp b/doc/snippets/MatrixBase_eigenvalues.cpp
new file mode 100644
index 0000000..039f887
--- /dev/null
+++ b/doc/snippets/MatrixBase_eigenvalues.cpp
@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+VectorXcd eivals = ones.eigenvalues();
+cout << "The eigenvalues of the 3x3 matrix of ones are:" << endl << eivals << endl;
diff --git a/doc/snippets/MatrixBase_end_int.cpp b/doc/snippets/MatrixBase_end_int.cpp
new file mode 100644
index 0000000..03c54a9
--- /dev/null
+++ b/doc/snippets/MatrixBase_end_int.cpp
@@ -0,0 +1,5 @@
+RowVector4i v = RowVector4i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "Here is v.tail(2):" << endl << v.tail(2) << endl;
+v.tail(2).setZero();
+cout << "Now the vector v is:" << endl << v << endl;
diff --git a/doc/snippets/MatrixBase_eval.cpp b/doc/snippets/MatrixBase_eval.cpp
new file mode 100644
index 0000000..1df3aa0
--- /dev/null
+++ b/doc/snippets/MatrixBase_eval.cpp
@@ -0,0 +1,12 @@
+Matrix2f M = Matrix2f::Random();
+Matrix2f m;
+m = M;
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Now we want to copy a column into a row." << endl;
+cout << "If we do m.col(1) = m.row(0), then m becomes:" << endl;
+m.col(1) = m.row(0);
+cout << m << endl << "which is wrong!" << endl;
+cout << "Now let us instead do m.col(1) = m.row(0).eval(). Then m becomes" << endl;
+m = M;
+m.col(1) = m.row(0).eval();
+cout << m << endl << "which is right." << endl;
diff --git a/doc/snippets/MatrixBase_extract.cpp b/doc/snippets/MatrixBase_extract.cpp
new file mode 100644
index 0000000..c96220f
--- /dev/null
+++ b/doc/snippets/MatrixBase_extract.cpp
@@ -0,0 +1,13 @@
+#ifndef _MSC_VER
+  #warning deprecated
+#endif
+/* deprecated
+Matrix3i m = Matrix3i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the upper-triangular matrix extracted from m:" << endl
+     << m.part<Eigen::UpperTriangular>() << endl;
+cout << "Here is the strictly-upper-triangular matrix extracted from m:" << endl
+     << m.part<Eigen::StrictlyUpperTriangular>() << endl;
+cout << "Here is the unit-lower-triangular matrix extracted from m:" << endl
+     << m.part<Eigen::UnitLowerTriangular>() << endl;
+*/
\ No newline at end of file
diff --git a/doc/snippets/MatrixBase_fixedBlock_int_int.cpp b/doc/snippets/MatrixBase_fixedBlock_int_int.cpp
new file mode 100644
index 0000000..3201127
--- /dev/null
+++ b/doc/snippets/MatrixBase_fixedBlock_int_int.cpp
@@ -0,0 +1,5 @@
+Matrix4d m = Vector4d(1,2,3,4).asDiagonal();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.fixed<2, 2>(2, 2):" << endl << m.block<2, 2>(2, 2) << endl;
+m.block<2, 2>(2, 0) = m.block<2, 2>(2, 2);
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_identity.cpp b/doc/snippets/MatrixBase_identity.cpp
new file mode 100644
index 0000000..b5c1e59
--- /dev/null
+++ b/doc/snippets/MatrixBase_identity.cpp
@@ -0,0 +1 @@
+cout << Matrix<double, 3, 4>::Identity() << endl;
diff --git a/doc/snippets/MatrixBase_identity_int_int.cpp b/doc/snippets/MatrixBase_identity_int_int.cpp
new file mode 100644
index 0000000..918649d
--- /dev/null
+++ b/doc/snippets/MatrixBase_identity_int_int.cpp
@@ -0,0 +1 @@
+cout << MatrixXd::Identity(4, 3) << endl;
diff --git a/doc/snippets/MatrixBase_inverse.cpp b/doc/snippets/MatrixBase_inverse.cpp
new file mode 100644
index 0000000..a56142e
--- /dev/null
+++ b/doc/snippets/MatrixBase_inverse.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Its inverse is:" << endl << m.inverse() << endl;
diff --git a/doc/snippets/MatrixBase_isDiagonal.cpp b/doc/snippets/MatrixBase_isDiagonal.cpp
new file mode 100644
index 0000000..5b1d599
--- /dev/null
+++ b/doc/snippets/MatrixBase_isDiagonal.cpp
@@ -0,0 +1,6 @@
+Matrix3d m = 10000 * Matrix3d::Identity();
+m(0,2) = 1;
+cout << "Here's the matrix m:" << endl << m << endl;
+cout << "m.isDiagonal() returns: " << m.isDiagonal() << endl;
+cout << "m.isDiagonal(1e-3) returns: " << m.isDiagonal(1e-3) << endl;
+
diff --git a/doc/snippets/MatrixBase_isIdentity.cpp b/doc/snippets/MatrixBase_isIdentity.cpp
new file mode 100644
index 0000000..17b756c
--- /dev/null
+++ b/doc/snippets/MatrixBase_isIdentity.cpp
@@ -0,0 +1,5 @@
+Matrix3d m = Matrix3d::Identity();
+m(0,2) = 1e-4;
+cout << "Here's the matrix m:" << endl << m << endl;
+cout << "m.isIdentity() returns: " << m.isIdentity() << endl;
+cout << "m.isIdentity(1e-3) returns: " << m.isIdentity(1e-3) << endl;
diff --git a/doc/snippets/MatrixBase_isOnes.cpp b/doc/snippets/MatrixBase_isOnes.cpp
new file mode 100644
index 0000000..f82f628
--- /dev/null
+++ b/doc/snippets/MatrixBase_isOnes.cpp
@@ -0,0 +1,5 @@
+Matrix3d m = Matrix3d::Ones();
+m(0,2) += 1e-4;
+cout << "Here's the matrix m:" << endl << m << endl;
+cout << "m.isOnes() returns: " << m.isOnes() << endl;
+cout << "m.isOnes(1e-3) returns: " << m.isOnes(1e-3) << endl;
diff --git a/doc/snippets/MatrixBase_isOrthogonal.cpp b/doc/snippets/MatrixBase_isOrthogonal.cpp
new file mode 100644
index 0000000..b22af06
--- /dev/null
+++ b/doc/snippets/MatrixBase_isOrthogonal.cpp
@@ -0,0 +1,6 @@
+Vector3d v(1,0,0);
+Vector3d w(1e-4,0,1);
+cout << "Here's the vector v:" << endl << v << endl;
+cout << "Here's the vector w:" << endl << w << endl;
+cout << "v.isOrthogonal(w) returns: " << v.isOrthogonal(w) << endl;
+cout << "v.isOrthogonal(w,1e-3) returns: " << v.isOrthogonal(w,1e-3) << endl;
diff --git a/doc/snippets/MatrixBase_isUnitary.cpp b/doc/snippets/MatrixBase_isUnitary.cpp
new file mode 100644
index 0000000..3877da3
--- /dev/null
+++ b/doc/snippets/MatrixBase_isUnitary.cpp
@@ -0,0 +1,5 @@
+Matrix3d m = Matrix3d::Identity();
+m(0,2) = 1e-4;
+cout << "Here's the matrix m:" << endl << m << endl;
+cout << "m.isUnitary() returns: " << m.isUnitary() << endl;
+cout << "m.isUnitary(1e-3) returns: " << m.isUnitary(1e-3) << endl;
diff --git a/doc/snippets/MatrixBase_isZero.cpp b/doc/snippets/MatrixBase_isZero.cpp
new file mode 100644
index 0000000..c2cfe22
--- /dev/null
+++ b/doc/snippets/MatrixBase_isZero.cpp
@@ -0,0 +1,5 @@
+Matrix3d m = Matrix3d::Zero();
+m(0,2) = 1e-4;
+cout << "Here's the matrix m:" << endl << m << endl;
+cout << "m.isZero() returns: " << m.isZero() << endl;
+cout << "m.isZero(1e-3) returns: " << m.isZero(1e-3) << endl;
diff --git a/doc/snippets/MatrixBase_leftCols_int.cpp b/doc/snippets/MatrixBase_leftCols_int.cpp
new file mode 100644
index 0000000..6ea984e
--- /dev/null
+++ b/doc/snippets/MatrixBase_leftCols_int.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.leftCols(2):" << endl;
+cout << a.leftCols(2) << endl;
+a.leftCols(2).setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_marked.cpp b/doc/snippets/MatrixBase_marked.cpp
new file mode 100644
index 0000000..f607121
--- /dev/null
+++ b/doc/snippets/MatrixBase_marked.cpp
@@ -0,0 +1,14 @@
+#ifndef _MSC_VER
+  #warning deprecated
+#endif
+/*
+Matrix3d m = Matrix3d::Zero();
+m.part<Eigen::UpperTriangular>().setOnes();
+cout << "Here is the matrix m:" << endl << m << endl;
+Matrix3d n = Matrix3d::Ones();
+n.part<Eigen::LowerTriangular>() *= 2;
+cout << "Here is the matrix n:" << endl << n << endl;
+cout << "And now here is m.inverse()*n, taking advantage of the fact that"
+        " m is upper-triangular:" << endl
+     << m.marked<Eigen::UpperTriangular>().solveTriangular(n);
+*/
\ No newline at end of file
diff --git a/doc/snippets/MatrixBase_noalias.cpp b/doc/snippets/MatrixBase_noalias.cpp
new file mode 100644
index 0000000..3b54a79
--- /dev/null
+++ b/doc/snippets/MatrixBase_noalias.cpp
@@ -0,0 +1,3 @@
+Matrix2d a, b, c; a << 1,2,3,4; b << 5,6,7,8;
+c.noalias() = a * b; // this computes the product directly to c
+cout << c << endl;
diff --git a/doc/snippets/MatrixBase_ones.cpp b/doc/snippets/MatrixBase_ones.cpp
new file mode 100644
index 0000000..02c767c
--- /dev/null
+++ b/doc/snippets/MatrixBase_ones.cpp
@@ -0,0 +1,2 @@
+cout << Matrix2d::Ones() << endl;
+cout << 6 * RowVector4i::Ones() << endl;
diff --git a/doc/snippets/MatrixBase_ones_int.cpp b/doc/snippets/MatrixBase_ones_int.cpp
new file mode 100644
index 0000000..2ef188e
--- /dev/null
+++ b/doc/snippets/MatrixBase_ones_int.cpp
@@ -0,0 +1,2 @@
+cout << 6 * RowVectorXi::Ones(4) << endl;
+cout << VectorXf::Ones(2) << endl;
diff --git a/doc/snippets/MatrixBase_ones_int_int.cpp b/doc/snippets/MatrixBase_ones_int_int.cpp
new file mode 100644
index 0000000..60f5a31
--- /dev/null
+++ b/doc/snippets/MatrixBase_ones_int_int.cpp
@@ -0,0 +1 @@
+cout << MatrixXi::Ones(2,3) << endl;
diff --git a/doc/snippets/MatrixBase_operatorNorm.cpp b/doc/snippets/MatrixBase_operatorNorm.cpp
new file mode 100644
index 0000000..355246f
--- /dev/null
+++ b/doc/snippets/MatrixBase_operatorNorm.cpp
@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+cout << "The operator norm of the 3x3 matrix of ones is "
+     << ones.operatorNorm() << endl;
diff --git a/doc/snippets/MatrixBase_part.cpp b/doc/snippets/MatrixBase_part.cpp
new file mode 100644
index 0000000..d3e7f48
--- /dev/null
+++ b/doc/snippets/MatrixBase_part.cpp
@@ -0,0 +1,13 @@
+#ifndef _MSC_VER
+  #warning deprecated
+#endif
+/*
+Matrix3d m = Matrix3d::Zero();
+m.part<Eigen::StrictlyUpperTriangular>().setOnes();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "And let us now compute m*m.adjoint() in a very optimized way" << endl
+     << "taking advantage of the symmetry." << endl;
+Matrix3d n;
+n.part<Eigen::SelfAdjoint>() = (m*m.adjoint()).lazy();
+cout << "The result is:" << endl << n << endl;
+*/
\ No newline at end of file
diff --git a/doc/snippets/MatrixBase_prod.cpp b/doc/snippets/MatrixBase_prod.cpp
new file mode 100644
index 0000000..d2f27bd
--- /dev/null
+++ b/doc/snippets/MatrixBase_prod.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the product of all the coefficients:" << endl << m.prod() << endl;
diff --git a/doc/snippets/MatrixBase_random.cpp b/doc/snippets/MatrixBase_random.cpp
new file mode 100644
index 0000000..65fc524
--- /dev/null
+++ b/doc/snippets/MatrixBase_random.cpp
@@ -0,0 +1 @@
+cout << 100 * Matrix2i::Random() << endl;
diff --git a/doc/snippets/MatrixBase_random_int.cpp b/doc/snippets/MatrixBase_random_int.cpp
new file mode 100644
index 0000000..f161d03
--- /dev/null
+++ b/doc/snippets/MatrixBase_random_int.cpp
@@ -0,0 +1 @@
+cout << VectorXi::Random(2) << endl;
diff --git a/doc/snippets/MatrixBase_random_int_int.cpp b/doc/snippets/MatrixBase_random_int_int.cpp
new file mode 100644
index 0000000..3f0f7dd
--- /dev/null
+++ b/doc/snippets/MatrixBase_random_int_int.cpp
@@ -0,0 +1 @@
+cout << MatrixXi::Random(2,3) << endl;
diff --git a/doc/snippets/MatrixBase_replicate.cpp b/doc/snippets/MatrixBase_replicate.cpp
new file mode 100644
index 0000000..3ce52bc
--- /dev/null
+++ b/doc/snippets/MatrixBase_replicate.cpp
@@ -0,0 +1,4 @@
+MatrixXi m = MatrixXi::Random(2,3);
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "m.replicate<3,2>() = ..." << endl;
+cout << m.replicate<3,2>() << endl;
diff --git a/doc/snippets/MatrixBase_replicate_int_int.cpp b/doc/snippets/MatrixBase_replicate_int_int.cpp
new file mode 100644
index 0000000..b1dbc70
--- /dev/null
+++ b/doc/snippets/MatrixBase_replicate_int_int.cpp
@@ -0,0 +1,4 @@
+Vector3i v = Vector3i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "v.replicate(2,5) = ..." << endl;
+cout << v.replicate(2,5) << endl;
diff --git a/doc/snippets/MatrixBase_reverse.cpp b/doc/snippets/MatrixBase_reverse.cpp
new file mode 100644
index 0000000..f545a28
--- /dev/null
+++ b/doc/snippets/MatrixBase_reverse.cpp
@@ -0,0 +1,8 @@
+MatrixXi m = MatrixXi::Random(3,4);
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the reverse of m:" << endl << m.reverse() << endl;
+cout << "Here is the coefficient (1,0) in the reverse of m:" << endl
+     << m.reverse()(1,0) << endl;
+cout << "Let us overwrite this coefficient with the value 4." << endl;
+m.reverse()(1,0) = 4;
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_rightCols_int.cpp b/doc/snippets/MatrixBase_rightCols_int.cpp
new file mode 100644
index 0000000..cb51340
--- /dev/null
+++ b/doc/snippets/MatrixBase_rightCols_int.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.rightCols(2):" << endl;
+cout << a.rightCols(2) << endl;
+a.rightCols(2).setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_row.cpp b/doc/snippets/MatrixBase_row.cpp
new file mode 100644
index 0000000..b15e626
--- /dev/null
+++ b/doc/snippets/MatrixBase_row.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Identity();
+m.row(1) = Vector3d(4,5,6);
+cout << m << endl;
diff --git a/doc/snippets/MatrixBase_rowwise.cpp b/doc/snippets/MatrixBase_rowwise.cpp
new file mode 100644
index 0000000..ae93964
--- /dev/null
+++ b/doc/snippets/MatrixBase_rowwise.cpp
@@ -0,0 +1,5 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the sum of each row:" << endl << m.rowwise().sum() << endl;
+cout << "Here is the maximum absolute value of each row:"
+     << endl << m.cwiseAbs().rowwise().maxCoeff() << endl;
diff --git a/doc/snippets/MatrixBase_segment_int_int.cpp b/doc/snippets/MatrixBase_segment_int_int.cpp
new file mode 100644
index 0000000..70cd6d2
--- /dev/null
+++ b/doc/snippets/MatrixBase_segment_int_int.cpp
@@ -0,0 +1,5 @@
+RowVector4i v = RowVector4i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "Here is v.segment(1, 2):" << endl << v.segment(1, 2) << endl;
+v.segment(1, 2).setZero();
+cout << "Now the vector v is:" << endl << v << endl;
diff --git a/doc/snippets/MatrixBase_select.cpp b/doc/snippets/MatrixBase_select.cpp
new file mode 100644
index 0000000..ae5477f
--- /dev/null
+++ b/doc/snippets/MatrixBase_select.cpp
@@ -0,0 +1,6 @@
+MatrixXi m(3, 3);
+m << 1, 2, 3,
+     4, 5, 6,
+     7, 8, 9;
+m = (m.array() >= 5).select(-m, m);
+cout << m << endl;
diff --git a/doc/snippets/MatrixBase_set.cpp b/doc/snippets/MatrixBase_set.cpp
new file mode 100644
index 0000000..50ecf5f
--- /dev/null
+++ b/doc/snippets/MatrixBase_set.cpp
@@ -0,0 +1,13 @@
+Matrix3i m1;
+m1 << 1, 2, 3,
+      4, 5, 6,
+      7, 8, 9;
+cout << m1 << endl << endl;
+Matrix3i m2 = Matrix3i::Identity();
+m2.block(0,0, 2,2) << 10, 11, 12, 13;
+cout << m2 << endl << endl;
+Vector2i v1;
+v1 << 14, 15;
+m2 << v1.transpose(), 16,
+      v1, m1.block(1,1,2,2);
+cout << m2 << endl;
diff --git a/doc/snippets/MatrixBase_setIdentity.cpp b/doc/snippets/MatrixBase_setIdentity.cpp
new file mode 100644
index 0000000..4fd0aa2
--- /dev/null
+++ b/doc/snippets/MatrixBase_setIdentity.cpp
@@ -0,0 +1,3 @@
+Matrix4i m = Matrix4i::Zero();
+m.block<3,3>(1,0).setIdentity();
+cout << m << endl;
diff --git a/doc/snippets/MatrixBase_setOnes.cpp b/doc/snippets/MatrixBase_setOnes.cpp
new file mode 100644
index 0000000..4cef9c1
--- /dev/null
+++ b/doc/snippets/MatrixBase_setOnes.cpp
@@ -0,0 +1,3 @@
+Matrix4i m = Matrix4i::Random();
+m.row(1).setOnes();
+cout << m << endl;
diff --git a/doc/snippets/MatrixBase_setRandom.cpp b/doc/snippets/MatrixBase_setRandom.cpp
new file mode 100644
index 0000000..e2c257d
--- /dev/null
+++ b/doc/snippets/MatrixBase_setRandom.cpp
@@ -0,0 +1,3 @@
+Matrix4i m = Matrix4i::Zero();
+m.col(1).setRandom();
+cout << m << endl;
diff --git a/doc/snippets/MatrixBase_setZero.cpp b/doc/snippets/MatrixBase_setZero.cpp
new file mode 100644
index 0000000..9b5b958
--- /dev/null
+++ b/doc/snippets/MatrixBase_setZero.cpp
@@ -0,0 +1,3 @@
+Matrix4i m = Matrix4i::Random();
+m.row(1).setZero();
+cout << m << endl;
diff --git a/doc/snippets/MatrixBase_start_int.cpp b/doc/snippets/MatrixBase_start_int.cpp
new file mode 100644
index 0000000..c261d2b
--- /dev/null
+++ b/doc/snippets/MatrixBase_start_int.cpp
@@ -0,0 +1,5 @@
+RowVector4i v = RowVector4i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "Here is v.head(2):" << endl << v.head(2) << endl;
+v.head(2).setZero();
+cout << "Now the vector v is:" << endl << v << endl;
diff --git a/doc/snippets/MatrixBase_template_int_bottomRows.cpp b/doc/snippets/MatrixBase_template_int_bottomRows.cpp
new file mode 100644
index 0000000..f9ea892
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_bottomRows.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.bottomRows<2>():" << endl;
+cout << a.bottomRows<2>() << endl;
+a.bottomRows<2>().setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_template_int_end.cpp b/doc/snippets/MatrixBase_template_int_end.cpp
new file mode 100644
index 0000000..f5ccb00
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_end.cpp
@@ -0,0 +1,5 @@
+RowVector4i v = RowVector4i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "Here is v.tail(2):" << endl << v.tail<2>() << endl;
+v.tail<2>().setZero();
+cout << "Now the vector v is:" << endl << v << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_block_int_int_int_int.cpp b/doc/snippets/MatrixBase_template_int_int_block_int_int_int_int.cpp
new file mode 100644
index 0000000..4dced03
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_block_int_int_int_int.cpp
@@ -0,0 +1,5 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the block:" << endl << m.block<2, Dynamic>(1, 1, 2, 3) << endl;
+m.block<2, Dynamic>(1, 1, 2, 3).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_bottomLeftCorner.cpp b/doc/snippets/MatrixBase_template_int_int_bottomLeftCorner.cpp
new file mode 100644
index 0000000..847892a
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_bottomLeftCorner.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.bottomLeftCorner<2,2>():" << endl;
+cout << m.bottomLeftCorner<2,2>() << endl;
+m.bottomLeftCorner<2,2>().setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_bottomLeftCorner_int_int.cpp b/doc/snippets/MatrixBase_template_int_int_bottomLeftCorner_int_int.cpp
new file mode 100644
index 0000000..a1edcc8
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_bottomLeftCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.bottomLeftCorner<2,Dynamic>(2,2):" << endl;
+cout << m.bottomLeftCorner<2,Dynamic>(2,2) << endl;
+m.bottomLeftCorner<2,Dynamic>(2,2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_bottomRightCorner.cpp b/doc/snippets/MatrixBase_template_int_int_bottomRightCorner.cpp
new file mode 100644
index 0000000..abacb01
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_bottomRightCorner.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.bottomRightCorner<2,2>():" << endl;
+cout << m.bottomRightCorner<2,2>() << endl;
+m.bottomRightCorner<2,2>().setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_bottomRightCorner_int_int.cpp b/doc/snippets/MatrixBase_template_int_int_bottomRightCorner_int_int.cpp
new file mode 100644
index 0000000..a65508f
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_bottomRightCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.bottomRightCorner<2,Dynamic>(2,2):" << endl;
+cout << m.bottomRightCorner<2,Dynamic>(2,2) << endl;
+m.bottomRightCorner<2,Dynamic>(2,2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_topLeftCorner.cpp b/doc/snippets/MatrixBase_template_int_int_topLeftCorner.cpp
new file mode 100644
index 0000000..1899d90
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_topLeftCorner.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.topLeftCorner<2,2>():" << endl;
+cout << m.topLeftCorner<2,2>() << endl;
+m.topLeftCorner<2,2>().setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_topLeftCorner_int_int.cpp b/doc/snippets/MatrixBase_template_int_int_topLeftCorner_int_int.cpp
new file mode 100644
index 0000000..fac761f
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_topLeftCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.topLeftCorner<2,Dynamic>(2,2):" << endl;
+cout << m.topLeftCorner<2,Dynamic>(2,2) << endl;
+m.topLeftCorner<2,Dynamic>(2,2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_topRightCorner.cpp b/doc/snippets/MatrixBase_template_int_int_topRightCorner.cpp
new file mode 100644
index 0000000..c3a1771
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_topRightCorner.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.topRightCorner<2,2>():" << endl;
+cout << m.topRightCorner<2,2>() << endl;
+m.topRightCorner<2,2>().setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_int_topRightCorner_int_int.cpp b/doc/snippets/MatrixBase_template_int_int_topRightCorner_int_int.cpp
new file mode 100644
index 0000000..a17acc0
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_int_topRightCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.topRightCorner<2,Dynamic>(2,2):" << endl;
+cout << m.topRightCorner<2,Dynamic>(2,2) << endl;
+m.topRightCorner<2,Dynamic>(2,2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_template_int_leftCols.cpp b/doc/snippets/MatrixBase_template_int_leftCols.cpp
new file mode 100644
index 0000000..1c425d9
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_leftCols.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.leftCols<2>():" << endl;
+cout << a.leftCols<2>() << endl;
+a.leftCols<2>().setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_template_int_rightCols.cpp b/doc/snippets/MatrixBase_template_int_rightCols.cpp
new file mode 100644
index 0000000..fc8c0d9
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_rightCols.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.rightCols<2>():" << endl;
+cout << a.rightCols<2>() << endl;
+a.rightCols<2>().setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_template_int_segment.cpp b/doc/snippets/MatrixBase_template_int_segment.cpp
new file mode 100644
index 0000000..e448b40
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_segment.cpp
@@ -0,0 +1,5 @@
+RowVector4i v = RowVector4i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "Here is v.segment<2>(1):" << endl << v.segment<2>(1) << endl;
+v.segment<2>(2).setZero();
+cout << "Now the vector v is:" << endl << v << endl;
diff --git a/doc/snippets/MatrixBase_template_int_start.cpp b/doc/snippets/MatrixBase_template_int_start.cpp
new file mode 100644
index 0000000..d336b37
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_start.cpp
@@ -0,0 +1,5 @@
+RowVector4i v = RowVector4i::Random();
+cout << "Here is the vector v:" << endl << v << endl;
+cout << "Here is v.head(2):" << endl << v.head<2>() << endl;
+v.head<2>().setZero();
+cout << "Now the vector v is:" << endl << v << endl;
diff --git a/doc/snippets/MatrixBase_template_int_topRows.cpp b/doc/snippets/MatrixBase_template_int_topRows.cpp
new file mode 100644
index 0000000..0110251
--- /dev/null
+++ b/doc/snippets/MatrixBase_template_int_topRows.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.topRows<2>():" << endl;
+cout << a.topRows<2>() << endl;
+a.topRows<2>().setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_topLeftCorner_int_int.cpp b/doc/snippets/MatrixBase_topLeftCorner_int_int.cpp
new file mode 100644
index 0000000..e52cb3b
--- /dev/null
+++ b/doc/snippets/MatrixBase_topLeftCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.topLeftCorner(2, 2):" << endl;
+cout << m.topLeftCorner(2, 2) << endl;
+m.topLeftCorner(2, 2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_topRightCorner_int_int.cpp b/doc/snippets/MatrixBase_topRightCorner_int_int.cpp
new file mode 100644
index 0000000..811fa56
--- /dev/null
+++ b/doc/snippets/MatrixBase_topRightCorner_int_int.cpp
@@ -0,0 +1,6 @@
+Matrix4i m = Matrix4i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is m.topRightCorner(2, 2):" << endl;
+cout << m.topRightCorner(2, 2) << endl;
+m.topRightCorner(2, 2).setZero();
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_topRows_int.cpp b/doc/snippets/MatrixBase_topRows_int.cpp
new file mode 100644
index 0000000..f2d75f1
--- /dev/null
+++ b/doc/snippets/MatrixBase_topRows_int.cpp
@@ -0,0 +1,6 @@
+Array44i a = Array44i::Random();
+cout << "Here is the array a:" << endl << a << endl;
+cout << "Here is a.topRows(2):" << endl;
+cout << a.topRows(2) << endl;
+a.topRows(2).setZero();
+cout << "Now the array a is:" << endl << a << endl;
diff --git a/doc/snippets/MatrixBase_transpose.cpp b/doc/snippets/MatrixBase_transpose.cpp
new file mode 100644
index 0000000..88eea83
--- /dev/null
+++ b/doc/snippets/MatrixBase_transpose.cpp
@@ -0,0 +1,8 @@
+Matrix2i m = Matrix2i::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the transpose of m:" << endl << m.transpose() << endl;
+cout << "Here is the coefficient (1,0) in the transpose of m:" << endl
+     << m.transpose()(1,0) << endl;
+cout << "Let us overwrite this coefficient with the value 0." << endl;
+m.transpose()(1,0) = 0;
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/MatrixBase_zero.cpp b/doc/snippets/MatrixBase_zero.cpp
new file mode 100644
index 0000000..6064936
--- /dev/null
+++ b/doc/snippets/MatrixBase_zero.cpp
@@ -0,0 +1,2 @@
+cout << Matrix2d::Zero() << endl;
+cout << RowVector4i::Zero() << endl;
diff --git a/doc/snippets/MatrixBase_zero_int.cpp b/doc/snippets/MatrixBase_zero_int.cpp
new file mode 100644
index 0000000..370a9ba
--- /dev/null
+++ b/doc/snippets/MatrixBase_zero_int.cpp
@@ -0,0 +1,2 @@
+cout << RowVectorXi::Zero(4) << endl;
+cout << VectorXf::Zero(2) << endl;
diff --git a/doc/snippets/MatrixBase_zero_int_int.cpp b/doc/snippets/MatrixBase_zero_int_int.cpp
new file mode 100644
index 0000000..4099c5d
--- /dev/null
+++ b/doc/snippets/MatrixBase_zero_int_int.cpp
@@ -0,0 +1 @@
+cout << MatrixXi::Zero(2,3) << endl;
diff --git a/doc/snippets/Matrix_resize_NoChange_int.cpp b/doc/snippets/Matrix_resize_NoChange_int.cpp
new file mode 100644
index 0000000..acdf18c
--- /dev/null
+++ b/doc/snippets/Matrix_resize_NoChange_int.cpp
@@ -0,0 +1,3 @@
+MatrixXd m(3,4);
+m.resize(NoChange, 5);
+cout << "m: " << m.rows() << " rows, " << m.cols() << " cols" << endl;
diff --git a/doc/snippets/Matrix_resize_int.cpp b/doc/snippets/Matrix_resize_int.cpp
new file mode 100644
index 0000000..044c789
--- /dev/null
+++ b/doc/snippets/Matrix_resize_int.cpp
@@ -0,0 +1,6 @@
+VectorXd v(10);
+v.resize(3);
+RowVector3d w;
+w.resize(3); // this is legal, but has no effect
+cout << "v: " << v.rows() << " rows, " << v.cols() << " cols" << endl;
+cout << "w: " << w.rows() << " rows, " << w.cols() << " cols" << endl;
diff --git a/doc/snippets/Matrix_resize_int_NoChange.cpp b/doc/snippets/Matrix_resize_int_NoChange.cpp
new file mode 100644
index 0000000..5c37c90
--- /dev/null
+++ b/doc/snippets/Matrix_resize_int_NoChange.cpp
@@ -0,0 +1,3 @@
+MatrixXd m(3,4);
+m.resize(5, NoChange);
+cout << "m: " << m.rows() << " rows, " << m.cols() << " cols" << endl;
diff --git a/doc/snippets/Matrix_resize_int_int.cpp b/doc/snippets/Matrix_resize_int_int.cpp
new file mode 100644
index 0000000..bfd4741
--- /dev/null
+++ b/doc/snippets/Matrix_resize_int_int.cpp
@@ -0,0 +1,9 @@
+MatrixXd m(2,3);
+m << 1,2,3,4,5,6;
+cout << "here's the 2x3 matrix m:" << endl << m << endl;
+cout << "let's resize m to 3x2. This is a conservative resizing because 2*3==3*2." << endl;
+m.resize(3,2);
+cout << "here's the 3x2 matrix m:" << endl << m << endl;
+cout << "now let's resize m to size 2x2. This is NOT a conservative resizing, so it becomes uninitialized:" << endl;
+m.resize(2,2);
+cout << m << endl;
diff --git a/doc/snippets/Matrix_setConstant_int.cpp b/doc/snippets/Matrix_setConstant_int.cpp
new file mode 100644
index 0000000..ff5a86c
--- /dev/null
+++ b/doc/snippets/Matrix_setConstant_int.cpp
@@ -0,0 +1,3 @@
+VectorXf v;
+v.setConstant(3, 5);
+cout << v << endl;
diff --git a/doc/snippets/Matrix_setConstant_int_int.cpp b/doc/snippets/Matrix_setConstant_int_int.cpp
new file mode 100644
index 0000000..32b950c
--- /dev/null
+++ b/doc/snippets/Matrix_setConstant_int_int.cpp
@@ -0,0 +1,3 @@
+MatrixXf m;
+m.setConstant(3, 3, 5);
+cout << m << endl;
diff --git a/doc/snippets/Matrix_setIdentity_int_int.cpp b/doc/snippets/Matrix_setIdentity_int_int.cpp
new file mode 100644
index 0000000..a659671
--- /dev/null
+++ b/doc/snippets/Matrix_setIdentity_int_int.cpp
@@ -0,0 +1,3 @@
+MatrixXf m;
+m.setIdentity(3, 3);
+cout << m << endl;
diff --git a/doc/snippets/Matrix_setOnes_int.cpp b/doc/snippets/Matrix_setOnes_int.cpp
new file mode 100644
index 0000000..752cb35
--- /dev/null
+++ b/doc/snippets/Matrix_setOnes_int.cpp
@@ -0,0 +1,3 @@
+VectorXf v;
+v.setOnes(3);
+cout << v << endl;
diff --git a/doc/snippets/Matrix_setOnes_int_int.cpp b/doc/snippets/Matrix_setOnes_int_int.cpp
new file mode 100644
index 0000000..1ffb66b
--- /dev/null
+++ b/doc/snippets/Matrix_setOnes_int_int.cpp
@@ -0,0 +1,3 @@
+MatrixXf m;
+m.setOnes(3, 3);
+cout << m << endl;
diff --git a/doc/snippets/Matrix_setRandom_int.cpp b/doc/snippets/Matrix_setRandom_int.cpp
new file mode 100644
index 0000000..e160dd7
--- /dev/null
+++ b/doc/snippets/Matrix_setRandom_int.cpp
@@ -0,0 +1,3 @@
+VectorXf v;
+v.setRandom(3);
+cout << v << endl;
diff --git a/doc/snippets/Matrix_setRandom_int_int.cpp b/doc/snippets/Matrix_setRandom_int_int.cpp
new file mode 100644
index 0000000..80cda11
--- /dev/null
+++ b/doc/snippets/Matrix_setRandom_int_int.cpp
@@ -0,0 +1,3 @@
+MatrixXf m;
+m.setRandom(3, 3);
+cout << m << endl;
diff --git a/doc/snippets/Matrix_setZero_int.cpp b/doc/snippets/Matrix_setZero_int.cpp
new file mode 100644
index 0000000..0fb16c1
--- /dev/null
+++ b/doc/snippets/Matrix_setZero_int.cpp
@@ -0,0 +1,3 @@
+VectorXf v;
+v.setZero(3);
+cout << v << endl;
diff --git a/doc/snippets/Matrix_setZero_int_int.cpp b/doc/snippets/Matrix_setZero_int_int.cpp
new file mode 100644
index 0000000..ad883b9
--- /dev/null
+++ b/doc/snippets/Matrix_setZero_int_int.cpp
@@ -0,0 +1,3 @@
+MatrixXf m;
+m.setZero(3, 3);
+cout << m << endl;
diff --git a/doc/snippets/PartialPivLU_solve.cpp b/doc/snippets/PartialPivLU_solve.cpp
new file mode 100644
index 0000000..fa3570a
--- /dev/null
+++ b/doc/snippets/PartialPivLU_solve.cpp
@@ -0,0 +1,7 @@
+MatrixXd A = MatrixXd::Random(3,3);
+MatrixXd B = MatrixXd::Random(3,2);
+cout << "Here is the invertible matrix A:" << endl << A << endl;
+cout << "Here is the matrix B:" << endl << B << endl;
+MatrixXd X = A.lu().solve(B);
+cout << "Here is the (unique) solution X to the equation AX=B:" << endl << X << endl;
+cout << "Relative error: " << (A*X-B).norm() / B.norm() << endl;
diff --git a/doc/snippets/PartialRedux_count.cpp b/doc/snippets/PartialRedux_count.cpp
new file mode 100644
index 0000000..c7b3097
--- /dev/null
+++ b/doc/snippets/PartialRedux_count.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the count of elements larger or equal than 0.5 of each row:" << endl << (m.array() >= 0.5).rowwise().count() << endl;
diff --git a/doc/snippets/PartialRedux_maxCoeff.cpp b/doc/snippets/PartialRedux_maxCoeff.cpp
new file mode 100644
index 0000000..e8fd382
--- /dev/null
+++ b/doc/snippets/PartialRedux_maxCoeff.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the maximum of each column:" << endl << m.colwise().maxCoeff() << endl;
diff --git a/doc/snippets/PartialRedux_minCoeff.cpp b/doc/snippets/PartialRedux_minCoeff.cpp
new file mode 100644
index 0000000..d717bc0
--- /dev/null
+++ b/doc/snippets/PartialRedux_minCoeff.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the minimum of each column:" << endl << m.colwise().minCoeff() << endl;
diff --git a/doc/snippets/PartialRedux_norm.cpp b/doc/snippets/PartialRedux_norm.cpp
new file mode 100644
index 0000000..dbcf290
--- /dev/null
+++ b/doc/snippets/PartialRedux_norm.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the norm of each column:" << endl << m.colwise().norm() << endl;
diff --git a/doc/snippets/PartialRedux_prod.cpp b/doc/snippets/PartialRedux_prod.cpp
new file mode 100644
index 0000000..aacf09c
--- /dev/null
+++ b/doc/snippets/PartialRedux_prod.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the product of each row:" << endl << m.rowwise().prod() << endl;
diff --git a/doc/snippets/PartialRedux_squaredNorm.cpp b/doc/snippets/PartialRedux_squaredNorm.cpp
new file mode 100644
index 0000000..9f3293e
--- /dev/null
+++ b/doc/snippets/PartialRedux_squaredNorm.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the square norm of each row:" << endl << m.rowwise().squaredNorm() << endl;
diff --git a/doc/snippets/PartialRedux_sum.cpp b/doc/snippets/PartialRedux_sum.cpp
new file mode 100644
index 0000000..ec82d3e
--- /dev/null
+++ b/doc/snippets/PartialRedux_sum.cpp
@@ -0,0 +1,3 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the sum of each row:" << endl << m.rowwise().sum() << endl;
diff --git a/doc/snippets/RealQZ_compute.cpp b/doc/snippets/RealQZ_compute.cpp
new file mode 100644
index 0000000..a18da42
--- /dev/null
+++ b/doc/snippets/RealQZ_compute.cpp
@@ -0,0 +1,17 @@
+MatrixXf A = MatrixXf::Random(4,4);
+MatrixXf B = MatrixXf::Random(4,4);
+RealQZ<MatrixXf> qz(4); // preallocate space for 4x4 matrices
+qz.compute(A,B);  // A = Q S Z,  B = Q T Z
+
+// print original matrices and result of decomposition
+cout << "A:\n" << A << "\n" << "B:\n" << B << "\n";
+cout << "S:\n" << qz.matrixS() << "\n" << "T:\n" << qz.matrixT() << "\n";
+cout << "Q:\n" << qz.matrixQ() << "\n" << "Z:\n" << qz.matrixZ() << "\n";
+
+// verify precision
+cout << "\nErrors:"
+  << "\n|A-QSZ|: " << (A-qz.matrixQ()*qz.matrixS()*qz.matrixZ()).norm()
+  << ", |B-QTZ|: " << (B-qz.matrixQ()*qz.matrixT()*qz.matrixZ()).norm()
+  << "\n|QQ* - I|: " << (qz.matrixQ()*qz.matrixQ().adjoint() - MatrixXf::Identity(4,4)).norm()
+  << ", |ZZ* - I|: " << (qz.matrixZ()*qz.matrixZ().adjoint() - MatrixXf::Identity(4,4)).norm()
+  << "\n";
diff --git a/doc/snippets/RealSchur_RealSchur_MatrixType.cpp b/doc/snippets/RealSchur_RealSchur_MatrixType.cpp
new file mode 100644
index 0000000..a5530dc
--- /dev/null
+++ b/doc/snippets/RealSchur_RealSchur_MatrixType.cpp
@@ -0,0 +1,10 @@
+MatrixXd A = MatrixXd::Random(6,6);
+cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
+
+RealSchur<MatrixXd> schur(A);
+cout << "The orthogonal matrix U is:" << endl << schur.matrixU() << endl;
+cout << "The quasi-triangular matrix T is:" << endl << schur.matrixT() << endl << endl;
+
+MatrixXd U = schur.matrixU();
+MatrixXd T = schur.matrixT();
+cout << "U * T * U^T = " << endl << U * T * U.transpose() << endl;
diff --git a/doc/snippets/RealSchur_compute.cpp b/doc/snippets/RealSchur_compute.cpp
new file mode 100644
index 0000000..20c2611
--- /dev/null
+++ b/doc/snippets/RealSchur_compute.cpp
@@ -0,0 +1,6 @@
+MatrixXf A = MatrixXf::Random(4,4);
+RealSchur<MatrixXf> schur(4);
+schur.compute(A, /* computeU = */ false);
+cout << "The matrix T in the decomposition of A is:" << endl << schur.matrixT() << endl;
+schur.compute(A.inverse(), /* computeU = */ false);
+cout << "The matrix T in the decomposition of A^(-1) is:" << endl << schur.matrixT() << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp b/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
new file mode 100644
index 0000000..73a7f62
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
@@ -0,0 +1,7 @@
+SelfAdjointEigenSolver<Matrix4f> es;
+Matrix4f X = Matrix4f::Random(4,4);
+Matrix4f A = X + X.transpose();
+es.compute(A);
+cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl;
+es.compute(A + Matrix4f::Identity(4,4)); // re-use es to compute eigenvalues of A+I
+cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp b/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp
new file mode 100644
index 0000000..3599b17
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp
@@ -0,0 +1,17 @@
+MatrixXd X = MatrixXd::Random(5,5);
+MatrixXd A = X + X.transpose();
+cout << "Here is a random symmetric 5x5 matrix, A:" << endl << A << endl << endl;
+
+SelfAdjointEigenSolver<MatrixXd> es(A);
+cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl;
+cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;
+
+double lambda = es.eigenvalues()[0];
+cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
+VectorXd v = es.eigenvectors().col(0);
+cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
+cout << "... and A * v = " << endl << A * v << endl << endl;
+
+MatrixXd D = es.eigenvalues().asDiagonal();
+MatrixXd V = es.eigenvectors();
+cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp b/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
new file mode 100644
index 0000000..bbb821e
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
@@ -0,0 +1,16 @@
+MatrixXd X = MatrixXd::Random(5,5);
+MatrixXd A = X + X.transpose();
+cout << "Here is a random symmetric matrix, A:" << endl << A << endl;
+X = MatrixXd::Random(5,5);
+MatrixXd B = X * X.transpose();
+cout << "and a random postive-definite matrix, B:" << endl << B << endl << endl;
+
+GeneralizedSelfAdjointEigenSolver<MatrixXd> es(A,B);
+cout << "The eigenvalues of the pencil (A,B) are:" << endl << es.eigenvalues() << endl;
+cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;
+
+double lambda = es.eigenvalues()[0];
+cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
+VectorXd v = es.eigenvectors().col(0);
+cout << "If v is the corresponding eigenvector, then A * v = " << endl << A * v << endl;
+cout << "... and lambda * B * v = " << endl << lambda * B * v << endl << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_compute_MatrixType.cpp b/doc/snippets/SelfAdjointEigenSolver_compute_MatrixType.cpp
new file mode 100644
index 0000000..2975cc3
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_compute_MatrixType.cpp
@@ -0,0 +1,7 @@
+SelfAdjointEigenSolver<MatrixXf> es(4);
+MatrixXf X = MatrixXf::Random(4,4);
+MatrixXf A = X + X.transpose();
+es.compute(A);
+cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl;
+es.compute(A + MatrixXf::Identity(4,4)); // re-use es to compute eigenvalues of A+I
+cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_compute_MatrixType2.cpp b/doc/snippets/SelfAdjointEigenSolver_compute_MatrixType2.cpp
new file mode 100644
index 0000000..07c92a1
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_compute_MatrixType2.cpp
@@ -0,0 +1,9 @@
+MatrixXd X = MatrixXd::Random(5,5);
+MatrixXd A = X * X.transpose();
+X = MatrixXd::Random(5,5);
+MatrixXd B = X * X.transpose();
+
+GeneralizedSelfAdjointEigenSolver<MatrixXd> es(A,B,EigenvaluesOnly);
+cout << "The eigenvalues of the pencil (A,B) are:" << endl << es.eigenvalues() << endl;
+es.compute(B,A,false);
+cout << "The eigenvalues of the pencil (B,A) are:" << endl << es.eigenvalues() << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_eigenvalues.cpp b/doc/snippets/SelfAdjointEigenSolver_eigenvalues.cpp
new file mode 100644
index 0000000..0ff33c6
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_eigenvalues.cpp
@@ -0,0 +1,4 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+SelfAdjointEigenSolver<MatrixXd> es(ones);
+cout << "The eigenvalues of the 3x3 matrix of ones are:" 
+     << endl << es.eigenvalues() << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_eigenvectors.cpp b/doc/snippets/SelfAdjointEigenSolver_eigenvectors.cpp
new file mode 100644
index 0000000..cfc8b0d
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_eigenvectors.cpp
@@ -0,0 +1,4 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+SelfAdjointEigenSolver<MatrixXd> es(ones);
+cout << "The first eigenvector of the 3x3 matrix of ones is:" 
+     << endl << es.eigenvectors().col(1) << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_operatorInverseSqrt.cpp b/doc/snippets/SelfAdjointEigenSolver_operatorInverseSqrt.cpp
new file mode 100644
index 0000000..114c65f
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_operatorInverseSqrt.cpp
@@ -0,0 +1,9 @@
+MatrixXd X = MatrixXd::Random(4,4);
+MatrixXd A = X * X.transpose();
+cout << "Here is a random positive-definite matrix, A:" << endl << A << endl << endl;
+
+SelfAdjointEigenSolver<MatrixXd> es(A);
+cout << "The inverse square root of A is: " << endl;
+cout << es.operatorInverseSqrt() << endl;
+cout << "We can also compute it with operatorSqrt() and inverse(). That yields: " << endl;
+cout << es.operatorSqrt().inverse() << endl;
diff --git a/doc/snippets/SelfAdjointEigenSolver_operatorSqrt.cpp b/doc/snippets/SelfAdjointEigenSolver_operatorSqrt.cpp
new file mode 100644
index 0000000..eeacca7
--- /dev/null
+++ b/doc/snippets/SelfAdjointEigenSolver_operatorSqrt.cpp
@@ -0,0 +1,8 @@
+MatrixXd X = MatrixXd::Random(4,4);
+MatrixXd A = X * X.transpose();
+cout << "Here is a random positive-definite matrix, A:" << endl << A << endl << endl;
+
+SelfAdjointEigenSolver<MatrixXd> es(A);
+MatrixXd sqrtA = es.operatorSqrt();
+cout << "The square root of A is: " << endl << sqrtA << endl;
+cout << "If we square this, we get: " << endl << sqrtA*sqrtA << endl;
diff --git a/doc/snippets/SelfAdjointView_eigenvalues.cpp b/doc/snippets/SelfAdjointView_eigenvalues.cpp
new file mode 100644
index 0000000..be19867
--- /dev/null
+++ b/doc/snippets/SelfAdjointView_eigenvalues.cpp
@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+VectorXd eivals = ones.selfadjointView<Lower>().eigenvalues();
+cout << "The eigenvalues of the 3x3 matrix of ones are:" << endl << eivals << endl;
diff --git a/doc/snippets/SelfAdjointView_operatorNorm.cpp b/doc/snippets/SelfAdjointView_operatorNorm.cpp
new file mode 100644
index 0000000..f380f55
--- /dev/null
+++ b/doc/snippets/SelfAdjointView_operatorNorm.cpp
@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+cout << "The operator norm of the 3x3 matrix of ones is "
+     << ones.selfadjointView<Lower>().operatorNorm() << endl;
diff --git a/doc/snippets/TopicAliasing_block.cpp b/doc/snippets/TopicAliasing_block.cpp
new file mode 100644
index 0000000..03282f4
--- /dev/null
+++ b/doc/snippets/TopicAliasing_block.cpp
@@ -0,0 +1,7 @@
+MatrixXi mat(3,3); 
+mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
+cout << "Here is the matrix mat:\n" << mat << endl;
+
+// This assignment shows the aliasing problem
+mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
+cout << "After the assignment, mat = \n" << mat << endl;
diff --git a/doc/snippets/TopicAliasing_block_correct.cpp b/doc/snippets/TopicAliasing_block_correct.cpp
new file mode 100644
index 0000000..6fee580
--- /dev/null
+++ b/doc/snippets/TopicAliasing_block_correct.cpp
@@ -0,0 +1,7 @@
+MatrixXi mat(3,3); 
+mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
+cout << "Here is the matrix mat:\n" << mat << endl;
+
+// The eval() solves the aliasing problem
+mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval();
+cout << "After the assignment, mat = \n" << mat << endl;
diff --git a/doc/snippets/TopicAliasing_cwise.cpp b/doc/snippets/TopicAliasing_cwise.cpp
new file mode 100644
index 0000000..7049f6c
--- /dev/null
+++ b/doc/snippets/TopicAliasing_cwise.cpp
@@ -0,0 +1,20 @@
+MatrixXf mat(2,2); 
+mat << 1, 2,  4, 7;
+cout << "Here is the matrix mat:\n" << mat << endl << endl;
+
+mat = 2 * mat;
+cout << "After 'mat = 2 * mat', mat = \n" << mat << endl << endl;
+
+
+mat = mat - MatrixXf::Identity(2,2);
+cout << "After the subtraction, it becomes\n" << mat << endl << endl;
+
+
+ArrayXXf arr = mat;
+arr = arr.square();
+cout << "After squaring, it becomes\n" << arr << endl << endl;
+
+// Combining all operations in one statement:
+mat << 1, 2,  4, 7;
+mat = (2 * mat - MatrixXf::Identity(2,2)).array().square();
+cout << "Doing everything at once yields\n" << mat << endl << endl;
diff --git a/doc/snippets/TopicAliasing_mult1.cpp b/doc/snippets/TopicAliasing_mult1.cpp
new file mode 100644
index 0000000..cd7e900
--- /dev/null
+++ b/doc/snippets/TopicAliasing_mult1.cpp
@@ -0,0 +1,4 @@
+MatrixXf matA(2,2); 
+matA << 2, 0,  0, 2;
+matA = matA * matA;
+cout << matA;
diff --git a/doc/snippets/TopicAliasing_mult2.cpp b/doc/snippets/TopicAliasing_mult2.cpp
new file mode 100644
index 0000000..a3ff568
--- /dev/null
+++ b/doc/snippets/TopicAliasing_mult2.cpp
@@ -0,0 +1,10 @@
+MatrixXf matA(2,2), matB(2,2); 
+matA << 2, 0,  0, 2;
+
+// Simple but not quite as efficient
+matB = matA * matA;
+cout << matB << endl << endl;
+
+// More complicated but also more efficient
+matB.noalias() = matA * matA;
+cout << matB;
diff --git a/doc/snippets/TopicAliasing_mult3.cpp b/doc/snippets/TopicAliasing_mult3.cpp
new file mode 100644
index 0000000..1d12a6c
--- /dev/null
+++ b/doc/snippets/TopicAliasing_mult3.cpp
@@ -0,0 +1,4 @@
+MatrixXf matA(2,2); 
+matA << 2, 0,  0, 2;
+matA.noalias() = matA * matA;
+cout << matA;
diff --git a/doc/snippets/TopicStorageOrders_example.cpp b/doc/snippets/TopicStorageOrders_example.cpp
new file mode 100644
index 0000000..0623ef0
--- /dev/null
+++ b/doc/snippets/TopicStorageOrders_example.cpp
@@ -0,0 +1,18 @@
+Matrix<int, 3, 4, ColMajor> Acolmajor;
+Acolmajor << 8, 2, 2, 9,
+             9, 1, 4, 4,
+	     3, 5, 4, 5;
+cout << "The matrix A:" << endl;
+cout << Acolmajor << endl << endl; 
+
+cout << "In memory (column-major):" << endl;
+for (int i = 0; i < Acolmajor.size(); i++)
+  cout << *(Acolmajor.data() + i) << "  ";
+cout << endl << endl;
+
+Matrix<int, 3, 4, RowMajor> Arowmajor = Acolmajor;
+cout << "In memory (row-major):" << endl;
+for (int i = 0; i < Arowmajor.size(); i++)
+  cout << *(Arowmajor.data() + i) << "  ";
+cout << endl;
+
diff --git a/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp b/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp
new file mode 100644
index 0000000..a260124
--- /dev/null
+++ b/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp
@@ -0,0 +1,9 @@
+MatrixXd X = MatrixXd::Random(5,5);
+MatrixXd A = X + X.transpose();
+cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl;
+Tridiagonalization<MatrixXd> triOfA(A);
+MatrixXd Q = triOfA.matrixQ();
+cout << "The orthogonal matrix Q is:" << endl << Q << endl;
+MatrixXd T = triOfA.matrixT();
+cout << "The tridiagonal matrix T is:" << endl << T << endl << endl;
+cout << "Q * T * Q^T = " << endl << Q * T * Q.transpose() << endl;
diff --git a/doc/snippets/Tridiagonalization_compute.cpp b/doc/snippets/Tridiagonalization_compute.cpp
new file mode 100644
index 0000000..0062a99
--- /dev/null
+++ b/doc/snippets/Tridiagonalization_compute.cpp
@@ -0,0 +1,9 @@
+Tridiagonalization<MatrixXf> tri;
+MatrixXf X = MatrixXf::Random(4,4);
+MatrixXf A = X + X.transpose();
+tri.compute(A);
+cout << "The matrix T in the tridiagonal decomposition of A is: " << endl;
+cout << tri.matrixT() << endl;
+tri.compute(2*A); // re-use tri to compute eigenvalues of 2A
+cout << "The matrix T in the tridiagonal decomposition of 2A is: " << endl;
+cout << tri.matrixT() << endl;
diff --git a/doc/snippets/Tridiagonalization_decomposeInPlace.cpp b/doc/snippets/Tridiagonalization_decomposeInPlace.cpp
new file mode 100644
index 0000000..93dcfca
--- /dev/null
+++ b/doc/snippets/Tridiagonalization_decomposeInPlace.cpp
@@ -0,0 +1,10 @@
+MatrixXd X = MatrixXd::Random(5,5);
+MatrixXd A = X + X.transpose();
+cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl;
+
+VectorXd diag(5);
+VectorXd subdiag(4);
+internal::tridiagonalization_inplace(A, diag, subdiag, true);
+cout << "The orthogonal matrix Q is:" << endl << A << endl;
+cout << "The diagonal of the tridiagonal matrix T is:" << endl << diag << endl;
+cout << "The subdiagonal of the tridiagonal matrix T is:" << endl << subdiag << endl;
diff --git a/doc/snippets/Tridiagonalization_diagonal.cpp b/doc/snippets/Tridiagonalization_diagonal.cpp
new file mode 100644
index 0000000..6eec821
--- /dev/null
+++ b/doc/snippets/Tridiagonalization_diagonal.cpp
@@ -0,0 +1,13 @@
+MatrixXcd X = MatrixXcd::Random(4,4);
+MatrixXcd A = X + X.adjoint();
+cout << "Here is a random self-adjoint 4x4 matrix:" << endl << A << endl << endl;
+
+Tridiagonalization<MatrixXcd> triOfA(A);
+MatrixXd T = triOfA.matrixT();
+cout << "The tridiagonal matrix T is:" << endl << T << endl << endl;
+
+cout << "We can also extract the diagonals of T directly ..." << endl;
+VectorXd diag = triOfA.diagonal();
+cout << "The diagonal is:" << endl << diag << endl; 
+VectorXd subdiag = triOfA.subDiagonal();
+cout << "The subdiagonal is:" << endl << subdiag << endl;
diff --git a/doc/snippets/Tridiagonalization_householderCoefficients.cpp b/doc/snippets/Tridiagonalization_householderCoefficients.cpp
new file mode 100644
index 0000000..e5d8728
--- /dev/null
+++ b/doc/snippets/Tridiagonalization_householderCoefficients.cpp
@@ -0,0 +1,6 @@
+Matrix4d X = Matrix4d::Random(4,4);
+Matrix4d A = X + X.transpose();
+cout << "Here is a random symmetric 4x4 matrix:" << endl << A << endl;
+Tridiagonalization<Matrix4d> triOfA(A);
+Vector3d hc = triOfA.householderCoefficients();
+cout << "The vector of Householder coefficients is:" << endl << hc << endl;
diff --git a/doc/snippets/Tridiagonalization_packedMatrix.cpp b/doc/snippets/Tridiagonalization_packedMatrix.cpp
new file mode 100644
index 0000000..0f55d0c
--- /dev/null
+++ b/doc/snippets/Tridiagonalization_packedMatrix.cpp
@@ -0,0 +1,8 @@
+Matrix4d X = Matrix4d::Random(4,4);
+Matrix4d A = X + X.transpose();
+cout << "Here is a random symmetric 4x4 matrix:" << endl << A << endl;
+Tridiagonalization<Matrix4d> triOfA(A);
+Matrix4d pm = triOfA.packedMatrix();
+cout << "The packed matrix M is:" << endl << pm << endl;
+cout << "The diagonal and subdiagonal corresponds to the matrix T, which is:" 
+     << endl << triOfA.matrixT() << endl;
diff --git a/doc/snippets/Tutorial_AdvancedInitialization_Block.cpp b/doc/snippets/Tutorial_AdvancedInitialization_Block.cpp
new file mode 100644
index 0000000..96e40ac
--- /dev/null
+++ b/doc/snippets/Tutorial_AdvancedInitialization_Block.cpp
@@ -0,0 +1,5 @@
+MatrixXf matA(2, 2);
+matA << 1, 2, 3, 4;
+MatrixXf matB(4, 4);
+matB << matA, matA/10, matA/10, matA;
+std::cout << matB << std::endl;
diff --git a/doc/snippets/Tutorial_AdvancedInitialization_CommaTemporary.cpp b/doc/snippets/Tutorial_AdvancedInitialization_CommaTemporary.cpp
new file mode 100644
index 0000000..50cff4c
--- /dev/null
+++ b/doc/snippets/Tutorial_AdvancedInitialization_CommaTemporary.cpp
@@ -0,0 +1,4 @@
+MatrixXf mat = MatrixXf::Random(2, 3);
+std::cout << mat << std::endl << std::endl;
+mat = (MatrixXf(2,2) << 0, 1, 1, 0).finished() * mat;
+std::cout << mat << std::endl;
diff --git a/doc/snippets/Tutorial_AdvancedInitialization_Join.cpp b/doc/snippets/Tutorial_AdvancedInitialization_Join.cpp
new file mode 100644
index 0000000..84e8715
--- /dev/null
+++ b/doc/snippets/Tutorial_AdvancedInitialization_Join.cpp
@@ -0,0 +1,11 @@
+RowVectorXd vec1(3);
+vec1 << 1, 2, 3;
+std::cout << "vec1 = " << vec1 << std::endl;
+
+RowVectorXd vec2(4);
+vec2 << 1, 4, 9, 16;;
+std::cout << "vec2 = " << vec2 << std::endl;
+
+RowVectorXd joined(7);
+joined << vec1, vec2;
+std::cout << "joined = " << joined << std::endl;
diff --git a/doc/snippets/Tutorial_AdvancedInitialization_LinSpaced.cpp b/doc/snippets/Tutorial_AdvancedInitialization_LinSpaced.cpp
new file mode 100644
index 0000000..c6a73ab
--- /dev/null
+++ b/doc/snippets/Tutorial_AdvancedInitialization_LinSpaced.cpp
@@ -0,0 +1,7 @@
+ArrayXXf table(10, 4);
+table.col(0) = ArrayXf::LinSpaced(10, 0, 90);
+table.col(1) = M_PI / 180 * table.col(0);
+table.col(2) = table.col(1).sin();
+table.col(3) = table.col(1).cos();
+std::cout << "  Degrees   Radians      Sine    Cosine\n";
+std::cout << table << std::endl;
diff --git a/doc/snippets/Tutorial_AdvancedInitialization_ThreeWays.cpp b/doc/snippets/Tutorial_AdvancedInitialization_ThreeWays.cpp
new file mode 100644
index 0000000..cb74576
--- /dev/null
+++ b/doc/snippets/Tutorial_AdvancedInitialization_ThreeWays.cpp
@@ -0,0 +1,20 @@
+const int size = 6;
+MatrixXd mat1(size, size);
+mat1.topLeftCorner(size/2, size/2)     = MatrixXd::Zero(size/2, size/2);
+mat1.topRightCorner(size/2, size/2)    = MatrixXd::Identity(size/2, size/2);
+mat1.bottomLeftCorner(size/2, size/2)  = MatrixXd::Identity(size/2, size/2);
+mat1.bottomRightCorner(size/2, size/2) = MatrixXd::Zero(size/2, size/2);
+std::cout << mat1 << std::endl << std::endl;
+
+MatrixXd mat2(size, size);
+mat2.topLeftCorner(size/2, size/2).setZero();
+mat2.topRightCorner(size/2, size/2).setIdentity();
+mat2.bottomLeftCorner(size/2, size/2).setIdentity();
+mat2.bottomRightCorner(size/2, size/2).setZero();
+std::cout << mat2 << std::endl << std::endl;
+
+MatrixXd mat3(size, size);
+mat3 << MatrixXd::Zero(size/2, size/2), MatrixXd::Identity(size/2, size/2),
+        MatrixXd::Identity(size/2, size/2), MatrixXd::Zero(size/2, size/2);
+std::cout << mat3 << std::endl;
+
diff --git a/doc/snippets/Tutorial_AdvancedInitialization_Zero.cpp b/doc/snippets/Tutorial_AdvancedInitialization_Zero.cpp
new file mode 100644
index 0000000..76a36a3
--- /dev/null
+++ b/doc/snippets/Tutorial_AdvancedInitialization_Zero.cpp
@@ -0,0 +1,13 @@
+std::cout << "A fixed-size array:\n";
+Array33f a1 = Array33f::Zero();
+std::cout << a1 << "\n\n";
+
+
+std::cout << "A one-dimensional dynamic-size array:\n";
+ArrayXf a2 = ArrayXf::Zero(3);
+std::cout << a2 << "\n\n";
+
+
+std::cout << "A two-dimensional dynamic-size array:\n";
+ArrayXXf a3 = ArrayXXf::Zero(3, 4);
+std::cout << a3 << "\n";
diff --git a/doc/snippets/Tutorial_Map_rowmajor.cpp b/doc/snippets/Tutorial_Map_rowmajor.cpp
new file mode 100644
index 0000000..fd45ace
--- /dev/null
+++ b/doc/snippets/Tutorial_Map_rowmajor.cpp
@@ -0,0 +1,7 @@
+int array[8];
+for(int i = 0; i < 8; ++i) array[i] = i;
+cout << "Column-major:\n" << Map<Matrix<int,2,4> >(array) << endl;
+cout << "Row-major:\n" << Map<Matrix<int,2,4,RowMajor> >(array) << endl;
+cout << "Row-major using stride:\n" <<
+  Map<Matrix<int,2,4>, Unaligned, Stride<1,4> >(array) << endl;
+
diff --git a/doc/snippets/Tutorial_Map_using.cpp b/doc/snippets/Tutorial_Map_using.cpp
new file mode 100644
index 0000000..e5e499f
--- /dev/null
+++ b/doc/snippets/Tutorial_Map_using.cpp
@@ -0,0 +1,21 @@
+typedef Matrix<float,1,Dynamic> MatrixType;
+typedef Map<MatrixType> MapType;
+typedef Map<const MatrixType> MapTypeConst;   // a read-only map
+const int n_dims = 5;
+  
+MatrixType m1(n_dims), m2(n_dims);
+m1.setRandom();
+m2.setRandom();
+float *p = &m2(0);  // get the address storing the data for m2
+MapType m2map(p,m2.size());   // m2map shares data with m2
+MapTypeConst m2mapconst(p,m2.size());  // a read-only accessor for m2
+
+cout << "m1: " << m1 << endl;
+cout << "m2: " << m2 << endl;
+cout << "Squared euclidean distance: " << (m1-m2).squaredNorm() << endl;
+cout << "Squared euclidean distance, using map: " <<
+  (m1-m2map).squaredNorm() << endl;
+m2map(3) = 7;   // this will change m2, since they share the same array
+cout << "Updated m2: " << m2 << endl;
+cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << endl;
+/* m2mapconst(2) = 5; */   // this yields a compile-time error
diff --git a/doc/snippets/Tutorial_commainit_01.cpp b/doc/snippets/Tutorial_commainit_01.cpp
new file mode 100644
index 0000000..47ba31d
--- /dev/null
+++ b/doc/snippets/Tutorial_commainit_01.cpp
@@ -0,0 +1,5 @@
+Matrix3f m;
+m << 1, 2, 3,
+     4, 5, 6,
+     7, 8, 9;
+std::cout << m;
diff --git a/doc/snippets/Tutorial_commainit_01b.cpp b/doc/snippets/Tutorial_commainit_01b.cpp
new file mode 100644
index 0000000..2adb2e2
--- /dev/null
+++ b/doc/snippets/Tutorial_commainit_01b.cpp
@@ -0,0 +1,5 @@
+Matrix3f m;
+m.row(0) << 1, 2, 3;
+m.block(1,0,2,2) << 4, 5, 7, 8;
+m.col(2).tail(2) << 6, 9;		    
+std::cout << m;
diff --git a/doc/snippets/Tutorial_commainit_02.cpp b/doc/snippets/Tutorial_commainit_02.cpp
new file mode 100644
index 0000000..c960d6a
--- /dev/null
+++ b/doc/snippets/Tutorial_commainit_02.cpp
@@ -0,0 +1,7 @@
+int rows=5, cols=5;
+MatrixXf m(rows,cols);
+m << (Matrix3f() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished(),
+     MatrixXf::Zero(3,cols-3),
+     MatrixXf::Zero(rows-3,3),
+     MatrixXf::Identity(rows-3,cols-3);
+cout << m;
diff --git a/doc/snippets/Tutorial_solve_matrix_inverse.cpp b/doc/snippets/Tutorial_solve_matrix_inverse.cpp
new file mode 100644
index 0000000..fff3244
--- /dev/null
+++ b/doc/snippets/Tutorial_solve_matrix_inverse.cpp
@@ -0,0 +1,6 @@
+Matrix3f A;
+Vector3f b;
+A << 1,2,3,  4,5,6,  7,8,10;
+b << 3, 3, 4;
+Vector3f x = A.inverse() * b;
+cout << "The solution is:" << endl << x << endl;
diff --git a/doc/snippets/Tutorial_solve_multiple_rhs.cpp b/doc/snippets/Tutorial_solve_multiple_rhs.cpp
new file mode 100644
index 0000000..5411a44
--- /dev/null
+++ b/doc/snippets/Tutorial_solve_multiple_rhs.cpp
@@ -0,0 +1,10 @@
+Matrix3f A(3,3);
+A << 1,2,3,  4,5,6,  7,8,10;
+Matrix<float,3,2> B;
+B << 3,1, 3,1, 4,1;
+Matrix<float,3,2> X;
+X = A.fullPivLu().solve(B);
+cout << "The solution with right-hand side (3,3,4) is:" << endl;
+cout << X.col(0) << endl;
+cout << "The solution with right-hand side (1,1,1) is:" << endl;
+cout << X.col(1) << endl;
diff --git a/doc/snippets/Tutorial_solve_reuse_decomposition.cpp b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp
new file mode 100644
index 0000000..3ca0645
--- /dev/null
+++ b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp
@@ -0,0 +1,13 @@
+Matrix3f A(3,3);
+A << 1,2,3,  4,5,6,  7,8,10;
+PartialPivLU<Matrix3f> luOfA(A); // compute LU decomposition of A
+Vector3f b;
+b << 3,3,4;
+Vector3f x;
+x = luOfA.solve(b);
+cout << "The solution with right-hand side (3,3,4) is:" << endl;
+cout << x << endl;
+b << 1,1,1;
+x = luOfA.solve(b);
+cout << "The solution with right-hand side (1,1,1) is:" << endl;
+cout << x << endl;
diff --git a/doc/snippets/Tutorial_solve_singular.cpp b/doc/snippets/Tutorial_solve_singular.cpp
new file mode 100644
index 0000000..abff1ef
--- /dev/null
+++ b/doc/snippets/Tutorial_solve_singular.cpp
@@ -0,0 +1,9 @@
+Matrix3f A;
+Vector3f b;
+A << 1,2,3,  4,5,6,  7,8,9;
+b << 3, 3, 4;
+cout << "Here is the matrix A:" << endl << A << endl;
+cout << "Here is the vector b:" << endl << b << endl;
+Vector3f x;
+x = A.lu().solve(b);
+cout << "The solution is:" << endl << x << endl;
diff --git a/doc/snippets/Tutorial_solve_triangular.cpp b/doc/snippets/Tutorial_solve_triangular.cpp
new file mode 100644
index 0000000..9d13f22
--- /dev/null
+++ b/doc/snippets/Tutorial_solve_triangular.cpp
@@ -0,0 +1,8 @@
+Matrix3f A;
+Vector3f b;
+A << 1,2,3,  0,5,6,  0,0,10;
+b << 3, 3, 4;
+cout << "Here is the matrix A:" << endl << A << endl;
+cout << "Here is the vector b:" << endl << b << endl;
+Vector3f x = A.triangularView<Upper>().solve(b);
+cout << "The solution is:" << endl << x << endl;
diff --git a/doc/snippets/Tutorial_solve_triangular_inplace.cpp b/doc/snippets/Tutorial_solve_triangular_inplace.cpp
new file mode 100644
index 0000000..16ae633
--- /dev/null
+++ b/doc/snippets/Tutorial_solve_triangular_inplace.cpp
@@ -0,0 +1,6 @@
+Matrix3f A;
+Vector3f b;
+A << 1,2,3,  0,5,6,  0,0,10;
+b << 3, 3, 4;
+A.triangularView<Upper>().solveInPlace(b);
+cout << "The solution is:" << endl << b << endl;
diff --git a/doc/snippets/Vectorwise_reverse.cpp b/doc/snippets/Vectorwise_reverse.cpp
new file mode 100644
index 0000000..2f6a350
--- /dev/null
+++ b/doc/snippets/Vectorwise_reverse.cpp
@@ -0,0 +1,10 @@
+MatrixXi m = MatrixXi::Random(3,4);
+cout << "Here is the matrix m:" << endl << m << endl;
+cout << "Here is the rowwise reverse of m:" << endl << m.rowwise().reverse() << endl;
+cout << "Here is the colwise reverse of m:" << endl << m.colwise().reverse() << endl;
+
+cout << "Here is the coefficient (1,0) in the rowise reverse of m:" << endl
+<< m.rowwise().reverse()(1,0) << endl;
+cout << "Let us overwrite this coefficient with the value 4." << endl;
+//m.colwise().reverse()(1,0) = 4;
+cout << "Now the matrix m is:" << endl << m << endl;
diff --git a/doc/snippets/class_FullPivLU.cpp b/doc/snippets/class_FullPivLU.cpp
new file mode 100644
index 0000000..fce7fac
--- /dev/null
+++ b/doc/snippets/class_FullPivLU.cpp
@@ -0,0 +1,16 @@
+typedef Matrix<double, 5, 3> Matrix5x3;
+typedef Matrix<double, 5, 5> Matrix5x5;
+Matrix5x3 m = Matrix5x3::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+Eigen::FullPivLU<Matrix5x3> lu(m);
+cout << "Here is, up to permutations, its LU decomposition matrix:"
+     << endl << lu.matrixLU() << endl;
+cout << "Here is the L part:" << endl;
+Matrix5x5 l = Matrix5x5::Identity();
+l.block<5,3>(0,0).triangularView<StrictlyLower>() = lu.matrixLU();
+cout << l << endl;
+cout << "Here is the U part:" << endl;
+Matrix5x3 u = lu.matrixLU().triangularView<Upper>();
+cout << u << endl;
+cout << "Let us now reconstruct the original matrix m:" << endl;
+cout << lu.permutationP().inverse() * l * u * lu.permutationQ().inverse() << endl;
diff --git a/doc/snippets/compile_snippet.cpp.in b/doc/snippets/compile_snippet.cpp.in
new file mode 100644
index 0000000..894cd52
--- /dev/null
+++ b/doc/snippets/compile_snippet.cpp.in
@@ -0,0 +1,12 @@
+#include <Eigen/Dense>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main(int, char**)
+{
+  cout.precision(3);
+  ${snippet_source_code}
+  return 0;
+}
diff --git a/doc/snippets/tut_arithmetic_redux_minmax.cpp b/doc/snippets/tut_arithmetic_redux_minmax.cpp
new file mode 100644
index 0000000..f4ae7f4
--- /dev/null
+++ b/doc/snippets/tut_arithmetic_redux_minmax.cpp
@@ -0,0 +1,12 @@
+  Matrix3f m = Matrix3f::Random();
+  std::ptrdiff_t i, j;
+  float minOfM = m.minCoeff(&i,&j);
+  cout << "Here is the matrix m:\n" << m << endl;
+  cout << "Its minimum coefficient (" << minOfM 
+       << ") is at position (" << i << "," << j << ")\n\n";
+
+  RowVector4i v = RowVector4i::Random();
+  int maxOfV = v.maxCoeff(&i);
+  cout << "Here is the vector v: " << v << endl;
+  cout << "Its maximum coefficient (" << maxOfV 
+       << ") is at position " << i << endl;
diff --git a/doc/snippets/tut_arithmetic_transpose_aliasing.cpp b/doc/snippets/tut_arithmetic_transpose_aliasing.cpp
new file mode 100644
index 0000000..c8e4746
--- /dev/null
+++ b/doc/snippets/tut_arithmetic_transpose_aliasing.cpp
@@ -0,0 +1,5 @@
+Matrix2i a; a << 1, 2, 3, 4;
+cout << "Here is the matrix a:\n" << a << endl;
+
+a = a.transpose(); // !!! do NOT do this !!!
+cout << "and the result of the aliasing effect:\n" << a << endl;
\ No newline at end of file
diff --git a/doc/snippets/tut_arithmetic_transpose_conjugate.cpp b/doc/snippets/tut_arithmetic_transpose_conjugate.cpp
new file mode 100644
index 0000000..88496b2
--- /dev/null
+++ b/doc/snippets/tut_arithmetic_transpose_conjugate.cpp
@@ -0,0 +1,12 @@
+MatrixXcf a = MatrixXcf::Random(2,2);
+cout << "Here is the matrix a\n" << a << endl;
+
+cout << "Here is the matrix a^T\n" << a.transpose() << endl;
+
+
+cout << "Here is the conjugate of a\n" << a.conjugate() << endl;
+
+
+cout << "Here is the matrix a^*\n" << a.adjoint() << endl;
+
+
diff --git a/doc/snippets/tut_arithmetic_transpose_inplace.cpp b/doc/snippets/tut_arithmetic_transpose_inplace.cpp
new file mode 100644
index 0000000..7a069ff
--- /dev/null
+++ b/doc/snippets/tut_arithmetic_transpose_inplace.cpp
@@ -0,0 +1,6 @@
+MatrixXf a(2,3); a << 1, 2, 3, 4, 5, 6;
+cout << "Here is the initial matrix a:\n" << a << endl;
+
+
+a.transposeInPlace();
+cout << "and after being transposed:\n" << a << endl;
\ No newline at end of file
diff --git a/doc/snippets/tut_matrix_assignment_resizing.cpp b/doc/snippets/tut_matrix_assignment_resizing.cpp
new file mode 100644
index 0000000..cf18998
--- /dev/null
+++ b/doc/snippets/tut_matrix_assignment_resizing.cpp
@@ -0,0 +1,5 @@
+MatrixXf a(2,2);
+std::cout << "a is of size " << a.rows() << "x" << a.cols() << std::endl;
+MatrixXf b(3,3);
+a = b;
+std::cout << "a is now of size " << a.rows() << "x" << a.cols() << std::endl;