Reading the code for eigvals!, it always tries to promote the matrix to some kind of “normal” float so it can work with LAPACK.geevx!… I think. But that clearly defeats the purpose of using something like DecFP.
julia> using DecFP, LinearAlgebra
julia> p = [d"4" d"7"; d"1" d"30"]
2×2 Matrix{Dec64}:
4.0 7.0
1.0 30.0
julia> eigvals!(p)
ERROR: MethodError: no method matching eigvals!(::Matrix{Dec64})
Why doesn’t Julia have a default function for eigvals that doesn’t rely on external libraries? Is there a package that solves this?
Computing eigenvalues inherently defeats the purpose of using a type like Dec64. Even if the input matrix consists of exact Dec64 decimals, the eigenvalues/vectors will generally not be decimal and will have roundoff errors. So you might as well convert the matrix to Float64 when calling eigvals.
Because eigenvalue computation is extremely nontrivial, so it’s not a matter of just writing a 60-line textbook-style fallback function like we can for LU factorization. That’s why eigenvalue computation for generic number types was left to packages like GenericLinearAlgebra.jl and GenericSchur.jl.
But in the case of Dec64, as I said above, you might as well convert to Float64 and use the high-performance LAPACK implementation. (There is more benefit to generic linear algebra if you have an extended-precision type like Dec128, but even then it is probably faster to go with a binary FP type like DoubleFloats than to use slow software-based decimal FP.)
Really, the only reason to use decimal FP is in situations where you need to exactly preserve human inputs and exact calculations thereon, e.g. for financial records. (Or for pedagogy.)
1 Like