From 4683dadc03890ad9d43ca38e755cb76b4638c828 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 25 Nov 2025 06:30:12 -0600 Subject: [PATCH] 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` --- src/tricg.jl | 53 ++++++++++++++++++++++++++++------------- test/test_warm_start.jl | 8 +++++++ 2 files changed, 44 insertions(+), 17 deletions(-) diff --git a/src/tricg.jl b/src/tricg.jl index 74bd12562..817ee401d 100644 --- a/src/tricg.jl +++ b/src/tricg.jl @@ -33,22 +33,25 @@ 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 `τ = 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 +73,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 +183,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 +219,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 diff --git a/test/test_warm_start.jl b/test/test_warm_start.jl index 3887faa71..55ff109bc 100644 --- a/test/test_warm_start.jl +++ b/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]