diff --git a/docs/src/examples/usymlqr.md b/docs/src/examples/usymlqr.md index 0074f1fd2..b3b1c168c 100644 --- a/docs/src/examples/usymlqr.md +++ b/docs/src/examples/usymlqr.md @@ -1,6 +1,6 @@ ```@example usymlqr using LinearAlgebra, Printf, SparseArrays -using Krylov +using Krylov, LinearOperators # Identity matrix. eye(n::Int) = sparse(1.0 * I, n, n) @@ -39,4 +39,14 @@ d = [0*b; c] r = d - K * [x; y] resid = norm(r) @printf("USYMLQR: Relative residual: %8.1e\n", resid) + +# [D A] [x] = [b] +# [Aᴴ 0] [y] [c] +opH = BlockDiagonalOperator(inv(D), eye(n)) +(x, y, stats) = usymlqr(A, b, c, M=inv(D)) +K = [D A; A' zeros(n,n)] +d = [b; c] +r = d - K * [x; y] +resid = sqrt(dot(r, opH * r)) +@printf("USYMLQR: Relative residual: %8.1e\n", resid) ``` diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md index a77d8d553..e0c2b6346 100644 --- a/docs/src/preconditioners.md +++ b/docs/src/preconditioners.md @@ -111,7 +111,7 @@ Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq) ### Saddle-point and symmetric quasi-definite systems -[`TriCG`](@ref tricg) and [`TriMR`](@ref trimr) can take advantage of the structure of Hermitian systems $Kz = d$ with the 2x2 block structure +[`TriCG`](@ref tricg), [`TriMR`](@ref trimr) and [`USYMLQR`](@ref usymlqr) can take advantage of the structure of Hermitian systems $Kz = d$ with the 2x2 block structure ```math \begin{bmatrix} \tau E & \phantom{-}A \\ A^H & \nu F \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} b \\ c \end{bmatrix}, ``` @@ -119,6 +119,8 @@ Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq) |:---------------:|:---------------------:|:--------------------:|:---------------------:|:--------------------:| | Arguments | `M` with `ldiv=false` | `M` with `ldiv=true` | `N` with `ldiv=false` | `N` with `ldiv=true` | +In the special case of [USYMLQR](@ref usymlqr), $\tau = 1$ and $\nu = 0$. + !!! warning The preconditioners `M` and `N` must be hermitian and positive definite. @@ -133,7 +135,7 @@ Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq) | Arguments | `C` and `E` with `ldiv=false` | `C` and `E` with `ldiv=true` | `D` and `F` with `ldiv=false` | `D` and `F` with `ldiv=true` | !!! note - Our implementations of [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr), [`USYMLQR`](@ref usymlqr) and [`TriLQR`](@ref trilqr) don't support preconditioning. + Our implementations of [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr) and [`TriLQR`](@ref trilqr) don't support preconditioning. ## Packages that provide preconditioners diff --git a/src/krylov_workspaces.jl b/src/krylov_workspaces.jl index 132c06bc6..a4df8c4f8 100644 --- a/src/krylov_workspaces.jl +++ b/src/krylov_workspaces.jl @@ -3239,10 +3239,10 @@ mutable struct UsymlqrWorkspace{T,FC,Sm,Sn} <: _KrylovWorkspace{T,FC,Sm,Sn} x :: Sm y :: Sn z :: Sn - vₖ₋₁ :: Sm - vₖ :: Sm - uₖ₋₁ :: Sn - uₖ :: Sn + M⁻¹vₖ₋₁ :: Sm + M⁻¹vₖ :: Sm + N⁻¹uₖ₋₁ :: Sn + N⁻¹uₖ :: Sn p :: Sn q :: Sm d̅ :: Sm @@ -3250,6 +3250,8 @@ mutable struct UsymlqrWorkspace{T,FC,Sm,Sn} <: _KrylovWorkspace{T,FC,Sm,Sn} wₖ₋₁ :: Sn Δx :: Sm Δy :: Sn + vₖ :: Sm + uₖ :: Sn warm_start :: Bool stats :: SimpleStats{T} end @@ -3264,10 +3266,10 @@ function UsymlqrWorkspace(kc::KrylovConstructor{Sm,Sn}) where {Sm,Sn} x = similar(kc.vm) y = similar(kc.vn) z = similar(kc.vn) - vₖ₋₁ = similar(kc.vm) - vₖ = similar(kc.vm) - uₖ₋₁ = similar(kc.vn) - uₖ = similar(kc.vn) + M⁻¹vₖ₋₁ = similar(kc.vm) + M⁻¹vₖ = similar(kc.vm) + N⁻¹uₖ₋₁ = similar(kc.vn) + N⁻¹uₖ = similar(kc.vn) p = similar(kc.vn) q = similar(kc.vm) d̅ = similar(kc.vm) @@ -3275,8 +3277,10 @@ function UsymlqrWorkspace(kc::KrylovConstructor{Sm,Sn}) where {Sm,Sn} wₖ₋₁ = similar(kc.vn) Δx = similar(kc.vm_empty) Δy = similar(kc.vn_empty) + vₖ = similar(kc.vm_empty) + uₖ = similar(kc.vn_empty) stats = SimpleStats(0, false, false, false, 0, T[], T[], T[], 0.0, 0.0, "unknown") - workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, vₖ₋₁, vₖ, uₖ₋₁, uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, false, stats) + workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, M⁻¹vₖ₋₁, M⁻¹vₖ, N⁻¹uₖ₋₁, N⁻¹uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, vₖ, uₖ, false, stats) workspace.stats.allocation_timer = start_allocation_time |> ktimer return workspace end @@ -3289,10 +3293,10 @@ function UsymlqrWorkspace(m::Integer, n::Integer, Sm::Type, Sn::Type) x = Sm(undef, m) y = Sn(undef, n) z = Sn(undef, n) - vₖ₋₁ = Sm(undef, m) - vₖ = Sm(undef, m) - uₖ₋₁ = Sn(undef, n) - uₖ = Sn(undef, n) + M⁻¹vₖ₋₁ = Sm(undef, m) + M⁻¹vₖ = Sm(undef, m) + N⁻¹uₖ₋₁ = Sn(undef, n) + N⁻¹uₖ = Sn(undef, n) p = Sn(undef, n) q = Sm(undef, m) d̅ = Sm(undef, m) @@ -3300,10 +3304,12 @@ function UsymlqrWorkspace(m::Integer, n::Integer, Sm::Type, Sn::Type) wₖ₋₁ = Sn(undef, n) Δx = Sm(undef, 0) Δy = Sn(undef, 0) + vₖ = Sm(undef, 0) + uₖ = Sn(undef, 0) Sm = isconcretetype(Sm) ? Sm : typeof(x) Sn = isconcretetype(Sn) ? Sn : typeof(y) stats = SimpleStats(0, false, false, false, 0, T[], T[], T[], 0.0, 0.0, "unknown") - workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, vₖ₋₁, vₖ, uₖ₋₁, uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, false, stats) + workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, M⁻¹vₖ₋₁, M⁻¹vₖ, N⁻¹uₖ₋₁, N⁻¹uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, vₖ, uₖ, false, stats) workspace.stats.allocation_timer = start_allocation_time |> ktimer return workspace end diff --git a/src/usymlqr.jl b/src/usymlqr.jl index 3fdb05ad7..e9d3edf0f 100644 --- a/src/usymlqr.jl +++ b/src/usymlqr.jl @@ -18,11 +18,11 @@ export usymlqr, usymlqr! """ (x, y, stats) = usymlqr(A, b::AbstractVector{FC}, c::AbstractVector{FC}; - ls::Bool=true, ln::Bool=true, ldiv::Bool=false, - atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0, - timemax::Float64=Inf, verbose::Int=0, - history::Bool=false, callback=workspace->false, - iostream::IO=kstdout) + M=I, N=I, ls::Bool=true, ln::Bool=true, + ldiv::Bool=false, atol::T=√eps(T), rtol::T=√eps(T), + itmax::Int=0, timemax::Float64=Inf, + verbose::Int=0, history::Bool=false, + callback=workspace->false, iostream::IO=kstdout) `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`. `FC` is `T` or `Complex{T}`. @@ -33,23 +33,29 @@ USYMLQR can be warm-started from initial guesses `x0` and `y0` where `kwargs` ar Solve the symmetric saddle-point system - [ I A ] [ x ] = [ b ] + [ E A ] [ x ] = [ b ] [ Aᴴ ] [ y ] [ c ] -by way of the Saunders-Simon-Yip tridiagonalization using USYMLQ and USYMQR methods. +where E = M⁻¹ ≻ 0 by way of the Saunders-Simon-Yip tridiagonalization using USYMLQ and USYMQR methods. The method solves the least-squares problem when `ls = true` - [ I A ] [ r ] = [ b ] + [ E A ] [ r ] = [ b ] [ Aᴴ ] [ s ] [ 0 ] and the least-norm problem when `ln = true` - [ I A ] [ w ] = [ 0 ] + [ E A ] [ w ] = [ 0 ] [ Aᴴ ] [ z ] [ c ] and simply adds the solutions. + [ M O ] + [ 0 N ] + +indicates the weighted norm in which residuals are measured. +It's the Euclidean norm when `M` and `N` are identity operators. + #### Interface To easily switch between Krylov methods, use the generic interface [`krylov_solve`](@ref) with `method = :usymlqr`. @@ -69,6 +75,8 @@ For an in-place variant that reuses memory across solves, see [`usymlqr!`](@ref) #### Keyword arguments +* `M`: linear operator that models a Hermitian positive-definite matrix of size `m` used for centered preconditioning of the partitioned system; +* `N`: linear operator that models a Hermitian positive-definite matrix of size `n` used for centered preconditioning of the partitioned system; * `ls`: define whether the least-squares problem is solved; * `ln`: define whether the least-norm problem is solved; * `ldiv`: define whether the preconditioners use `ldiv!` or `mul!`; @@ -114,7 +122,9 @@ def_args_usymlqr = (:(A ), def_optargs_usymlqr = (:(x0::AbstractVector), :(y0::AbstractVector)) -def_kwargs_usymlqr = (:(; ls::Bool = true ), +def_kwargs_usymlqr = (:(; M = I ), + :(; N = I ), + :(; ls::Bool = true ), :(; ln::Bool = true ), :(; ldiv::Bool = false ), :(; atol::T = √eps(T) ), @@ -130,7 +140,7 @@ def_kwargs_usymlqr = extract_parameters.(def_kwargs_usymlqr) args_usymlqr = (:A, :b, :c) optargs_usymlqr = (:x0, :y0) -kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) +kwargs_usymlqr = (:M, :N, :ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) @eval begin function usymlqr!(workspace :: UsymlqrWorkspace{T,FC,Sm,Sn}, $(def_args_usymlqr...); $(def_kwargs_usymlqr...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, Sm <: AbstractVector{FC}, Sn <: AbstractVector{FC}} @@ -145,6 +155,10 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi (ls || ln) || error("The keyword arguments `ls` and `ln` can't be both `false`.") (verbose > 0) && @printf(iostream, "USYMLQR: system of %d equations in %d variables\n", m+n, m+n) + # Check M = Iₘ and N = Iₙ + MisI = (M === I) + NisI = (N === I) + # Check type consistency eltype(A) == FC || @warn "eltype(A) ≠ $FC. This could lead to errors or additional allocations in operator-vector products." ktypeof(b) == Sm || error("ktypeof(b) is not a subtype of $Sm") @@ -154,10 +168,14 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi Aᴴ = A' # Set up workspace. + allocate_if(!MisI, workspace, :vₖ, Sm, workspace.x) # The length of vₖ is m + allocate_if(!NisI, workspace, :uₖ, Sn, workspace.y) # The length of uₖ is n Δx, Δy = workspace.Δx, workspace.Δy - vₖ₋₁, vₖ, q, xₖ, rₖ = workspace.vₖ₋₁, workspace.vₖ, workspace.q, workspace.x, workspace.r - uₖ₋₁, uₖ, p, yₖ, zₖ = workspace.uₖ₋₁, workspace.uₖ, workspace.p, workspace.y, workspace.z + M⁻¹vₖ₋₁, M⁻¹vₖ, q, xₖ, rₖ = workspace.M⁻¹vₖ₋₁, workspace.M⁻¹vₖ, workspace.q, workspace.x, workspace.r + N⁻¹uₖ₋₁, N⁻¹uₖ, p, yₖ, zₖ = workspace.N⁻¹uₖ₋₁, workspace.N⁻¹uₖ, workspace.p, workspace.y, workspace.z d̅, wₖ₋₂, wₖ₋₁ = workspace.d̅, workspace.wₖ₋₂, workspace.wₖ₋₁ + vₖ = MisI ? M⁻¹vₖ : workspace.vₖ + uₖ = NisI ? N⁻¹uₖ : workspace.uₖ warm_start = workspace.warm_start b₀ = warm_start ? q : b c₀ = warm_start ? p : c @@ -169,9 +187,9 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi iter = 0 itmax == 0 && (itmax = n+m) - # Initialize orthogonal tridiagonalization process. - kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 - kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 + # Initialize orthogonal tridiagonalization process with respect to elliptic norms. + kfill!(M⁻¹vₖ₋₁, zero(FC)) # v₀ = 0 + kfill!(N⁻¹uₖ₋₁, zero(FC)) # u₀ = 0 # [ I A ] [ xₖ ] = [ b - Δx - AΔy ] = [ b₀ ] # [ Aᴴ ] [ yₖ ] [ c - AᴴΔx ] [ c₀ ] @@ -189,24 +207,30 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi kfill!(yₖ, zero(FC)) kfill!(zₖ, zero(FC)) - # β₁v₁ = b ↔ β₁v₁ = b - kcopy!(m, vₖ, b₀) - βₖ = knorm(m, vₖ) # β₁ = ‖v₁‖ + # β₁Ev₁ = b ↔ β₁v₁ = Mb + kcopy!(m, M⁻¹vₖ, b₀) + MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv) + βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ) # β₁ = ‖v₁‖_E if βₖ ≠ 0 - kdiv!(m, vₖ, βₖ) + kdiv!(m, M⁻¹vₖ, βₖ) + MisI || kdiv!(m, vₖ, βₖ) else # v₁ = 0 such that v₁ ⊥ Span{v₁, ..., vₖ} - kfill!(vₖ, zero(FC)) + kfill!(M⁻¹vₖ, zero(FC)) + MisI || kfill!(vₖ, zero(FC)) end - # γ₁u₁ = c ↔ γ₁u₁ = c - kcopy!(n, uₖ, c₀) - γₖ = knorm(n, uₖ) # γ₁ = ‖u₁‖ + # γ₁Fu₁ = c ↔ γ₁u₁ = Nc + kcopy!(n, N⁻¹uₖ, c₀) + NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv) + γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ) # γ₁ = ‖u₁‖_F if γₖ ≠ 0 - kdiv!(n, uₖ, γₖ) + kdiv!(n, N⁻¹uₖ, γₖ) + NisI || kdiv!(n, uₖ, γₖ) else # u₁ = 0 such that u₁ ⊥ Span{u₁, ..., uₖ} - kfill!(uₖ, zero(FC)) + kfill!(N⁻¹uₖ, zero(FC)) + NisI || kfill!(uₖ, zero(FC)) end cₖ₋₂ = cₖ₋₁ = cₖ = -one(T) # Givens cosines used for the QR factorization of Tₖ₊₁.ₖ @@ -245,30 +269,33 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi # Update iteration index. iter = iter + 1 - # Continue the orthogonal tridiagonalization process. - # AUₖ = VₖTₖ + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ - # AᴴVₖ = Uₖ(Tₖ)ᴴ + γₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ + # Continue the orthogonal tridiagonalization process with respect to elliptic norms. + # AUₖ = EVₖTₖ + βₖ₊₁Evₖ₊₁(eₖ)ᵀ = EVₖ₊₁Tₖ₊₁.ₖ + # AᴴVₖ = FUₖ(Tₖ)ᴴ + γₖ₊₁Fuₖ₊₁(eₖ)ᵀ = FUₖ₊₁(Tₖ.ₖ₊₁)ᴴ - mul!(q, A , uₖ) # Forms vₖ₊₁ : q ← Auₖ - mul!(p, Aᴴ, vₖ) # Forms uₖ₊₁ : p ← Aᴴvₖ + mul!(q, A , uₖ) # Forms Evₖ₊₁ : q ← Auₖ + mul!(p, Aᴴ, vₖ) # Forms Fuₖ₊₁ : p ← Aᴴvₖ if iter ≥ 2 - kaxpy!(m, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁ - kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - βₖ * uₖ₋₁ + kaxpy!(m, -γₖ, M⁻¹vₖ₋₁, q) # q ← q - γₖ * M⁻¹vₖ₋₁ + kaxpy!(n, -βₖ, N⁻¹uₖ₋₁, p) # p ← p - βₖ * N⁻¹uₖ₋₁ end αₖ = kdot(m, vₖ, q) # αₖ = ⟨uₖ,q⟩ - kaxpy!(m, - αₖ , vₖ, q) # q ← q - αₖ * vₖ - kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ + kaxpy!(m, - αₖ , M⁻¹vₖ, q) # q ← q - αₖ * M⁻¹vₖ + kaxpy!(n, -conj(αₖ), N⁻¹uₖ, p) # p ← p - ᾱₖ * N⁻¹uₖ + + # Update M⁻¹vₖ₋₁ and N⁻¹uₖ₋₁ + kcopy!(m, M⁻¹vₖ₋₁, M⁻¹vₖ) + kcopy!(n, N⁻¹uₖ₋₁, N⁻¹uₖ) - # Update vₖ₋₁ and uₖ₋₁ - kcopy!(m, vₖ₋₁, vₖ) - kcopy!(n, uₖ₋₁, uₖ) + # Compute M⁻¹vₖ and N⁻¹uₖ + MisI || mulorldiv!(M⁻¹vₖ, M, q, ldiv) # βₖ₊₁vₖ₊₁ = MAuₖ - γₖvₖ₋₁ - αₖvₖ + NisI || mulorldiv!(N⁻¹uₖ, N, p, ldiv) # γₖ₊₁uₖ₊₁ = NAᴴvₖ - βₖuₖ₋₁ - ᾱₖuₖ - # Compute βₖ₊₁ and γₖ₊₁ - βₖ₊₁ = knorm(m, q) # βₖ₊₁ = ‖vₖ₊₁‖ - γₖ₊₁ = knorm(n, p) # γₖ₊₁ = ‖uₖ₊₁‖ + βₖ₊₁ = MisI ? knorm(m, q) : knorm_elliptic(m, M⁻¹vₖ, q) # βₖ₊₁ = ‖vₖ₊₁‖_E + γₖ₊₁ = NisI ? knorm(n, p) : knorm_elliptic(n, N⁻¹uₖ, p) # γₖ₊₁ = ‖uₖ₊₁‖_F # Update the QR factorization of Tₖ₊₁.ₖ = Qₖ [ Rₖ ]. # [ Oᵀ ] @@ -340,6 +367,10 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi ϕₖ = cₖ * ϕbarₖ ϕbarₖ₊₁ = conj(sₖ) * ϕbarₖ + # println("iter = $iter") + # println("Exact ‖Aᴴrₖ₋₁‖ = ", norm(A' * (b - A * yₖ))) + # println("Exact ‖Aᴴrₖ₋₁‖ = ", norm(A' * rₖ)) + # Update the solution yₖ = Wₖfₖ. # yₖ ← yₖ₋₁ + ϕₖ * wₖ kaxpy!(n, ϕₖ, wₖ, yₖ) @@ -360,6 +391,12 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi # Update ϕbarₖ₊₁ ϕbarₖ = ϕbarₖ₊₁ + # println("Estimate ‖Aᴴrₖ₋₁‖ = ", AᴴrNorm) + # tₖ = b - A * yₖ + # println("Exact ‖rₖ‖ = ", norm(tₖ)) + # println("Exact ‖rₖ‖ = ", norm(rₖ)) + # println("Estimate ‖rₖ‖_M = ", norm(rNorm_ls)) + # Update solved_ls solved_ls = rNorm_ls ≤ ε_ls iter == 1 && (κ = atol + rtol * AᴴrNorm) @@ -433,24 +470,36 @@ kwargs_usymlqr = (:ls, :ln, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hi # Update ηₖ₋₁ ηₖ₋₁ = ηₖ + # println("iter = $iter") + # println("Exact ‖rₖ‖_LQ = ", norm(c - A' * xₖ)) + # println("Exact ‖rₖ‖_LQ = ", norm(c + A' * A * zₖ)) + # println("Estimate ‖rₖ‖_LQ = ", norm(rNorm_ln)) + # println("Exact ‖xₖ + Azₖ‖_LQ = ", norm(xₖ + A * zₖ)) + # Update solved_ln solved_ln = rNorm_ln ≤ ε_ln end # Compute vₖ₊₁ and uₖ₊₁. if βₖ₊₁ ≠ zero(T) + MisI || kdiv!(m, M⁻¹vₖ, βₖ₊₁) kdivcopy!(m, vₖ, q, βₖ₊₁) # vₖ₊₁ = q / βₖ₊₁ + # println("vₖᵀMvₖ: ", dot(vₖ, M⁻¹vₖ)) else # If βₖ₊₁ == 0 then vₖ₊₁ = 0 and Auₖ ∈ Span{v₁, ..., vₖ} # We can keep vₖ₊₁ = 0 such that vₖ₊₁ ⊥ Span{v₁, ..., vₖ} + MisI || kfill!(M⁻¹vₖ, zero(FC)) kfill!(vₖ, zero(FC)) end if γₖ₊₁ ≠ zero(T) + NisI || kdiv!(n, N⁻¹uₖ, γₖ₊₁) kdivcopy!(n, uₖ, p, γₖ₊₁) # uₖ₊₁ = p / γₖ₊₁ + # println("uₖᵀNuₖ: ", dot(uₖ, N⁻¹uₖ)) else # If γₖ₊₁ == 0 then uₖ₊₁ = 0 and Aᴴvₖ ∈ Span{u₁, ..., uₖ} # We can keep uₖ₊₁ = 0 such that uₖ₊₁ ⊥ Span{u₁, ..., uₖ} + NisI || kfill!(N⁻¹uₖ, zero(FC)) kfill!(uₖ, zero(FC)) end diff --git a/test/test_usymlqr.jl b/test/test_usymlqr.jl index 7d56f9505..316ad0bac 100644 --- a/test/test_usymlqr.jl +++ b/test/test_usymlqr.jl @@ -73,5 +73,29 @@ @test typeof(workspace.x) === typeof(b) @test typeof(workspace.y) === typeof(d) @test workspace.stats.solved + + (x, y, stats) = usymlqr(A, b, c, M=D⁻¹, ls=true, ln=false) + K = [D A; A' zeros(n, n)] + d = [b; zeros(FC, n)] + r = d - K * [x; y] + resid = sqrt(dot(r, H⁻¹ * r) |> real) / sqrt(dot(d, H⁻¹ * d) |> real) + @printf("USYMLQR: Relative residual: %8.1e\n", resid) + @test(resid ≤ usymlqr_tol) + + (x, y, stats) = usymlqr(A, b, c, M=D⁻¹, ls=false, ln=true) + K = [D A; A' zeros(n, n)] + d = [zeros(FC, m); c] + r = d - K * [x; y] + resid = sqrt(dot(r, H⁻¹ * r) |> real) / sqrt(dot(d, H⁻¹ * d) |> real) + @printf("USYMLQR: Relative residual: %8.1e\n", resid) + @test(resid ≤ usymlqr_tol) + + (x, y, stats) = usymlqr(A, b, c, M=D⁻¹) + K = [D A; A' zeros(n, n)] + d = [b; c] + r = d - K * [x; y] + resid = sqrt(dot(r, H⁻¹ * r) |> real) / sqrt(dot(d, H⁻¹ * d) |> real) + @printf("USYMLQR: Relative residual: %8.1e\n", resid) + @test(resid ≤ usymlqr_tol) end end