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: 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/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/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..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) @@ -87,6 +82,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), @@ -120,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]) diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index 3aafc07..ca0f7bd 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} @@ -238,7 +239,9 @@ include("matrixfun/exponentiate.jl") include("matrixfun/expintegrator.jl") # custom vector types -include("recursivevec.jl") include("innerproductvec.jl") +# deprecations +include("deprecated.jl") + end 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...)) diff --git a/src/eigsolve/eigsolve.jl b/src/eigsolve/eigsolve.jl index fb1adf1..d687f8e 100644 --- a/src/eigsolve/eigsolve.jl +++ b/src/eigsolve/eigsolve.jl @@ -182,7 +182,8 @@ function eigsolve(A::AbstractMatrix, 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/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/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) diff --git a/test/recursivevec.jl b/test/nestedtuple.jl similarity index 95% rename from test/recursivevec.jl rename to test/nestedtuple.jl index 08d5a4c..ec3973a 100644 --- a/test/recursivevec.jl +++ b/test/nestedtuple.jl @@ -1,3 +1,4 @@ +# 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) 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")