VectorInterface

.gitignore

@@ -2,3 +2,4 @@
 *.jl.*.cov
 *.jl.mem
 .DS_Store
+Manifest.toml

Project.toml

@@ -1,20 +1,22 @@
 name = "KrylovKit"
 uuid = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 authors = ["Jutho Haegeman"]
-version = "0.6.1"
+version = "0.7.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [compat]
 Aqua = "0.6, 0.7, 0.8"
 ChainRulesCore = "1"
 ChainRulesTestUtils = "1"
 FiniteDifferences = "0.12"
 GPUArraysCore = "0.1"
+VectorInterface = "0.4"
 LinearAlgebra = "1"
 Random = "1"
 Printf = "1"

docs/src/index.md

@@ -61,20 +61,8 @@ 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
-    following functions (with `v` and `w` two instances of this type and `α, β` scalars
-    (i.e. `Number`)):
-    *   `Base.:*(α, v)`: multiply `v` with a scalar `α`, which can be of a different scalar
-        type; in particular this method is used to create vectors similar to `v` but with a
-        different type of underlying scalars.
-    *   `Base.similar(v)`: a way to construct vectors which are exactly similar to `v`
-    *   `LinearAlgebra.mul!(w, v, α)`: out of place scalar multiplication; multiply
-        vector `v` with scalar `α` and store the result in `w`
-    *   `LinearAlgebra.rmul!(v, α)`: in-place scalar multiplication of `v` with `α`; in
-        particular with `α = false`, `v` is the corresponding zero vector
-    *   `LinearAlgebra.axpy!(α, v, w)`: store in `w` the result of `α*v + w`
-    *   `LinearAlgebra.axpby!(α, v, β, w)`: store in `w` the result of `α*v + β*w`
-    *   `LinearAlgebra.dot(v,w)`: compute the inner product of two vectors
-    *   `LinearAlgebra.norm(v)`: compute the 2-norm of a vector
+    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,
@@ -84,14 +72,14 @@ However, KrylovKit.jl distinguishes itself from the previous packages in the fol
     *   [`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. `dot`)
+    *   [`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.
 
 ## Current functionality
 
 The following algorithms are currently implemented
-*   `linsolve`: [`CG`](@ref), [`GMRES`](@ref)
+*   `linsolve`: [`CG`](@ref), [`GMRES`](@ref), [`BiCGStab`](@ref)
 *   `eigsolve`: a Krylov-Schur algorithm (i.e. with tick restarts) for extremal eigenvalues
     of normal (i.e. not generalized) eigenvalue problems, corresponding to
     [`Lanczos`](@ref) for real symmetric or complex hermitian linear maps, and to

docs/src/man/implementation.md

@@ -22,11 +22,11 @@ KrylovKit.orthonormalize
 
 The expansion coefficients of a general vector in terms of a given orthonormal basis can be obtained as
 ```@docs
-KrylovKit.project!
+KrylovKit.project!!
 ```
 whereas the inverse calculation is obtained as
 ```@docs
-KrylovKit.unproject!
+KrylovKit.unproject!!
 ```
 
 An orthonormal basis can be transformed using a rank-1 update using

src/KrylovKit.jl

@@ -20,15 +20,16 @@ The high level interface of KrylovKit is provided by the following functions:
     computes a linear combination of the `ϕⱼ` functions which generalize `ϕ₀(z) = exp(z)`.
 """
 module KrylovKit
-
+using VectorInterface
+using VectorInterface: add!!
 using LinearAlgebra
 using Printf
 using ChainRulesCore
 using GPUArraysCore
 const IndexRange = AbstractRange{Int}
 
 export linsolve, eigsolve, geneigsolve, svdsolve, schursolve, exponentiate, expintegrator
-export orthogonalize, orthogonalize!, orthonormalize, orthonormalize!
+export orthogonalize, orthogonalize!!, orthonormalize, orthonormalize!!
 export basis, rayleighquotient, residual, normres, rayleighextension
 export initialize, initialize!, expand!, shrink!
 export ClassicalGramSchmidt, ClassicalGramSchmidt2, ClassicalGramSchmidtIR

src/adrules/linsolve.jl

@@ -66,18 +66,20 @@
         ∂self = NoTangent()
         ∂x₀ = ZeroTangent()
         ∂algorithm = NoTangent()
-        (∂b, reverse_info) = linsolve(fᴴ, x̄, (zero(a₀) * zero(a₁)) * x̄, algorithm,
-                                      conj(a₀), conj(a₁))
+        T = VectorInterface.promote_scale(VectorInterface.promote_scale(x̄, a₀),
+                                          scalartype(a₁))
+        ∂b, reverse_info = linsolve(fᴴ, x̄, zerovector(x̄, T), algorithm, conj(a₀),
+                                    conj(a₁))
         if reverse_info.converged == 0
             @warn "Linear problem for reverse rule did not converge." reverse_info
         end
-        ∂f = @thunk(f_pullback(-conj(a₁) * ∂b)[1])
-        ∂a₀ = @thunk(-dot(x, ∂b))
+        ∂f = @thunk(f_pullback(scale(∂b, -conj(a₁)))[1])
+        ∂a₀ = @thunk(-inner(x, ∂b))
         # ∂a₁ = @thunk(-dot(f(x), ∂b))
         if a₀ == zero(a₀) && a₁ == one(a₁)
-            ∂a₁ = @thunk(-dot(b, ∂b))
+            ∂a₁ = @thunk(-inner(b, ∂b))
         else
-            ∂a₁ = @thunk(-dot((b - a₀ * x) / a₁, ∂b))
+            ∂a₁ = @thunk(-inner(scale!!(add(b, x, -a₀), inv(a₁)), ∂b))
         end
         return ∂self, ∂f, ∂b, ∂x₀, ∂algorithm, ∂a₀, ∂a₁
     end
@@ -91,17 +93,17 @@
     (x, info) = linsolve(A, b, x₀, algorithm, a₀, a₁)
 
     if Δb isa ChainRulesCore.AbstractZero
-        rhs = zero(b)
+        rhs = zerovector(b)
     else
-        rhs = (1 - Δa₁) * Δb
+        rhs = scale(Δb, (1 - Δa₁))
     end
     if !iszero(Δa₀)
-        rhs = axpy!(-Δa₀, x, rhs)
+        rhs = add!!(rhs, x, -Δa₀)
     end
     if !iszero(ΔA)
         rhs = mul!(rhs, ΔA, x, -a₁, true)
     end
-    (Δx, forward_info) = linsolve(A, rhs, zero(rhs), algorithm, a₀, a₁)
+    (Δx, forward_info) = linsolve(A, rhs, zerovector(rhs), algorithm, a₀, a₁)
     if info.converged > 0 && forward_info.converged == 0
         @warn "The tangent linear problem did not converge, whereas the primal linear problem did."
     end
@@ -121,17 +123,17 @@
     (x, info) = linsolve(f, b, x₀, algorithm, a₀, a₁)
 
     if Δb isa AbstractZero
-        rhs = false * b
+        rhs = zerovector(b)
     else
-        rhs = (1 - Δa₁) * Δb
+        rhs = scale(Δb, (1 - Δa₁))
     end
     if !iszero(Δa₀)
-        rhs = axpy!(-Δa₀, x, rhs)
+        rhs = add!!(rhs, x, -Δa₀)
     end
     if !(Δf isa AbstractZero)
-        rhs = axpy!(-a₁, frule_via_ad(config, (Δf, ZeroTangent()), f, x), rhs)
+        rhs = add!!(rhs, frule_via_ad(config, (Δf, ZeroTangent()), f, x), -a₀)
     end
-    (Δx, forward_info) = linsolve(f, rhs, false * rhs, algorithm, a₀, a₁)
+    (Δx, forward_info) = linsolve(f, rhs, zerovector(rhs), algorithm, a₀, a₁)
     if info.converged > 0 && forward_info.converged == 0
         @warn "The tangent linear problem did not converge, whereas the primal linear problem did."
     end

src/apply.jl

@@ -4,7 +4,7 @@ apply(f, x) = f(x)
 function apply(operator, x, α₀, α₁)
     y = apply(operator, x)
     if α₀ != zero(α₀) || α₁ != one(α₁)
-        axpby!(α₀, x, α₁, y)
+        y = add!!(y, x, α₀, α₁)
     end
 
     return y

src/dense/givens.jl

@@ -29,9 +29,9 @@ end
 
 function _rmul!(b::OrthonormalBasis, G::Givens)
     q1, q2 = b[G.i1], b[G.i2]
-    q1old = mul!(similar(q1), q1, true)
-    q1 = axpby!(-conj(G.s), q2, G.c, q1)
-    q2 = axpby!(G.s, q1old, G.c, q2)
+    q1′ = add(q1, q2, -conj(G.s), G.c)
+    q2′ = add!!(q2, q1, G.s, G.c)
+    b[G.i1], b[G.i2] = q1′, q2′
     return b
 end
 

src/dense/reflector.jl

@@ -145,10 +145,10 @@ function LinearAlgebra.rmul!(b::OrthonormalBasis, H::Householder)
     r = H.r
     β = H.β
     iszero(β) && return b
-    w = similar(b[first(r)])
+    w = zerovector(b[first(r)])
     @inbounds begin
-        unproject!(w, b, v, 1, 0, r)
-        rank1update!(b, w, v, -β, 1, r)
+        w = unproject!!(w, b, v, 1, 0, r)
+        b = rank1update!(b, w, v, -β, 1, r)
     end
     return b
 end

src/eigsolve/arnoldi.jl

@@ -106,7 +106,7 @@ function schursolve(A, x₀, howmany::Int, which::Selector, alg::Arnoldi)
         [B * u for u in cols(U, 1:howmany)]
     end
     residuals = let r = residual(fact)
-        [last(u) * r for u in cols(U, 1:howmany)]
+        [scale(r, last(u)) for u in cols(U, 1:howmany)]
     end
     normresiduals = [normres(fact) * abs(last(u)) for u in cols(U, 1:howmany)]
 
@@ -145,7 +145,7 @@ function eigsolve(A, x₀, howmany::Int, which::Selector, alg::Arnoldi)
         [B * v for v in cols(V)]
     end
     residuals = let r = residual(fact)
-        [last(v) * r for v in cols(V)]
+        [scale(r, last(v)) for v in cols(V)]
     end
     normresiduals = [normres(fact) * abs(last(v)) for v in cols(V)]
 
@@ -271,7 +271,7 @@ function _schursolve(A, x₀, howmany::Int, which::Selector, alg::Arnoldi)
             B = basis(fact)
             basistransform!(B, view(U, :, 1:keep))
             r = residual(fact)
-            B[keep + 1] = rmul!(r, 1 / normres(fact))
+            B[keep + 1] = scale!!(r, 1 / normres(fact))
 
             # Shrink Arnoldi factorization
             fact = shrink!(fact, keep)