How do I use eigen! with my input matrix?

Specific Domains Numerics
1

Hello

I have an idea, but i don’t know how to implemented it. I would like to pass a temporary variables to store the eigenvectors of eigen! - by default eigen! uses the original matrix.

I am confident that this is possible to do, because inside eigen function, a new matrix to store the eigenvectors must be created. I need some help to point out where this happen, then my naive thoughs would be to create a new function, and pass my matrix as argument.

The rest of algorithms does not need to change, therefore, less probability of issues.

Below I show my expected behavior, if, somehow, it helps.

using LinearAlgebra
A = rand(4,4)
λ₁, ψ₁ = eigen(A)
λ₂, ψ₂ = eigen!(A)

@show all(λ₁ .≈ λ₂)
@show all(ψ₁ .≈ ψ₂)

A_mine = zeros(A)
λ₃ = eigen!(A, A_mine) ### what I need
@show all(ψ₁ .≈ A_mine) ### of course should be true
2

you might as well just do this

julia> a= rand(3,3); b=similar(a);

julia> my_eigen(a,b) = b.=eigen(a).vectors

julia> my_eigen(a,b);

julia> x,y = eigen(a);

julia> all(b .≈ y)
true
1 Like
3

But the function eigen will allocate memory somewhere inside it, before we put the content inside b.
I would like to avoid the eigen to create a new matrix temporarily.

4

I mean it’s about the same where you do the copying:
https://github.com/JuliaLang/julia/blob/af3cd59e13c93fcb3e3788a79b0168127473d18a/stdlib/LinearAlgebra/src/eigen.jl#L235-L237

or

https://github.com/JuliaLang/julia/blob/af3cd59e13c93fcb3e3788a79b0168127473d18a/stdlib/LinearAlgebra/src/eigen.jl#L175

because what you want is a vector (for eigenvalue) and a matrix (for eigenvectors), which needs additional memory at some point

5

Hummmmm

I thought something more complex was happning behind the scenes, but AA is just a copy of A and then the eigen! function is called.

Looking here :

 eval, evec = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]] return Eigen(sorteig!(eval, evec, sortby)...)

Seems that evec comes from LAPACK directly, and I cannot offer an empty matrix as input. Am I right ?

1 Like