From 6856bc972464d8c3dec8746d497c32533114d1f2 Mon Sep 17 00:00:00 2001 From: Jutho Date: Tue, 5 Nov 2024 22:41:59 +0100 Subject: [PATCH 1/9] cleanup and deprecate recursivevec --- .../KrylovKitChainRulesCoreExt.jl | 1 - ext/KrylovKitChainRulesCoreExt/constructor.jl | 6 - ext/KrylovKitChainRulesCoreExt/eigsolve.jl | 1 + src/KrylovKit.jl | 2 +- src/recursivevec.jl | 117 ------------------ 5 files changed, 2 insertions(+), 125 deletions(-) delete mode 100644 ext/KrylovKitChainRulesCoreExt/constructor.jl delete mode 100644 src/recursivevec.jl diff --git a/ext/KrylovKitChainRulesCoreExt/KrylovKitChainRulesCoreExt.jl b/ext/KrylovKitChainRulesCoreExt/KrylovKitChainRulesCoreExt.jl index 64a67b6..bf31bf0 100644 --- a/ext/KrylovKitChainRulesCoreExt/KrylovKitChainRulesCoreExt.jl +++ b/ext/KrylovKitChainRulesCoreExt/KrylovKitChainRulesCoreExt.jl @@ -11,6 +11,5 @@ include("utilities.jl") include("linsolve.jl") include("eigsolve.jl") include("svdsolve.jl") -include("constructor.jl") end # module diff --git a/ext/KrylovKitChainRulesCoreExt/constructor.jl b/ext/KrylovKitChainRulesCoreExt/constructor.jl deleted file mode 100644 index 04c13ee..0000000 --- a/ext/KrylovKitChainRulesCoreExt/constructor.jl +++ /dev/null @@ -1,6 +0,0 @@ -function ChainRulesCore.rrule(::Type{RecursiveVec}, A) - function RecursiveVec_pullback(ΔA) - return NoTangent(), ΔA.vecs - end - return RecursiveVec(A), RecursiveVec_pullback -end diff --git a/ext/KrylovKitChainRulesCoreExt/eigsolve.jl b/ext/KrylovKitChainRulesCoreExt/eigsolve.jl index e1c3131..b87f6f2 100644 --- a/ext/KrylovKitChainRulesCoreExt/eigsolve.jl +++ b/ext/KrylovKitChainRulesCoreExt/eigsolve.jl @@ -87,6 +87,7 @@ function make_eigsolve_pullback(config, f, fᴴ, x₀, howmany, which, alg_prima end end end + # Compute actual pullback data: #------------------------------ ws = compute_eigsolve_pullback_data(Δvals, Δvecs, view(vals, 1:n), view(vecs, 1:n), diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index 3aafc07..06e5f2e 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -238,7 +238,7 @@ include("matrixfun/exponentiate.jl") include("matrixfun/expintegrator.jl") # custom vector types -include("recursivevec.jl") include("innerproductvec.jl") +Base.@deprecate(RecursiveVec(args...), tuple(args...)) end diff --git a/src/recursivevec.jl b/src/recursivevec.jl deleted file mode 100644 index cc7379c..0000000 --- a/src/recursivevec.jl +++ /dev/null @@ -1,117 +0,0 @@ -""" - v = RecursiveVec(vecs) - -Create a new vector `v` from an existing (homogeneous or heterogeneous) list of vectors -`vecs` with one or more elements, represented as a `Tuple` or `AbstractVector`. The elements -of `vecs` can be any type of vectors that are supported by KrylovKit. For a heterogeneous -list, it is best to use a tuple for reasons of type stability, while for a homogeneous list, -either a `Tuple` or a `Vector` can be used. From a mathematical perspectve, `v` represents -the direct sum of the vectors in `vecs`. Scalar multiplication and addition of vectors `v` -acts simultaneously on all elements of `v.vecs`. The inner product corresponds to the sum -of the inner products of the individual vectors in the list `v.vecs`. - -The vector `v` also adheres to the iteration syntax, but where it will just produce the -individual vectors in `v.vecs`. Hence, `length(v) = length(v.vecs)`. It can also be indexed, -so that `v[i] = v.vecs[i]`, which can be useful in writing a linear map that acts on `v`. -""" -struct RecursiveVec{T<:Union{Tuple,AbstractVector}} - vecs::T -end -function RecursiveVec(arg1::AbstractVector{T}) where {T} - if isbitstype(T) - return RecursiveVec((arg1,)) - else - return RecursiveVec{typeof(arg1)}(arg1) - end -end -RecursiveVec(arg1, args...) = RecursiveVec((arg1, args...)) - -Base.getindex(v::RecursiveVec, i) = v.vecs[i] - -Base.iterate(v::RecursiveVec, args...) = iterate(v.vecs, args...) - -Base.IteratorEltype(::Type{RecursiveVec{T}}) where {T} = Base.IteratorEltype(T) -Base.IteratorSize(::Type{RecursiveVec{T}}) where {T} = Base.IteratorSize(T) - -Base.eltype(v::RecursiveVec) = eltype(v.vecs) -Base.size(v::RecursiveVec) = size(v.vecs) -Base.length(v::RecursiveVec) = length(v.vecs) - -Base.first(v::RecursiveVec) = first(v.vecs) -Base.last(v::RecursiveVec) = last(v.vecs) - -Base.:-(v::RecursiveVec) = RecursiveVec(map(-, v.vecs)) -Base.:+(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(+, v.vecs, w.vecs)) -Base.:-(v::RecursiveVec, w::RecursiveVec) = RecursiveVec(map(-, v.vecs, w.vecs)) -Base.:*(v::RecursiveVec, a::Number) = RecursiveVec(map(x -> x * a, v.vecs)) -Base.:*(a::Number, v::RecursiveVec) = RecursiveVec(map(x -> a * x, v.vecs)) -Base.:/(v::RecursiveVec, a::Number) = RecursiveVec(map(x -> x / a, v.vecs)) -Base.:\(a::Number, v::RecursiveVec) = RecursiveVec(map(x -> a \ x, v.vecs)) - -function Base.similar(v::RecursiveVec) - return RecursiveVec(similar.(v.vecs)) -end - -function Base.copy!(w::RecursiveVec, v::RecursiveVec) - @assert length(w) == length(v) - @inbounds for i in 1:length(w) - copyto!(w[i], v[i]) - end - return w -end - -function LinearAlgebra.dot(v::RecursiveVec{T}, w::RecursiveVec{T}) where {T} - return sum(dot.(v.vecs, w.vecs)) -end - -VectorInterface.scalartype(::Type{RecursiveVec{T}}) where {T} = scalartype(T) - -function VectorInterface.zerovector(v::RecursiveVec, T::Type{<:Number}) - return RecursiveVec(zerovector(v.vecs, T)) -end - -function VectorInterface.scale(v::RecursiveVec, a::Number) - return RecursiveVec(scale(v.vecs, a)) -end - -function VectorInterface.scale!(v::RecursiveVec, a::Number) - scale!(v.vecs, a) - return v -end - -function VectorInterface.scale!(w::RecursiveVec, v::RecursiveVec, a::Number) - scale!(w.vecs, v.vecs, a) - return w -end - -function VectorInterface.scale!!(x::RecursiveVec, a::Number) - return RecursiveVec(scale!!(x.vecs, a)) -end - -function VectorInterface.scale!!(w::RecursiveVec, - v::RecursiveVec, a::Number) - return RecursiveVec(scale!!(w.vecs, v.vecs, a)) -end - -function VectorInterface.add(w::RecursiveVec, v::RecursiveVec, a::Number=One(), - b::Number=One()) - return RecursiveVec(add(w.vecs, v.vecs, a, b)) -end - -function VectorInterface.add!(w::RecursiveVec, v::RecursiveVec, a::Number=One(), - b::Number=One()) - add!(w.vecs, v.vecs, a, b) - return w -end - -function VectorInterface.add!!(w::RecursiveVec, v::RecursiveVec, - a::Number=One(), - b::Number=One()) - return RecursiveVec(add!!(w.vecs, v.vecs, a, b)) -end - -function VectorInterface.inner(v::RecursiveVec{T}, w::RecursiveVec{T}) where {T} - return inner(v.vecs, w.vecs) -end - -VectorInterface.norm(v::RecursiveVec) = VectorInterface.norm(v.vecs) From e3957d3657caaf4b5b46932c85d79e8dc146186a Mon Sep 17 00:00:00 2001 From: Jutho Date: Tue, 5 Nov 2024 22:42:47 +0100 Subject: [PATCH 2/9] make x0 initialisation GPU compatible --- Project.toml | 6 +- src/KrylovKit.jl | 1 + src/eigsolve/eigsolve.jl | 4 +- src/eigsolve/geneigsolve.jl | 9 +- src/eigsolve/svdsolve.jl | 3 +- src/leastsquares/lsmr.jl | 1 + src/linsolve/lgmres.jl | 145 ++++++++++++++++++++ src/matrixfun/geometricseries.jl | 224 +++++++++++++++++++++++++++++++ 8 files changed, 385 insertions(+), 8 deletions(-) create mode 100644 src/leastsquares/lsmr.jl create mode 100644 src/linsolve/lgmres.jl create mode 100644 src/matrixfun/geometricseries.jl diff --git a/Project.toml b/Project.toml index 7caab5f..d625b5c 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PackageExtensionCompat = "65ce6f38-6b18-4e1d-a461-8949797d7930" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [weakdeps] @@ -27,7 +28,7 @@ PackageExtensionCompat = "1" Printf = "1" Random = "1" Test = "1" -TestExtras = "0.2" +TestExtras = "0.2,0.3" VectorInterface = "0.4" Zygote = "0.6" julia = "1.6" @@ -37,10 +38,9 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Test", "Aqua", "Random", "TestExtras", "ChainRulesTestUtils", "ChainRulesCore", "FiniteDifferences", "Zygote"] +test = ["Test", "Aqua", "TestExtras", "ChainRulesTestUtils", "ChainRulesCore", "FiniteDifferences", "Zygote"] diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index 06e5f2e..1f903c0 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -24,6 +24,7 @@ using VectorInterface using VectorInterface: add!! using LinearAlgebra using Printf +using Random using GPUArraysCore using PackageExtensionCompat const IndexRange = AbstractRange{Int} diff --git a/src/eigsolve/eigsolve.jl b/src/eigsolve/eigsolve.jl index fb1adf1..744467c 100644 --- a/src/eigsolve/eigsolve.jl +++ b/src/eigsolve/eigsolve.jl @@ -177,12 +177,14 @@ EigSorter(f::F; rev=false) where {F} = EigSorter{F}(f, rev) const Selector = Union{Symbol,EigSorter} +# TODO: initialize with `Random.rand!(similar(A, size(A, 1)))` for GPU arrays function eigsolve(A::AbstractMatrix, howmany::Int=1, which::Selector=:LM, T::Type=eltype(A); kwargs...) - return eigsolve(A, rand(T, size(A, 1)), howmany, which; kwargs...) + x₀ = Random.rand!(similar(A, T, size(A, 1))) + return eigsolve(A, x₀, howmany, which; kwargs...) end function eigsolve(f, n::Int, howmany::Int=1, which::Selector=:LM, T::Type=Float64; diff --git a/src/eigsolve/geneigsolve.jl b/src/eigsolve/geneigsolve.jl index ddddbe8..10a8cb1 100644 --- a/src/eigsolve/geneigsolve.jl +++ b/src/eigsolve/geneigsolve.jl @@ -150,21 +150,24 @@ function geneigsolve(AB::Tuple{AbstractMatrix,AbstractMatrix}, if !(size(AB[1], 1) == size(AB[1], 2) == size(AB[2], 1) == size(AB[2], 2)) throw(DimensionMismatch("Matrices `A` and `B` should be square and have matching size")) end - return geneigsolve(AB, rand(T, size(AB[1], 1)), howmany::Int, which; kwargs...) + x₀ = Random.rand!(similar(AB[1], T, size(AB[1], 1))) + return geneigsolve(AB, x₀, howmany::Int, which; kwargs...) end function geneigsolve(AB::Tuple{Any,AbstractMatrix}, howmany::Int=1, which::Selector=:LM, T=eltype(AB[2]); kwargs...) - return geneigsolve(AB, rand(T, size(AB[2], 1)), howmany, which; kwargs...) + x₀ = Random.rand!(similar(AB[2], T, size(AB[2], 1))) + return geneigsolve(AB, x₀, howmany, which; kwargs...) end function geneigsolve(AB::Tuple{AbstractMatrix,Any}, howmany::Int=1, which::Selector=:LM, T=eltype(AB[1]); kwargs...) - return geneigsolve(AB, rand(T, size(AB[1], 1)), howmany, which; kwargs...) + x₀ = Random.rand!(similar(AB[1], T, size(AB[1], 1))) + return geneigsolve(AB, x₀, howmany, which; kwargs...) end function geneigsolve(f, diff --git a/src/eigsolve/svdsolve.jl b/src/eigsolve/svdsolve.jl index 2043fab..2e8d5e3 100644 --- a/src/eigsolve/svdsolve.jl +++ b/src/eigsolve/svdsolve.jl @@ -121,7 +121,8 @@ function svdsolve(A::AbstractMatrix, which::Selector=:LR, T::Type=eltype(A); kwargs...) - return svdsolve(A, rand(T, size(A, 1)), howmany, which; kwargs...) + x₀ = Random.rand!(similar(A, T, size(A, 1))) + return svdsolve(A, x₀, howmany, which; kwargs...) end function svdsolve(f, n::Int, howmany::Int=1, which::Selector=:LR, T::Type=Float64; kwargs...) diff --git a/src/leastsquares/lsmr.jl b/src/leastsquares/lsmr.jl new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/leastsquares/lsmr.jl @@ -0,0 +1 @@ + diff --git a/src/linsolve/lgmres.jl b/src/linsolve/lgmres.jl new file mode 100644 index 0000000..ba26155 --- /dev/null +++ b/src/linsolve/lgmres.jl @@ -0,0 +1,145 @@ +function linsolve(operator, b, x₀, alg::GMRES, a₀::Number=0, a₁::Number=1) + # Initial function operation and division defines number type + y₀ = apply(operator, x₀) + T = typeof(dot(b, y₀) / norm(b) * one(a₀) * one(a₁)) + α₀ = convert(T, a₀)::T + α₁ = convert(T, a₁)::T + # Continue computing r = b - a₀ * x₀ - a₁ * operator(x₀) + r = one(T) * b # mul!(similar(b, T), b, 1) + r = iszero(α₀) ? r : axpy!(-α₀, x₀, r) + r = axpy!(-α₁, y₀, r) + x = mul!(similar(r), x₀, 1) + β = norm(r) + S = typeof(β) + + # Algorithm parameters + maxiter = alg.maxiter + krylovdim = alg.krylovdim + tol::S = alg.tol + + # Check for early return + if β < tol + if alg.verbosity > 0 + @info """GMRES linsolve converged without any iterations: + * norm of residual = $β + * number of operations = 1""" + end + return (x, ConvergenceInfo(1, r, β, 0, 1)) + end + + # Initialize data structures + y = Vector{T}(undef, krylovdim + 1) + gs = Vector{Givens{T}}(undef, krylovdim) + R = fill(zero(T), (krylovdim, krylovdim)) + numiter = 0 + numops = 1 # operator has been applied once to determine T + + iter = ArnoldiIterator(operator, r, alg.orth) + fact = initialize(iter) + numops += 1 # start applies operator once + + while numiter < maxiter # restart loop + numiter += 1 + y[1] = β + k = 1 + H = rayleighquotient(fact) + R[1, 1] = α₀ + α₁ * H[1, 1] + gs[1], R[1, 1] = givens(R[1, 1], α₁ * normres(fact), 1, 2) + y[2] = zero(T) + lmul!(gs[1], y) + β = convert(S, abs(y[2])) + if alg.verbosity > 2 + msg = "GMRES linsolve in iter $numiter; step $k: " + msg *= "normres = " + msg *= @sprintf("%.12e", β) + @info msg + end + + while (β > tol && length(fact) < krylovdim) # inner arnoldi loop + fact = expand!(iter, fact) + numops += 1 # expand! applies the operator once + k = length(fact) + H = rayleighquotient(fact) + + # copy Arnoldi Hessenberg matrix into R + @inbounds begin + for i in 1:(k - 1) + R[i, k] = α₁ * H[i, k] + end + R[k, k] = α₀ + α₁ * H[k, k] + end + + # Apply Givens rotations + Rk = view(R, :, k) + @inbounds for i in 1:(k - 1) + lmul!(gs[i], Rk) + end + gs[k], R[k, k] = givens(R[k, k], α₁ * normres(fact), k, k + 1) + + # Apply Givens rotations to right hand side + y[k + 1] = zero(T) + lmul!(gs[k], y) + + # New error + β = convert(S, abs(y[k + 1])) + if alg.verbosity > 2 + msg = "GMRES linsolve in iter $numiter; step $k: " + msg *= "normres = " + msg *= @sprintf("%.12e", β) + @info msg + end + end + if alg.verbosity > 1 + msg = "GMRES linsolve in iter $numiter; finished at step $k: " + msg *= "normres = " + msg *= @sprintf("%.12e", β) + @info msg + end + + # Solve upper triangular system + y2 = copy(y) + ldiv!(UpperTriangular(R), y, 1:k) + + # Update x + V = basis(fact) + @inbounds for i in 1:k + x = axpy!(y[i], V[i], x) + end + + if β > tol + # Recompute residual without reevaluating operator + w = residual(fact) + push!(V, rmul!(w, 1 / normres(fact))) + for i in 1:k + rmul!(V, gs[i]') + end + r = mul!(r, V[k + 1], y[k + 1]) + else + # Recompute residual and its norm explicitly, to ensure that no + # numerical errors have accumulated + r = mul!(r, b, 1) + r = axpy!(-1, apply(operator, x, α₀, α₁), r) + numops += 1 + β = norm(r) + if β < tol + if alg.verbosity > 0 + @info """GMRES linsolve converged at iteration $numiter, step $k: + * norm of residual = $β + * number of operations = $numops""" + end + return (x, ConvergenceInfo(1, r, β, numiter, numops)) + end + end + + # Restart Arnoldi factorization with new r + iter = ArnoldiIterator(operator, r, alg.orth) + fact = initialize!(iter, fact) + end + + if alg.verbosity > 0 + @warn """GMRES linsolve finished without converging after $numiter iterations: + * norm of residual = $β + * number of operations = $numops""" + end + return (x, ConvergenceInfo(0, r, β, numiter, numops)) +end diff --git a/src/matrixfun/geometricseries.jl b/src/matrixfun/geometricseries.jl new file mode 100644 index 0000000..4d87f2d --- /dev/null +++ b/src/matrixfun/geometricseries.jl @@ -0,0 +1,224 @@ +# """ +# function exponentiate(A, t::Number, x; kwargs...) +# function exponentiate(A, t::Number, x, algorithm) +# +# Compute ``y = exp(t*A) x``, where `A` is a general linear map, i.e. a `AbstractMatrix` or +# just a general function or callable object and `x` is of any Julia type with vector like +# behavior. +# +# ### Arguments: +# The linear map `A` can be an `AbstractMatrix` (dense or sparse) or a general function or +# callable object that implements the action of the linear map on a vector. If `A` is an +# `AbstractMatrix`, `x` is expected to be an `AbstractVector`, otherwise `x` can be of any +# type that behaves as a vector and supports the required methods (see KrylovKit docs). +# +# The time parameter `t` can be real or complex, and it is better to choose `t` e.g. imaginary +# and `A` hermitian, then to absorb the imaginary unit in an antihermitian `A`. For the +# former, the Lanczos scheme is used to built a Krylov subspace, in which an approximation to +# the exponential action of the linear map is obtained. The argument `x` can be of any type +# and should be in the domain of `A`. +# +# +# ### Return values: +# The return value is always of the form `y, info = eigsolve(...)` with +# * `y`: the result of the computation, i.e. `y = exp(t*A)*x` +# * `info`: an object of type [`ConvergenceInfo`], which has the following fields +# - `info.converged::Int`: 0 or 1 if the solution `y` was approximated up to the +# requested tolerance `tol`. +# - `info.residual::Nothing`: value `nothing`, there is no concept of a residual in +# this case +# - `info.normres::Real`: an estimate (upper bound) of the error between the +# approximate and exact solution +# - `info.numops::Int`: number of times the linear map was applied, i.e. number of times +# `f` was called, or a vector was multiplied with `A` +# - `info.numiter::Int`: number of times the Krylov subspace was restarted (see below) +# !!! warning "Check for convergence" +# No warning is printed if not all requested eigenvalues were converged, so always check +# if `info.converged >= howmany`. +# +# ### Keyword arguments: +# Keyword arguments and their default values are given by: +# * `krylovdim = 30`: the maximum dimension of the Krylov subspace that will be constructed. +# Note that the dimension of the vector space is not known or checked, e.g. `x₀` should +# not necessarily support the `Base.length` function. If you know the actual problem +# dimension is smaller than the default value, it is useful to reduce the value of +# `krylovdim`, though in principle this should be detected. +# * `tol = 1e-12`: the requested accuracy (corresponding to the 2-norm of the residual for +# Schur vectors, not the eigenvectors). If you work in e.g. single precision (`Float32`), +# you should definitely change the default value. +# * `maxiter::Int = 100`: the number of times the Krylov subspace can be rebuilt; see below +# for further details on the algorithms. +# * `info::Int = 0`: the level of verbosity, default is zero (no output) +# * `issymmetric`: if the linear map is symmetric, only meaningful if `T<:Real` +# * `ishermitian`: if the linear map is hermitian +# The default value for the last two depends on the method. If an `AbstractMatrix` is used, +# `issymmetric` and `ishermitian` are checked for that matrix, ortherwise the default values +# are `issymmetric = false` and `ishermitian = T <: Real && issymmetric`. +# +# ### Algorithm +# The last method, without default values and keyword arguments, is the one that is finally +# called, and can also be used directly. Here, one specifies the algorithm explicitly as +# either [`Lanczos`](@ref), for real symmetric or complex hermitian problems, or +# [`Arnoldi`](@ref), for general problems. Note that these names refer to the process for +# building the Krylov subspace. +# +# !!! warning "`Arnoldi` not yet implented" +# """ +function geometricseries end + +function geometricseries(A, b, x₀=rmul!(similar(b), false); kwargs...) + alg = eigselector(A, promote_type(eltype(b), eltype(x₀)); kwargs...) + return geometricseries(A, b, x₀, alg) +end + +function geometricseries(A, b, x₀, alg::Lanczos) + # Initial function operation and division defines number type + y₀ = apply(A, x₀) + numops = 1 + T = typeof(dot(b, y₀) / norm(b)) + + r = copyto!(similar(b, T), b) + r = axpy!(-1, x₀, r) + r = axpy!(+1, y₀, r) + β = norm(r) + S = typeof(β) + tol::S = alg.tol + x = copyto!(similar(r), x₀) + + # krylovdim and related allocations + krylovdim = alg.krylovdim + PP1 = Matrix{S}(undef, (krylovdim, krylovdim)) + PP2 = Matrix{S}(undef, (krylovdim, krylovdim)) + GG = Matrix{S}(undef, (krylovdim, krylovdim)) + + # initialize iterator + iter = LanczosIterator(A, r, alg.orth) + fact = initialize(iter; info=alg.info - 2) + numops += 1 + sizehint!(fact, krylovdim) + maxiter = alg.maxiter + numiter = 0 + while true + while normres(fact) > tol && length(fact) < krylovdim + fact = expand!(iter, fact; info=alg.info - 2) + numops += 1 + end + K = fact.k # current Krylov dimension + V = basis(fact) + + # Small matrix exponential and error estimation + H = rayleighquotient(fact) + P1 = copyto!(view(PP1, 1:K, 1:K), rayleighquotient(fact)) + P2 = copyto!(view(PP2, 1:K, 1:K), rayleighquotient(fact)) + GG = copyto!(view(GG, 1:K, 1:K), I) + + while (GG[K, 1] + P1[K, 1]) * normres(fact) < tol + GG .+= P1 + mul!(P2, P1, H) + P1, P2 = P2, P1 + end + unproject!(x, V, view(GG, :, 1), β, 1) + r = apply(A, x) + numops += 1 + r = axpy!(+1, b, r) + r = axpy!(-1, x, r) + β = norm(r) + numiter += 1 + + if alg.info > 1 + @info "Geometric series in iteration $numiter: norm residual $β" + end + if β <= tol + if alg.info > 0 + @info "Geometric series finished in iteration $numiter: norm residual $β" + end + return x, ConvergenceInfo(1, r, β, numiter, numops) + end + if numiter == maxiter + if alg.info > 0 + @warn "Geometric series finished without convergence after $numiter iterations" + end + return x, ConvergenceInfo(0, r, β, numiter, numops) + end + r = rmul!(r, 1 / β) + iter = LanczosIterator(A, r, alg.orth) + fact = initialize!(iter, fact; info=alg.info - 2) + numops += 1 + end +end + +function geometricseries(A, b, x₀, alg::Arnoldi) + # Initial function operation and division defines number type + y₀ = apply(A, x₀) + numops = 1 + T = typeof(dot(b, y₀) / norm(b)) + + r = copyto!(similar(b, T), b) + r = axpy!(-1, x₀, r) + r = axpy!(+1, y₀, r) + β = norm(r) + tol::typeof(β) = alg.tol + x = copyto!(similar(r), x₀) + + # krylovdim and related allocations + krylovdim = alg.krylovdim + HH = Matrix{T}(undef, (krylovdim, krylovdim)) + PP1 = Matrix{T}(undef, (krylovdim, krylovdim)) + PP2 = Matrix{T}(undef, (krylovdim, krylovdim)) + GG = Matrix{T}(undef, (krylovdim, krylovdim)) + + # initialize iterator + iter = ArnoldiIterator(A, r, alg.orth) + fact = initialize(iter; info=alg.info - 2) + numops += 1 + sizehint!(fact, krylovdim) + maxiter = alg.maxiter + numiter = 0 + while true + while normres(fact) > tol && length(fact) < krylovdim + fact = expand!(iter, fact; info=alg.info - 2) + numops += 1 + end + K = fact.k # current Krylov dimension + V = basis(fact) + + # Small matrix exponential and error estimation + H = copyto!(view(HH, 1:K, 1:K), rayleighquotient(fact)) + P1 = copyto!(view(PP1, 1:K, 1:K), rayleighquotient(fact)) + P2 = copyto!(view(PP2, 1:K, 1:K), rayleighquotient(fact)) + GG = copyto!(view(GG, 1:K, 1:K), I) + + while (GG[K, 1] + P1[K, 1]) * normres(fact) < tol + GG .+= P1 + mul!(P2, P1, H) + P1, P2 = P2, P1 + end + unproject!(x, V, view(GG, :, 1), β, 1) + r = apply(A, x) + numops += 1 + r = axpy!(+1, b, r) + r = axpy!(-1, x, r) + β = norm(r) + numiter += 1 + + if alg.info > 1 + @info "Geometric series in iteration $numiter: norm residual $β" + end + if β <= tol + if alg.info > 0 + @info "Geometric series finished in iteration $numiter: norm residual $β" + end + return x, ConvergenceInfo(1, r, β, numiter, numops) + end + if numiter == maxiter + if alg.info > 0 + @warn "Geometric series finished without convergence after $numiter iterations" + end + return x, ConvergenceInfo(0, r, β, numiter, numops) + end + r = rmul!(r, 1 / β) + iter = ArnoldiIterator(A, r, alg.orth) + fact = initialize!(iter, fact; info=alg.info - 2) + numops += 1 + end +end From 6e1512e90160bb75500ee0d6edeb8cb210fb4d79 Mon Sep 17 00:00:00 2001 From: Jutho Date: Tue, 5 Nov 2024 22:48:25 +0100 Subject: [PATCH 3/9] update CI and reduce number of combinations --- .github/workflows/ci-nightly.yml | 6 ++---- .github/workflows/ci.yml | 14 +++++++------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/.github/workflows/ci-nightly.yml b/.github/workflows/ci-nightly.yml index 679f455..5955468 100644 --- a/.github/workflows/ci-nightly.yml +++ b/.github/workflows/ci-nightly.yml @@ -24,16 +24,14 @@ jobs: - 'nightly' os: - ubuntu-latest - - macOS-latest - - windows-latest arch: - x64 steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 96c2f40..d118bc8 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -21,7 +21,8 @@ jobs: fail-fast: false matrix: version: - - '1.6' # LTS version + - '1.6' + - 'lts' - '1' # automatically expands to the latest stable 1.x release of Julia os: - ubuntu-latest @@ -31,11 +32,11 @@ jobs: - x64 steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest - uses: julia-actions/julia-processcoverage@v1 @@ -52,20 +53,19 @@ jobs: fail-fast: false matrix: version: - - '1.6' # LTS version + - 'lts' - '1' # automatically expands to the latest stable 1.x release of Julia os: - ubuntu-latest - - macOS-latest arch: - x64 steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest env: From 4dfe2bcfe1eaeabd0cc65755e55e18054b255b9e Mon Sep 17 00:00:00 2001 From: Jutho Date: Tue, 5 Nov 2024 22:51:26 +0100 Subject: [PATCH 4/9] remove erroneous files --- src/leastsquares/lsmr.jl | 1 - src/linsolve/lgmres.jl | 145 -------------------- src/matrixfun/geometricseries.jl | 224 ------------------------------- 3 files changed, 370 deletions(-) delete mode 100644 src/leastsquares/lsmr.jl delete mode 100644 src/linsolve/lgmres.jl delete mode 100644 src/matrixfun/geometricseries.jl diff --git a/src/leastsquares/lsmr.jl b/src/leastsquares/lsmr.jl deleted file mode 100644 index 8b13789..0000000 --- a/src/leastsquares/lsmr.jl +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/linsolve/lgmres.jl b/src/linsolve/lgmres.jl deleted file mode 100644 index ba26155..0000000 --- a/src/linsolve/lgmres.jl +++ /dev/null @@ -1,145 +0,0 @@ -function linsolve(operator, b, x₀, alg::GMRES, a₀::Number=0, a₁::Number=1) - # Initial function operation and division defines number type - y₀ = apply(operator, x₀) - T = typeof(dot(b, y₀) / norm(b) * one(a₀) * one(a₁)) - α₀ = convert(T, a₀)::T - α₁ = convert(T, a₁)::T - # Continue computing r = b - a₀ * x₀ - a₁ * operator(x₀) - r = one(T) * b # mul!(similar(b, T), b, 1) - r = iszero(α₀) ? r : axpy!(-α₀, x₀, r) - r = axpy!(-α₁, y₀, r) - x = mul!(similar(r), x₀, 1) - β = norm(r) - S = typeof(β) - - # Algorithm parameters - maxiter = alg.maxiter - krylovdim = alg.krylovdim - tol::S = alg.tol - - # Check for early return - if β < tol - if alg.verbosity > 0 - @info """GMRES linsolve converged without any iterations: - * norm of residual = $β - * number of operations = 1""" - end - return (x, ConvergenceInfo(1, r, β, 0, 1)) - end - - # Initialize data structures - y = Vector{T}(undef, krylovdim + 1) - gs = Vector{Givens{T}}(undef, krylovdim) - R = fill(zero(T), (krylovdim, krylovdim)) - numiter = 0 - numops = 1 # operator has been applied once to determine T - - iter = ArnoldiIterator(operator, r, alg.orth) - fact = initialize(iter) - numops += 1 # start applies operator once - - while numiter < maxiter # restart loop - numiter += 1 - y[1] = β - k = 1 - H = rayleighquotient(fact) - R[1, 1] = α₀ + α₁ * H[1, 1] - gs[1], R[1, 1] = givens(R[1, 1], α₁ * normres(fact), 1, 2) - y[2] = zero(T) - lmul!(gs[1], y) - β = convert(S, abs(y[2])) - if alg.verbosity > 2 - msg = "GMRES linsolve in iter $numiter; step $k: " - msg *= "normres = " - msg *= @sprintf("%.12e", β) - @info msg - end - - while (β > tol && length(fact) < krylovdim) # inner arnoldi loop - fact = expand!(iter, fact) - numops += 1 # expand! applies the operator once - k = length(fact) - H = rayleighquotient(fact) - - # copy Arnoldi Hessenberg matrix into R - @inbounds begin - for i in 1:(k - 1) - R[i, k] = α₁ * H[i, k] - end - R[k, k] = α₀ + α₁ * H[k, k] - end - - # Apply Givens rotations - Rk = view(R, :, k) - @inbounds for i in 1:(k - 1) - lmul!(gs[i], Rk) - end - gs[k], R[k, k] = givens(R[k, k], α₁ * normres(fact), k, k + 1) - - # Apply Givens rotations to right hand side - y[k + 1] = zero(T) - lmul!(gs[k], y) - - # New error - β = convert(S, abs(y[k + 1])) - if alg.verbosity > 2 - msg = "GMRES linsolve in iter $numiter; step $k: " - msg *= "normres = " - msg *= @sprintf("%.12e", β) - @info msg - end - end - if alg.verbosity > 1 - msg = "GMRES linsolve in iter $numiter; finished at step $k: " - msg *= "normres = " - msg *= @sprintf("%.12e", β) - @info msg - end - - # Solve upper triangular system - y2 = copy(y) - ldiv!(UpperTriangular(R), y, 1:k) - - # Update x - V = basis(fact) - @inbounds for i in 1:k - x = axpy!(y[i], V[i], x) - end - - if β > tol - # Recompute residual without reevaluating operator - w = residual(fact) - push!(V, rmul!(w, 1 / normres(fact))) - for i in 1:k - rmul!(V, gs[i]') - end - r = mul!(r, V[k + 1], y[k + 1]) - else - # Recompute residual and its norm explicitly, to ensure that no - # numerical errors have accumulated - r = mul!(r, b, 1) - r = axpy!(-1, apply(operator, x, α₀, α₁), r) - numops += 1 - β = norm(r) - if β < tol - if alg.verbosity > 0 - @info """GMRES linsolve converged at iteration $numiter, step $k: - * norm of residual = $β - * number of operations = $numops""" - end - return (x, ConvergenceInfo(1, r, β, numiter, numops)) - end - end - - # Restart Arnoldi factorization with new r - iter = ArnoldiIterator(operator, r, alg.orth) - fact = initialize!(iter, fact) - end - - if alg.verbosity > 0 - @warn """GMRES linsolve finished without converging after $numiter iterations: - * norm of residual = $β - * number of operations = $numops""" - end - return (x, ConvergenceInfo(0, r, β, numiter, numops)) -end diff --git a/src/matrixfun/geometricseries.jl b/src/matrixfun/geometricseries.jl deleted file mode 100644 index 4d87f2d..0000000 --- a/src/matrixfun/geometricseries.jl +++ /dev/null @@ -1,224 +0,0 @@ -# """ -# function exponentiate(A, t::Number, x; kwargs...) -# function exponentiate(A, t::Number, x, algorithm) -# -# Compute ``y = exp(t*A) x``, where `A` is a general linear map, i.e. a `AbstractMatrix` or -# just a general function or callable object and `x` is of any Julia type with vector like -# behavior. -# -# ### Arguments: -# The linear map `A` can be an `AbstractMatrix` (dense or sparse) or a general function or -# callable object that implements the action of the linear map on a vector. If `A` is an -# `AbstractMatrix`, `x` is expected to be an `AbstractVector`, otherwise `x` can be of any -# type that behaves as a vector and supports the required methods (see KrylovKit docs). -# -# The time parameter `t` can be real or complex, and it is better to choose `t` e.g. imaginary -# and `A` hermitian, then to absorb the imaginary unit in an antihermitian `A`. For the -# former, the Lanczos scheme is used to built a Krylov subspace, in which an approximation to -# the exponential action of the linear map is obtained. The argument `x` can be of any type -# and should be in the domain of `A`. -# -# -# ### Return values: -# The return value is always of the form `y, info = eigsolve(...)` with -# * `y`: the result of the computation, i.e. `y = exp(t*A)*x` -# * `info`: an object of type [`ConvergenceInfo`], which has the following fields -# - `info.converged::Int`: 0 or 1 if the solution `y` was approximated up to the -# requested tolerance `tol`. -# - `info.residual::Nothing`: value `nothing`, there is no concept of a residual in -# this case -# - `info.normres::Real`: an estimate (upper bound) of the error between the -# approximate and exact solution -# - `info.numops::Int`: number of times the linear map was applied, i.e. number of times -# `f` was called, or a vector was multiplied with `A` -# - `info.numiter::Int`: number of times the Krylov subspace was restarted (see below) -# !!! warning "Check for convergence" -# No warning is printed if not all requested eigenvalues were converged, so always check -# if `info.converged >= howmany`. -# -# ### Keyword arguments: -# Keyword arguments and their default values are given by: -# * `krylovdim = 30`: the maximum dimension of the Krylov subspace that will be constructed. -# Note that the dimension of the vector space is not known or checked, e.g. `x₀` should -# not necessarily support the `Base.length` function. If you know the actual problem -# dimension is smaller than the default value, it is useful to reduce the value of -# `krylovdim`, though in principle this should be detected. -# * `tol = 1e-12`: the requested accuracy (corresponding to the 2-norm of the residual for -# Schur vectors, not the eigenvectors). If you work in e.g. single precision (`Float32`), -# you should definitely change the default value. -# * `maxiter::Int = 100`: the number of times the Krylov subspace can be rebuilt; see below -# for further details on the algorithms. -# * `info::Int = 0`: the level of verbosity, default is zero (no output) -# * `issymmetric`: if the linear map is symmetric, only meaningful if `T<:Real` -# * `ishermitian`: if the linear map is hermitian -# The default value for the last two depends on the method. If an `AbstractMatrix` is used, -# `issymmetric` and `ishermitian` are checked for that matrix, ortherwise the default values -# are `issymmetric = false` and `ishermitian = T <: Real && issymmetric`. -# -# ### Algorithm -# The last method, without default values and keyword arguments, is the one that is finally -# called, and can also be used directly. Here, one specifies the algorithm explicitly as -# either [`Lanczos`](@ref), for real symmetric or complex hermitian problems, or -# [`Arnoldi`](@ref), for general problems. Note that these names refer to the process for -# building the Krylov subspace. -# -# !!! warning "`Arnoldi` not yet implented" -# """ -function geometricseries end - -function geometricseries(A, b, x₀=rmul!(similar(b), false); kwargs...) - alg = eigselector(A, promote_type(eltype(b), eltype(x₀)); kwargs...) - return geometricseries(A, b, x₀, alg) -end - -function geometricseries(A, b, x₀, alg::Lanczos) - # Initial function operation and division defines number type - y₀ = apply(A, x₀) - numops = 1 - T = typeof(dot(b, y₀) / norm(b)) - - r = copyto!(similar(b, T), b) - r = axpy!(-1, x₀, r) - r = axpy!(+1, y₀, r) - β = norm(r) - S = typeof(β) - tol::S = alg.tol - x = copyto!(similar(r), x₀) - - # krylovdim and related allocations - krylovdim = alg.krylovdim - PP1 = Matrix{S}(undef, (krylovdim, krylovdim)) - PP2 = Matrix{S}(undef, (krylovdim, krylovdim)) - GG = Matrix{S}(undef, (krylovdim, krylovdim)) - - # initialize iterator - iter = LanczosIterator(A, r, alg.orth) - fact = initialize(iter; info=alg.info - 2) - numops += 1 - sizehint!(fact, krylovdim) - maxiter = alg.maxiter - numiter = 0 - while true - while normres(fact) > tol && length(fact) < krylovdim - fact = expand!(iter, fact; info=alg.info - 2) - numops += 1 - end - K = fact.k # current Krylov dimension - V = basis(fact) - - # Small matrix exponential and error estimation - H = rayleighquotient(fact) - P1 = copyto!(view(PP1, 1:K, 1:K), rayleighquotient(fact)) - P2 = copyto!(view(PP2, 1:K, 1:K), rayleighquotient(fact)) - GG = copyto!(view(GG, 1:K, 1:K), I) - - while (GG[K, 1] + P1[K, 1]) * normres(fact) < tol - GG .+= P1 - mul!(P2, P1, H) - P1, P2 = P2, P1 - end - unproject!(x, V, view(GG, :, 1), β, 1) - r = apply(A, x) - numops += 1 - r = axpy!(+1, b, r) - r = axpy!(-1, x, r) - β = norm(r) - numiter += 1 - - if alg.info > 1 - @info "Geometric series in iteration $numiter: norm residual $β" - end - if β <= tol - if alg.info > 0 - @info "Geometric series finished in iteration $numiter: norm residual $β" - end - return x, ConvergenceInfo(1, r, β, numiter, numops) - end - if numiter == maxiter - if alg.info > 0 - @warn "Geometric series finished without convergence after $numiter iterations" - end - return x, ConvergenceInfo(0, r, β, numiter, numops) - end - r = rmul!(r, 1 / β) - iter = LanczosIterator(A, r, alg.orth) - fact = initialize!(iter, fact; info=alg.info - 2) - numops += 1 - end -end - -function geometricseries(A, b, x₀, alg::Arnoldi) - # Initial function operation and division defines number type - y₀ = apply(A, x₀) - numops = 1 - T = typeof(dot(b, y₀) / norm(b)) - - r = copyto!(similar(b, T), b) - r = axpy!(-1, x₀, r) - r = axpy!(+1, y₀, r) - β = norm(r) - tol::typeof(β) = alg.tol - x = copyto!(similar(r), x₀) - - # krylovdim and related allocations - krylovdim = alg.krylovdim - HH = Matrix{T}(undef, (krylovdim, krylovdim)) - PP1 = Matrix{T}(undef, (krylovdim, krylovdim)) - PP2 = Matrix{T}(undef, (krylovdim, krylovdim)) - GG = Matrix{T}(undef, (krylovdim, krylovdim)) - - # initialize iterator - iter = ArnoldiIterator(A, r, alg.orth) - fact = initialize(iter; info=alg.info - 2) - numops += 1 - sizehint!(fact, krylovdim) - maxiter = alg.maxiter - numiter = 0 - while true - while normres(fact) > tol && length(fact) < krylovdim - fact = expand!(iter, fact; info=alg.info - 2) - numops += 1 - end - K = fact.k # current Krylov dimension - V = basis(fact) - - # Small matrix exponential and error estimation - H = copyto!(view(HH, 1:K, 1:K), rayleighquotient(fact)) - P1 = copyto!(view(PP1, 1:K, 1:K), rayleighquotient(fact)) - P2 = copyto!(view(PP2, 1:K, 1:K), rayleighquotient(fact)) - GG = copyto!(view(GG, 1:K, 1:K), I) - - while (GG[K, 1] + P1[K, 1]) * normres(fact) < tol - GG .+= P1 - mul!(P2, P1, H) - P1, P2 = P2, P1 - end - unproject!(x, V, view(GG, :, 1), β, 1) - r = apply(A, x) - numops += 1 - r = axpy!(+1, b, r) - r = axpy!(-1, x, r) - β = norm(r) - numiter += 1 - - if alg.info > 1 - @info "Geometric series in iteration $numiter: norm residual $β" - end - if β <= tol - if alg.info > 0 - @info "Geometric series finished in iteration $numiter: norm residual $β" - end - return x, ConvergenceInfo(1, r, β, numiter, numops) - end - if numiter == maxiter - if alg.info > 0 - @warn "Geometric series finished without convergence after $numiter iterations" - end - return x, ConvergenceInfo(0, r, β, numiter, numops) - end - r = rmul!(r, 1 / β) - iter = ArnoldiIterator(A, r, alg.orth) - fact = initialize!(iter, fact; info=alg.info - 2) - numops += 1 - end -end From 6d3bcbf6a593aee5560f054708de741281047afc Mon Sep 17 00:00:00 2001 From: Jutho Date: Tue, 5 Nov 2024 22:58:16 +0100 Subject: [PATCH 5/9] remove todo Co-authored-by: Lukas Devos --- src/eigsolve/eigsolve.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/eigsolve/eigsolve.jl b/src/eigsolve/eigsolve.jl index 744467c..d687f8e 100644 --- a/src/eigsolve/eigsolve.jl +++ b/src/eigsolve/eigsolve.jl @@ -177,7 +177,6 @@ EigSorter(f::F; rev=false) where {F} = EigSorter{F}(f, rev) const Selector = Union{Symbol,EigSorter} -# TODO: initialize with `Random.rand!(similar(A, size(A, 1)))` for GPU arrays function eigsolve(A::AbstractMatrix, howmany::Int=1, which::Selector=:LM, From 97795a1a0542e3e73c4b478fd1bcaa8b14b0fea8 Mon Sep 17 00:00:00 2001 From: Jutho Date: Tue, 5 Nov 2024 23:38:27 +0100 Subject: [PATCH 6/9] fix docs and tests for recursivevec removal --- docs/src/index.md | 24 ++++++++---------------- docs/src/man/intro.md | 6 +++--- src/orthonormal.jl | 2 +- test/{recursivevec.jl => nestedtuple.jl} | 8 ++++---- test/runtests.jl | 2 +- 5 files changed, 17 insertions(+), 25 deletions(-) rename test/{recursivevec.jl => nestedtuple.jl} (90%) diff --git a/docs/src/index.md b/docs/src/index.md index 392fd9a..8419068 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -49,7 +49,7 @@ There are already a fair number of packages with Krylov-based or other iterative contains implementations of [high order exponential integrators](https://docs.juliadiffeq.org/latest/solvers/split_ode_solve/#OrdinaryDiffEq.jl-2) with adaptive Krylov-subspace calculations for solving semilinear and nonlinear ODEs. -These packages have certainly inspired and influenced the development of KrylovKit.jl. +Some of these packages have certainly inspired and influenced the development of KrylovKit.jl. However, KrylovKit.jl distinguishes itself from the previous packages in the following ways: 1. KrylovKit accepts general functions to represent the linear map or operator that defines @@ -63,21 +63,13 @@ However, KrylovKit.jl distinguishes itself from the previous packages in the fol 2. KrylovKit does not assume that the vectors involved in the problem are actual subtypes of `AbstractVector`. Any Julia object that behaves as a vector is supported, so in particular higher-dimensional arrays or any custom user type that supports the - interface as defined in - [`VectorInterface.jl`](https://github.com/Jutho/VectorInterface.jl) - - Algorithms in KrylovKit.jl are tested against such a minimal implementation (named - `MinimalVec`) in the test suite. This type is only defined in the tests. However, - KrylovKit provides two types implementing this interface and slightly more, to make - them behave more like `AbstractArrays` (e.g. also `Base.:+` etc), which can facilitate - certain applications: - * [`RecursiveVec`](@ref) can be used for grouping a set of vectors into a single - vector like structure (can be used recursively). This is more robust than trying to - use nested `Vector{<:Vector}` types. - * [`InnerProductVec`](@ref) can be used to redefine the inner product (i.e. `inner`) - and corresponding norm (`norm`) of an already existing vector like object. The - latter should help with implementing certain type of preconditioners. - + interface as defined in [`VectorInterface.jl`](https://github.com/Jutho/VectorInterface.jl). + Aside from arrays filled with scalar entries, this includes tuples, named tuples, and + arbitrarily nested combinations of tuples and arrays. Furthermore, `CuArray` objects + are fully supported as vectors, so that the application of the linear operator on the + vector can be executed on a GPU. The computations performed within the Krylov subspace, + such as diagonalising the projected matrix, are however always performed on the CPU. + 3. Since version 0.8, KrylovKit.jl supports reverse-mode AD by defining `ChainRulesCore.rrule` definitions for the most common functionality (`linsolve`, `eigsolve`, `svdsolve`). Hence, reverse mode AD engines that are compatible with the [ChainRules](https://juliadiff.org/ChainRulesCore.jl/dev/) diff --git a/docs/src/man/intro.md b/docs/src/man/intro.md index 1f16685..b60bb8c 100644 --- a/docs/src/man/intro.md +++ b/docs/src/man/intro.md @@ -76,9 +76,9 @@ results..., info = problemsolver(A, args..., algorithm(; kwargs...)) Most `algorithm` constructions take the same keyword arguments (`tol`, `krylovdim`, `maxiter` and `verbosity`) discussed above. -As mentioned before, there are two auxiliary structs that can be used to define new vectors, -namely +While KrylovKit.jl does currently not provide a general interface for including +preconditioners, it is possible to e.g. use a modified inner product. KrylovKit.jl provides +a specific type for this purpose: ```@docs -RecursiveVec InnerProductVec ``` diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 5c3ec97..ac40048 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -375,7 +375,7 @@ end orthogonalize(v, args...) = orthogonalize!(true * v, args...) function orthogonalize!!(v::T, b::OrthonormalBasis{T}, alg::Orthogonalizer) where {T} - S = promote_type(eltype(v), eltype(T)) + S = promote_type(scalartype(v), scalartype(T)) c = Vector{S}(undef, length(b)) return orthogonalize!!(v, b, c, alg) end diff --git a/test/recursivevec.jl b/test/nestedtuple.jl similarity index 90% rename from test/recursivevec.jl rename to test/nestedtuple.jl index 08d5a4c..7ff1fc9 100644 --- a/test/recursivevec.jl +++ b/test/nestedtuple.jl @@ -3,13 +3,13 @@ @testset for orth in (cgs2, mgs2, cgsr, mgsr) A = rand(T, (n, n)) v = rand(T, (n,)) - v2 = RecursiveVec(v, zero(v)) + v2 = (v, zero(v)) alg = Lanczos(; orth=orth, krylovdim=2 * n, maxiter=1, tol=tolerance(T)) D, V, info = eigsolve(v2, n, :LR, alg) do x x1, x2 = x y1 = A * x2 y2 = A' * x1 - return RecursiveVec(y1, y2) + return (y1, y2) end @test info.converged >= n S = D[1:n] @@ -27,14 +27,14 @@ end A = rand(T, (N, 2 * N)) v = rand(T, (N,)) w = rand(T, (2 * N,)) - v2 = RecursiveVec(v, w) + v2 = (v, w) alg = Lanczos(; orth=orth, krylovdim=n, maxiter=300, tol=tolerance(T)) n1 = div(n, 2) D, V, info = eigsolve(v2, n1, :LR, alg) do x x1, x2 = x y1 = A * x2 y2 = A' * x1 - return RecursiveVec(y1, y2) + return (y1, y2) end @test info.converged >= n1 S = D[1:n1] diff --git a/test/runtests.jl b/test/runtests.jl index 16a35c3..aa402e5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,7 +36,7 @@ include("svdsolve.jl") include("expintegrator.jl") include("linalg.jl") -include("recursivevec.jl") +include("nestedtuple.jl") include("ad.jl") From 371a0e288d4bc8e32e1743d9a55a4c9daacaff3b Mon Sep 17 00:00:00 2001 From: Jutho Date: Thu, 7 Nov 2024 11:07:34 +0100 Subject: [PATCH 7/9] remove unused code block --- ext/KrylovKitChainRulesCoreExt/eigsolve.jl | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/ext/KrylovKitChainRulesCoreExt/eigsolve.jl b/ext/KrylovKitChainRulesCoreExt/eigsolve.jl index b87f6f2..bb14a0b 100644 --- a/ext/KrylovKitChainRulesCoreExt/eigsolve.jl +++ b/ext/KrylovKitChainRulesCoreExt/eigsolve.jl @@ -65,11 +65,6 @@ function make_eigsolve_pullback(config, f, fᴴ, x₀, howmany, which, alg_prima return ∂self, ∂f, ∂x₀, ∂howmany, ∂which, ∂alg end end - if n < length(vals) && vals[n + 1] ≈ conj(vals[n]) - # this can probably only happen for real problems, where it would be problematic - # to split complex conjugate pairs in solving the tangent problem - n += 1 - end Δvals = fill(zero(vals[1]), n) if n_vals > 0 Δvals[1:n_vals] .= view(_Δvals, 1:n_vals) @@ -121,12 +116,15 @@ function compute_eigsolve_pullback_data(Δvals, Δvecs, vals, vecs, info, which, continue end - # General case : - - # for the case where `f` is a real matrix, we can expect the following simplication - # TODO: can we implement this within our general approach, or generalise this to also - # cover the case where `f` is a function? - # if i > 1 && eltype(A) <: Real && + # TODO: Is the following useful and correct? + # (given that Δvecs might contain weird tangent types) + # The following only holds if `f` represents a real linear operator, which we cannot + # check explicitly, unless `f isa AbstractMatrix`.` + # However, exact equality between conjugate pairs of eigenvalues and eigenvectors + # seems sufficient to guarantee this + # Also, we can only be sure to know how to apply complex conjugation when the + # vectors are of type `AbstractArray{T}` with `T` the scalar type + # if i > 1 && ws[i - 1] isa AbstractArray{T} && # vals[i] == conj(vals[i - 1]) && Δvals[i] == conj(Δvals[i - 1]) && # vecs[i] == conj(vecs[i - 1]) && Δvecs[i] == conj(Δvecs[i - 1]) # ws[i] = conj(ws[i - 1]) From d164897cf0f68ff929ff9966cc2ae2bacf0171a4 Mon Sep 17 00:00:00 2001 From: Jutho Date: Thu, 7 Nov 2024 11:07:56 +0100 Subject: [PATCH 8/9] keep RecursiveVec in tests for deprecations --- src/KrylovKit.jl | 4 +++- test/nestedtuple.jl | 9 +++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index 1f903c0..ca0f7bd 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -240,6 +240,8 @@ include("matrixfun/expintegrator.jl") # custom vector types include("innerproductvec.jl") -Base.@deprecate(RecursiveVec(args...), tuple(args...)) + +# deprecations +include("deprecated.jl") end diff --git a/test/nestedtuple.jl b/test/nestedtuple.jl index 7ff1fc9..ec3973a 100644 --- a/test/nestedtuple.jl +++ b/test/nestedtuple.jl @@ -1,15 +1,16 @@ +# TODO: Remove RecursiveVec and simply use tuple when RecursiveVec is removed. @testset "RecursiveVec - singular values full" begin @testset for T in (Float32, Float64, ComplexF32, ComplexF64) @testset for orth in (cgs2, mgs2, cgsr, mgsr) A = rand(T, (n, n)) v = rand(T, (n,)) - v2 = (v, zero(v)) + v2 = RecursiveVec(v, zero(v)) alg = Lanczos(; orth=orth, krylovdim=2 * n, maxiter=1, tol=tolerance(T)) D, V, info = eigsolve(v2, n, :LR, alg) do x x1, x2 = x y1 = A * x2 y2 = A' * x1 - return (y1, y2) + return RecursiveVec(y1, y2) end @test info.converged >= n S = D[1:n] @@ -27,14 +28,14 @@ end A = rand(T, (N, 2 * N)) v = rand(T, (N,)) w = rand(T, (2 * N,)) - v2 = (v, w) + v2 = RecursiveVec(v, w) alg = Lanczos(; orth=orth, krylovdim=n, maxiter=300, tol=tolerance(T)) n1 = div(n, 2) D, V, info = eigsolve(v2, n1, :LR, alg) do x x1, x2 = x y1 = A * x2 y2 = A' * x1 - return (y1, y2) + return RecursiveVec(y1, y2) end @test info.converged >= n1 S = D[1:n1] From c7fb66dadad6512d778a7fe8f6fb3348c47680e4 Mon Sep 17 00:00:00 2001 From: Jutho Date: Thu, 7 Nov 2024 11:11:12 +0100 Subject: [PATCH 9/9] actually add deprecated.jl --- src/deprecated.jl | 1 + 1 file changed, 1 insertion(+) create mode 100644 src/deprecated.jl diff --git a/src/deprecated.jl b/src/deprecated.jl new file mode 100644 index 0000000..96bb856 --- /dev/null +++ b/src/deprecated.jl @@ -0,0 +1 @@ +Base.@deprecate(RecursiveVec(args...), tuple(args...))