linalg: Eigenvalues and Eigenvectors

[perazz]

Computing eigenvalues and eigenvectors of a square matrix: $A \cdot \bar{v} - \lambda \cdot \bar{v} = 0$ . Based on LAPACK General (*GEEV) or real symmetric (*SYEV) / complex Hermitian (*HEEV) operations.

• [x] base implementation
• [x] tests
• [x] exclude unsupported xdp
• [x] documentation
• [x] submodule
• [x] examples
• [x] subroutine interface

General matrix

Prior art

• Numpy: eigenvalues, eigenvectors = eig(a)
• Numpy: eigenvalues = eigvals(a)
• Scipy: eig(a, b=None, left=False, right=True, overwrite_a=False, overwrite_b=False, check_finite=True, homogeneous_eigvals=False)

Proposed implementation

• Eigenvalues only: two function interfaces lambda = eigvals(A [, err])

• Eigendecomposition: subroutine interface call eig(a,lambda [,right] [,left] [,overwrite_a] [,err])

  • Eigenvalues of real matrices can be complex-conjugate. To avoid confusion, right, left are postprocessed and always returned as complex arrays. LAPACK instead stores (re,im) of conjugates as consecutive real values [j, j+1] in a real array, this would be very confusing and hard to track down.

  • The eigendecomposition can be used either for eigenvalues only, or to retrieve left or right eigenvectors.

Real Symmetric / Complex Hermitian

Prior art

• Numpy: eigenvalues, eigenvectors = eigh(a, uplo='L')
• Numpy: eigenvalues = eigvalsh(a)
• Scipy: eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=<object object>, eigvals=<object object>, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None)

Proposed implementation

• Eigenvalues only: two function interfaces lambda = eigvalsh(A [, err])

• Eigendecomposition: subroutine interface call eigh(a,lambda [,vectors] [,upper_a] [,overwrite_a] [,err] )

  • option to access lower/upper triangle part only
  • Avoid allocation when possible e.g. when providing vectors that can be used as temporary storage for a.

Note: the LAPACK backends are not pure due to functions with intent(out) arguments. The current proposed interfaces are ready to be made pure (e.g., function interface with all intent(in) arguments) to make it easy when the backend will be pure.

I believe this is ready for consideration, thank you @jvdp1 @jalvesz @fortran-lang/stdlib