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

src/tricg.jl

@@ -33,22 +33,26 @@ 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`. `tricg` is applicable when `E` and `F` are represented through their
+inverses or have an efficient `ldiv!` implementation.
 
 By default, TriCG solves Hermitian and quasi-definite linear systems with `τ = 1` and `ν = -1`.
+TriCG could breakdown if the choice of `(τ, ν)` leads to matrices that are not SQD. An example is
+saddle-point systems where`τ = 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, 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 +74,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`, in which case `M` and `N` must support `mul!` as well as `ldiv!` (PDMats.jl may be helpful); or
+* `M` and `N` are either `I` or the corresponding coefficient (`τ` or `ν`) is zero.
 
 #### Keyword arguments
 
@@ -178,8 +184,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'
@@ -205,28 +213,58 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax
     kfill!(xₖ, zero(FC))
     kfill!(yₖ, zero(FC))
 
+    # Initialization, far enough to compute convergence tolerance
+    # Some of this may have to be repeated if we are warm-starting, but
+    # the convergence tolerance should be independent of whether we are warm-starting.
+    # β₁Ev₁ = b
+    kcopy!(m, M⁻¹vₖ, b)  # M⁻¹vₖ ← b
+    MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv)
+    βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ)  # β₁ = ‖v₁‖_E
+    # γ₁Fu₁ = c
+    kcopy!(n, N⁻¹uₖ, c)  # N⁻¹uₖ ← c
+    NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv)
+    γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ)  # γ₁ = ‖u₁‖_F
+    # Convergence tolerance
+    rNorm = sqrt(γₖ^2 + βₖ^2)
+    ε = atol + rtol * rNorm
+
     iter = 0
     itmax == 0 && (itmax = m+n)
 
     # Initialize preconditioned orthogonal tridiagonalization process.
     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
+      # Repeat the initialization with the updated right-hand side
+      kcopy!(m, M⁻¹vₖ, b₀)  # M⁻¹vₖ ← b₀
+      MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv)
+      βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ)  # β₁ = ‖v₁‖_E
+      kcopy!(n, N⁻¹uₖ, c₀)  # N⁻¹uₖ ← c₀
+      NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv)
+      γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ)  # γ₁ = ‖u₁‖_F
     end
 
-    # β₁Ev₁ = b ↔ β₁v₁ = Mb
-    kcopy!(m, M⁻¹vₖ, b₀)  # M⁻¹vₖ ← b₀
-    MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv)
-    βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ)  # β₁ = ‖v₁‖_E
     if βₖ ≠ 0
       kdiv!(m, M⁻¹vₖ, βₖ)
       MisI || kdiv!(m, vₖ, βₖ)
@@ -236,10 +274,6 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax
       MisI || kfill!(vₖ, zero(FC))
     end
 
-    # γ₁Fu₁ = c ↔ γ₁u₁ = Nc
-    kcopy!(n, N⁻¹uₖ, c₀)  # M⁻¹uₖ ← c₀
-    NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv)
-    γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ)  # γ₁ = ‖u₁‖_F
     if γₖ ≠ 0
       kdiv!(n, N⁻¹uₖ, γₖ)
       NisI || kdiv!(n, uₖ, γₖ)
@@ -258,7 +292,6 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax
     # Compute ‖r₀‖² = (γ₁)² + (β₁)²
     rNorm = sqrt(γₖ^2 + βₖ^2)
     history && push!(rNorms, rNorm)
-    ε = atol + rtol * rNorm
 
     (verbose > 0) && @printf(iostream, "%5s  %7s  %7s  %7s  %5s\n", "k", "‖rₖ‖", "βₖ₊₁", "γₖ₊₁", "timer")
     kdisplay(iter, verbose) && @printf(iostream, "%5d  %7.1e  %7.1e  %7.1e  %.2fs\n", iter, rNorm, βₖ, γₖ, start_time |> ktimer)

src/trimr.jl

@@ -69,7 +69,9 @@ For an in-place variant that reuses memory across solves, see [`trimr!`](@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`, in which case `M` and `N` must support `mul!` as well as `ldiv!` (PDMats.jl may be helpful); or
+* `M` and `N` are either `I` or the corresponding coefficient (`τ` or `ν`) is zero.
 
 #### Keyword arguments
 
