Best algorithm for computing eigenvalue decomposition of a $3 \times 3$ symmetrix matrix In one of my applications, I need to compute the eigenvalue decomposition of a $3 \times 3$ symmetric matrix. What algorithms can I use? Which is the most efficient one? 
More specifically, the matrix is a Hessian matrix of a 3d function $f(x_1,x_2,x_3)$. 
$$\begin{bmatrix}
\dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1\,\partial x_2} & \dfrac{\partial^2 f}{\partial x_1\,\partial x_3} \\[2.2ex]
\dfrac{\partial^2 f}{\partial x_2\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \dfrac{\partial^2 f}{\partial x_2\,\partial x_3} \\[2.2ex]
\dfrac{\partial^2 f}{\partial x_3\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_3\,\partial x_2} & \dfrac{\partial^2 f}{\partial x_3^2}
\end{bmatrix}
$$
 A: If the matrix is positive definite, then I find that the following is very effective. I will provide justification also.
The largest eigenvalue is greater than the largest diagonal. So largest diagonal entry is a good approximation to the largest eigenvalue. Pick a random initial guess for eigenvector. 
If $e$ is a guessed eigenvector, and $\mu$ an approximation to an eigenvalue, then the following inverse power sequence
$$
\hat e = \left(A - \mu I\right)^{-1} e \\
\hat e = \frac{\hat e}{\left|| e \right||} ~\text{any norm is okay here.}\\
\hat \mu = \frac{\hat e^T A \hat e} {\hat e^T \hat e} 
$$
is a better approximation. When inverting, you can avoid division by zero by adding a small quantity to the pivots (I use 10^{-10}). 
So starting with the largest diagonal, you can get the largest eigenvalue in a 2 to 3 iterations (often less).
The smallest eigenvalue is less than the smallest diagonal entry. So starting with that as your initial guess, you can get the smallest eigenvalue and eigenvector pair.
Since sum of the diagonals equals the sum of the eigenvalues, you can get the third eigenvalue without too much calculations. You can iterate to get better accuracy if you want.
Every time you have eigenvectors, the remaining eigenvector(s) have to be perpendicular to the ones you have. So a simple Gramm-Schmidt process can be used to find an initial guess for eigenvector but I find that a random vector is just as effective.
Just saw the update
If the matrix is not positive definite, then you have to have some way to guess an eigenvalue. Here is what I do in this case.


*

*Pick a random vector and let initial guess of eigenvalue as $\mu = (e^T A e)/(e^Te)$ and proceed as before.

*Once you have isolated an eigenvalue and eigenvector, pick a new $e$ that is orthogonal to the known eigenvector and proceed as in 1 above

*Once you have two eigenvalues and eigenectors, you can find the third easily. The eigenvector will be orthogonal to the other two eigenvectors and sum of the eigenvalues must be equal to sum of the diagonals.


If you want a worked out example, let me know.
