Document SelfAdjointEigenSolver and add examples.
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..f05d67d
--- /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;
+
+SelfAdjointEigenSolver<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..4b0f110
--- /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();
+
+SelfAdjointEigenSolver<MatrixXd> es(A,B,false);
+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;