Open awvwgk opened 3 years ago
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.
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.
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).
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.
:+1:
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.