diff --git a/Project.toml b/Project.toml index c1b7be4..0338712 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "InfinitesimalGenerators" uuid = "2fce0c6f-5f0b-5c85-85c9-2ffe1d5ee30d" -version = "0.5.1" +version = "0.5.2" [deps] Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" diff --git a/src/principal_eigenvalue.jl b/src/principal_eigenvalue.jl index 2a52c9f..9bcb9ea 100644 --- a/src/principal_eigenvalue.jl +++ b/src/principal_eigenvalue.jl @@ -18,14 +18,19 @@ In other words, all eigenvalues of 𝕋 have real part <= 0. This means that """ function principal_eigenvalue(𝕋; r0 = ones(size(𝕋, 1))) a, η, r = 0.0, 0.0, r0 - try - # faster in certain cases. Check that all sum up to zero - @assert maximum(abs.(sum(𝕋, dims = 1))) < 1e-9 - # standard way of solving Ax = 0 is to do inverse iteration https://stackoverflow.com/questions/33563401/lapack-routines-for-solving-a-x-0 - vals, vecs = Arpack.eigs(𝕋; v0 = collect(r0), nev = 1, which = :LM, sigma = 0.0) - η = vals[1] - r = vecs[:, 1] - catch + if maximum(abs.(sum(𝕋, dims = 1))) < 1e-9 + # if columns sum up to zero + # we know principal is asssociated with zero + if 𝕋 isa Tridiagonal + η = 0.0 + r = [1.0 ; - Tridiagonal(𝕋.dl[2:end], 𝕋.d[2:end], 𝕋.du[2:end]) \ vec(𝕋[2:end, 1])] + else + # standard way of solving Ax = 0 is to do inverse iteration https://stackoverflow.com/questions/33563401/lapack-routines-for-solving-a-x-0 + vals, vecs = Arpack.eigs(𝕋; v0 = collect(r0), nev = 1, which = :LM, sigma = 0.0) + η = vals[1] + r = vecs[:, 1] + end + else a = - minimum(diag(𝕋)) try vals, vecs = Arpack.eigs(𝕋 + a * I; v0 = collect(r0), nev = 1, which = :LM)