From 31d546423bd6158f093286e00fea6da00bf09125 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 17 Oct 2024 17:00:00 -0500 Subject: [PATCH] An implementation of USYMLQR --- docs/make.jl | 1 + docs/src/api.md | 1 + docs/src/examples/usymlqr.md | 25 ++ docs/src/matrix_free.md | 2 +- docs/src/preconditioners.md | 2 +- docs/src/processes.md | 2 +- docs/src/solvers/sp_sqd.md | 7 + docs/src/storage.md | 6 +- src/Krylov.jl | 1 + src/block_krylov_solvers.jl | 6 +- src/krylov_solve.jl | 1 + src/krylov_solvers.jl | 281 +++++++++++-------- src/minres_qlp.jl | 4 +- src/trilqr.jl | 6 +- src/usymlq.jl | 8 +- src/usymlqr.jl | 516 +++++++++++++++++++++++++++++++++++ src/usymqr.jl | 8 +- test/runtests.jl | 1 + test/test_allocations.jl | 18 ++ test/test_mp.jl | 6 +- test/test_solvers.jl | 7 +- test/test_usymlqr.jl | 29 ++ test/test_verbose.jl | 4 +- test/test_warm_start.jl | 12 + 24 files changed, 822 insertions(+), 132 deletions(-) create mode 100644 docs/src/examples/usymlqr.md create mode 100644 src/usymlqr.jl create mode 100644 test/test_usymlqr.jl diff --git a/docs/make.jl b/docs/make.jl index e62aad428..8062e3e68 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -38,6 +38,7 @@ makedocs( "MINARES" => "examples/minares.md", "TriCG" => "examples/tricg.md", "TriMR" => "examples/trimr.md", + "USYMLQR" => "examples/usymlqr.md", "BICGSTAB" => "examples/bicgstab.md", "DQGMRES" => "examples/dqgmres.md", "BLOCK-GMRES" => "examples/block_gmres.md", diff --git a/docs/src/api.md b/docs/src/api.md index b3050c16e..ac1b6a918 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -32,6 +32,7 @@ DqgmresSolver GmresSolver UsymlqSolver UsymqrSolver +UsymlqrSolver TricgSolver TrimrSolver TrilqrSolver diff --git a/docs/src/examples/usymlqr.md b/docs/src/examples/usymlqr.md new file mode 100644 index 000000000..903daa8aa --- /dev/null +++ b/docs/src/examples/usymlqr.md @@ -0,0 +1,25 @@ +```@example usymlqr +using Krylov, LinearOperators +using LinearAlgebra, Printf, SparseArrays + +# Identity matrix. +eye(n::Int) = sparse(1.0 * I, n, n) + +# Saddle-point systems +n = m = 5 +A = [2^(i/j)*j + (-1)^(i-j) * n*(i-1) for i = 1:n, j = 1:n] +b = ones(n) +D = diagm(0 => [2.0 * i for i = 1:n]) +m, n = size(A) +c = -b + +# [D A] [x] = [b] +# [Aᴴ 0] [y] [c] +opH = BlockDiagonalOperator(inv(D), eye(n)) +(x, y, stats) = usymlqr(A, b, c, M=D, ldiv=true) +K = [D A; A' zeros(n,n)] +B = [b; c] +r = B - K * [x; y] +resid = sqrt(dot(r, opH * r)) +@printf("USYMLQR: Relative residual: %8.1e\n", resid) +``` diff --git a/docs/src/matrix_free.md b/docs/src/matrix_free.md index 95b87bc86..135cffb76 100644 --- a/docs/src/matrix_free.md +++ b/docs/src/matrix_free.md @@ -40,7 +40,7 @@ Some methods only require `A * v` products, whereas other ones also require `A' | CG, CR, CAR | CGLS, CRLS, CGNE, CRMR | | SYMMLQ, CG-LANCZOS, MINRES, MINRES-QLP, MINARES | LSLQ, LSQR, LSMR, LNLQ, CRAIG, CRAIGMR | | DIOM, FOM, DQGMRES, GMRES, FGMRES, BLOCK-GMRES | BiLQ, QMR, BiLQR, USYMLQ, USYMQR, TriLQR | -| CGS, BICGSTAB | TriCG, TriMR | +| CGS, BICGSTAB | TriCG, TriMR, USYMLQR | | CG-LANCZOS-SHIFT | CGLS-LANCZOS-SHIFT | !!! info diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md index 983bc51c7..00b460827 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}, ``` diff --git a/docs/src/processes.md b/docs/src/processes.md index 4583f13e5..681e6623f 100644 --- a/docs/src/processes.md +++ b/docs/src/processes.md @@ -294,7 +294,7 @@ T_{k,k+1} = The function [`saunders_simon_yip`](@ref saunders_simon_yip(::Any, ::AbstractVector{FC}, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\beta_1$, $T_{k+1,k}$, $U_{k+1}$, $\bar{\gamma}_1$ and $T_{k,k+1}^H$. -Related methods: [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr), [`TriLQR`](@ref trilqr), [`TriCG`](@ref tricg) and [`TriMR`](@ref trimr). +Related methods: [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr), [`USYMLQR`](@ref usymlqr), [`TriLQR`](@ref trilqr), [`TriCG`](@ref tricg) and [`TriMR`](@ref trimr). !!! note The Saunders-Simon-Yip is equivalent to the block-Lanczos process applied to $\begin{bmatrix} 0 & A \\ A^H & 0 \end{bmatrix}$ with initial matrix $\begin{bmatrix} b & 0 \\ 0 & c \end{bmatrix}$. diff --git a/docs/src/solvers/sp_sqd.md b/docs/src/solvers/sp_sqd.md index 4ee4ab09b..4d8e0b8c5 100644 --- a/docs/src/solvers/sp_sqd.md +++ b/docs/src/solvers/sp_sqd.md @@ -15,3 +15,10 @@ tricg! trimr trimr! ``` + +## USYMLQR + +```@docs +usymlqr! +usymlqr +``` diff --git a/docs/src/storage.md b/docs/src/storage.md index db51f3661..2639f5aac 100644 --- a/docs/src/storage.md +++ b/docs/src/storage.md @@ -93,9 +93,9 @@ Each table summarizes the storage requirements of Krylov methods recommended to #### Saddle-point and Hermitian quasi-definite systems -| Methods | [`TriCG`](@ref tricg) | [`TriMR`](@ref trimr) | -|:--------:|:---------------------:|:---------------------:| -| Storage | $6n + 6m$ | $8n + 8m$ | +| Methods | [`TriCG`](@ref tricg) | [`TriMR`](@ref trimr) | [`USYMLQR`](@ref usymlqr) | +|:--------:|:---------------------:|:---------------------:|:-------------------------:| +| Storage | $6n + 6m$ | $8n + 8m$ | $7n + 6m$ | #### Generalized saddle-point and non-Hermitian partitioned systems diff --git a/src/Krylov.jl b/src/Krylov.jl index 2578d308b..b429e65e7 100644 --- a/src/Krylov.jl +++ b/src/Krylov.jl @@ -34,6 +34,7 @@ include("gpmr.jl") include("usymlq.jl") include("usymqr.jl") +include("usymlqr.jl") include("tricg.jl") include("trimr.jl") include("trilqr.jl") diff --git a/src/block_krylov_solvers.jl b/src/block_krylov_solvers.jl index a1e8db13a..7944dec37 100644 --- a/src/block_krylov_solvers.jl +++ b/src/block_krylov_solvers.jl @@ -8,14 +8,14 @@ const BLOCK_KRYLOV_SOLVERS = Dict(:block_gmres => :BlockGmresSolver) abstract type BlockKrylovSolver{T,FC,SV,SM} end """ -Type for storing the vectors required by the in-place version of BLOCK-GMRES. +Workspace for the in-place version of BLOCK-GMRES. -The outer constructors +The outer constructors: solver = BlockGmresSolver(m, n, p, memory, SV, SM) solver = BlockGmresSolver(A, B, memory = 5) -may be used in order to create these vectors. +can be used to initialize this workspace. `memory` is set to `div(n,p)` if the value given is larger than `div(n,p)`. `memory` is an optional argument in the second constructor. """ diff --git a/src/krylov_solve.jl b/src/krylov_solve.jl index 69f3d40ae..7ec7f2a1b 100644 --- a/src/krylov_solve.jl +++ b/src/krylov_solve.jl @@ -57,6 +57,7 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs (:FgmresSolver , :fgmres , args_fgmres , def_args_fgmres , optargs_fgmres , def_optargs_fgmres , kwargs_fgmres , def_kwargs_fgmres ) (:FomSolver , :fom , args_fom , def_args_fom , optargs_fom , def_optargs_fom , kwargs_fom , def_kwargs_fom ) (:GpmrSolver , :gpmr , args_gpmr , def_args_gpmr , optargs_gpmr , def_optargs_gpmr , kwargs_gpmr , def_kwargs_gpmr ) + (:UsymlqrSolver , :usymlqr , args_usymlqr , def_args_usymlqr , optargs_usymlqr , def_optargs_usymlqr , kwargs_usymlqr , def_kwargs_usymlqr ) (:CgLanczosShiftSolver , :cg_lanczos_shift , args_cg_lanczos_shift , def_args_cg_lanczos_shift , (), (), kwargs_cg_lanczos_shift , def_kwargs_cg_lanczos_shift ) (:CglsLanczosShiftSolver, :cgls_lanczos_shift, args_cgls_lanczos_shift, def_args_cgls_lanczos_shift, (), (), kwargs_cgls_lanczos_shift, def_kwargs_cgls_lanczos_shift) ] diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl index 2793699f5..0642a0395 100644 --- a/src/krylov_solvers.jl +++ b/src/krylov_solvers.jl @@ -3,7 +3,7 @@ CgLanczosShiftSolver, MinresQlpSolver, DqgmresSolver, DiomSolver, UsymlqSolver, UsymqrSolver, TricgSolver, TrimrSolver, TrilqrSolver, CgsSolver, BicgstabSolver, BilqSolver, QmrSolver, BilqrSolver, CglsSolver, CglsLanczosShiftSolver, CrlsSolver, CgneSolver, CrmrSolver, LslqSolver, LsqrSolver, LsmrSolver, LnlqSolver, CraigSolver, CraigmrSolver, -GmresSolver, FomSolver, GpmrSolver, FgmresSolver, CarSolver, MinaresSolver +GmresSolver, FomSolver, GpmrSolver, UsymlqrSolver, FgmresSolver, CarSolver, MinaresSolver export solve!, solution, nsolution, statistics, issolved, issolved_primal, issolved_dual, niterations, Aprod, Atprod, Bprod, warm_start! @@ -45,6 +45,7 @@ const KRYLOV_SOLVERS = Dict( :lnlq => :LnlqSolver , :craig => :CraigSolver , :craigmr => :CraigmrSolver , + :usymlqr => :UsymlqrSolver , :cg_lanczos_shift => :CgLanczosShiftSolver , :cgls_lanczos_shift => :CglsLanczosShiftSolver, ) @@ -53,14 +54,14 @@ const KRYLOV_SOLVERS = Dict( abstract type KrylovSolver{T,FC,S} end """ -Type for storing the vectors required by the in-place version of MINRES. +Workspace for the in-place version of MINRES. -The outer constructors +The outer constructors: solver = MinresSolver(m, n, S; window :: Int=5) solver = MinresSolver(A, b; window :: Int=5) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct MinresSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -102,14 +103,14 @@ function MinresSolver(A, b; window :: Int=5) end """ -Type for storing the vectors required by the in-place version of MINARES. +Workspace for the in-place version of MINARES. -The outer constructors +The outer constructors: solver = MinaresSolver(m, n, S) solver = MinaresSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct MinaresSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -151,14 +152,14 @@ function MinaresSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CG. +Workspace for the in-place version of CG. -The outer constructors +The outer constructors: solver = CgSolver(m, n, S) solver = CgSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CgSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -194,14 +195,14 @@ function CgSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CR. +Workspace for the in-place version of CR. -The outer constructors +The outer constructors: solver = CrSolver(m, n, S) solver = CrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -239,14 +240,14 @@ function CrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CAR. +Workspace for the in-place version of CAR. -The outer constructors +The outer constructors: solver = CarSolver(m, n, S) solver = CarSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CarSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -288,14 +289,14 @@ function CarSolver(A, b) end """ -Type for storing the vectors required by the in-place version of SYMMLQ. +Workspace for the in-place version of SYMMLQ. -The outer constructors +The outer constructors: solver = SymmlqSolver(m, n, S) solver = SymmlqSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct SymmlqSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -339,14 +340,14 @@ function SymmlqSolver(A, b; window :: Int=5) end """ -Type for storing the vectors required by the in-place version of CG-LANCZOS. +Workspace for the in-place version of CG-LANCZOS. -The outer constructors +The outer constructors: solver = CgLanczosSolver(m, n, S) solver = CgLanczosSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CgLanczosSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -384,14 +385,14 @@ function CgLanczosSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CG-LANCZOS-SHIFT. +Workspace for the in-place version of CG-LANCZOS-SHIFT. -The outer constructors +The outer constructors: solver = CgLanczosShiftSolver(m, n, nshifts, S) solver = CgLanczosShiftSolver(A, b, nshifts) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CgLanczosShiftSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -442,14 +443,14 @@ function CgLanczosShiftSolver(A, b, nshifts) end """ -Type for storing the vectors required by the in-place version of MINRES-QLP. +Workspace for the in-place version of MINRES-QLP. -The outer constructors +The outer constructors: solver = MinresQlpSolver(m, n, S) solver = MinresQlpSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct MinresQlpSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -489,14 +490,14 @@ function MinresQlpSolver(A, b) end """ -Type for storing the vectors required by the in-place version of DQGMRES. +Workspace for the in-place version of DQGMRES. -The outer constructors +The outer constructors: solver = DqgmresSolver(m, n, memory, S) solver = DqgmresSolver(A, b, memory = 20) -may be used in order to create these vectors. +can be used to initialize this workspace. `memory` is set to `n` if the value given is larger than `n`. `memory` is an optional argument in the second constructor. """ @@ -543,14 +544,14 @@ function DqgmresSolver(A, b, memory = 20) end """ -Type for storing the vectors required by the in-place version of DIOM. +Workspace for the in-place version of DIOM. -The outer constructors +The outer constructors: solver = DiomSolver(m, n, memory, S) solver = DiomSolver(A, b, memory = 20) -may be used in order to create these vectors. +can be used to initialize this workspace. `memory` is set to `n` if the value given is larger than `n`. `memory` is an optional argument in the second constructor. """ @@ -595,14 +596,14 @@ function DiomSolver(A, b, memory = 20) end """ -Type for storing the vectors required by the in-place version of USYMLQ. +Workspace for the in-place version of USYMLQ. -The outer constructors +The outer constructors: solver = UsymlqSolver(m, n, S) solver = UsymlqSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct UsymlqSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -644,14 +645,14 @@ function UsymlqSolver(A, b) end """ -Type for storing the vectors required by the in-place version of USYMQR. +Workspace for the in-place version of USYMQR. -The outer constructors +The outer constructors: solver = UsymqrSolver(m, n, S) solver = UsymqrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct UsymqrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -695,14 +696,14 @@ function UsymqrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of TRICG. +Workspace for the in-place version of TRICG. -The outer constructors +The outer constructors: solver = TricgSolver(m, n, S) solver = TricgSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct TricgSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -758,14 +759,14 @@ function TricgSolver(A, b) end """ -Type for storing the vectors required by the in-place version of TRIMR. +Workspace for the in-place version of TRIMR. -The outer constructors +The outer constructors: solver = TrimrSolver(m, n, S) solver = TrimrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct TrimrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -829,14 +830,14 @@ function TrimrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of TRILQR. +Workspace for the in-place version of TRILQR. -The outer constructors +The outer constructors: solver = TrilqrSolver(m, n, S) solver = TrilqrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct TrilqrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -886,14 +887,14 @@ function TrilqrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CGS. +Workspace for the in-place version of CGS. -The outer constructorss +The outer constructors:s solver = CgsSolver(m, n, S) solver = CgsSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CgsSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -935,14 +936,14 @@ function CgsSolver(A, b) end """ -Type for storing the vectors required by the in-place version of BICGSTAB. +Workspace for the in-place version of BICGSTAB. -The outer constructors +The outer constructors: solver = BicgstabSolver(m, n, S) solver = BicgstabSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct BicgstabSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -984,14 +985,14 @@ function BicgstabSolver(A, b) end """ -Type for storing the vectors required by the in-place version of BILQ. +Workspace for the in-place version of BILQ. -The outer constructors +The outer constructors: solver = BilqSolver(m, n, S) solver = BilqSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct BilqSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1037,14 +1038,14 @@ function BilqSolver(A, b) end """ -Type for storing the vectors required by the in-place version of QMR. +Workspace for the in-place version of QMR. -The outer constructors +The outer constructors: solver = QmrSolver(m, n, S) solver = QmrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct QmrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1092,14 +1093,14 @@ function QmrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of BILQR. +Workspace for the in-place version of BILQR. -The outer constructors +The outer constructors: solver = BilqrSolver(m, n, S) solver = BilqrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct BilqrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1149,14 +1150,14 @@ function BilqrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CGLS. +Workspace for the in-place version of CGLS. -The outer constructors +The outer constructors: solver = CglsSolver(m, n, S) solver = CglsSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CglsSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1251,14 +1252,14 @@ function CglsLanczosShiftSolver(A, b, nshifts) end """ -Type for storing the vectors required by the in-place version of CRLS. +Workspace for the in-place version of CRLS. -The outer constructors +The outer constructors: solver = CrlsSolver(m, n, S) solver = CrlsSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CrlsSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1297,14 +1298,14 @@ function CrlsSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CGNE. +Workspace for the in-place version of CGNE. -The outer constructors +The outer constructors: solver = CgneSolver(m, n, S) solver = CgneSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CgneSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1341,14 +1342,14 @@ function CgneSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CRMR. +Workspace for the in-place version of CRMR. -The outer constructors +The outer constructors: solver = CrmrSolver(m, n, S) solver = CrmrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CrmrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1385,14 +1386,14 @@ function CrmrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of LSLQ. +Workspace for the in-place version of LSLQ. -The outer constructors +The outer constructors: solver = LslqSolver(m, n, S) solver = LslqSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct LslqSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1433,14 +1434,14 @@ function LslqSolver(A, b; window :: Int=5) end """ -Type for storing the vectors required by the in-place version of LSQR. +Workspace for the in-place version of LSQR. -The outer constructors +The outer constructors: solver = LsqrSolver(m, n, S) solver = LsqrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct LsqrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1481,14 +1482,14 @@ function LsqrSolver(A, b; window :: Int=5) end """ -Type for storing the vectors required by the in-place version of LSMR. +Workspace for the in-place version of LSMR. -The outer constructors +The outer constructors: solver = LsmrSolver(m, n, S) solver = LsmrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct LsmrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1531,14 +1532,14 @@ function LsmrSolver(A, b; window :: Int=5) end """ -Type for storing the vectors required by the in-place version of LNLQ. +Workspace for the in-place version of LNLQ. -The outer constructors +The outer constructors: solver = LnlqSolver(m, n, S) solver = LnlqSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct LnlqSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1581,14 +1582,14 @@ function LnlqSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CRAIG. +Workspace for the in-place version of CRAIG. -The outer constructors +The outer constructors: solver = CraigSolver(m, n, S) solver = CraigSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CraigSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1631,14 +1632,14 @@ function CraigSolver(A, b) end """ -Type for storing the vectors required by the in-place version of CRAIGMR. +Workspace for the in-place version of CRAIGMR. -The outer constructors +The outer constructors: solver = CraigmrSolver(m, n, S) solver = CraigmrSolver(A, b) -may be used in order to create these vectors. +can be used to initialize this workspace. """ mutable struct CraigmrSolver{T,FC,S} <: KrylovSolver{T,FC,S} m :: Int @@ -1685,14 +1686,14 @@ function CraigmrSolver(A, b) end """ -Type for storing the vectors required by the in-place version of GMRES. +Workspace for the in-place version of GMRES. -The outer constructors +The outer constructors: solver = GmresSolver(m, n, memory, S) solver = GmresSolver(A, b, memory = 20) -may be used in order to create these vectors. +can be used to initialize this workspace. `memory` is set to `n` if the value given is larger than `n`. `memory` is an optional argument in the second constructor. """ @@ -1740,14 +1741,14 @@ function GmresSolver(A, b, memory = 20) end """ -Type for storing the vectors required by the in-place version of FGMRES. +Workspace for the in-place version of FGMRES. -The outer constructors +The outer constructors: solver = FgmresSolver(m, n, memory, S) solver = FgmresSolver(A, b, memory = 20) -may be used in order to create these vectors. +can be used to initialize this workspace. `memory` is set to `n` if the value given is larger than `n`. `memory` is an optional argument in the second constructor. """ @@ -1795,14 +1796,14 @@ function FgmresSolver(A, b, memory = 20) end """ -Type for storing the vectors required by the in-place version of FOM. +Workspace for the in-place version of FOM. -The outer constructors +The outer constructors: solver = FomSolver(m, n, memory, S) solver = FomSolver(A, b, memory = 20) -may be used in order to create these vectors. +can be used to initialize this workspace. `memory` is set to `n` if the value given is larger than `n`. `memory` is an optional argument in the second constructor. """ @@ -1847,14 +1848,14 @@ function FomSolver(A, b, memory = 20) end """ -Type for storing the vectors required by the in-place version of GPMR. +Workspace for the in-place version of GPMR. -The outer constructors +The outer constructors: solver = GpmrSolver(m, n, memory, S) solver = GpmrSolver(A, b, memory = 20) -may be used in order to create these vectors. +can be used to initialize this workspace. `memory` is set to `n + m` if the value given is larger than `n + m`. `memory` is an optional argument in the second constructor. """ @@ -1912,6 +1913,71 @@ function GpmrSolver(A, b, memory = 20) GpmrSolver(m, n, memory, S) end +""" +Workspace for the in-place version of USYMLQR. + +The outer constructors: + + solver = UsymlqrSolver(m, n, S) + solver = UsymlqrSolver(A, b) + +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 + z :: S + M⁻¹uₖ₋₁ :: S + M⁻¹uₖ :: S + N⁻¹vₖ₋₁ :: S + N⁻¹vₖ :: S + p :: S + q :: 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) + 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) + 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, 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 + +function UsymlqrSolver(A, b) + m, n = size(A) + S = ktypeof(b) + UsymlqrSolver(m, n, S) +end + """ solution(solver) @@ -2013,6 +2079,7 @@ for (KS, fun, nsol, nA, nAt, warm_start) in [ (:FgmresSolver , :fgmres! , 1, 1, 0, true ) (:FomSolver , :fom! , 1, 1, 0, true ) (:GpmrSolver , :gpmr! , 2, 1, 0, true ) + (:UsymlqrSolver , :usymlqr! , 2, 1, 1, true ) (:CgLanczosShiftSolver , :cg_lanczos_shift! , 1, 1, 0, false) (:CglsLanczosShiftSolver, :cgls_lanczos_shift!, 1, 1, 1, false) ] @@ -2044,7 +2111,7 @@ for (KS, fun, nsol, nA, nAt, warm_start) in [ issolved(solver :: $KS) = solver.stats.solved end if $warm_start - if $KS in (BilqrSolver, TrilqrSolver, TricgSolver, TrimrSolver, GpmrSolver) + if $KS in (BilqrSolver, TrilqrSolver, UsymlqrSolver, TricgSolver, TrimrSolver, GpmrSolver) function warm_start!(solver :: $KS, x0, y0) n = length(solver.x) m = length(solver.y) diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl index f981da740..9c96346e4 100644 --- a/src/minres_qlp.jl +++ b/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ₖ₊₁] diff --git a/src/trilqr.jl b/src/trilqr.jl index 307fc0b67..4ea3685a2 100644 --- a/src/trilqr.jl +++ b/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⟩ diff --git a/src/usymlq.jl b/src/usymlq.jl index 18b859fa1..9474c1c59 100644 --- a/src/usymlq.jl +++ b/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 βₖ₊₁ = knorm(m, q) # βₖ₊₁ = ‖q‖ γₖ₊₁ = knorm(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ₖ diff --git a/src/usymlqr.jl b/src/usymlqr.jl new file mode 100644 index 000000000..c1b71ca2f --- /dev/null +++ b/src/usymlqr.jl @@ -0,0 +1,516 @@ +# An implementation of USYMLQR for the solution of symmetric saddle-point systems. +# +# This method is described in +# +# M. A. Saunders, H. D. Simon, and E. L. Yip +# Two Conjugate-Gradient-Type Methods for Unsymmetric Linear Equations. +# SIAM Journal on Numerical Analysis, 25(4), pp. 927--940, 1988. +# +# A. Buttari, D. Orban, D. Ruiz and D. Titley-Peloquin +# A tridiagonalization method for symmetric saddle-point systems. +# SIAM Journal on Scientific Computing, 41(5), pp. 409--432, 2019 +# +# Dominique Orban, +# Alexis Montoison, -- +# Montréal, November 2021 -- Chicago, October 2024. + +export usymlqr, usymlqr! + +""" + (x, y, stats) = usymlqr(A, b::AbstractVector{FC}, c::AbstractVector{FC}; + M=I, N=I, ldiv::Bool=false, atol::T=√eps(T), + rtol::T=√eps(T), itmax::Int=0, timemax::Float64=Inf, + verbose::Int=0, history::Bool=false, + callback=solver->false, iostream::IO=kstdout) + +`T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`. +`FC` is `T` or `Complex{T}`. + + (x, y, stats) = usymlqr(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...) + +USYMLQR can be warm-started from initial guesses `x0` and `y0` where `kwargs` are the same keyword arguments as above. + +Solve the symmetric saddle-point system + + [ E A ] [ x ] = [ b ] + [ Aᴴ ] [ y ] [ c ] + +where E = M⁻¹ ≻ 0 by way of the Saunders-Simon-Yip tridiagonalization using USYMLQ and USYMQR methods. +The method solves the least-squares problem + + [ E A ] [ s ] = [ b ] + [ Aᴴ ] [ t ] [ 0 ] + +and the least-norm problem + + [ 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. + +#### Input arguments + +* `A`: a linear operator that models a matrix of dimension m × n; +* `b`: a vector of length m; +* `c`: a vector of length n. + +#### Optional arguments + +* `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. + +#### 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; +* `ldiv`: define whether the preconditioners use `ldiv!` or `mul!`; +* `atol`: absolute stopping tolerance based on the residual norm; +* `rtol`: relative stopping tolerance based on the residual norm; +* `itmax`: the maximum number of iterations. If `itmax=0`, the default number of iterations is set to `m+n`; +* `timemax`: the time limit in seconds; +* `verbose`: additional details can be kdisplayed if verbose mode is enabled (verbose > 0). Information will be kdisplayed every `verbose` iterations; +* `history`: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms; +* `callback`: function or functor called as `callback(solver)` that returns `true` if the Krylov method should terminate, and `false` otherwise; +* `iostream`: stream to which output is logged. + +#### Output arguments + +* `x`: a dense vector of length m; +* `y`: a dense vector of length n; +* `stats`: statistics collected on the run in a [`SimpleStats`](@ref) structure. + +#### References + +* M. A. Saunders, H. D. Simon, and E. L. Yip, [*Two Conjugate-Gradient-Type Methods for Unsymmetric Linear Equations*](https://doi.org/10.1137/0725052), SIAM Journal on Numerical Analysis, 25(4), pp. 927--940, 1988. +* A. Buttari, D. Orban, D. Ruiz and D. Titley-Peloquin, [*A tridiagonalization method for symmetric saddle-point and quasi-definite systems*](https://doi.org/10.1137/18M1194900), SIAM Journal on Scientific Computing, 41(5), pp. 409--432, 2019. +""" +function usymlqr end + +""" + solver = usymlqr!(solver::UsymlqrSolver, A, b, c; kwargs...) + solver = usymlqr!(solver::UsymlqrSolver, A, b, c, x0, y0; kwargs...) + +where `kwargs` are keyword arguments of [`usymlqr`](@ref). + +See [`UsymlqrSolver`](@ref) for more details about the `solver`. +""" +function usymlqr! end + +def_args_usymlqr = (:(A ), + :(b::AbstractVector{FC}), + :(c::AbstractVector{FC})) + +def_optargs_usymlqr = (:(x0::AbstractVector), + :(y0::AbstractVector)) + +def_kwargs_usymlqr = (:(; transfer_to_usymcg::Bool = true), + :(; M = I ), + :(; N = I ), + :(; ldiv::Bool = false ), + :(; atol::T = √eps(T) ), + :(; rtol::T = √eps(T) ), + :(; itmax::Int = 0 ), + :(; timemax::Float64 = Inf ), + :(; verbose::Int = 0 ), + :(; history::Bool = false ), + :(; callback = solver -> false ), + :(; iostream::IO = kstdout )) + +def_kwargs_usymlqr = extract_parameters.(def_kwargs_usymlqr) + +args_usymlqr = (:A, :b, :c) +optargs_usymlqr = (:x0, :y0) +kwargs_usymlqr = (:transfer_to_usymcg, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) + +@eval begin + function usymlqr!(solver :: UsymlqrSolver{T,FC,S}, $(def_args_usymlqr...); $(def_kwargs_usymlqr...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: AbstractVector{FC}} + + # Timer + start_time = time_ns() + timemax_ns = 1e9 * timemax + + m, n = size(A) + (m == solver.m && n == solver.n) || error("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($m, $n)") + length(b) == m || error("Inconsistent problem size") + length(c) == n || error("Inconsistent problem size") + (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) <: S || error("ktypeof(b) is not a subtype of $S") + ktypeof(c) <: S || error("ktypeof(c) is not a subtype of $S") + + # Compute the adjoint of A + Aᴴ = A' + + # Set up workspace. + allocate_if(!MisI, solver, :uₖ, S, m) + allocate_if(!NisI, solver, :vₖ, S, n) + Δxz, Δry = solver.Δx, solver.Δy + M⁻¹uₖ₋₁, M⁻¹uₖ, N⁻¹vₖ₋₁, N⁻¹vₖ = solver.M⁻¹uₖ₋₁, solver.M⁻¹uₖ, solver.N⁻¹vₖ₋₁, solver.N⁻¹vₖ + rₖ, xₖ, yₖ, zₖ, p, q = solver.r, solver.x, solver.y, solver.z, solver.p, solver.q + d̅, wₖ₋₂, wₖ₋₁ = solver.d̅, solver.wₖ₋₂, solver.wₖ₋₁ + uₖ = MisI ? M⁻¹uₖ : solver.uₖ + vₖ = NisI ? N⁻¹vₖ : solver.vₖ + warm_start = solver.warm_start + b₀ = warm_start ? q : b + c₀ = warm_start ? p : c + + stats = solver.stats + rNorms = stats.residuals + reset!(stats) + + iter = 0 + itmax == 0 && (itmax = n+m) + + # Initial solutions r₀, x₀, y₀ and z₀. + @kfill!(rₖ, zero(FC)) + @kfill!(xₖ, zero(FC)) + @kfill!(yₖ, zero(FC)) + @kfill!(zₖ, zero(FC)) + + # Initialize preconditioned orthogonal tridiagonalization process. + @kfill!(M⁻¹uₖ₋₁, zero(FC)) # u₀ = 0 + @kfill!(N⁻¹vₖ₋₁, zero(FC)) # v₀ = 0 + + # [ I A ] [ xₖ ] = [ b - Δx - AΔy ] = [ b₀ ] + # [ Aᴴ ] [ yₖ ] [ c - AᴴΔx ] [ c₀ ] + if warm_start + mul!(b₀, A, Δy) + @kaxpy!(m, one(T), Δx, b₀) + @kaxpby!(m, one(T), b, -one(T), b₀) + mul!(c₀, Aᴴ, Δx) + @kaxpby!(n, one(T), c, -one(T), c₀) + end + + # β₁Eu₁ = b ↔ β₁u₁ = Mb + M⁻¹uₖ .= b₀ + MisI || mul!(uₖ, M, M⁻¹uₖ) + βₖ = sqrt(@kdot(m, uₖ, M⁻¹uₖ)) # β₁ = ‖u₁‖_E + if βₖ ≠ 0 + @kscal!(m, 1 / βₖ, M⁻¹uₖ) + MisI || @kscal!(m, 1 / βₖ, uₖ) + else + error("b must be nonzero") + end + + # γ₁Fv₁ = c ↔ γ₁v₁ = Nc + N⁻¹vₖ .= c₀ + NisI || mul!(vₖ, N, N⁻¹vₖ) + γₖ = sqrt(@kdot(n, vₖ, N⁻¹vₖ)) # γ₁ = ‖v₁‖_F + if γₖ ≠ 0 + @kscal!(n, 1 / γₖ, N⁻¹vₖ) + NisI || @kscal!(n, 1 / γₖ, vₖ) + else + error("c must be nonzero") + end + + ε = atol + rtol * rNorm + κ = zero(T) + (verbose > 0) && @printf(iostream, "%4s %7s %7s %7s\n", "k", "αₖ", "βₖ", "γₖ") + kdisplay(iter, verbose) && @printf(iostream, "%4d %7.1e %7.1e %7.1e\n", iter, αₖ, βₖ, γₖ) + + cₖ₋₂ = cₖ₋₁ = cₖ = one(T) # Givens cosines used for the QR factorization of Tₖ₊₁.ₖ + sₖ₋₂ = sₖ₋₁ = sₖ = zero(FC) # Givens sines used for the QR factorization of Tₖ₊₁.ₖ + @kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Vₖ(Rₖ)⁻¹ + @kfill!(wₖ₋₁, zero(FC)) # Column k-1 of Wₖ = Vₖ(Rₖ)⁻¹ + ϕbarₖ = βₖ # ϕbarₖ is the last component of f̄ₖ = (Qₖ)ᴴβ₁e₁ + @kfill!(d̅, zero(FC)) # Last column of D̅ₖ = UₖQₖ + ηₖ₋₁ = ηbarₖ = zero(FC) # ηₖ₋₁ and ηbarₖ are the last components of h̄ₖ = (Rₖ)⁻ᵀγ₁e₁ + ηₖ₋₂ = θₖ = zero(FC) # ζₖ₋₂ and θₖ are used to update ηₖ₋₁ and ηbarₖ + δbarₖ₋₁ = δbarₖ = zero(FC) # Coefficients of Rₖ₋₁ and Rₖ modified over the course of two iterations + + # Stopping criterion. + rNorm_LS = β₁ + rNorm_LN = γ₁ + solved_LS = rNorm_LS ≤ ε + solved_LN = rNorm_LN ≤ ε + solved = solved_LS && solved_LN + inconsistent = false + tired = iter ≥ itmax + status = "unknown" + ill_cond = false + user_requested_exit = false + overtimed = false + + while !(solved || tired || ill_cond || inconsistent || user_requested_exit || overtimed) + # Update iteration index. + iter = iter + 1 + + # Continue the SSY tridiagonalization process. + # AVₖ = EUₖTₖ + βₖ₊₁Euₖ₊₁(eₖ)ᵀ = EUₖ₊₁Tₖ₊₁.ₖ + # AᴴUₖ = FVₖ(Tₖ)ᴴ + γₖ₊₁Fvₖ₊₁(eₖ)ᵀ = FVₖ₊₁(Tₖ.ₖ₊₁)ᴴ + + mul!(q, A , vₖ) # Forms Euₖ₊₁ : q ← Avₖ + mul!(p, Aᴴ, uₖ) # Forms Fvₖ₊₁ : p ← Aᴴuₖ + + if iter ≥ 2 + @kaxpy!(m, -γₖ, M⁻¹uₖ₋₁, q) # q ← q - γₖ * M⁻¹uₖ₋₁ + @kaxpy!(n, -βₖ, N⁻¹vₖ₋₁, p) # p ← p - βₖ * N⁻¹vₖ₋₁ + end + + αₖ = @kdot(m, uₖ, q) # αₖ = ⟨uₖ,q⟩ + + @kaxpy!(m, - αₖ , M⁻¹uₖ, q) # q ← q - αₖ * M⁻¹uₖ + @kaxpy!(n, -conj(αₖ), N⁻¹vₖ, p) # p ← p - ᾱₖ * N⁻¹vₖ + + # Compute vₖ₊₁ and uₖ₊₁ + MisI || mulorldiv!(uₖ₊₁, M, q, ldiv) # βₖ₊₁uₖ₊₁ = MAvₖ - γₖuₖ₋₁ - αₖuₖ + NisI || mulorldiv!(vₖ₊₁, N, p, ldiv) # γₖ₊₁vₖ₊₁ = NAᴴuₖ - βₖvₖ₋₁ - ᾱₖvₖ + + βₖ₊₁ = sqrt(@kdotr(m, uₖ₊₁, q)) # βₖ₊₁ = ‖uₖ₊₁‖_E + γₖ₊₁ = sqrt(@kdotr(n, vₖ₊₁, p)) # γₖ₊₁ = ‖vₖ₊₁‖_F + + # Update M⁻¹uₖ₋₁ and N⁻¹vₖ₋₁ + M⁻¹uₖ₋₁ .= M⁻¹uₖ + N⁻¹vₖ₋₁ .= N⁻¹vₖ + + # Update the QR factorization of Tₖ₊₁.ₖ = Qₖ [ Rₖ ]. + # [ Oᵀ ] + # [ α₁ γ₂ 0 • • • 0 ] [ δ₁ λ₁ ϵ₁ 0 • • 0 ] + # [ β₂ α₂ γ₃ • • ] [ 0 δ₂ λ₂ • • • ] + # [ 0 • • • • • ] [ • • δ₃ • • • • ] + # [ • • • • • • • ] = Qₖ [ • • • • • 0 ] + # [ • • • • • 0 ] [ • • • • ϵₖ₋₂] + # [ • • • • γₖ ] [ • • • λₖ₋₁] + # [ • • βₖ αₖ ] [ 0 • • • • 0 δₖ ] + # [ 0 • • • • 0 βₖ₊₁] [ 0 • • • • • 0 ] + # + # If k = 1, we don't have any previous reflexion. + # If k = 2, we apply the last reflexion. + # If k ≥ 3, we only apply the two previous reflexions. + + # Apply previous Givens reflections Qₖ₋₂.ₖ₋₁ + if iter ≥ 3 + # [cₖ₋₂ sₖ₋₂] [0 ] = [ ϵₖ₋₂ ] + # [s̄ₖ₋₂ -cₖ₋₂] [γₖ] [λbarₖ₋₁] + ϵₖ₋₂ = sₖ₋₂ * γₖ + λbarₖ₋₁ = -cₖ₋₂ * γₖ + end + + # Apply previous Givens reflections Qₖ₋₁.ₖ + if iter ≥ 2 + iter == 2 && (λbarₖ₋₁ = γₖ) + # [cₖ₋₁ sₖ₋₁] [λbarₖ₋₁] = [λₖ₋₁ ] + # [s̄ₖ₋₁ -cₖ₋₁] [ αₖ ] [δbarₖ] + λₖ₋₁ = cₖ₋₁ * λbarₖ₋₁ + sₖ₋₁ * αₖ + δbarₖ = conj(sₖ₋₁) * λbarₖ₋₁ - cₖ₋₁ * αₖ + end + + # Compute and apply current Givens reflection Qₖ.ₖ₊₁ + iter == 1 && (δbarₖ = αₖ) + # [cₖ sₖ] [δbarₖ] = [δₖ] + # [s̄ₖ -cₖ] [βₖ₊₁ ] [0 ] + (cₖ, sₖ, δₖ) = sym_givens(δbarₖ, βₖ₊₁) + + # Compute f̄ₖ₊₁ = [ fₖ ] = (Qₖ)ᴴβ₁e₁ + # [ϕbarₖ₊₁] + # + # [cₖ sₖ] [ϕbarₖ] = [ ϕₖ ] + # [s̄ₖ -cₖ] [ 0 ] [ϕbarₖ₊₁] + ϕₖ = cₖ * ϕbarₖ + ϕbarₖ₊₁ = conj(sₖ) * ϕbarₖ + + # Compute the direction wₖ, the last column of Wₖ = Vₖ(Rₖ)⁻¹ ⟷ (Rₖ)ᵀ(Wₖ)ᵀ = (Vₖ)ᵀ. + # w₁ = v₁ / δ₁ + if iter == 1 + wₖ = wₖ₋₁ + @kaxpy!(n, one(FC), vₖ, wₖ) + wₖ .= wₖ ./ δₖ + end + # w₂ = (v₂ - λ₁w₁) / δ₂ + if iter == 2 + wₖ = wₖ₋₂ + @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) + @kaxpy!(n, one(FC), vₖ, wₖ) + wₖ .= wₖ ./ δₖ + end + # wₖ = (vₖ - λₖ₋₁wₖ₋₁ - ϵₖ₋₂wₖ₋₂) / δₖ + if iter ≥ 3 + @kscal!(n, -ϵₖ₋₂, wₖ₋₂) + wₖ = wₖ₋₂ + @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) + @kaxpy!(n, one(FC), vₖ, wₖ) + wₖ .= wₖ ./ δₖ + end + + # Update the solution xₖ. + # xₖ ← xₖ₋₁ + ϕₖ * wₖ + @kaxpy!(n, ϕₖ, wₖ, xₖ) + + # Update the residual rₖ. + # rₖ ← |sₖ|² * rₖ₋₁ - cₖ * ϕbarₖ₊₁ * uₖ₊₁ + @kaxpby!(n, cₖ * ϕbarₖ₊₁, q, abs2(sₖ), rₖ) + + # Compute ‖rₖ‖ = |ϕbarₖ₊₁|. + rNorm = abs(ϕbarₖ₊₁) + history && push!(rNorms, rNorm) + + # Compute ‖Aᴴrₖ₋₁‖ = |ϕbarₖ| * √(|δbarₖ|² + |λbarₖ|²). + AᴴrNorm = abs(ϕbarₖ) * √(abs2(δbarₖ) + abs2(cₖ₋₁ * γₖ₊₁)) + history && push!(AᴴrNorms, AᴴrNorm) + + # Compute ηₖ₋₁ and ηbarₖ, last components of the solution of (Rₖ)ᵀh̄ₖ = γ₁e₁ + # [δbar₁] [ηbar₁] = [γ₁] + if iter == 1 + ωₖ = γₖ + end + # [δ₁ 0 ] [ η₁ ] = [γ₁] + # [λ₁ δbar₂] [ηbar₂] [0 ] + if iter == 2 + ωₖ₋₁ = ηₖ + ηₖ₋₁ = ηₖ₋₁ / δₖ₋₁ + ωₖ = -λₖ₋₁ * ζₖ₋₁ + end + # [λₖ₋₂ δₖ₋₁ 0 ] [ηₖ₋₂ ] = [0] + # [ϵₖ₋₂ λₖ₋₁ δbarₖ] [ηₖ₋₁ ] [0] + # [ηbarₖ] + if iter ≥ 3 + ηₖ₋₂ = ζₖ₋₁ + ωₖ₋₁ = ηₖ + ηₖ₋₁ = ηₖ₋₁ / δₖ₋₁ + ωₖ = -ϵₖ₋₂ * ηₖ₋₂ - λₖ₋₁ * ηₖ₋₁ + end + + # Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = UₖQₖ. + # Note: D̄ₖ represents the matrix P̄ₖ in the paper of USYMLQR. + # [d̅ₖ₋₁ uₖ] [cₖ s̄ₖ] = [dₖ₋₁ d̅ₖ] ⟷ dₖ₋₁ = cₖ * d̅ₖ₋₁ + sₖ * uₖ + # [sₖ -cₖ] ⟷ d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ + if iter ≥ 2 + # Compute solution yₖ. + # (yᴸ)ₖ₋₁ ← (yᴸ)ₖ₋₂ + ηₖ₋₁ * dₖ₋₁ + @kaxpy!(n, ηₖ₋₁ * cₖ, d̅, x) + @kaxpy!(n, ηₖ₋₁ * sₖ, uₖ, x) + end + + # Compute d̅ₖ. + if iter == 1 + # d̅₁ = u₁ + @kcopy!(n, uₖ, d̅) # d̅ ← vₖ + else + # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ + @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅) + end + + # Compute USYMLQ residual norm + # ‖rₖ‖ = √(|μₖ|² + |ωₖ|²) + if iter == 1 + rNorm_lq = bNorm + else + μₖ = βₖ * (sₖ₋₁ * ζₖ₋₂ - cₖ₋₁ * cₖ * ζₖ₋₁) + αₖ * sₖ * ζₖ₋₁ + ωₖ = βₖ₊₁ * sₖ * ζₖ₋₁ + rNorm_lq = sqrt(abs2(μₖ) + abs2(ωₖ)) + end + history && push!(rNorms, rNorm_lq) + + # Compute USYMCG residual norm + # ‖rₖ‖ = |ρₖ| + if transfer_to_usymcg && (abs(δbarₖ) > eps(T)) + ζbarₖ = ηₖ / δbarₖ + ρₖ = βₖ₊₁ * (sₖ * ζₖ₋₁ - cₖ * ζbarₖ) + rNorm_cg = abs(ρₖ) + end + + # Compute zₖ. + @kaxpy!(n, -ηₖ, wₖ, zₖ) + + # Compute uₖ₊₁ and vₖ₊₁. + @kcopy!(m, uₖ, uₖ₋₁) # uₖ₋₁ ← uₖ + @kcopy!(n, vₖ, vₖ₋₁) # vₖ₋₁ ← vₖ + + if βₖ₊₁ ≠ zero(T) + uₖ .= q ./ βₖ₊₁ # βₖ₊₁uₖ₊₁ = q + end + if γₖ₊₁ ≠ zero(T) + vₖ .= p ./ γₖ₊₁ # γₖ₊₁vₖ₊₁ = p + end + + # Update directions for x. + if iter ≥ 2 + @kswap(wₖ₋₂, wₖ₋₁) + end + + # Update sₖ₋₂, cₖ₋₂, sₖ₋₁, cₖ₋₁, ϕbarₖ, γₖ, βₖ. + if iter ≥ 2 + sₖ₋₂ = sₖ₋₁ + cₖ₋₂ = cₖ₋₁ + end + sₖ₋₁ = sₖ + cₖ₋₁ = cₖ + ϕbarₖ = ϕbarₖ₊₁ + γₖ = γₖ₊₁ + βₖ = βₖ₊₁ + + # Update δbarₖ₋₁. + # δbarₖ₋₁ = δbarₖ + + # Update stopping criterion. + ill_cond_lim = one(T) / Acond ≤ ctol + ill_cond_mach = one(T) + one(T) / Acond ≤ one(T) + ill_cond = ill_cond_mach || ill_cond_lim + tired = iter ≥ itmax + solved = solved_LS && solved_LN + + iter == 1 && (κ = atol + rtol * AᴴrNorm) + user_requested_exit = callback(solver) :: Bool + solved = rNorm ≤ ε + solved_lq = rNorm_lq ≤ ε + solved_cg = transfer_to_usymcg && (abs(δbarₖ) > eps(T)) && (rNorm_cg ≤ ε) + inconsistent = !solved && AᴴrNorm ≤ κ + tired = iter ≥ itmax + timer = time_ns() - start_time + overtimed = timer > timemax_ns + kdisplay(iter, verbose) && @printf(iostream, "%7.1e\n", rNorm_lq) + kdisplay(iter, verbose) && @printf(iostream, "%4d %8.1e %7.1e %7.1e %7.1e %7.1e %7.1e ", iter, αₖ, βₖ, γₖ, Anorm, Acond, rNorm_qr) + kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %8.1e %.2fs\n", iter, rNorm, AᴴrNorm, ktimer(start_time)) + end + (verbose > 0) && @printf(iostream, "\n") + + # Compute USYMCG point + # (yᶜ)ₖ ← (yᴸ)ₖ₋₁ + ηbarₖ * d̅ₖ + # (zᶜ)ₖ ← (zᴸ)ₖ₋₁ - ηbarₖ * w̄ₖ + if solved_cg + @kaxpy!(n, ηbarₖ, d̅, yₖ) + @kaxpy!(m, -ηbarₖ, w̄, zₖ) + end + + # Termination status + tired && (status = "maximum number of iterations exceeded") + # solved_lq && (status = "solution xᴸ good enough given atol and rtol") + # solved_cg && (status = "solution xᶜ good enough given atol and rtol") + solved && (status = "solution good enough given atol and rtol") + ill_cond_mach && (status = "condition number seems too large for this machine") + ill_cond_lim && (status = "condition number exceeds tolerance") + user_requested_exit && (status = "user-requested exit") + overtimed && (status = "time limit exceeded") + + # Compute the solution the saddle point system + # xₖ ← xₖ + zₖ + # yₖ ← yₖ + rₖ + @kaxpy!(n, one(FC), zₖ, xₖ) + @kaxpy!(m, one(FC), rₖ, yₖ) + + # Update xₖ and yₖ + warm_start && @kaxpy!(n, one(FC), Δxz, xₖ) + warm_start && @kaxpy!(m, one(FC), Δyr, yₖ) + solver.warm_start = false + + # Update stats + stats.niter = iter + stats.solved = solved + stats.inconsistent = false + stats.timer = ktimer(start_time) + stats.status = status + return solver + end +end diff --git a/src/usymqr.jl b/src/usymqr.jl index e385ef43f..ad44326af 100644 --- a/src/usymqr.jl +++ b/src/usymqr.jl @@ -201,8 +201,10 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, 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⟩ @@ -293,7 +295,7 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, AᴴrNorm = abs(ζbarₖ) * √(abs2(δbarₖ) + abs2(cₖ₋₁ * γₖ₊₁)) history && push!(AᴴrNorms, AᴴrNorm) - # Compute uₖ₊₁ and uₖ₊₁. + # Compute vₖ₊₁ and uₖ₊₁. kcopy!(m, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ diff --git a/test/runtests.jl b/test/runtests.jl index d342af4f8..b157b399e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ include("test_stats.jl") include("test_block_processes.jl") include("test_processes.jl") +include("test_usymlqr.jl") include("test_minares.jl") include("test_car.jl") include("test_fgmres.jl") diff --git a/test/test_allocations.jl b/test/test_allocations.jl index e017a75ae..5dd793ff2 100644 --- a/test/test_allocations.jl +++ b/test/test_allocations.jl @@ -615,6 +615,24 @@ @test inplace_trimr_bytes == 0 end + @testset "USYMLQR" begin + # USYMLQR needs: + # - 7 n-vectors: x, z, N⁻¹vₖ₋₁, N⁻¹vₖ, p, wₖ₋₂, wₖ₋₁ + # - 6 m-vectors: r, y, M⁻¹uₖ₋₁, M⁻¹uₖ, q, d̅ + storage_usymlqr(m, n) = 7 * n + 6 * m + storage_usymlqr_bytes(m, n) = nbits_FC * storage_usymlqr(m, n) + + expected_usymlqr_bytes = storage_usymlqr_bytes(k, n) + usymlqr(Au, c, b) # warmup + actual_usymlqr_bytes = @allocated usymlqr(Au, c, b) + @test expected_usymlqr_bytes ≤ actual_usymlqr_bytes ≤ 1.02 * expected_usymlqr_bytes + + solver = UsymlqrSolver(Au, c) + usymlqr!(solver, Au, c, b) # warmup + inplace_usymlqr_bytes = @allocated usymlqr!(solver, Au, c, b) + @test inplace_usymlqr_bytes == 0 + end + @testset "GPMR" begin # GPMR needs: # - 2 m-vectors: x, q diff --git a/test/test_mp.jl b/test/test_mp.jl index 3cf9e01d5..6d5d0b1bb 100644 --- a/test/test_mp.jl +++ b/test/test_mp.jl @@ -16,7 +16,7 @@ x, _ = @eval $fn($A, $b, $c) elseif fn in (:trilqr, :bilqr) x, t, _ = @eval $fn($A, $b, $c) - elseif fn in (:tricg, :trimr) + elseif fn in (:tricg, :trimr, :usymlqr) x, y, _ = @eval $fn($A, $b, $c) elseif fn == :gpmr x, y, _ = @eval $fn($A, $B, $b, $c) @@ -34,6 +34,10 @@ @test norm(x + A * y - b) ≤ Κ * (atol + norm([b; c]) * rtol) @test norm(A' * x - y - c) ≤ Κ * (atol + norm([b; c]) * rtol) @test eltype(y) == FC + if fn == :usymlqr + @test norm(x + A * y - b) ≤ Κ * (atol + norm([b; c]) * rtol) + @test norm(A' * x - c) ≤ Κ * (atol + norm([b; c]) * rtol) + @test eltype(y) == FC elseif fn == :gpmr @test norm(x + A * y - b) ≤ Κ * (atol + norm([b; c]) * rtol) @test norm(B * x + y - c) ≤ Κ * (atol + norm([b; c]) * rtol) diff --git a/test/test_solvers.jl b/test/test_solvers.jl index 6c49502be..64ab17715 100644 --- a/test/test_solvers.jl +++ b/test/test_solvers.jl @@ -49,6 +49,7 @@ function test_solvers(FC) $solvers[:usymlq] = $(KRYLOV_SOLVERS[:usymlq])($m, $n, $S) $solvers[:tricg] = $(KRYLOV_SOLVERS[:tricg])($m, $n, $S) $solvers[:trimr] = $(KRYLOV_SOLVERS[:trimr])($m, $n, $S) + $solvers[:usymlqr] = $(KRYLOV_SOLVERS[:usymlqr])($m, $n, $S) $solvers[:gpmr] = $(KRYLOV_SOLVERS[:gpmr])($n, $m, $mem, $S) $solvers[:cg_lanczos_shift] = $(KRYLOV_SOLVERS[:cg_lanczos_shift])($n, $n, $nshifts, $S) end @@ -74,7 +75,7 @@ function test_solvers(FC) method == :cgls_lanczos_shift && @test_throws ErrorException("solver.nshifts = $(solver.nshifts) is inconsistent with length(shifts) = $(length(shifts2))") solve!(solver, Ao, b, shifts2) method ∈ (:bilqr, :trilqr) && @test_throws ErrorException("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($n2, $n2)") solve!(solver, A2, b2, b2) method == :gpmr && @test_throws ErrorException("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($n2, $m2)") solve!(solver, Ao2, Au2, b2, c2) - method ∈ (:tricg, :trimr) && @test_throws ErrorException("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($n2, $m2)") solve!(solver, Ao2, b2, c2) + method ∈ (:tricg, :trimr, :usymlqr) && @test_throws ErrorException("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($n2, $m2)") solve!(solver, Ao2, b2, c2) method == :usymlq && @test_throws ErrorException("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($m2, $n2)") solve!(solver, Au2, c2, b2) method == :usymqr && @test_throws ErrorException("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($n2, $m2)") solve!(solver, Ao2, b2, c2) end @@ -90,7 +91,7 @@ function test_solvers(FC) method == :cgls_lanczos_shift && solve!(solver, Ao, b, shifts, timemax=timemax) method ∈ (:bilqr, :trilqr) && solve!(solver, A, b, b, timemax=timemax) method == :gpmr && solve!(solver, Ao, Au, b, c, timemax=timemax) - method ∈ (:tricg, :trimr) && solve!(solver, Au, c, b, timemax=timemax) + method ∈ (:tricg, :trimr, :usymlqr) && solve!(solver, Au, c, b, timemax=timemax) method == :usymlq && solve!(solver, Au, c, b, timemax=timemax) method == :usymqr && solve!(solver, Ao, b, c, timemax=timemax) @test solver.stats.status == "time limit exceeded" @@ -147,7 +148,7 @@ function test_solvers(FC) @test issolved_dual(solver) end - if method ∈ (:tricg, :trimr, :gpmr) + if method ∈ (:tricg, :trimr, :usymlqr, :gpmr) method == :gpmr ? solve!(solver, Ao, Au, b, c) : solve!(solver, Au, c, b) niter = niterations(solver) @test Aprod(solver) == niter diff --git a/test/test_usymlqr.jl b/test/test_usymlqr.jl new file mode 100644 index 000000000..0fad97ef2 --- /dev/null +++ b/test/test_usymlqr.jl @@ -0,0 +1,29 @@ +@testset "usymlqr" begin + usymlqr_tol = 1.0e-6 + + # Test saddle-point systems + A, b, D = saddle_point() + m, n = size(A) + c = -b + D⁻¹ = sparse(inv(D)) + N⁻¹ = eye(n) + H⁻¹ = blockdiag(D⁻¹, N⁻¹) + + (x, y, stats) = usymlqr(A, b, c, M=D⁻¹) + K = [D A; A' zeros(n, n)] + B = [b; c] + r = B - K * [x; y] + resid = sqrt(dot(r, H⁻¹ * r)) / sqrt(dot(B, H⁻¹ * B)) + @printf("USYMLQR: Relative residual: %8.1e\n", resid) + @test(resid ≤ usymlqr_tol) + + (x, y, stats) = usymlqr(A, b, c) + K = [eye(m) A; A' zeros(n, n)] + B = [b; c] + r = B - K * [x; y] + resid = norm(r) / norm(B) + @printf("USYMLQR: Relative residual: %8.1e\n", resid) + @test(resid ≤ usymlqr_tol) +end + +test_usymlqr() diff --git a/test/test_verbose.jl b/test/test_verbose.jl index 06856d22a..84909ab87 100644 --- a/test/test_verbose.jl +++ b/test/test_verbose.jl @@ -15,13 +15,13 @@ function test_verbose(FC) for fn in (:cg, :cgls, :usymqr, :cgne, :cgs, :crmr, :cg_lanczos, :dqgmres, :diom, :cr, :gpmr, :lslq, :lsqr, :lsmr, :lnlq, :craig, :bicgstab, :craigmr, :crls, :symmlq, :minres, :bilq, :minres_qlp, :qmr, :usymlq, :tricg, :trimr, :trilqr, :bilqr, :gmres, :fom, - :car, :minares, :fgmres, :cg_lanczos_shift, :cgls_lanczos_shift) + :car, :minares, :fgmres, :usymlqr :cg_lanczos_shift, :cgls_lanczos_shift) @testset "$fn" begin io = IOBuffer() if fn in (:trilqr, :bilqr) @eval $fn($A, $b, $b, verbose=1, iostream=$io) - elseif fn in (:tricg, :trimr) + elseif fn in (:tricg, :trimr, :usymlqr) @eval $fn($Au, $c, $b, verbose=1, iostream=$io) elseif fn in (:lnlq, :craig, :craigmr, :cgne, :crmr) @eval $fn($Au, $c, verbose=1, iostream=$io) diff --git a/test/test_warm_start.jl b/test/test_warm_start.jl index 47508a63a..e79ba3f61 100644 --- a/test/test_warm_start.jl +++ b/test/test_warm_start.jl @@ -68,6 +68,18 @@ function test_warm_start(FC) resid = norm(r) / norm([b; b]) @test(resid ≤ tol) + # USYMLQR + x, y, stats = usymlqr(A, b, b, x0, y0) + r = [b - x - A * y; b - A' * x] + resid = norm(r) / norm([b; b]) + @test(resid ≤ tol) + + solver = UsymlqrSolver(A, b) + solve!(solver, A, b, b, x0, y0) + r = [b - solver.x - A * solver.y; b - A' * solver.x] + resid = norm(r) / norm([b; b]) + @test(resid ≤ tol) + # GPMR x, y, stats = gpmr(A, A', b, b, x0, y0) r = [b - x - A * y; b - A' * x - y]