Add in-place Krylov methods

Author: amontoison

docs/src/api.md

@@ -11,8 +11,36 @@ Krylov.AdjointStats
 ## Solver Types
 
 ```@docs
-Krylov.KrylovSolver
-Krylov.MinresSolver
+KrylovSolver
+MinresSolver
+CgSolver
+CrSolver
+SymmlqSolver
+CgLanczosSolver
+CgLanczosShiftSolver
+MinresQlpSolver
+DqgmresSolver
+DiomSolver
+UsymlqSolver
+UsymqrSolver
+TricgSolver
+TrimrSolver
+TrilqrSolver
+CgsSolver
+BicgstabSolver
+BilqSolver
+QmrSolver
+BilqrSolver
+CglsSolver
+CrlsSolver
+CgneSolver
+CrmrSolver
+LslqSolver
+LsqrSolver
+LsmrSolver
+LnlqSolver
+CraigSolver
+CraigmrSolver
 ```
 
 ## Utilities
@@ -22,6 +50,7 @@ Krylov.roots_quadratic
 Krylov.sym_givens
 Krylov.to_boundary
 Krylov.vec2str
+Krylov.ktypeof
 Krylov.kzeros
 Krylov.kones
 ```

src/bicgstab.jl

@@ -13,7 +13,7 @@
 # Alexis Montoison, <alexis.montoison@polymtl.ca>
 # Montréal, October 2020.
 
-export bicgstab
+export bicgstab, bicgstab!
 
 """
     (x, stats) = bicgstab(A, b::AbstractVector{T}; c::AbstractVector{T}=b,
@@ -41,9 +41,14 @@ This implementation allows a left preconditioner `M` and a right preconditioner
 * H. A. van der Vorst, *Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems*, SIAM Journal on Scientific and Statistical Computing, 13(2), pp. 631--644, 1992.
 * G. L.G. Sleijpen and D. R. Fokkema, *BiCGstab(ℓ) for linear equations involving unsymmetric matrices with complex spectrum*, Electronic Transactions on Numerical Analysis, 1, pp. 11--32, 1993.
 """
-function bicgstab(A, b :: AbstractVector{T}; c :: AbstractVector{T}=b,
-                  M=opEye(), N=opEye(), atol :: T=√eps(T), rtol :: T=√eps(T),
-                  itmax :: Int=0, verbose :: Int=0, history :: Bool=false) where T <: AbstractFloat
+function bicgstab(A, b :: AbstractVector{T}; kwargs...) where T <: AbstractFloat
+  solver = BicgstabSolver(A, b)
+  bicgstab!(solver, A, b; kwargs...)
+end
+
+function bicgstab!(solver :: BicgstabSolver{T,S}, A, b :: AbstractVector{T}; c :: AbstractVector{T}=b,
+                   M=opEye(), N=opEye(), atol :: T=√eps(T), rtol :: T=√eps(T),
+                   itmax :: Int=0, verbose :: Int=0, history :: Bool=false) where {S, T <: AbstractFloat}
 
   n, m = size(A)
   m == n || error("System must be square")
@@ -56,18 +61,19 @@ function bicgstab(A, b :: AbstractVector{T}; c :: AbstractVector{T}=b,
 
   # Check type consistency
   eltype(A) == T || error("eltype(A) ≠ $T")
+  ktypeof(b) == S || error("ktypeof(b) ≠ $S")
+  ktypeof(c) == S || error("ktypeof(c) ≠ $S")
   MisI || (eltype(M) == T) || error("eltype(M) ≠ $T")
   NisI || (eltype(N) == T) || error("eltype(N) ≠ $T")
 
-  # Determine the storage type of b
-  S = typeof(b)
-
   # Set up workspace.
-  x = kzeros(S, n)  # x₀
-  s = kzeros(S, n)  # s₀
-  v = kzeros(S, n)  # v₀
-  r = copy(M * b)   # r₀
-  p = copy(r)       # p₁
+  x, s, v, r, p = solver.x, solver.s, solver.v, solver.r, solver.p
+
+  x .= zero(T) # x₀
+  s .= zero(T) # s₀
+  v .= zero(T) # v₀
+  r .= (M * b) # r₀
+  p .= r       # p₁
 
   α = one(T) # α₀
   ω = one(T) # ω₀

src/bilq.jl

@@ -10,7 +10,7 @@
 # Alexis Montoison, <alexis.montoison@polymtl.ca>
 # Montreal, February 2019.
 
-export bilq
+export bilq, bilq!
 
 """
     (x, stats) = bilq(A, b::AbstractVector{T}; c::AbstractVector{T}=b,
@@ -29,9 +29,14 @@ when it exists. The transfer is based on the residual norm.
 
 * A. Montoison and D. Orban, *BiLQ: An Iterative Method for Nonsymmetric Linear Systems with a Quasi-Minimum Error Property*, SIAM Journal on Matrix Analysis and Applications, 41(3), pp. 1145--1166, 2020.
 """
-function bilq(A, b :: AbstractVector{T}; c :: AbstractVector{T}=b,
-              atol :: T=√eps(T), rtol :: T=√eps(T), transfer_to_bicg :: Bool=true,
-              itmax :: Int=0, verbose :: Int=0, history :: Bool=false) where T <: AbstractFloat
+function bilq(A, b :: AbstractVector{T}; kwargs...) where T <: AbstractFloat
+  solver = BilqSolver(A, b)
+  bilq!(solver, A, b; kwargs...)
+end
+
+function bilq!(solver :: BilqSolver{T,S}, A, b :: AbstractVector{T}; c :: AbstractVector{T}=b,
+               atol :: T=√eps(T), rtol :: T=√eps(T), transfer_to_bicg :: Bool=true,
+               itmax :: Int=0, verbose :: Int=0, history :: Bool=false) where {S, T <: AbstractFloat}
 
   n, m = size(A)
   m == n || error("System must be square")
@@ -40,15 +45,17 @@ function bilq(A, b :: AbstractVector{T}; c :: AbstractVector{T}=b,
 
   # Check type consistency
   eltype(A) == T || error("eltype(A) ≠ $T")
+  ktypeof(b) == S || error("ktypeof(b) ≠ $S")
+  ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
   Aᵀ = A'
 
-  # Determine the storage type of b
-  S = typeof(b)
+  # Set up workspace.
+  uₖ₋₁, uₖ, vₖ₋₁, vₖ, x, d̅ = solver.uₖ₋₁, solver.uₖ, solver.vₖ₋₁, solver.vₖ, solver.x, solver.d̅
 
   # Initial solution x₀ and residual norm ‖r₀‖.
-  x = kzeros(S, n)
+  x .= zero(T)
   bNorm = @knrm2(n, b)  # ‖r₀‖
   bNorm == 0 && return (x, SimpleStats(true, false, [bNorm], T[], "x = 0 is a zero-residual solution"))
 
@@ -64,16 +71,15 @@ function bilq(A, b :: AbstractVector{T}; c :: AbstractVector{T}=b,
   bᵗc = @kdot(n, b, c)  # ⟨b,c⟩
   bᵗc == 0 && return (x, SimpleStats(false, false, [bNorm], T[], "Breakdown bᵀc = 0"))
 
-  # Set up workspace.
   βₖ = √(abs(bᵗc))           # β₁γ₁ = bᵀc
   γₖ = bᵗc / βₖ              # β₁γ₁ = bᵀc
-  vₖ₋₁ = kzeros(S, n)        # v₀ = 0
-  uₖ₋₁ = kzeros(S, n)        # u₀ = 0
-  vₖ = b / βₖ                # v₁ = b / β₁
-  uₖ = c / γₖ                # u₁ = c / γ₁
+  vₖ₋₁ .= zero(T)            # v₀ = 0
+  uₖ₋₁ .= zero(T)            # u₀ = 0
+  vₖ .= b ./ βₖ              # v₁ = b / β₁
+  uₖ .= c ./ γₖ              # u₁ = c / γ₁
   cₖ₋₁ = cₖ = -one(T)        # Givens cosines used for the LQ factorization of Tₖ
   sₖ₋₁ = sₖ = zero(T)        # Givens sines used for the LQ factorization of Tₖ
-  d̅ = kzeros(S, n)           # Last column of D̅ₖ = Vₖ(Qₖ)ᵀ
+  d̅ .= zero(T)               # Last column of D̅ₖ = Vₖ(Qₖ)ᵀ
   ζₖ₋₁ = ζbarₖ = zero(T)     # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
   ζₖ₋₂ = ηₖ = zero(T)        # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
   δbarₖ₋₁ = δbarₖ = zero(T)  # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations

src/bilqr.jl

@@ -10,7 +10,7 @@
 # Alexis Montoison, <alexis.montoison@polymtl.ca>
 # Montreal, July 2019.
 
-export bilqr
+export bilqr, bilqr!
 
 """
     (x, t, stats) = bilqr(A, b::AbstractVector{T}, c::AbstractVector{T};
@@ -32,9 +32,14 @@ BiCG point, when it exists. The transfer is based on the residual norm.
 
 * A. Montoison and D. Orban, *BiLQ: An Iterative Method for Nonsymmetric Linear Systems with a Quasi-Minimum Error Property*, SIAM Journal on Matrix Analysis and Applications, 41(3), pp. 1145--1166, 2020.
 """
-function bilqr(A, b :: AbstractVector{T}, c :: AbstractVector{T};
-               atol :: T=√eps(T), rtol :: T=√eps(T), transfer_to_bicg :: Bool=true,
-               itmax :: Int=0, verbose :: Int=0, history :: Bool=false) where T <: AbstractFloat
+function bilqr(A, b :: AbstractVector{T}, c :: AbstractVector{T}; kwargs...) where T <: AbstractFloat
+  solver = BilqrSolver(A, b)
+  bilqr!(solver, A, b, c; kwargs...)
+end
+
+function bilqr!(solver :: BilqrSolver{T,S}, A, b :: AbstractVector{T}, c :: AbstractVector{T};
+                atol :: T=√eps(T), rtol :: T=√eps(T), transfer_to_bicg :: Bool=true,
+                itmax :: Int=0, verbose :: Int=0, history :: Bool=false) where {S, T <: AbstractFloat}
 
   n, m = size(A)
   m == n || error("Systems must be square")
@@ -44,19 +49,21 @@ function bilqr(A, b :: AbstractVector{T}, c :: AbstractVector{T};
 
   # Check type consistency
   eltype(A) == T || error("eltype(A) ≠ $T")
+  ktypeof(b) == S || error("ktypeof(b) ≠ $S")
+  ktypeof(c) == S || error("ktypeof(c) ≠ $S")
 
   # Compute the adjoint of A
   Aᵀ = A'
 
-  # Determine the storage type of b
-  S = typeof(b)
+  # Set up workspace.
+  uₖ₋₁, uₖ, vₖ₋₁, vₖ, x, t, d̅, wₖ₋₃, wₖ₋₂ = solver.uₖ₋₁, solver.uₖ, solver.vₖ₋₁, solver.vₖ, solver.x, solver.t, solver.d̅, solver.wₖ₋₃, solver.wₖ₋₂
 
   # Initial solution x₀ and residual norm ‖r₀‖ = ‖b - Ax₀‖.
-  x = kzeros(S, n)      # x₀
+  x .= zero(T)          # x₀
   bNorm = @knrm2(n, b)  # rNorm = ‖r₀‖
 
   # Initial solution t₀ and residual norm ‖s₀‖ = ‖c - Aᵀt₀‖.
-  t = kzeros(S, n)      # t₀
+  t .= zero(T)          # t₀
   cNorm = @knrm2(n, c)  # sNorm = ‖s₀‖
 
   iter = 0
@@ -76,21 +83,21 @@ function bilqr(A, b :: AbstractVector{T}, c :: AbstractVector{T};
   # Set up workspace.
   βₖ = √(abs(bᵗc))           # β₁γ₁ = bᵀc
   γₖ = bᵗc / βₖ              # β₁γ₁ = bᵀc
-  vₖ₋₁ = kzeros(S, n)        # v₀ = 0
-  uₖ₋₁ = kzeros(S, n)        # u₀ = 0
-  vₖ = b / βₖ                # v₁ = b / β₁
-  uₖ = c / γₖ                # u₁ = c / γ₁
+  vₖ₋₁ .= zero(T)            # v₀ = 0
+  uₖ₋₁ .= zero(T)            # u₀ = 0
+  vₖ .= b ./ βₖ              # v₁ = b / β₁
+  uₖ .= c ./ γₖ              # u₁ = c / γ₁
   cₖ₋₁ = cₖ = -one(T)        # Givens cosines used for the LQ factorization of Tₖ
   sₖ₋₁ = sₖ = zero(T)        # Givens sines used for the LQ factorization of Tₖ
-  d̅ = kzeros(S, n)           # Last column of D̅ₖ = Vₖ(Qₖ)ᵀ
+  d̅ .= zero(T)               # Last column of D̅ₖ = Vₖ(Qₖ)ᵀ
   ζₖ₋₁ = ζbarₖ = zero(T)     # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
   ζₖ₋₂ = ηₖ = zero(T)        # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
   δbarₖ₋₁ = δbarₖ = zero(T)  # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations
   ψbarₖ₋₁ = ψₖ₋₁ = zero(T)   # ψₖ₋₁ and ψbarₖ are the last components of h̅ₖ = Qₖγ₁e₁
   norm_vₖ = bNorm / βₖ       # ‖vₖ‖ is used for residual norm estimates
   ϵₖ₋₃ = λₖ₋₂ = zero(T)      # Components of Lₖ₋₁
-  wₖ₋₃ = kzeros(S, n)        # Column k-3 of Wₖ = Uₖ(Lₖ)⁻ᵀ
-  wₖ₋₂ = kzeros(S, n)        # Column k-2 of Wₖ = Uₖ(Lₖ)⁻ᵀ
+  wₖ₋₃ .= zero(T)            # Column k-3 of Wₖ = Uₖ(Lₖ)⁻ᵀ
+  wₖ₋₂ .= zero(T)            # Column k-2 of Wₖ = Uₖ(Lₖ)⁻ᵀ
   τₖ = zero(T)               # τₖ is used for the dual residual norm estimate
 
   # Stopping criterion.

src/cg.jl

@@ -13,7 +13,7 @@
 # Dominique Orban, <dominique.orban@gerad.ca>
 # Salt Lake City, UT, March 2015.
 
-export cg
+export cg, cg!
 
 
 """
@@ -37,10 +37,15 @@ with `n = length(b)`.
 
 * M. R. Hestenes and E. Stiefel, *Methods of conjugate gradients for solving linear systems*, Journal of Research of the National Bureau of Standards, 49(6), pp. 409--436, 1952.
 """
-function cg(A, b :: AbstractVector{T};
-            M=opEye(), atol :: T=√eps(T), rtol :: T=√eps(T),
-            itmax :: Int=0, radius :: T=zero(T), linesearch :: Bool=false,
-            verbose :: Int=0, history :: Bool=false) where T <: AbstractFloat
+function cg(A, b :: AbstractVector{T}; kwargs...) where T <: AbstractFloat
+  solver = CgSolver(A, b)
+  cg!(solver, A, b; kwargs...)
+end
+
+function cg!(solver :: CgSolver{T,S}, A, b :: AbstractVector{T};
+             M=opEye(), atol :: T=√eps(T), rtol :: T=√eps(T),
+             itmax :: Int=0, radius :: T=zero(T), linesearch :: Bool=false,
+             verbose :: Int=0, history :: Bool=false) where {S, T <: AbstractFloat}
 
   linesearch && (radius > 0) && error("`linesearch` set to `true` but trust-region radius > 0")
 
@@ -50,16 +55,16 @@ function cg(A, b :: AbstractVector{T};
 
   # Check type consistency
   eltype(A) == T || error("eltype(A) ≠ $T")
+  ktypeof(b) == S || error("ktypeof(b) ≠ $S")
   isa(M, opEye) || (eltype(M) == T) || error("eltype(M) ≠ $T")
 
-  # Determine the storage type of b
-  S = typeof(b)
+  # Set up workspace.
+  x, r, p = solver.x, solver.r, solver.p
 
-  # Initial state.
-  x = kzeros(S, n)
-  r = copy(b)
+  x .= zero(T)
+  r .= b
   z = M * r
-  p = copy(z)
+  p .= z
   γ = @kdot(n, r, z)
   γ == 0 && return x, SimpleStats(true, false, [zero(T)], T[], "x = 0 is a zero-residual solution")