GenericArpack.jl Documentation
This library is a julia translation of the Arpack library. I strongly recommend browsing the documentation on the functions in that library:
The major portion of the Julia code is a translation of these functions.
However, we also provide our own user-friendly interface to these codes. This is inspired by the eigs and svds interface in the Julia wrappers Arpack.jl, however it has a number of simplifications.
Eigenvalue, eigenvectors, standard, generalized.
The documentation needs to talk about matrices and eigenvectors, of course.
A standard eigenvalue and eigenvector is a pair:
- $Ax = \lambda x$ where $x$ is not-zero and $||x|| = 1$ by convention
If $A$is symmetic, then all of the $\lambda$ are real. Note that the vector $x$ is only determined up to sign, as well as rotation if there are multliple linearly independent eigenvectors with the same eigenvalue.
Likewise if $A$ is Hermitian, then all of the $\lambda$ are still real whereas the vectors are complex-valued. Also, the vectors $x$ are only unique up to scaling by complex-value of magnitude 1. (i.e. $e^{i\theta}$ for an angle $\theta$).
A generalized eigenvalue and eigenvector is a pair:
- $Ax = \lambda Bx$ where $x$ is not-zero and $||Bx|| = 1$ by convention.
The vectors are just as unique as in the other cases.
There are many complexities with eigenvalues and subtleties in the defintion. Please see the ARPACK manual for exact detail on what we mean here as this is a translation of ARPACK.
ARPACK Notes and Key Parameters
The three key parameters of the Arpack methods are:
kwhichncv
The goal of Arpack is to compute a partial eigenspace of a matrix.
The size of this eigenspace is the number of dimension k requested.
In order to decide which eigenspace we are targeting, we need to tell it! Arpack uses the natural variable: which in order to provide this information.
which: which part of the spectrum to compute. There are a few choices::LM(the default) find the largest magnitude eigenspace:SMfind the smallest magnitude eigenspace:LAfind the largest algebraic value eigenspace:SAfind the smallest algebraic value eigenspace:BEfind both smallest and largest algebraic eigenspace
It turns out that ARPACK uses a strategy based on looking for a larger eigenspace than only k vectors. The parameter ncv (number of compute vectors) determines this size. These parameters must satisfy $k < ncv \le =n$ (where $n$ is the dimension of the matrix).
If your eigenproblem isn't converging, it is often useful to increase ncv instead as well as increasing maxiter. This is because if the eigenspace you are trying to find is poorly separated, it's easier to find if you can find any separation of a larger eigenspace.
Examples of real symmetric, complex Hermitian, and singular value decomposition
The best way to understand the library is just to see a variety of examples.
There is more than one way to do something in GenericArpack.jl. This is by design and there are various limits to the interfaces. Also, there exist multiple synonymous calls. This is make it easier for someone reading your code to understand what you meant, and often because there is a small hint about what a natural default type would be. For instance: hermeigs causes us to use ComplexF64 as a default type, whereas symeigs will use Float64 instead. Of course, if you all symeigs and specify ComplexF64, it'll work. Likewise, if you call hermeigs and specify Float64, but on the other hand, your readers will probably be annoyed by that confusion, so use sparingly.
Find the largest eigenvalues of a real symmetric matrix.
Let A be any type in Julia that implements Base.size and LineareAlgebra.mul!(y,A,x). For instance, A can be a SparseMatrixCSC, a Matrix, a Tridiagonal, a Diagonal Moreover, we assume that A is symmetric (and we won't check). This can be a common error.
julia> n = 100100julia> A = Tridiagonal(-ones(n-1), 2*ones(n), -ones(n-1))100×100 Tridiagonal{Float64, Vector{Float64}}: 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋮ ⋱ ⋮ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0julia> symeigs(A, 6)GenericArpack.ArpackEigen{Float64, Float64, ArpackSimpleOp{Tridiagonal{Float64, Vector{Float64}}}, ArpackState{Float64}, :I} eigenspace: LM values: 6-element Vector{Float64}: 3.9652704964444903 3.9758608794815036 3.984539744726552 3.991298695938015 3.9961311942671656 3.999032564583961 vectors: 100×6 Matrix{Float64}: -0.0261102 0.0217972 -0.0174631 0.0131121 -0.00874848 0.00437636 0.0513136 -0.0430682 0.0346562 -0.0261102 0.0174631 -0.00874848 -0.0747349 0.0632996 -0.0513136 0.038881 -0.0261102 0.0131121 0.0955607 -0.082003 0.0671776 -0.0513136 0.0346562 -0.0174631 -0.113068 0.098727 -0.082003 0.0632996 -0.0430682 0.0217972 0.126648 -0.113068 0.0955607 -0.0747349 0.0513136 -0.0261102 -0.13583 0.124679 -0.107641 0.0855198 -0.0593604 0.0303979 0.140294 -0.133281 0.118057 -0.0955607 0.0671776 -0.0346562 -0.139886 0.138665 -0.126648 0.10477 -0.0747349 0.038881 0.13462 -0.140702 0.133281 -0.113068 0.082003 -0.0430682 ⋮ ⋮ -0.139886 -0.138665 -0.126648 -0.10477 -0.0747349 -0.038881 0.140294 0.133281 0.118057 0.0955607 0.0671776 0.0346562 -0.13583 -0.124679 -0.107641 -0.0855198 -0.0593604 -0.0303979 0.126648 0.113068 0.0955607 0.0747349 0.0513136 0.0261102 -0.113068 -0.098727 -0.082003 -0.0632996 -0.0430682 -0.0217972 0.0955607 0.082003 0.0671776 0.0513136 0.0346562 0.0174631 -0.0747349 -0.0632996 -0.0513136 -0.038881 -0.0261102 -0.0131121 0.0513136 0.0430682 0.0346562 0.0261102 0.0174631 0.00874848 -0.0261102 -0.0217972 -0.0174631 -0.0131121 -0.00874848 -0.00437636
We can also use any of the following equivalent calls
julia> eigs(Symmetric(A), 6)GenericArpack.ArpackEigen{Float64, Float64, ArpackSimpleOp{Symmetric{Float64, Tridiagonal{Float64, Vector{Float64}}}}, ArpackState{Float64}, :I} eigenspace: LM values: 6-element Vector{Float64}: 3.9652704964444903 3.9758608794815036 3.984539744726552 3.991298695938015 3.9961311942671656 3.999032564583961 vectors: 100×6 Matrix{Float64}: -0.0261102 0.0217972 -0.0174631 0.0131121 -0.00874848 0.00437636 0.0513136 -0.0430682 0.0346562 -0.0261102 0.0174631 -0.00874848 -0.0747349 0.0632996 -0.0513136 0.038881 -0.0261102 0.0131121 0.0955607 -0.082003 0.0671776 -0.0513136 0.0346562 -0.0174631 -0.113068 0.098727 -0.082003 0.0632996 -0.0430682 0.0217972 0.126648 -0.113068 0.0955607 -0.0747349 0.0513136 -0.0261102 -0.13583 0.124679 -0.107641 0.0855198 -0.0593604 0.0303979 0.140294 -0.133281 0.118057 -0.0955607 0.0671776 -0.0346562 -0.139886 0.138665 -0.126648 0.10477 -0.0747349 0.038881 0.13462 -0.140702 0.133281 -0.113068 0.082003 -0.0430682 ⋮ ⋮ -0.139886 -0.138665 -0.126648 -0.10477 -0.0747349 -0.038881 0.140294 0.133281 0.118057 0.0955607 0.0671776 0.0346562 -0.13583 -0.124679 -0.107641 -0.0855198 -0.0593604 -0.0303979 0.126648 0.113068 0.0955607 0.0747349 0.0513136 0.0261102 -0.113068 -0.098727 -0.082003 -0.0632996 -0.0430682 -0.0217972 0.0955607 0.082003 0.0671776 0.0513136 0.0346562 0.0174631 -0.0747349 -0.0632996 -0.0513136 -0.038881 -0.0261102 -0.0131121 0.0513136 0.0430682 0.0346562 0.0261102 0.0174631 0.00874848 -0.0261102 -0.0217972 -0.0174631 -0.0131121 -0.00874848 -0.00437636julia> hermeigs(Symmetric(A), 6)GenericArpack.ArpackEigen{Float64, ComplexF64, ArpackSimpleOp{Symmetric{Float64, Tridiagonal{Float64, Vector{Float64}}}}, ArpackState{Float64}, :I} eigenspace: LM values: 6-element Vector{Float64}: 3.965270496444535 3.9758608794815165 3.9845397447265563 3.991298695938034 3.9961311942671838 3.999032564583959 vectors: 100×6 Matrix{ComplexF64}: 0.000180037-0.0261096im … 0.0040213-0.00172675im -0.000353822+0.0513124im -0.0080387+0.00345184im 0.000515318-0.0747331im 0.0120483-0.00517358im -0.000658918+0.0955584im -0.0160463+0.00689032im 0.000779634-0.113065im 0.0200288-0.00860039im -0.000873273+0.126645im … -0.0239918+0.0103021im 0.000936585-0.135827im 0.0279317-0.0119939im -0.000967369+0.140291im -0.0318445+0.0136741im 0.000964557-0.139883im 0.0357266-0.0153411im -0.000928246+0.134617im -0.039574+0.0169932im ⋮ ⋱ ⋮ 0.000964557-0.139883im -0.0357266+0.0153411im -0.000967369+0.140291im 0.0318445-0.0136741im 0.000936585-0.135827im -0.0279317+0.0119939im -0.000873273+0.126645im 0.0239918-0.0103021im 0.000779634-0.113065im … -0.0200288+0.00860039im -0.000658918+0.0955584im 0.0160463-0.00689032im 0.000515318-0.0747331im -0.0120483+0.00517358im -0.000353822+0.0513124im 0.0080387-0.00345184im 0.000180037-0.0261096im -0.0040213+0.00172675im
Finding the largest eigenvalues with an operator
All GenericArpack calls are mapped to a single computational interface with an ArpackOp type.
In the previous block, we can't use eigs and must use symeigs. That's because right now, we only have the symmetric eigensolvers ported. However, since the op type has no way to tell eigs that it should pick the symmetric or non-symmetric version, we need to tell it.
So if you see
MethodError: no method matching eigs(::ArpackSimpleOp{ ... })That means you just need to call symeigs instead! This will probaby get fixed at some point.
Complex Hermitian
We also extend Arpack's symmetric solvers with the ability to solve Hermitian problems.
julia> n = 100100julia> A = Tridiagonal(-ones(n-1)*1im, 2*ones(n).+0im, ones(n-1)*1im)100×100 Tridiagonal{ComplexF64, Vector{ComplexF64}}: 2.0+0.0im 0.0+1.0im ⋅ … ⋅ ⋅ ⋅ -0.0-1.0im 2.0+0.0im 0.0+1.0im ⋅ ⋅ ⋅ ⋅ -0.0-1.0im 2.0+0.0im ⋅ ⋅ ⋅ ⋅ ⋅ -0.0-1.0im ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0+1.0im ⋅ ⋅ ⋅ ⋅ ⋅ 2.0+0.0im 0.0+1.0im ⋅ ⋅ ⋅ ⋅ -0.0-1.0im 2.0+0.0im 0.0+1.0im ⋅ ⋅ ⋅ ⋅ -0.0-1.0im 2.0+0.0imjulia> hermeigs(A, 6)GenericArpack.ArpackEigen{Float64, ComplexF64, ArpackSimpleOp{Tridiagonal{ComplexF64, Vector{ComplexF64}}}, ArpackState{Float64}, :I} eigenspace: LM values: 6-element Vector{Float64}: 3.965270496444513 3.9758608794815284 3.9845397447265345 3.9912986959380317 3.996131194267193 3.9990325645839797 vectors: 100×6 Matrix{ComplexF64}: 0.0255806+0.00523195im … -0.0043763+2.21686e-5im 0.0102822-0.0502729im 4.43158e-5+0.00874837im -0.0732191-0.0149754im 0.013112-6.64201e-5im -0.0191484+0.0936225im -8.84602e-5-0.0174629im 0.110774+0.0226565im -0.0217969+0.000110415im 0.0253777-0.124079im … 0.000132262+0.0261099im -0.133075-0.0272175im 0.0303975-0.000153982im -0.0281121+0.137449im -0.000175553-0.0346558im 0.137049+0.0280304im -0.0388805+0.000196954im 0.0269752-0.13189im 0.000218164+0.0430677im ⋮ ⋱ ⋮ 0.0280304-0.137049im -0.000196954-0.0388805im -0.137449-0.0281121im -0.0346558+0.000175553im -0.0272175+0.133075im 0.000153982+0.0303975im 0.124079+0.0253777im 0.0261099-0.000132262im 0.0226565-0.110774im … -0.000110415-0.0217969im -0.0936225-0.0191484im -0.0174629+8.84602e-5im -0.0149754+0.0732191im 6.64201e-5+0.013112im 0.0502729+0.0102822im 0.00874837-4.43158e-5im 0.00523195-0.0255806im -2.21686e-5-0.0043763im
SVD (Singular Value Decomposition) via the Normal Equations
Gene Golub often railed against the use of the normal equations for SVD computations. However, he would also use them when appropriate. For computing singular values via ARPACK, the normal equations offer a few advantages.
- smallest singular values
- they reduce computation as the dimension of the vectors is smaller
The downside is that they also reduce maximum obtainable accuracy. So please consider using the high-precision types we enable in GenericArpack.jl, if accuracy is desired.
(I need a better SVD example here, but alas.)
julia> using GenericArpack, SparseArrays, StableRNGsjulia> A = sprand(StableRNG(1), 200, 100, 5/200)200×100 SparseMatrixCSC{Float64, Int64} with 472 stored entries: ⠀⠠⠅⠀⠀⠀⠀⠰⠤⠐⠀⠠⠠⡀⠠⢀⠉⢰⠀⠀ ⠈⡈⠀⡅⡁⠠⠃⠅⢀⠐⠈⠠⡁⠂⠃⠠⠀⠈⠀⠀ ⢄⠦⠈⠠⠁⠀⠈⠊⠈⠀⢌⠁⠀⠁⠀⠁⠁⠀⠉⠀ ⢂⡈⡂⠐⠠⠀⠀⠂⠊⠁⢀⠉⢀⠠⠂⠀⠈⠀⠠⠀ ⠐⠲⠠⠀⠀⠂⡄⠀⠠⣉⠐⢀⠡⠀⠂⠁⠀⠀⠌⠁ ⢁⡐⢄⠓⠢⢨⡐⠐⠂⠀⡀⠂⠀⠂⡀⠂⠂⠀⢀⠁ ⡀⠢⠤⠀⢈⠈⠀⠁⢀⠀⠀⠒⠀⣀⠂⠈⣀⠀⠀⢠ ⣐⣀⡀⡐⠒⠀⢈⠠⠈⠊⠀⠀⠉⠀⢐⠄⣐⠀⠐⢁ ⠀⠃⢂⡂⠘⠀⠀⢈⠂⠀⡄⠁⡂⠂⠐⠀⠁⡒⠈⠂ ⠈⡀⠀⠀⠁⠀⠐⠀⠐⠀⠀⢀⠀⠆⠐⠁⠀⢀⠐⠀ ⠀⢐⠁⢀⠠⡂⢀⠂⠈⠤⠁⠈⠀⡅⠁⠂⠔⠐⠂⠌ ⠈⠀⠀⡀⢀⠀⠈⠈⠀⠂⠀⠂⠄⠊⠆⠌⡄⠫⠤⡜ ⠁⠠⠔⢠⠀⠀⠀⠰⠀⠀⡀⠁⠆⠁⠀⠁⠁⡄⢀⠈ ⠠⡆⠀⡄⠄⠍⠀⠀⢐⠠⠀⠀⢀⠠⠍⠀⠠⠀⡀⠄ ⠀⠀⠠⢀⢀⠄⢉⢤⠀⡤⠀⠩⠀⠀⠉⡁⠄⠠⠀⠀ ⠠⠀⠀⠀⠤⠈⢀⠀⡀⠣⠀⡈⡄⠤⠊⠈⡀⠀⠠⢀ ⠀⡋⠀⠃⠂⡈⠘⡈⡐⠀⠓⡀⠀⢓⠐⡀⡓⠠⡎⠉ ⡘⢂⠀⠖⠁⠂⡘⠉⠀⠐⡀⡂⠃⡠⠀⡄⡀⣀⡀⠠ ⠀⡀⢀⠀⡀⢑⠐⠀⠀⢀⠀⠀⡠⠀⠀⡂⠈⡐⠑⢀ ⠈⠀⠀⠀⢀⠂⠠⢁⡈⠀⠁⠂⠈⢒⠠⠀⠀⠰⡀⠀julia> svds(A, 3)SVD{Float64, Float64, Matrix{Float64}} U factor: 200×3 Matrix{Float64}: -0.00366776 -0.0067998 0.00625225 -0.0368311 -0.0160961 -0.0131541 -0.0229894 0.00929866 -0.0209747 -0.00117525 -0.00190433 -0.00369561 -0.0305999 -0.011878 -0.0595654 -0.0165446 0.00449102 -0.017125 0.0 0.0 0.0 0.0 0.0 0.0 -0.0590323 -0.087415 -0.148514 -0.0694489 -0.111795 0.0807277 ⋮ -0.0346454 -0.0150924 0.0240361 -0.0708986 0.0193406 -0.0635748 -0.0364264 0.0199232 -0.0125269 -0.0501192 -0.066657 -0.0880075 -0.0666061 -0.109053 -0.0877816 -0.0064085 -0.00913391 -0.0220183 -0.00253531 -0.00422375 -0.00814783 -0.0714716 0.0137463 -0.109677 -0.0415146 -0.0741379 -0.153432 singular values: 3-element Vector{Float64}: 2.6944729446552693 2.4968736254169444 2.3785967629735754 Vt factor: 3×100 Matrix{Float64}: -0.0178603 -0.0281502 -0.0186638 … -0.0158545 -0.18666 -0.0380778 -0.0189334 -0.0177954 -0.0111649 0.0137575 0.177704 -0.0229597 -0.0228844 -0.02005 0.0240699 0.00480795 0.101246 0.0060408
We can also find the smallest subspace
julia> svds(A, 3; which=:SM) # smallestSVD{Float64, Float64, Matrix{Float64}} U factor: 200×3 Matrix{Float64}: 0.582817 -0.0134732 0.00607081 -0.155068 -0.0256579 -0.00537595 -0.00418642 -0.062982 0.0739621 -0.0107715 -0.18135 -0.0441208 0.0246127 0.0455897 -0.0964834 -0.0614142 0.192247 0.369308 0.0 0.0 0.0 0.0 0.0 0.0 0.0242838 0.11556 -0.0388435 -0.0164395 0.0235335 0.00997486 ⋮ -0.0288299 -0.0486256 -0.00814136 0.0146921 -0.0276828 -0.0537498 0.039303 0.0911307 -0.0604381 0.0549959 0.142183 -0.0632323 -0.0117026 -0.0344057 0.0675823 -0.0227149 0.177344 0.00728057 0.00362815 0.0593515 -0.000656137 0.0048518 0.108275 -0.000904173 -0.0124722 -0.0741284 0.057463 singular values: 3-element Vector{Float64}: 0.10806971046978665 0.21758284541859946 0.2303303405180739 Vt factor: 3×100 Matrix{Float64}: -0.0194694 0.00119836 0.00181863 … 0.0184082 0.000268922 -0.0689222 -0.0466259 0.00542729 -0.00773754 -0.00103489 -0.00309843 0.0261328 0.0399371 -0.00529339 -0.0470376
Or only the singular values (although this doesn't really save as much work as it might sound like) and use them to compute the matrix 2-norm.
julia> vals = svds(A, 4; which=:BE, ritzvec=false, ncv=12).S # both smallest and largest4-element Vector{Float64}: 0.10806971046977838 0.21758284541859638 2.496873625416946 2.6944729446552684julia> [vals[end]/vals[1] cond(Matrix(A))]1×2 Matrix{Float64}: 24.9327 24.9327
So we can estimate the matrix 2-norm with some reasonable accuracy (in a well-conditioned case).
Complex SVD (Singular Value Decomposition)
This operation is not in the pure ARPACK library.
julia> using GenericArpack, SparseArrays, StableRNGsjulia> A = sprand(StableRNG(1), ComplexF64, 200, 100, 5/200)200×100 SparseMatrixCSC{ComplexF64, Int64} with 472 stored entries: ⠀⠠⠅⠀⠀⠀⠀⠰⠤⠐⠀⠠⠠⡀⠠⢀⠉⢰⠀⠀ ⠈⡈⠀⡅⡁⠠⠃⠅⢀⠐⠈⠠⡁⠂⠃⠠⠀⠈⠀⠀ ⢄⠦⠈⠠⠁⠀⠈⠊⠈⠀⢌⠁⠀⠁⠀⠁⠁⠀⠉⠀ ⢂⡈⡂⠐⠠⠀⠀⠂⠊⠁⢀⠉⢀⠠⠂⠀⠈⠀⠠⠀ ⠐⠲⠠⠀⠀⠂⡄⠀⠠⣉⠐⢀⠡⠀⠂⠁⠀⠀⠌⠁ ⢁⡐⢄⠓⠢⢨⡐⠐⠂⠀⡀⠂⠀⠂⡀⠂⠂⠀⢀⠁ ⡀⠢⠤⠀⢈⠈⠀⠁⢀⠀⠀⠒⠀⣀⠂⠈⣀⠀⠀⢠ ⣐⣀⡀⡐⠒⠀⢈⠠⠈⠊⠀⠀⠉⠀⢐⠄⣐⠀⠐⢁ ⠀⠃⢂⡂⠘⠀⠀⢈⠂⠀⡄⠁⡂⠂⠐⠀⠁⡒⠈⠂ ⠈⡀⠀⠀⠁⠀⠐⠀⠐⠀⠀⢀⠀⠆⠐⠁⠀⢀⠐⠀ ⠀⢐⠁⢀⠠⡂⢀⠂⠈⠤⠁⠈⠀⡅⠁⠂⠔⠐⠂⠌ ⠈⠀⠀⡀⢀⠀⠈⠈⠀⠂⠀⠂⠄⠊⠆⠌⡄⠫⠤⡜ ⠁⠠⠔⢠⠀⠀⠀⠰⠀⠀⡀⠁⠆⠁⠀⠁⠁⡄⢀⠈ ⠠⡆⠀⡄⠄⠍⠀⠀⢐⠠⠀⠀⢀⠠⠍⠀⠠⠀⡀⠄ ⠀⠀⠠⢀⢀⠄⢉⢤⠀⡤⠀⠩⠀⠀⠉⡁⠄⠠⠀⠀ ⠠⠀⠀⠀⠤⠈⢀⠀⡀⠣⠀⡈⡄⠤⠊⠈⡀⠀⠠⢀ ⠀⡋⠀⠃⠂⡈⠘⡈⡐⠀⠓⡀⠀⢓⠐⡀⡓⠠⡎⠉ ⡘⢂⠀⠖⠁⠂⡘⠉⠀⠐⡀⡂⠃⡠⠀⡄⡀⣀⡀⠠ ⠀⡀⢀⠀⡀⢑⠐⠀⠀⢀⠀⠀⡠⠀⠀⡂⠈⡐⠑⢀ ⠈⠀⠀⠀⢀⠂⠠⢁⡈⠀⠁⠂⠈⢒⠠⠀⠀⠰⡀⠀julia> svds(A, 3)SVD{ComplexF64, Float64, Matrix{ComplexF64}} U factor: 200×3 Matrix{ComplexF64}: 0.000622283+0.000650681im … 0.00051443+0.000543965im -0.00927321+0.117742im 0.0113616+0.0428172im -0.0195419+0.0418711im 0.0269183-0.0133982im -0.000861276+0.00347481im -0.00645723-0.00473884im -0.0195111+0.0563448im -0.0461299-0.0302614im 0.00417568+0.0334827im … -0.0744713-0.0309124im 0.0+0.0im -0.0+0.0im 0.0+0.0im -0.0+0.0im -0.00268413+0.0210515im -0.024052+0.000259037im -0.00524457+0.0504315im -0.0472448+0.00713175im ⋮ ⋱ -0.00894584+0.0287871im -0.0400742-0.0365695im -0.0186+0.0510772im -0.0363485-0.0364106im 5.17672e-5+0.00578616im 0.0041762+0.00170046im -0.0161525+0.096482im -0.057708+0.0398135im -0.00852922+0.0182682im … -0.0252699-0.0242831im 0.00268336+0.00765789im -0.00310276+0.00644434im 0.0138138+0.03483im -0.0202434+0.00419133im 0.00237519+0.0176187im -0.0514271-0.023305im -0.0418854+0.0458577im -0.0314791-0.12014im singular values: 3-element Vector{Float64}: 3.578494653659692 3.287380399538958 3.1650809676802596 Vt factor: 3×100 Matrix{ComplexF64}: 0.0085072-0.0152773im 0.0227414-0.0697884im … 0.0319019-0.0196392im -0.0128495-0.0195603im 0.0120002+0.0599222im 0.00884387+0.00317647im -0.0153662+0.0154273im -0.000249934-0.0193798im -0.0534621-0.0252439im
Real and Complex SVD via an ArpackNormalOp
All of the complex SVD calls for any Arpack matrix type are mapped to a single call with an ArpackNormalOp
julia> svds(ArpackNormalOp(Float64, A), 3; which=:SM) # smallestSVD{Float64, Float64, Matrix{Float64}} U factor: 200×3 Matrix{Float64}: 0.582817 -0.0134732 0.00607081 -0.155068 -0.0256579 -0.00537595 -0.00418642 -0.062982 0.0739621 -0.0107715 -0.18135 -0.0441208 0.0246127 0.0455897 -0.0964834 -0.0614142 0.192247 0.369308 0.0 0.0 0.0 0.0 0.0 0.0 0.0242838 0.11556 -0.0388435 -0.0164395 0.0235335 0.00997486 ⋮ -0.0288299 -0.0486256 -0.00814136 0.0146921 -0.0276828 -0.0537498 0.039303 0.0911307 -0.0604381 0.0549959 0.142183 -0.0632323 -0.0117026 -0.0344057 0.0675823 -0.0227149 0.177344 0.00728057 0.00362815 0.0593515 -0.000656137 0.0048518 0.108275 -0.000904173 -0.0124722 -0.0741284 0.057463 singular values: 3-element Vector{Float64}: 0.10806971046978665 0.21758284541859946 0.2303303405180739 Vt factor: 3×100 Matrix{Float64}: -0.0194694 0.00119836 0.00181863 … 0.0184082 0.000268922 -0.0689222 -0.0466259 0.00542729 -0.00773754 -0.00103489 -0.00309843 0.0261328 0.0399371 -0.00529339 -0.0470376
Examples of Generalized Eigenvalue Problems
Generalized Eigenvalue Example (from Arpack dsdrv3 sample).
julia> n = 100100julia> A = Tridiagonal(-ones(n-1),2*ones(n),-ones(n-1)).*(n+1)100×100 Tridiagonal{Float64, Vector{Float64}}: 202.0 -101.0 ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ -101.0 202.0 -101.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -101.0 202.0 -101.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -101.0 202.0 -101.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -101.0 202.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -101.0 … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … -101.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 202.0 -101.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -101.0 202.0 -101.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -101.0 202.0 -101.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -101.0 202.0julia> B = SymTridiagonal(4*ones(n),ones(n-1)).*(1/(6*(n+1)))100×100 Matrix{Float64}: 0.00660066 0.00165017 0.0 … 0.0 0.0 0.0 0.00165017 0.00660066 0.00165017 0.0 0.0 0.0 0.0 0.00165017 0.00660066 0.0 0.0 0.0 0.0 0.0 0.00165017 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋮ ⋱ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.00165017 0.0 0.0 0.0 0.0 0.0 0.00660066 0.00165017 0.0 0.0 0.0 0.0 0.00165017 0.00660066 0.00165017 0.0 0.0 0.0 0.0 0.00165017 0.00660066julia> vals, vecs = eigs(Symmetric(A), Symmetric(B), 4)GenericArpack.ArpackEigen{Float64, Float64, ArpackSymmetricGeneralizedOp{Symmetric{Float64, Tridiagonal{Float64, Vector{Float64}}}, BunchKaufman{Float64, Matrix{Float64}}, Symmetric{Float64, Matrix{Float64}}}, ArpackState{Float64}, :G} eigenspace: LM values: 4-element Vector{Float64}: 121003.49732902284 121616.60247324086 122057.49457079486 122323.22366457571 vectors: 100×4 Matrix{Float64}: -0.30281 -0.227747 -0.152137 0.0761604 0.600939 0.453512 0.303685 -0.152247 -0.889778 -0.675331 -0.454058 0.228186 1.16486 0.891273 0.602675 -0.303905 -1.42193 -1.09946 -0.74896 0.37933 1.65702 1.29808 0.892348 -0.454387 -1.86649 -1.48541 -1.03228 0.529006 2.04711 1.65981 1.16822 -0.603112 -2.19608 -1.81977 -1.29965 0.676635 2.31109 1.96389 1.42604 -0.749503 ⋮ -2.19608 1.81977 -1.29965 -0.676635 2.04711 -1.65981 1.16822 0.603112 -1.86649 1.48541 -1.03228 -0.529006 1.65702 -1.29808 0.892348 0.454387 -1.42193 1.09946 -0.74896 -0.37933 1.16486 -0.891273 0.602675 0.303905 -0.889778 0.675331 -0.454058 -0.228186 0.600939 -0.453512 0.303685 0.152247 -0.30281 0.227747 -0.152137 -0.0761604
Complex Hermitian Generalized Eigenvalue Example
Not here yet! Send me one if you have one!
Shift-Invert Example
Not here yet! Send me one if you have one! (Need to test and debug the Shift Invert mode driver.)
Buckling Mode Example
Not here yet! Send me one if you have one! (Need to test and debug the Buckling mode driver.)
Examples with high-precision (or low-precision) types
Both svds and symeigs/hermeigs take in a set of types.
With just one type
Consider our initial tridiagonal matrix $A$. By default, we use a float-type based on the element-type of the input matrix $A$. But when the input includes a specific type, we use that:
julia> eigs(Float32, Symmetric(A), 6) # use Float32 for all computationsGenericArpack.ArpackEigen{Float32, Float32, ArpackSimpleOp{Symmetric{Float64, Tridiagonal{Float64, Vector{Float64}}}}, ArpackState{Float32}, :I} eigenspace: LM values: 6-element Vector{Float32}: 3.9652638 3.975855 3.9845366 3.9912934 3.9961293 3.9990337 vectors: 100×6 Matrix{Float32}: -0.0261102 0.0217993 -0.0174604 0.0131113 -0.00875244 0.00437105 0.0513135 -0.0430726 0.0346508 -0.0261086 0.017471 -0.00873789 -0.0747347 0.0633062 -0.0513053 0.0388788 -0.0261219 0.0130963 0.0955601 -0.0820118 0.0671664 -0.0513108 0.0346717 -0.0174421 -0.113067 0.098738 -0.0819889 0.0632965 -0.0430874 0.0217712 0.126647 -0.113081 0.0955433 -0.0747318 0.0513363 -0.0260794 -0.135828 0.124695 -0.10762 0.0855172 -0.0593866 0.0303624 0.140292 -0.133299 0.118033 -0.095559 0.0672069 -0.0346163 -0.139884 0.138685 -0.12662 0.10477 -0.0747671 0.0388367 0.134619 -0.140723 0.13325 -0.113069 0.082038 -0.0430197 ⋮ ⋮ -0.13987 -0.138678 -0.12671 -0.104741 -0.0746934 -0.0389015 0.140281 0.133294 0.118116 0.0955331 0.0671401 0.0346748 -0.135819 -0.124692 -0.107695 -0.0854945 -0.0593271 -0.0304143 0.12664 0.113079 0.0956095 0.0747122 0.0512845 0.0261243 -0.113061 -0.0987375 -0.0820456 -0.06328 -0.0430434 -0.0218088 0.0955561 0.0820121 0.0672129 0.0512973 0.034636 0.0174724 -0.0747315 -0.0633071 -0.0513407 -0.0388684 -0.0260947 -0.0131192 0.0513114 0.0430735 0.0346747 0.0261015 0.0174527 0.00875327 -0.0261091 -0.0217998 -0.0174724 -0.0131077 -0.00874317 -0.00437874
If A has a different float type, then we may use that.
julia> Af32 = Float32.(A) # make a Float32 copy100×100 Tridiagonal{Float32, Vector{Float32}}: 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋮ ⋱ ⋮ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0 -1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -1.0 2.0julia> symeigs(Float32.(A), 6) # will switch default to Float32.GenericArpack.ArpackEigen{Float64, Float64, ArpackSimpleOp{Tridiagonal{Float32, Vector{Float32}}}, ArpackState{Float64}, :I} eigenspace: LM values: 6-element Vector{Float64}: 3.9652704964444903 3.9758608794815036 3.984539744726552 3.991298695938015 3.9961311942671656 3.999032564583961 vectors: 100×6 Matrix{Float64}: -0.0261102 0.0217972 -0.0174631 0.0131121 -0.00874848 0.00437636 0.0513136 -0.0430682 0.0346562 -0.0261102 0.0174631 -0.00874848 -0.0747349 0.0632996 -0.0513136 0.038881 -0.0261102 0.0131121 0.0955607 -0.082003 0.0671776 -0.0513136 0.0346562 -0.0174631 -0.113068 0.098727 -0.082003 0.0632996 -0.0430682 0.0217972 0.126648 -0.113068 0.0955607 -0.0747349 0.0513136 -0.0261102 -0.13583 0.124679 -0.107641 0.0855198 -0.0593604 0.0303979 0.140294 -0.133281 0.118057 -0.0955607 0.0671776 -0.0346562 -0.139886 0.138665 -0.126648 0.10477 -0.0747349 0.038881 0.13462 -0.140702 0.133281 -0.113068 0.082003 -0.0430682 ⋮ ⋮ -0.139886 -0.138665 -0.126648 -0.10477 -0.0747349 -0.038881 0.140294 0.133281 0.118057 0.0955607 0.0671776 0.0346562 -0.13583 -0.124679 -0.107641 -0.0855198 -0.0593604 -0.0303979 0.126648 0.113068 0.0955607 0.0747349 0.0513136 0.0261102 -0.113068 -0.098727 -0.082003 -0.0632996 -0.0430682 -0.0217972 0.0955607 0.082003 0.0671776 0.0513136 0.0346562 0.0174631 -0.0747349 -0.0632996 -0.0513136 -0.038881 -0.0261102 -0.0131121 0.0513136 0.0430682 0.0346562 0.0261102 0.0174631 0.00874848 -0.0261102 -0.0217972 -0.0174631 -0.0131121 -0.00874848 -0.00437636
But if you give an op, that has no natural type, so we default to Float64
julia> op = ArpackSimpleOp(Float32.(A))ArpackSimpleOp{Tridiagonal{Float32, Vector{Float32}}}(Float32[2.0 -1.0 … 0.0 0.0; -1.0 2.0 … 0.0 0.0; … ; 0.0 0.0 … 2.0 -1.0; 0.0 0.0 … -1.0 2.0])julia> symeigs(op, 6)GenericArpack.ArpackEigen{Float64, Float64, ArpackSimpleOp{Tridiagonal{Float32, Vector{Float32}}}, ArpackState{Float64}, :I} eigenspace: LM values: 6-element Vector{Float64}: 3.9652704964444903 3.9758608794815036 3.984539744726552 3.991298695938015 3.9961311942671656 3.999032564583961 vectors: 100×6 Matrix{Float64}: -0.0261102 0.0217972 -0.0174631 0.0131121 -0.00874848 0.00437636 0.0513136 -0.0430682 0.0346562 -0.0261102 0.0174631 -0.00874848 -0.0747349 0.0632996 -0.0513136 0.038881 -0.0261102 0.0131121 0.0955607 -0.082003 0.0671776 -0.0513136 0.0346562 -0.0174631 -0.113068 0.098727 -0.082003 0.0632996 -0.0430682 0.0217972 0.126648 -0.113068 0.0955607 -0.0747349 0.0513136 -0.0261102 -0.13583 0.124679 -0.107641 0.0855198 -0.0593604 0.0303979 0.140294 -0.133281 0.118057 -0.0955607 0.0671776 -0.0346562 -0.139886 0.138665 -0.126648 0.10477 -0.0747349 0.038881 0.13462 -0.140702 0.133281 -0.113068 0.082003 -0.0430682 ⋮ ⋮ -0.139886 -0.138665 -0.126648 -0.10477 -0.0747349 -0.038881 0.140294 0.133281 0.118057 0.0955607 0.0671776 0.0346562 -0.13583 -0.124679 -0.107641 -0.0855198 -0.0593604 -0.0303979 0.126648 0.113068 0.0955607 0.0747349 0.0513136 0.0261102 -0.113068 -0.098727 -0.082003 -0.0632996 -0.0430682 -0.0217972 0.0955607 0.082003 0.0671776 0.0513136 0.0346562 0.0174631 -0.0747349 -0.0632996 -0.0513136 -0.038881 -0.0261102 -0.0131121 0.0513136 0.0430682 0.0346562 0.0261102 0.0174631 0.00874848 -0.0261102 -0.0217972 -0.0174631 -0.0131121 -0.00874848 -0.00437636
(This is because an ArpackOp may wrap a variety of information and we don't need the type information.)
julia> op = ArpackSimpleOp(Float32.(A))ArpackSimpleOp{Tridiagonal{Float32, Vector{Float32}}}(Float32[2.0 -1.0 … 0.0 0.0; -1.0 2.0 … 0.0 0.0; … ; 0.0 0.0 … 2.0 -1.0; 0.0 0.0 … -1.0 2.0])julia> symeigs(Float32, op, 6)GenericArpack.ArpackEigen{Float32, Float32, ArpackSimpleOp{Tridiagonal{Float32, Vector{Float32}}}, ArpackState{Float32}, :I} eigenspace: LM values: 6-element Vector{Float32}: 3.9652638 3.975855 3.9845366 3.9912934 3.9961293 3.9990337 vectors: 100×6 Matrix{Float32}: -0.0261102 0.0217993 -0.0174604 0.0131113 -0.00875244 0.00437105 0.0513135 -0.0430726 0.0346508 -0.0261086 0.017471 -0.00873789 -0.0747347 0.0633062 -0.0513053 0.0388788 -0.0261219 0.0130963 0.0955601 -0.0820118 0.0671664 -0.0513108 0.0346717 -0.0174421 -0.113067 0.098738 -0.0819889 0.0632965 -0.0430874 0.0217712 0.126647 -0.113081 0.0955433 -0.0747318 0.0513363 -0.0260794 -0.135828 0.124695 -0.10762 0.0855172 -0.0593866 0.0303624 0.140292 -0.133299 0.118033 -0.095559 0.0672069 -0.0346163 -0.139884 0.138685 -0.12662 0.10477 -0.0747671 0.0388367 0.134619 -0.140723 0.13325 -0.113069 0.082038 -0.0430197 ⋮ ⋮ -0.13987 -0.138678 -0.12671 -0.104741 -0.0746934 -0.0389015 0.140281 0.133294 0.118116 0.0955331 0.0671401 0.0346748 -0.135819 -0.124692 -0.107695 -0.0854945 -0.0593271 -0.0304143 0.12664 0.113079 0.0956095 0.0747122 0.0512845 0.0261243 -0.113061 -0.0987375 -0.0820456 -0.06328 -0.0430434 -0.0218088 0.0955561 0.0820121 0.0672129 0.0512973 0.034636 0.0174724 -0.0747315 -0.0633071 -0.0513407 -0.0388684 -0.0260947 -0.0131192 0.0513114 0.0430735 0.0346747 0.0261015 0.0174527 0.00875327 -0.0261091 -0.0217998 -0.0174724 -0.0131077 -0.00874317 -0.00437874
With two types
In fact, the symeigs functions all take in two types: the type of the Eigeninformation (which must be be non-complex) and the type of the vector information (which can be complex).
julia> op = ArpackSimpleOp(Float16.(A))ArpackSimpleOp{Tridiagonal{Float16, Vector{Float16}}}(Float16[2.0 -1.0 … 0.0 0.0; -1.0 2.0 … 0.0 0.0; … ; 0.0 0.0 … 2.0 -1.0; 0.0 0.0 … -1.0 2.0])julia> symeigs(ComplexF64, Float32, op, 6)GenericArpack.ArpackEigen{Float32, ComplexF64, ArpackSimpleOp{Tridiagonal{Float16, Vector{Float16}}}, ArpackState{Float32}, :I} eigenspace: LM values: 6-element Vector{Float32}: 3.965262 3.975854 3.9845324 3.991294 3.996123 3.9990332 vectors: 100×6 Matrix{ComplexF64}: 0.000182627-0.0261136im … 0.00402896-0.00172913im -0.000358975+0.0513203im -0.00805422+0.0034565im 0.000523109-0.0747447im 0.0120717-0.0051804im -0.000669184+0.0955731im -0.0160774+0.00689916im 0.000792235-0.113082im 0.0200673-0.00861124im -0.000887953+0.126664im … -0.0240375+0.0103152im 0.000953151-0.135846im 0.0279838-0.0120093im -0.000985599+0.140311im -0.0319028+0.013692im 0.000984304-0.139902im 0.0357907-0.0153616im -0.000949314+0.134635im -0.0396439+0.0170163im ⋮ ⋱ ⋮ 0.000934804-0.139877im -0.0356781+0.0152971im -0.000940488+0.140284im 0.0318012-0.0136336im 0.000912888-0.13582im -0.0278937+0.0119574im -0.000852898+0.126639im 0.0239592-0.0102699im 0.000762638-0.11306im … -0.0200016+0.00857282im -0.000645268+0.0955539im 0.0160245-0.00686782im 0.000505069-0.0747297im -0.012032+0.00515648im -0.000347007+0.05131im 0.00802774-0.00344032im 0.000176683-0.0261084im -0.00401576+0.00172094im
These eigenvectors are correct, even though they aren't real-valued. It's just that they are less unique in the complex plane compared to the real-plane.
This allows us to use higher or lower-precision in the vector compared with the eigeninformation (which also is used to represent the Arnoldi factorization).
This uses BigFloat for the Eigenvalue information and Float64 for the vectors. This will realistically limit you to Float64 accuracy but might handle some edge cases better. (And BigFloat should really be a different type, but this was handy to write as it's built into Julia.)
julia> symeigs(Float64, BigFloat, A, 3; ncv=36)GenericArpack.ArpackEigen{BigFloat, Float64, ArpackSimpleOp{Tridiagonal{Float64, Vector{Float64}}}, ArpackState{BigFloat}, :I} eigenspace: LM values: 3-element Vector{BigFloat}: 3.991298695938037159508597021220321966797436385660381715852115264051838787433121 3.996131194267188652353734193440331996350294707921790683980865441165101665783157 3.999032564583976105529006234362963633374988313375663912614167292763251549371105 vectors: 100×3 Matrix{Float64}: 0.0131121 -0.00874848 -0.00437636 -0.0261102 0.0174631 0.00874848 0.038881 -0.0261102 -0.0131121 -0.0513136 0.0346562 0.0174631 0.0632996 -0.0430682 -0.0217972 -0.0747349 0.0513136 0.0261102 0.0855198 -0.0593604 -0.0303979 -0.0955607 0.0671776 0.0346562 0.10477 -0.0747349 -0.038881 -0.113068 0.082003 0.0430682 ⋮ -0.10477 -0.0747349 0.038881 0.0955607 0.0671776 -0.0346562 -0.0855198 -0.0593604 0.0303979 0.0747349 0.0513136 -0.0261102 -0.0632996 -0.0430682 0.0217972 0.0513136 0.0346562 -0.0174631 -0.038881 -0.0261102 0.0131121 0.0261102 0.0174631 -0.00874848 -0.0131121 -0.00874848 0.00437636
Advanced usage
Writing your own ArpackOp
It is "easy" (says the developer) to write your own ArpackOp type.
Here is the code to wrap a function F in an ArpackOp. This is the complete implementation of the ArpackFunctionOp
struct ArpackSimpleFunctionOp <: ArpackOp
F::Function
n::Int
end
arpack_mode(::ArpackSimpleFunctionOp) = 1
Base.size(op::ArpackSimpleFunctionOp) = op.n
bmat(::ArpackSimpleFunctionOp) = Val(:I)
opx!(y,op::ArpackSimpleFunctionOp,x) = op.F(y,x)
is_arpack_mode_valid_for_op(mode::Int, ::ArpackSimpleFunctionOp) = mode == 1 The Arpack Drivers
All of the Arpack interface is through the functions dsaupd and dseupd. (Yes, Lapack-heads, I know that it should be `saupd/seupd` in proper nomenclature, but I'll get around to fixing that at some point.)
The key differences from the standard interface are four new types of arguments for Julia.
ArpackDebugcontrols the Arpack debug messages (fromdebug.h)ArpackStatscontrols the Arpack statistics collected (fromstats.h)idonowis the Julia information for avoiding the reverse communication interface. It's slightly faster. But likely to make compile times longer. I'm leaving it in even though it'll be a source of bugs, ugh.ArpackStatetracks the computation specific state that was in the Fortransavevariables. This includes things like the random number seeds.
For stats if the type is nothing, then we don't track stats. Likewise for debug.
The other key difference is that all of the functions return the info value instead of setting the info parameter. (I didn't want a needed Ref type hanging around in the Julia code.)
Here is how a dsaupd call maps between the two libraries. This is one we use in the test code.
ido = Ref{Int}(0); bmat=:I; which=:LM;
state = ArpackState{Float64}() # (or whatever the Eigenvector info is...)
stats = nothing # okay to not track...
ierr, state = GenericArpack.dsaupd!(ido, Val(bmat), n, which, nev, tol, resid, ncv, V, ldv, iparam,
ipntr, workd, workl, lworkl, info_initv;
state, stats, debug # these are the new parameters, which are Julia keywoard params, so any order is okay!
)
arierr = arpack_dsaupd!(arido, bmat, n, which, nev, tol, arresid, ncv, arV, ldv, ariparam,
aripntr, arworkd, arworkl, lworkl, info_initv)If tol=0, the tol parameter is typically initalized by Arpack to dlamch("E") which is eps(Float64)/2. This parameter is then propagated to future calls since everything in Fortran is pass-by-reference.
Instead of worrying about this, we just set tol on each call to GenericArpack.dsaupd!; however this can cause some small issues if you try and compare our calls directly. So in those cases, I recommend setting tol to be anything non-zero so Arpack won't touch it.
Getting debug information
You can also pass ArpackStats and ArpackDebug to any of the higher-level drivers.
julia> eigs(Symmetric(A), 3; debug=ArpackDebug(maupd=2, maup2=1), stats=ArpackStats()) # lots of convergence info_saup2: **** Start of major iteration number **** 1 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 2 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 3 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 4 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 5 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 6 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 7 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 8 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 9 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 10 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 11 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 12 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 13 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 14 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 15 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 16 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 17 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 18 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 19 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 20 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 21 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 22 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 23 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 24 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 25 _saup2: no. of "converged" Ritz values at this iter: 0 _saup2: **** Start of major iteration number **** 26 _saup2: no. of "converged" Ritz values at this iter: 1 _saup2: **** Start of major iteration number **** 27 _saup2: no. of "converged" Ritz values at this iter: 1 _saup2: **** Start of major iteration number **** 28 _saup2: no. of "converged" Ritz values at this iter: 2 _saup2: **** Start of major iteration number **** 29 _saup2: no. of "converged" Ritz values at this iter: 2 _saup2: **** Start of major iteration number **** 30 _saup2: no. of "converged" Ritz values at this iter: 2 _saup2: **** Start of major iteration number **** 31 _saupd: number of update iterations taken 31 _saupd: number of "converged" Ritz values 3 _saupd: final Ritz values ------------------------- 3.991298695938041 , 3.9961311942671642 , 3.9990325645839695 , # 1-3 _saupd: corresponding error bounds ---------------------------------- 7.757113726812292e-17 , 4.814945065525502e-22 , 2.5019471633814614e-24 , # 1-3 ========================================== = Symmetric implicit Arnoldi update code = = Version Number: 2.4 = = Version Date: 07/31/96 = ========================================== = Julia Port Version: 1.0 = = Julia Port Version Date: 2022-04-16 = ========================================== = Summary of timing statistics = ========================================== Total number update iterations = 31 Total number of OP*x operations = 522 Total number of B*x operations = 0 Total number of reorthogonalization steps = 521 Total number of iterative refinement steps = 0 Total number of restart steps = 0 Total time in user OP*x operation = 0.013231515884399414 Total time in user B*x operation = 0.0 Total time in Arnoldi update routine = 0.05061793327331543 Total time in saup2 routine = 0.038987159729003906 Total time in basic Arnoldi iteration loop = 0.014371633529663086 Total time in reorthogonalization phase = 0.0005669593811035156 Total time in (re)start vector generation = 7.152557373046875e-6 Total time in trid eigenvalue subproblem = 0.0005733966827392578 Total time in getting the shifts = 2.3603439331054688e-5 Total time in applying the shifts = 0.0004820823669433594 Total time in convergence testing = 0.0 GenericArpack.ArpackEigen{Float64, Float64, ArpackSimpleOp{Symmetric{Float64, Tridiagonal{Float64, Vector{Float64}}}}, ArpackState{Float64}, :I} eigenspace: LM values: 3-element Vector{Float64}: 3.991298695938041 3.9961311942671642 3.9990325645839695 vectors: 100×3 Matrix{Float64}: 0.0131121 -0.00874848 -0.00437636 -0.0261102 0.0174631 0.00874848 0.038881 -0.0261102 -0.0131121 -0.0513136 0.0346562 0.0174631 0.0632996 -0.0430682 -0.0217972 -0.0747349 0.0513136 0.0261102 0.0855198 -0.0593604 -0.0303979 -0.0955607 0.0671776 0.0346562 0.10477 -0.0747349 -0.038881 -0.113068 0.082003 0.0430682 ⋮ -0.10477 -0.0747349 0.038881 0.0955607 0.0671776 -0.0346562 -0.0855198 -0.0593604 0.0303979 0.0747349 0.0513136 -0.0261102 -0.0632996 -0.0430682 0.0217972 0.0513136 0.0346562 -0.0174631 -0.038881 -0.0261102 0.0131121 0.0261102 0.0174631 -0.00874848 -0.0131121 -0.00874848 0.00437636