expv error estimate: hermitian only

src/krylov_phiv.jl

@@ -43,27 +43,27 @@ function expv(t, A, b; mode = :happy_breakdown, kwargs...)
     end
 end
 function _expv_hb(t::Tt, A, b;
-    expmethod = ExpMethodHigham2005Base(),
-    cache = nothing, kwargs_arnoldi...) where {Tt}
+        expmethod = ExpMethodHigham2005Base(),
+        cache = nothing, kwargs_arnoldi...) where {Tt}
     # Happy-breakdown mode: first construct Krylov subspace then expv!
     Ks = arnoldi(A, b; kwargs_arnoldi...)
     w = similar(b, promote_type(Tt, eltype(A), eltype(b)))
     expv!(w, t, Ks; cache = cache, expmethod = expmethod)
 end
 function _expv_ee(t::Tt, A, b; m = min(30, size(A, 1)), tol = 1e-7, rtol = √(tol),
-    ishermitian::Bool = LinearAlgebra.ishermitian(A),
-    expmethod = ExpMethodHigham2005Base()) where {Tt}
+        ishermitian::Bool = LinearAlgebra.ishermitian(A),
+        expmethod = ExpMethodHigham2005Base()) where {Tt}
     # Error-estimate mode: construction of Krylov subspace and expv! at the same time
     n = size(A, 1)
     T = promote_type(typeof(t), eltype(A), eltype(b))
     U = ishermitian ? real(T) : T
     Ks = KrylovSubspace{T, U}(n, m)
     w = similar(b, promote_type(Tt, eltype(A), eltype(b)))
     expv!(w, t, A, b, Ks, get_subspace_cache(Ks); atol = tol, rtol = rtol,
-        expmethod = expmethod)
+        ishermitian = ishermitian, expmethod = expmethod)
 end
 function expv(t::Tt, Ks::KrylovSubspace{T, U}; expmethod = ExpMethodHigham2005(),
-    kwargs...) where {Tt, T, U}
+        kwargs...) where {Tt, T, U}
     n = size(getV(Ks), 1)
     w = Vector{promote_type(Tt, T)}(undef, n)
     expv!(w, t, Ks; kwargs...)
@@ -74,7 +74,7 @@ end
 Non-allocating version of `expv` that uses precomputed Krylov subspace `Ks`.
 """
 function expv!(w::AbstractVector{Tw}, t::Real, Ks::KrylovSubspace{T, U};
-    cache = nothing, expmethod = ExpMethodHigham2005Base()) where {Tw, T, U}
+        cache = nothing, expmethod = ExpMethodHigham2005Base()) where {Tw, T, U}
     m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
     @assert length(w)==size(V, 1) "Dimension mismatch"
     if cache == nothing
@@ -102,7 +102,7 @@ end
 #       or Tw can be Float64, while t is ComplexF32 and T is Float64
 #       thus they can not share the same TypeVar.
 function expv!(w::AbstractVector{Complex{Tw}}, t::Complex{Tt}, Ks::KrylovSubspace{T, U};
-    cache = nothing, expmethod = ExpMethodHigham2005Base()) where {Tw, Tt, T, U}
+        cache = nothing, expmethod = ExpMethodHigham2005Base()) where {Tw, Tt, T, U}
     m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
     @assert length(w)==size(V, 1) "Dimension mismatch"
     if cache === nothing
@@ -130,9 +130,9 @@ function expv!(w::AbstractVector{Complex{Tw}}, t::Complex{Tt}, Ks::KrylovSubspac
 end
 
 function ExponentialUtilities.expv!(w::GPUArraysCore.AbstractGPUVector{Tw},
-    t::Real, Ks::KrylovSubspace{T, U};
-    cache = nothing,
-    expmethod = ExpMethodHigham2005Base()) where {Tw, T, U}
+        t::Real, Ks::KrylovSubspace{T, U};
+        cache = nothing,
+        expmethod = ExpMethodHigham2005Base()) where {Tw, T, U}
     m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
     @assert length(w)==size(V, 1) "Dimension mismatch"
     if cache === nothing
@@ -238,7 +238,7 @@ Compute the matrix-phi-vector products using a pre-constructed Krylov subspace.
     Formula (10).
 """
 function phiv(t, A, b, k; cache = nothing, correct = false, errest = false,
-    kwargs_arnoldi...)
+        kwargs_arnoldi...)
     Ks = arnoldi(A, b; kwargs_arnoldi...)
     w = Matrix{eltype(b)}(undef, length(b), k + 1)
     phiv!(w, t, Ks, k; cache = cache, correct = correct, errest = errest)
