Closed perazz closed 6 days ago
From Fortran Monthly call:
real
array of singular values (it is complex
by default). But, stop on errors if the matrix had complex
eigenvalue pairs. For example, NumPy returns real
eigenvalues when all imaginary parts are zero:
if not isComplexType(t) and all(w.imag == 0.0):
w = w.real
vt = vt.real
result_t = _realType(result_t)
else:
result_t = _complexType(result_t)
Thank you @jvdp1 for the reviews. As pointed out during our last Monthly call with @everythingfunctional @jalvesz, I've added the option to return the real
part of the eigenvalues only. In this latest case, we must check that there are no meaningful imaginary components, and raise an error if so:
Thanks a lot @jvdp1 @jalvesz for the reviews - this PR looks well polished now. You will find references for each of your comments to the related commit.
LGTM @perazz, I think this is ready to go :)
Thank you both @jvdp1 @jalvesz, I would say let's wait a couple more days, and then merge if there are no further comments.
Thank you, I will merge this one.
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.xdp
subroutine
interfaceGeneral matrix
Prior art
eigenvalues, eigenvectors = eig(a)
eigenvalues = eigvals(a)
eig(a, b=None, left=False, right=True, overwrite_a=False, overwrite_b=False, check_finite=True, homogeneous_eigvals=False)
Proposed implementation
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 ascomplex
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
eigenvalues, eigenvectors = eigh(a, uplo='L')
eigenvalues = eigvalsh(a)
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
lambda = eigvalsh(A [, err])
Eigendecomposition: subroutine interface
call eigh(a,lambda [,vectors] [,upper_a] [,overwrite_a] [,err] )
vectors
that can be used as temporary storage fora
.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 allintent(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