### Install ArnoldiMethod.jl Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Use the Julia package manager to add the ArnoldiMethod package. ```julia (v1.6) pkg> add ArnoldiMethod ``` -------------------------------- ### Allocate ArnoldiWorkspace examples Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Demonstrates various ways to initialize the Arnoldi workspace, including automatic and manual memory allocation. ```julia # allocates workspace for 20-dimensional Krylov subspace arnoldi = ArnoldiWorkspace(100, 20) # allocate workspace for 20-dimensional Krylov subspace, with initial vector ones(100) copied into # the first column of V arnoldi = ArnoldiWorkspace(ones(100), 20) # manually allocate workspace with V, H V = Matrix{Float64}(undef, 100, 21) H = Matrix{Float64}(undef, 21, 20) arnoldi = ArnoldiWorkspace(V, H) # manually allocate all arrays, including temporaries V, tmp = Matrix{Float64}(undef, 100, 21), Matrix{Float64}(undef, 100, 21) H, Q = Matrix{Float64}(undef, 21, 20), Matrix{Float64}(undef, 20, 20) arnoldi = ArnoldiWorkspace(V, H, V_tmp = tmp, Q = Q) ``` -------------------------------- ### Load ArnoldiMethod.jl Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Import the package into your Julia session. ```julia using ArnoldiMethod ``` -------------------------------- ### Resuming partial Schur decomposition Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Demonstrates initializing an ArnoldiWorkspace and performing iterative partial Schur decomposition updates. ```julia A = rand(100, 100) # Pre-allocate the relevant Krylov subspace arrays V, H = rand(100, 21), rand(21, 20) arnoldi = ArnoldiWorkspace(V, H) # Find a few eigenvalues _, info_1 = partialschur!(A, arnoldi, nev = 3, tol = 1e-12) # Then continue to find a couple more. Notice: 5 in total, so 2 more. Allow larger errors by # changin `tol`. F, info_2 = partialschur!(A, arnoldi, nev = 5, start_from = info_1.nconverged + 1 , tol = 1e-8) @show norm(A * F.Q - F.Q * F.R) ``` -------------------------------- ### ArnoldiWorkspace Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Constructor for pre-allocating memory for the Arnoldi relation matrices. ```APIDOC ## ArnoldiWorkspace(n, k) / ArnoldiWorkspace(v1, k) / ArnoldiWorkspace(V, H; V_tmp, Q) ### Description Holds the large arrays for the Arnoldi relation: Vₖ₊₁ and Hₖ. Used for pre-allocating memory or working with custom matrix types. ### Parameters - **n** (Int) - Required - Dimension of the Krylov subspace. - **k** (Int) - Required - Size of the subspace. - **v1** (Vector) - Optional - Initial vector. - **V, H, V_tmp, Q** (Matrix) - Optional - Manually allocated matrices for the Arnoldi relation. ``` -------------------------------- ### partialschur(A; v1, workspace, nev, which, tol, mindim, maxdim, restarts) Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Computes an approximate, partial Schur decomposition of a matrix A. It runs iteratively until eigenpairs are approximated to the prescribed tolerance or restarts are exhausted. ```APIDOC ## partialschur(A; v1, workspace, nev, which, tol, mindim, maxdim, restarts) ### Description Finds `nev` approximate eigenpairs of `A` with eigenvalues near a specified target. The matrix `A` can be any linear map that implements `mul!(y, A, x)`, `eltype` and `size`. ### Parameters - **A** (Any) - Required - The linear map/matrix to decompose. - **nev** (Int) - Optional - Number of eigenvalues (Default: min(6, size(A, 1))). - **which** (Symbol or Target) - Optional - Target eigenvalue selection (Default: :LM). - **tol** (Real) - Optional - Tolerance for convergence (Default: √eps). - **v1** (AbstractVector) - Optional - Starting vector for the Krylov subspace. - **mindim** (Int) - Optional - Minimum Krylov dimension. - **maxdim** (Int) - Optional - Maximum Krylov dimension. - **restarts** (Int) - Optional - Maximum number of restarts (Default: 200). ### Returns - **decomp** (PartialSchur) - The partial Schur decomposition struct. - **history** (History) - Information about the convergence of the method. ``` -------------------------------- ### Define History struct Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Tracks convergence status and computational effort of the Arnoldi method. ```julia History(mvproducts, nconverged, converged, nev) ``` -------------------------------- ### Compute Eigenvalues and Eigenvectors Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Compute the first ten eigenvalues and eigenvectors of a tridiagonal sparse matrix using partialschur and partialeigen. ```julia julia> using ArnoldiMethod, LinearAlgebra, SparseArrays julia> A = spdiagm( -1 => fill(-1.0, 99), 0 => fill(2.0, 100), 1 => fill(-1.0, 99) ); julia> decomp, history = partialschur(A, nev=10, tol=1e-6, which=:SR); julia> decomp PartialSchur decomposition (Float64) of dimension 10 eigenvalues: 10-element Array{Complex{Float64},1}: 0.0009674354160236865 + 0.0im 0.003868805732811139 + 0.0im 0.008701304061962657 + 0.0im 0.01546025527344699 + 0.0im 0.024139120518486677 + 0.0im 0.0347295035554728 + 0.0im 0.04722115887278571 + 0.0im 0.06160200160067088 + 0.0im 0.0778581192025522 + 0.0im 0.09597378493453936 + 0.0im julia> history Converged: 10 of 10 eigenvalues in 174 matrix-vector products julia> norm(A * decomp.Q - decomp.Q * decomp.R) 6.39386920955869e-8 julia> λs, X = partialeigen(decomp); julia> norm(A * X - X * Diagonal(λs)) 6.393869211477937e-8 ``` -------------------------------- ### Check Convergence History Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Access the convergence status and information from the history object. ```julia > history.converged true > @show history Converged after 359 matrix-vector products ``` -------------------------------- ### partialschur!(A, arnoldi; start_from, initialize, nev, which, tol, mindim, maxdim, restarts) Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable A variant of partialschur that operates on a pre-allocated ArnoldiWorkspace instance to avoid allocations or continue from an existing decomposition. ```APIDOC ## partialschur!(A, arnoldi; start_from, initialize, nev, which, tol, mindim, maxdim, restarts) ### Description Operates on a pre-allocated `ArnoldiWorkspace` instance. Useful for continuing from an existing decomposition, using custom array types, or avoiding repeated allocations. ### Parameters - **A** (Any) - Required - The linear map/matrix. - **arnoldi** (ArnoldiWorkspace) - Required - Pre-allocated workspace. - **start_from** (Int) - Optional - Number of Schur vectors to continue from. - **initialize** (Bool) - Optional - Whether to initialize the Krylov subspace (Default: true). - **nev** (Int) - Optional - Number of eigenvalues. - **which** (Symbol or Target) - Optional - Target eigenvalue selection. - **tol** (Real) - Optional - Tolerance for convergence. - **mindim** (Int) - Optional - Minimum Krylov dimension. - **maxdim** (Int) - Optional - Maximum Krylov dimension. - **restarts** (Int) - Optional - Maximum number of restarts. ### Returns - **decomp** (PartialSchur) - The partial Schur decomposition struct containing views of the workspace. - **history** (History) - Information about the convergence of the method. ``` -------------------------------- ### Shift-and-invert for smallest eigenvalues Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Uses LinearMaps.jl to invert matrix A, allowing the solver to find the largest eigenvalues of A⁻¹, which correspond to the smallest eigenvalues of A. ```julia using ArnoldiMethod, LinearAlgebra, LinearMaps # Define a matrix whose eigenvalues you want A = rand(100,100) # Factorizes A and builds a linear map that applies inv(A) to a vector. function construct_linear_map(A) F = factorize(A) LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, x), size(A,1), ismutating=true) end # Target the largest eigenvalues of the inverted problem decomp, = partialschur(construct_linear_map(A), nev=4, tol=1e-5, restarts=100, which=:LM) λs_inv, X = partialeigen(decomp) # Eigenvalues have to be inverted to find the smallest eigenvalues of the non-inverted problem. λs = 1 ./ λs_inv # Show that Ax = xλ @show norm(A * X - X * Diagonal(λs)) # 7.38473677258669e-6 ``` -------------------------------- ### Shift-and-invert for generalized eigenvalue problems Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Recasts Ax=Bxλ into A⁻¹Bx=xθ to target eigenvalues closest to the origin. ```julia using ArnoldiMethod, LinearAlgebra, LinearMaps # Define the matrices of the generalized eigenvalue problem A, B = rand(100, 100), rand(100, 100) struct ShiftAndInvert{TA,TB,TT} A_lu::TA B::TB temp::TT end function (M::ShiftAndInvert)(y, x) mul!(M.temp, M.B, x) ldiv!(y, M.A_lu, M.temp) end function construct_linear_map(A, B) a = ShiftAndInvert(factorize(A), B, Vector{eltype(A)}(undef, size(A, 1))) LinearMap{eltype(A)}(a, size(A, 1), ismutating = true) end # Target the largest eigenvalues of the inverted problem decomp, = partialschur( construct_linear_map(A, B), nev = 4, tol = 1e-5, restarts = 100, which = :LM, ) λs_inv, X = partialeigen(decomp) # Eigenvalues have to be inverted to find the smallest eigenvalues of the non-inverted problem. λs = 1 ./ λs_inv # Show that Ax = Bxλ @show norm(A * X - B * X * Diagonal(λs)) # 2.8043149027575927e-6 ``` -------------------------------- ### partialschur Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Computes the partial Schur decomposition of a matrix A. Users can provide an initial guess vector v1 to speed up convergence. ```APIDOC ## partialschur(A, v1=my_initial_vector) ### Description Computes the partial Schur decomposition of a matrix A. An optional initial guess vector v1 can be provided to build the Krylov subspace. ### Parameters - **A** (Matrix) - Required - The matrix to decompose. - **v1** (Vector) - Optional - Initial guess for the target eigenvector. ``` -------------------------------- ### partialeigen(P::PartialSchur) Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Transforms a partial Schur decomposition into an eigendecomposition by computing the eigendecomposition of the upper triangular matrix R. ```APIDOC ## partialeigen(P::PartialSchur) ### Description Transforms a partial Schur decomposition into an eigendecomposition. This is a post-processing step that computes the eigendecomposition of the upper triangular matrix R from the partial Schur decomposition. ### Signature `partialeigen(P::PartialSchur) → (Vector{<:Union{Real,Complex}}, Matrix{<:Union{Real,Complex}})` ### Parameters - **P** (PartialSchur) - Required - The partial Schur decomposition object to be transformed. ### Returns - **Vector{<:Union{Real,Complex}}** - The eigenvalues. - **Matrix{<:Union{Real,Complex}}** - The eigenvectors. ``` -------------------------------- ### Verify Decomposition Norm Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Check the residual norm of the computed partial Schur decomposition. ```julia > norm(A * decomp.Q - decomp.Q * decomp.R) ``` -------------------------------- ### Partial Schur Decomposition Function Signature Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable The primary interface for computing a partial Schur decomposition of a linear map A. ```julia partialschur(A; v1, workspace, nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History ``` -------------------------------- ### Partial Schur Return Values Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable The function returns a tuple containing the decomposition structure and the convergence history. ```julia decomp, history = partialschur(A, ...) ``` -------------------------------- ### Define PartialSchur struct Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Represents the orthonormal basis and quasi-upper triangular matrix resulting from the partial Schur decomposition. ```julia PartialSchur(Q, R, eigenvalues) ``` -------------------------------- ### Generalized eigenvalue problem with symmetric positive-definite B Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Uses Cholesky decomposition of B to transform the generalized problem into a standard Schur decomposition. ```julia using ArnoldiMethod, LinearAlgebra, LinearMaps struct WithBInnerProduct{TA,TL} A::TA L::TL end function (M::WithBInnerProduct)(y, x) # Julia unfortunately does not have in-place CHOLMOD solve, so this is far from optimal. tmp = M.L \ x mul!(y, M.A, tmp) y .= M.L' \ y return y end # Define the matrices of the generalized eigenvalue problem A = rand(100, 100) B = Diagonal(range(1.0, 2.0, length = 100)) # Reformulate the problem as standard Schur decomposition F = cholesky(B) linear_map = LinearMap{eltype(A)}(WithBInnerProduct(A, F.L), size(A, 1), ismutating = true) decomp, info = partialschur(linear_map, nev = 4, which = :LM, tol = 1e-10) # Translate back to the original problem Q = F.L' \ decomp.Q @show norm(Q' * A * Q - decomp.R) # 3.883933945390996e-14 @show norm(Q' * B * Q - I) # 3.1672155003480583e-15 ``` -------------------------------- ### Compute partial eigendecomposition Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Transforms a PartialSchur object into a tuple containing eigenvalues and eigenvectors. Note that this is unnecessary for real-symmetric or complex-Hermitian matrices. ```julia partialeigen(P::PartialSchur) → (Vector{<:Union{Real,Complex}}, Matrix{<:Union{Real,Complex}}) ``` -------------------------------- ### In-place Partial Schur Decomposition Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable A variant of partialschur that uses a pre-allocated ArnoldiWorkspace to avoid allocations and support custom array types. ```julia partialschur!(A, arnoldi; start_from, initialize, nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History ``` -------------------------------- ### Define ArnoldiWorkspace struct Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable Constructs the workspace for storing large arrays used in the Arnoldi relation. ```julia ArnoldiWorkspace(n, k) → ArnoldiWorkspace ArnoldiWorkspace(v1, k) → ArnoldiWorkspace ArnoldiWorkspace(V, H; V_tmp, Q) → ArnoldiWorkspace ``` === COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.