@@ -254,8 +254,8 @@ end
 Non-allocating version of 'phiv' that uses precomputed Krylov subspace `Ks`.
 """
 function phiv!(w::AbstractMatrix, t::Number, Ks::KrylovSubspace{T, U}, k::Integer;
-    cache = nothing, correct = false,
-    errest = false) where {T <: Number, U <: Number}
+        cache = nothing, correct = false,
+        errest = false) where {T <: Number, U <: Number}
     m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
     @assert size(w, 1)==size(V, 1) "Dimension mismatch"
     @assert size(w, 2)==k + 1 "Dimension mismatch"

src/krylov_phiv_error_estimate.jl

@@ -25,7 +25,7 @@ subspace matrix `T` with `α` on the diagonal and `β` on the
 super-/subdiagonal, diagonalizing via `stegr!`.
 """
 function expT!(α::AbstractVector{R}, β::AbstractVector{R}, t::Number,
-    cache::StegrCache{T, R}) where {T, R <: Real}
+        cache::StegrCache{T, R}) where {T, R <: Real}
     LAPACK.stegr!(α, β, cache.sw)
     sel = 1:length(α)
     @inbounds for i in sel
@@ -55,12 +55,19 @@ exact type decides which algorithm is used to compute the subspace
 exponential.
 """
 function expv!(w::AbstractVector{T}, t::Number, A, b::AbstractVector{T},
-    Ks::KrylovSubspace{T, B, B}, cache::HSC;
-    atol::B = 1.0e-8, rtol::B = 1.0e-4,
-    m = min(Ks.maxiter, size(A, 1)),
-    verbose::Bool = false,
-    expmethod = ExpMethodHigham2005Base()) where {B, T <: Number,
-    HSC <: HermitianSubspaceCache}
+        Ks::KrylovSubspace{T, B, B}, cache::HSC;
+        atol::B = 1.0e-8, rtol::B = 1.0e-4,
+        m = min(Ks.maxiter, size(A, 1)),
+        ishermitian::Bool = LinearAlgebra.ishermitian(A),
+        verbose::Bool = false,
+        expmethod = ExpMethodHigham2005Base()) where {B, T <: Number,
+        HSC <: HermitianSubspaceCache}
+    # TODO: this only implements the Lanczos algorithm for Hermitian matrices
+    # ks.H is tridiagonal, required for the expT! function above to call stegr!()
+    if !ishermitian
+        error("Error estimation not yet available for non-Hermitian matrices.")
+    end
+
     if m > Ks.maxiter
         resize!(Ks, m)
     else
@@ -77,7 +84,7 @@ function expv!(w::AbstractVector{T}, t::Number, A, b::AbstractVector{T},
     @. V[:, 1] = b / Ks.beta
 
     ε = atol + rtol * Ks.beta
-    verbose && @printf("Initial norm: β₀ %e, stopping threshold: %e\n", Ks.beta, ε)
+    verbose && @printf("Initial norm: β₀ %e, stopping threshold: %e\n", Ks.beta,ε)
 
     α = @diagview(H)
     β = @diagview(H, -1)
@@ -92,7 +99,7 @@ function expv!(w::AbstractVector{T}, t::Number, A, b::AbstractVector{T},
         #   Saad, Y. (1992). Analysis of some Krylov subspace
         #   approximations. SIAM Journal on Numerical Analysis.
         σ = β[j] * Ks.beta * abs(cache.v[j])
-        verbose && @printf("iter %d, α[%d] %e, β[%d] %e, σ %e\n", j, j, α[j], j, β[j], σ)
+        verbose && @printf("iter %d, α[%d] %e, β[%d] %e, σ %e\n", j, j, α[j], j, β[j],σ)
         if σ < ε
             Ks.m = j
             break