TRICG & TRIMR: support warm-start with ldiv=true

This addresses some potential inconsistencies in `tricg` and `trimr` in
the case `ldiv=true`:

- docs: fixes the statement about the relationship between `E`/`F` and
  `M`/`N`
- docs: describes the residual-norm in terms of `E`/`F` rather than
  `M`/`N` (only the former are unambiguous, given the problem statement)
- adds support for warm-starting when `ldiv=true` for arbitrary `M`/`N`

Author: timholy

src/tricg.jl

@@ -33,22 +33,24 @@ Given a matrix `A` of dimension m × n, TriCG solves the Hermitian linear system
     [ τE    A ] [ x ] = [ b ]
     [  Aᴴ  νF ] [ y ]   [ c ],
 
-of size (n+m) × (n+m) where `τ` and `ν` are real numbers, `E` = `M⁻¹` ≻ 0 and `F` = `N⁻¹` ≻ 0.
-TriCG could breakdown if `τ = 0` or `ν = 0`.
-It's recommended to use TriMR in these cases.
+of size (n+m) × (n+m) where `τ` and `ν` are real numbers, `E` ≻ 0 and `F` ≻ 0.
+`E` and `F` are related to the preconditioners `M` and `N` as `E = M⁻¹` and `F =
+N⁻¹` when `ldiv = false` (the default), and `E = M` and `F = N` when `ldiv =
+true`.
 
 By default, TriCG solves Hermitian and quasi-definite linear systems with `τ = 1` and `ν = -1`.
+TriCG could breakdown if `τ = 0` or `ν = 0`. It's recommended to use TriMR in these cases.
 
 TriCG is based on the preconditioned orthogonal tridiagonalization process
 and its relation with the preconditioned block-Lanczos process.
 
 The matrix
 
-    [ M   0 ]
-    [ 0   N ]
+    [ E⁻¹   0 ]
+    [ 0   F⁻¹ ]
 
 indicates the weighted norm in which residuals are measured, here denoted `‖·‖`.
-It's the Euclidean norm when `M` and `N` are identity operators.
+When `M` and `N` are identity operators and `τ` and `ν` are ±1, it's the Euclidean norm.
 
 TriCG stops when `itmax` iterations are reached or when `‖rₖ‖ ≤ atol + ‖r₀‖ * rtol`.
 `atol` is an absolute tolerance and `rtol` is a relative tolerance.
@@ -70,7 +72,9 @@ For an in-place variant that reuses memory across solves, see [`tricg!`](@ref).
 * `x0`: a vector of length `m` that represents an initial guess of the solution `x`;
 * `y0`: a vector of length `n` that represents an initial guess of the solution `y`.
 
-Warm-starting is supported only when `M` and `N` are either `I` or the corresponding coefficient (`τ` or `ν`) is zero.
+Warm-starting is supported only when:
+* ldiv = true; or
+* `M` and `N` are either `I` or the corresponding coefficient (`τ` or `ν`) is zero.
 
 #### Keyword arguments
 
@@ -178,8 +182,10 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax
     snd  && (τ = -one(T) ; ν = -one(T))
 
     warm_start = workspace.warm_start
-    warm_start && (τ ≠ 0) && !MisI && error("Warm-start with preconditioners is not supported.")
-    warm_start && (ν ≠ 0) && !NisI && error("Warm-start with preconditioners is not supported.")
+    if !ldiv
+      warm_start && (τ ≠ 0) && !MisI && error("Warm-start with preconditioners is not supported.")
+      warm_start && (ν ≠ 0) && !NisI && error("Warm-start with preconditioners is not supported.")
+    end
 
     # Compute the adjoint of A
     Aᴴ = A'
@@ -212,15 +218,27 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax
     kfill!(M⁻¹vₖ₋₁, zero(FC))  # v₀ = 0
     kfill!(N⁻¹uₖ₋₁, zero(FC))  # u₀ = 0
 
-    # [ τI    A ] [ xₖ ] = [ b -  τΔx - AΔy ] = [ b₀ ]
-    # [  Aᴴ  νI ] [ yₖ ]   [ c - AᴴΔx - νΔy ]   [ c₀ ]
     if warm_start
-      kmul!(b₀, A, Δy)
-      (τ ≠ 0) && kaxpy!(m, τ, Δx, b₀)
-      kaxpby!(m, one(FC), b, -one(FC), b₀)
-      kmul!(c₀, Aᴴ, Δx)
-      (ν ≠ 0) && kaxpy!(n, ν, Δy, c₀)
-      kaxpby!(n, one(FC), c, -one(FC), c₀)
+      if ldiv
+        # [ τM    A ] [ xₖ ] = [ b -  τMΔx - AΔy ] = [ b₀ ]
+        # [  Aᴴ  νN ] [ yₖ ]   [ c - AᴴΔx - νNΔy ]   [ c₀ ]
+        kmul!(b₀, A, Δy)
+        (τ ≠ 0) && mul!(b₀, M, Δx, τ, one(FC))  # b₀ ← τMΔx + AΔy
+        kaxpby!(m, one(FC), b, -one(FC), b₀)
+        kmul!(c₀, Aᴴ, Δx)
+        (ν ≠ 0) && mul!(c₀, N, Δy, ν, one(FC))  # c₀ ← νNΔy + AᴴΔx
+        kaxpby!(n, one(FC), c, -one(FC), c₀)
+      else
+        # Only supported for M = I or τ = 0, and N = I or ν = 0
+        # [ τI    A ] [ xₖ ] = [ b -  τΔx - AΔy ] = [ b₀ ]
+        # [  Aᴴ  νI ] [ yₖ ]   [ c - AᴴΔx - νΔy ]   [ c₀ ]
+        kmul!(b₀, A, Δy)
+        (τ ≠ 0) && kaxpy!(m, τ, Δx, b₀)
+        kaxpby!(m, one(FC), b, -one(FC), b₀)
+        kmul!(c₀, Aᴴ, Δx)
+        (ν ≠ 0) && kaxpy!(n, ν, Δy, c₀)
+        kaxpby!(n, one(FC), c, -one(FC), c₀)
+      end
     end
 
     # β₁Ev₁ = b ↔ β₁v₁ = Mb

test/test_warm_start.jl

@@ -94,6 +94,14 @@ function test_warm_start(FC)
     @test stats.niter <= 1
     @test_throws "x0 should have size $(length(x30))" tricg(A3x2, b3, c2, [1.0], y20)
     @test_throws "y0 should have size $(length(y20))" tricg(A3x2, b3, c2, x30, [1.0])
+
+    N = Diagonal([0.8, 1.2])
+    x30, y20, _ = tricg(A3x2, b3, c2; N, ldiv=true)
+    x3, y2, stats = tricg(A3x2, b3, c2, x30, y20; N, ldiv=true)
+    @test x3 ≈ x30 atol=1e-12
+    @test y2 ≈ y20 atol=1e-12
+    @test stats.niter <= 1
+
     workspace = TricgWorkspace(A, b)
     krylov_solve!(workspace, A, b, b, x0, y0)
     r = [b - workspace.x - A * workspace.y; b - A' * workspace.x + workspace.y]