@@ -183,8 +185,10 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, :
     sp   && (τ =  one(T) ; ν = zero(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'
@@ -218,21 +222,48 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, :
     kfill!(M⁻¹vₖ₋₁, zero(FC))  # v₀ = 0
     kfill!(N⁻¹uₖ₋₁, zero(FC))  # u₀ = 0
 
+    # Initialization, far enough to compute convergence tolerance
+    # Some of this may have to be repeated if we are warm-starting, but
+    # the convergence tolerance should be independent of whether we are warm-starting.
+    # β₁Ev₁ = b
+    kcopy!(m, M⁻¹vₖ, b)  # M⁻¹vₖ ← b
+    MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv)
+    βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ)  # β₁ = ‖v₁‖_E
+    # γ₁Fu₁ = c
+    kcopy!(n, N⁻¹uₖ, c)  # N⁻¹uₖ ← c
+    NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv)
+    γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ)  # γ₁ = ‖u₁‖_F
+    # Convergence tolerance
+    rNorm = sqrt(γₖ^2 + βₖ^2)
+    ε = atol + rtol * rNorm
+
     # [ τ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
+        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
+        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
+      # Repeat the initialization with the updated right-hand side
+      kcopy!(m, M⁻¹vₖ, b₀)  # M⁻¹vₖ ← b₀
+      MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv)
+      βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ)  # β₁ = ‖v₁‖_E
+      kcopy!(n, N⁻¹uₖ, c₀)  # N⁻¹uₖ ← c₀
+      NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv)
+      γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ)  # γ₁ = ‖u₁‖_F
     end
 
-    # β₁Ev₁ = b ↔ β₁v₁ = Mb
-    kcopy!(m, M⁻¹vₖ, b₀)  # M⁻¹vₖ ← b₀
-    MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv)
-    βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ)  # β₁ = ‖v₁‖_E
     if βₖ ≠ 0
       kdiv!(m, M⁻¹vₖ, βₖ)
       MisI || kdiv!(m, vₖ, βₖ)
@@ -242,10 +273,6 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, :
       MisI || kfill!(vₖ, zero(FC))
     end
 
-    # γ₁Fu₁ = c ↔ γ₁u₁ = Nc
-    kcopy!(n, N⁻¹uₖ, c₀)  # N⁻¹uₖ ← c₀
-    NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv)
-    γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ)  # γ₁ = ‖u₁‖_F
     if γₖ ≠ 0
       # u₁ = 0 such that u₁ ⊥ Span{u₁, ..., uₖ}
       kdiv!(n, N⁻¹uₖ, γₖ)
@@ -268,7 +295,6 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, :
     # Compute ‖r₀‖² = (γ₁)² + (β₁)²
     rNorm = sqrt(γₖ^2 + βₖ^2)
     history && push!(rNorms, rNorm)
-    ε = atol + rtol * rNorm
 
     (verbose > 0) && @printf(iostream, "%5s  %7s  %7s  %7s  %5s\n", "k", "‖rₖ‖", "βₖ₊₁", "γₖ₊₁", "timer")
     kdisplay(iter, verbose) && @printf(iostream, "%5d  %7.1e  %7.1e  %7.1e  %.2fs\n", iter, rNorm, βₖ, γₖ, start_time |> ktimer)

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]
@@ -122,6 +130,14 @@ function test_warm_start(FC)
     @test stats.niter <= 1
     @test_throws "x0 should have size $(length(x30))" trimr(A3x2, b3, c2, [1.0], y20)
     @test_throws "y0 should have size $(length(y20))" trimr(A3x2, b3, c2, x30, [1.0])
+
+    N = Diagonal([0.8, 1.2])
+    x30, y20, _ = trimr(A3x2, b3, c2; N, ldiv=true)
+    x3, y2, stats = trimr(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 = TrimrWorkspace(A, b)
     krylov_solve!(workspace, A, b, b, x0, y0)
     r = [b - workspace.x - A * workspace.y; b - A' * workspace.x + workspace.y]