deprecate RecursiveVec, improve x0 initialisation for GPU

.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

.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:

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"]

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/)

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
 ```

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

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])

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

src/deprecated.jl

@@ -0,0 +1 @@
+Base.@deprecate(RecursiveVec(args...), tuple(args...))

src/eigsolve/eigsolve.jl

@@ -182,7 +182,8 @@
                   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;

src/eigsolve/geneigsolve.jl

@@ -150,21 +150,24 @@
     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,

src/eigsolve/svdsolve.jl

@@ -121,7 +121,8 @@
                   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...)

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