Clean UsymlqrSolver

Author: amontoison

docs/src/storage.md

@@ -95,7 +95,7 @@ Each table summarizes the storage requirements of Krylov methods recommended to
 
 | Methods  | [`TriCG`](@ref tricg) | [`TriMR`](@ref trimr) | [`USYMLQR`](@ref usymlqr) |
 |:--------:|:---------------------:|:---------------------:|:-------------------------:|
-| Storage  | $6n + 6m$             | $8n + 8m$             | $Xn + Ym$                 |
+| Storage  | $6n + 6m$             | $8n + 8m$             | $7n + 6m$                 |
 
 #### Generalized saddle-point and non-Hermitian partitioned systems
 

src/krylov_solvers.jl

@@ -1920,35 +1920,49 @@ can be used to initialize this workspace.
 mutable struct UsymlqrSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   m          :: Int
   n          :: Int
+  r          :: S
   x          :: S
   y          :: S
-  M⁻¹vₖ₋₁    :: S
-  M⁻¹vₖ      :: S
-  N⁻¹uₖ₋₁    :: S
-  N⁻¹uₖ      :: S
+  z          :: S
+  M⁻¹uₖ₋₁    :: S
+  M⁻¹uₖ      :: S
+  N⁻¹vₖ₋₁    :: S
+  N⁻¹vₖ      :: S
   p          :: S
   q          :: S
-  vₖ         :: S
+  d̅          :: S
+  wₖ₋₂       :: S
+  wₖ₋₁       :: S
+  Δx         :: S
+  Δy         :: S
   uₖ         :: S
+  vₖ         :: S
   warm_start :: Bool
   stats      :: SimpleStats{T}
 end
 
 function UsymlqrSolver(m, n, S)
   FC      = eltype(S)
   T       = real(FC)
-  x       = S(undef, m)
-  y       = S(undef, n)
-  M⁻¹vₖ₋₁ = S(undef, m)
-  M⁻¹vₖ   = S(undef, m)
-  N⁻¹uₖ₋₁ = S(undef, n)
-  N⁻¹uₖ   = S(undef, n)
-  p       = S(undef, n)
+  r       = S(undef, m)
+  x       = S(undef, n)
+  y       = S(undef, m)
+  z       = S(undef, n)
+  M⁻¹uₖ₋₁ = S(undef, m)
+  M⁻¹uₖ   = S(undef, m)
+  N⁻¹vₖ₋₁ = S(undef, n)
+  N⁻¹vₖ   = S(undef, n)
   q       = S(undef, m)
-  vₖ      = S(undef, 0)
+  p       = S(undef, n)
+  d̅       = S(undef, m)
+  wₖ₋₂    = S(undef, n)
+  wₖ₋₁    = S(undef, n)
+  Δx      = S(undef, 0)
+  Δy      = S(undef, 0)
   uₖ      = S(undef, 0)
+  vₖ      = S(undef, 0)
   stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown")
-  solver = UsymlqrSolver{T,FC,S}(m, n, x, y, M⁻¹vₖ₋₁, M⁻¹vₖ, N⁻¹uₖ₋₁, N⁻¹uₖ, p, q, vₖ, uₖ, false, stats)
+  solver = UsymlqrSolver{T,FC,S}(m, n, r, x, y, z, M⁻¹uₖ₋₁, M⁻¹uₖ, N⁻¹vₖ₋₁, N⁻¹vₖ, q, p, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, uₖ, vₖ, false, stats)
   return solver
 end
 

src/minres_qlp.jl

@@ -279,8 +279,8 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
       # [sₖ -cₖ] [βₖ₊₁ ]   [0 ]
       (cₖ, sₖ, λₖ) = sym_givens(λbarₖ, βₖ₊₁)
 
-      # Compute [   zₖ  ] = (Qₖ)ᴴβ₁e₁
-      #         [ζbarₖ₊₁]
+      # Compute z̅ₖ₊₁ = [   zₖ  ] = (Qₖ)ᴴβ₁e₁
+      #                [ζbarₖ₊₁]
       #
       # [cₖ  sₖ] [ζbarₖ] = [   ζₖ  ]
       # [sₖ -cₖ] [  0  ]   [ζbarₖ₊₁]

src/trilqr.jl

@@ -201,8 +201,10 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
       mul!(q, A , uₖ)  # Forms vₖ₊₁ : q ← Auₖ
       mul!(p, Aᴴ, vₖ)  # Forms uₖ₊₁ : p ← Aᴴvₖ
 
-      @kaxpy!(m, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
-      @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - βₖ * uₖ₋₁
+      if iter ≥ 2
+        @kaxpy!(m, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
+        @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - βₖ * uₖ₋₁
+      end
 
       αₖ = @kdot(m, vₖ, q)      # αₖ = ⟨vₖ,q⟩
 

src/usymlq.jl

@@ -203,8 +203,10 @@ kwargs_usymlq = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
 
       αₖ = @kdot(m, vₖ, q)      # αₖ = ⟨vₖ,q⟩
 
-      @kaxpy!(m, -     αₖ , vₖ, q)    # q ← q - αₖ * vₖ
-      @kaxpy!(n, -conj(αₖ), uₖ, p)    # p ← p - ᾱₖ * uₖ
+      if iter ≥ 2
+        @kaxpy!(m, -     αₖ , vₖ, q)    # q ← q - αₖ * vₖ
+        @kaxpy!(n, -conj(αₖ), uₖ, p)    # p ← p - ᾱₖ * uₖ
+      end
 
       βₖ₊₁ = @knrm2(m, q)       # βₖ₊₁ = ‖q‖
       γₖ₊₁ = @knrm2(n, p)       # γₖ₊₁ = ‖p‖
@@ -281,7 +283,7 @@ kwargs_usymlq = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
         @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅)
       end
 
-      # Compute uₖ₊₁ and uₖ₊₁.
+      # Compute vₖ₊₁ and uₖ₊₁.
       @kcopy!(m, vₖ₋₁, vₖ)  # vₖ₋₁ ← vₖ
       @kcopy!(n, uₖ₋₁, uₖ)  # uₖ₋₁ ← uₖ