# How can I give a function to eigs() in julia like matlab

**URL:** <https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699>\
**Category:** General Usage\
**Tags:** question, package, matlab, functions, arpack\
**Created:** [November 2, 2023, 1:37pm UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699 "2023-11-02T13:37:42Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Meteoooor](https://avatars.discourse-cdn.com/v4/letter/m/f6c823/32.png) [@Meteoooor](https://discourse.julialang.org/u/Meteoooor)\
**Post date:** [November 2, 2023, 1:37pm UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699/1 "2023-11-02T13:37:42Z")

</div>

I have the following matlab code:

```julia
function [DomRVec, DomEig, Error] = MyEigs( Tup, Tdn, NumofValue, SpeString, Dbond, Dspin, InitVec )

VecDim = Dbond^2; 
Tup = reshape( permute( Tup, [1, 3, 2] ), [Dbond * Dspin, Dbond] ); % A(lrm) => A(lm,r), placed at upper
Tdn = reshape( permute( Tdn, [3, 2, 1] ), [Dspin * Dbond, Dbond] ); % B(lrm) => B(mr,l), placed at lower
ActionFun = @( x ) Action( Tup, Tdn, x, Dbond, Dspin );
[DomRVec, DomEig] = eigs( ActionFun, VecDim, NumofValue, SpeString, 'IsFunctionSymmetric', 0, 'Tolerance', 1E-14, 'MaxIterations', 300, 'StartVector', InitVec, ...
    'Display', 0 );

Error = Maxdiff( Action( Tup, Tdn, DomRVec, Dbond, Dspin ), DomRVec * DomEig );

end

function [y] = Action( Tup, Tdn, x, Dbond, Dspin )

d = size( x, 2 );
x = reshape( x, [Dbond, Dbond * d] );

if( d == 1 )
    y = reshape( Tup * x, [Dbond, Dspin * Dbond] );
    y = reshape( y * Tdn, [Dbond^2, 1] );
elseif( d > 1) 
    y = reshape( Tup * x, [Dbond, Dspin * Dbond, d] );
    y = reshape( permute( y, [1, 3, 2] ), [Dbond * d, Dspin * Dbond] );
    y = reshape( y * Tdn, [Dbond, d, Dbond] );
    y = reshape( permute( y, [1, 3, 2] ), [Dbond * Dbond, d] );    
end

end

```

After referring to the documentation and source code of _Arpack.jl_, I wrote the following code. It runs without errors, but the results are significantly different from those obtained with MATLAB. Moreover, the results vary each time I run the code. I can’t figure out what’s causing this inconsistency. Could anyone point out where the mistake might be, or suggest alternative methods to achieve the same functionality?

```julia
using LinearAlgebra
using Arpack

struct MyMatrix{T,S} <: AbstractArray{T, 2}
    Tup::S
    Tdn::S
    Dbond::Int
    Dspin::Int
end

function MyMatrix(Tup, Tdn, Dbond, Dspin)
    Tnew = typeof(zero(eltype(Tup))/sqrt(one(eltype(Tup))))
    Tup_new = convert(AbstractMatrix{Tnew}, Tup)
    Tdn_new = convert(AbstractMatrix{Tnew}, Tdn)
    MyMatrix{Tnew,typeof(Tup_new)}(Tup_new, Tdn_new, Dbond, Dspin)
end

function LinearAlgebra.mul!(y::StridedVector{T}, A::MyMatrix{T}, x::StridedVector{T}) where T
    d = 1
    x = reshape(x, (A.Dbond, A.Dbond * d))
    if d == 1
        y = reshape(A.Tup * x, (A.Dbond, A.Dspin * A.Dbond))
        y = reshape(y * A.Tdn, (A.Dbond^2, 1))

    return y
end

Base.size(A::MyMatrix) = (A.Dbond^2, A.Dbond^2)
LinearAlgebra.ishermitian(A::MyMatrix) = true

Tup, Tdn, Dbond, Dspin = reshape(1:18,3,3,2),reshape(2:19,3,3,2),3,2
Tup = reshape( permutedims( Tup, [1, 3, 2] ), Dbond * Dspin, Dbond )
Tdn = reshape( permutedims( Tdn, [3, 2, 1] ), Dspin * Dbond, Dbond )
A = MyMatrix(Tup, Tdn, Dbond, Dspin)
λ, ϕ = eigs(A; nev=1,which=:LR)

```

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [November 2, 2023, 2:38pm UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699/2 "2023-11-02T14:38:50Z")

</div>

Did you try [LinearMaps](https://julialinearalgebra.github.io/LinearMaps.jl/stable)?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 2, 2023, 2:54pm UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699/3 "2023-11-02T14:54:09Z")

</div>

> [@Meteoooor](#):
>
> After referring to the documentation and source code of _Arpack.jl_

Note that nowadays there are pure-julia implementations of restarted Arnoldi, which may be more convenient than Arpack, such as [KrylovKit.jl](https://jutho.github.io/KrylovKit.jl/stable/man/eig/).

Note that if your linear operator is Hermitian, you should make sure the iterative solver knows this since better methods are available in this case. (In that case, you could also try [LOBPCG in IterativeSolvers.jl](https://iterativesolvers.julialinearalgebra.org/dev/eigenproblems/lobpcg/).)

---

<div class="post-metadata">

**Author:** ![Meteoooor](https://avatars.discourse-cdn.com/v4/letter/m/f6c823/32.png) [@Meteoooor](https://discourse.julialang.org/u/Meteoooor)\
**Post date:** [November 3, 2023, 3:15am UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699/4 "2023-11-03T03:15:16Z")

</div>

thank u! It’s solved by KrylovKit.jl, but there are still some minor issues, the keyword parameter _howany_ cannot accurately control the number of output eigenvalues/eigenvectors. Take the following code as an example, the number of the output eigenvalues vary each time though I set _howany==1_:

```julia
function Action(Tup, Tdn, x, Dbond, Dspin)
    d = size(x)[2]
    x = reshape(x, Dbond, Dbond * d)
    if d == 1
        y = reshape(Tup * x, Dbond, Dspin * Dbond)
        y = reshape(y * Tdn, Dbond^2,1)
    return y
end
Tup, Tdn, Dbond, Dspin = rand(10,10,2),rand(10,10,2),10,2
Tup = reshape( permutedims( Tup, [1, 3, 2] ), Dbond * Dspin, Dbond )
Tdn = reshape( permutedims( Tdn, [3, 2, 1] ), Dspin * Dbond, Dbond );
x0=ones(100,1);

ActionFun = x -> Action(Tup, Tdn, x, Dbond, Dspin);
DomEig, DomRVec = eigsolve(ActionFun,x0,1,:LR);

```

I don’t know why this is happening 🤔

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 3, 2023, 7:17am UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699/5 "2023-11-03T07:17:41Z")

</div>

from the docs (at least):

> Compute at least `howmany` eigenvalues from the linear map encoded in the matrix `A`

Basically it will stop when at least 1 ev has been found

---

<div class="post-metadata">

**Author:** ![Meteoooor](https://avatars.discourse-cdn.com/v4/letter/m/f6c823/32.png) [@Meteoooor](https://discourse.julialang.org/u/Meteoooor)\
**Post date:** [November 3, 2023, 7:55am UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699/6 "2023-11-03T07:55:30Z")

</div>

Alright, it looks like I asked a stupid question 🫢, thank u.

---

<div class="post-metadata">

**Author:** ![Meteoooor](https://avatars.discourse-cdn.com/v4/letter/m/f6c823/32.png) [@Meteoooor](https://discourse.julialang.org/u/Meteoooor)\
**Post date:** [November 3, 2023, 7:57am UTC](https://discourse.julialang.org/t/how-can-i-give-a-function-to-eigs-in-julia-like-matlab/105699/7 "2023-11-03T07:57:05Z")

</div>

No, but it’s solved by KrylovKit.jl, thanks for your suggestion