Eigenvalues/eigenvectors for small matrix sizes

[awvwgk]

Sometimes I have need for calculating eigenvalues or eigenvectors for small symmetric matrices, usually 3×3, which feels a bit overkill to depend on LAPACK for such a project, especially given the complications of correctly linking tuned libraries like MKL.

For this purpose I implemented the analytical formula of the eigenvalues and an algorithm based on cross products to determine the eigenvectors as well and released the implemented to the public domain: https://github.com/awvwgk/diag3x3.

If there is interest I can contribute this version to stdlib. I'm using this in quite a few projects and also have some tests for the implementation rather than just the single file.

[ivan-pi]

I think would this beneficial, even if I am more interested in the simpler solution to Ax=b. But I have at least one problem, which could use the eig3by3 routine. In the long run, I guess the routines could be generalized for complex numbers and non-symmetric problems?

Would it to be responsibility of the user to assure the eigenvalue problem is consistent? (Meaning that rank(A) == 3.)

I once wrote my own Broyden's method to solve a small system of 3 nonlinear equations. The iteration involved solving a 3x3 system of linear equation for the update vector which involves the Jacobian matrix. At the time I was still new to Fortran, and a procedure like this would have saved me some trouble.

[ivan-pi]

One could also use this as quick solver for the (real) roots of a third order polynomial.

E.g. for the polynomial 6 - 5*x - 2*x**2 + x**3 you just need to find the eigenvalues of the matrix:

real :: a(3,3), roots(3)
a = diag(3,1) ! ones above diagonal
a(3,:) = [-6, 5, 2]

call eigval3x3(a,roots)

Edit: just noticed this requires a non-symmetric eigenvalue solver.

[awvwgk]

Would it to be responsibility of the user to assure the eigenvalue problem is consistent? (Meaning that rank(A) == 3.)

I feel like the answer should be yes, passing a general 4×4 matrix to a syev_3x3 routine might be out of our responsibility.

I would this beneficial, even if I am more interested in the simpler solution to Ax=b.

Symmetric linear equations would be in-scope. So far I was using those routines to find principal axes for Cartesian geometries, which requires the eigenvectors of a 3×3 symmetric eigenvalue problem (preferably in pure way).

[ivan-pi]

Following your usage description, I see this can be used in (solid) mechanics to calculate the principal stresses and stress invariants. Given a displacement vector field (and its gradient), one could use this routine to compute the scalar von Mises stress.

👍

[perazz]

Addressed by #816