diff --git a/docs/make.jl b/docs/make.jl index 9668d815b..8b1bddcd1 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 a256d4b0f..c4e63eae5 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..5fcbb19d6 --- /dev/null +++ b/docs/src/examples/usymlqr.md @@ -0,0 +1,27 @@ +```@example usymlqr +using Krylov, LinearOperators, LDLFactorizations +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] +llt_D = cholesky(D) +opD⁻¹ = LinearOperator(Float64, 5, 5, true, true, (y, v) -> ldiv!(y, llt_D, v)) +opH⁻¹ = BlockDiagonalOperator(opD⁻¹, eye(n)) +(x, y, stats) = usymlqr(A, b, c, M=opD⁻¹, sp=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 b1c639595..1efa17ffe 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 | !!! info GPMR is the only method that requires `A * v` and `B * w` products. diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md index 8030f9697..df9b47893 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 7452b0577..1ff338468 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 fcd74abbb..87143e7e3 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$ | $Xn + Ym$ | #### Generalized saddle-point and non-Hermitian partitioned systems diff --git a/src/Krylov.jl b/src/Krylov.jl index 760ad6074..3b03e793c 100644 --- a/src/Krylov.jl +++ b/src/Krylov.jl @@ -35,6 +35,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/krylov_solve.jl b/src/krylov_solve.jl index 2940bc7e8..28aa46b08 100644 --- a/src/krylov_solve.jl +++ b/src/krylov_solve.jl @@ -41,6 +41,7 @@ for (KS, fun, args, def_args, optargs, def_optargs, kwargs, def_kwargs) in [ (: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 ) ] @eval begin solve!(solver :: $KS{T,FC,S}, $(def_args...); $(def_kwargs...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}, S <: AbstractVector{FC}} = $(fun)(solver, $(args...); $(kwargs...)) diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl index d9f9fd9b1..ae9412132 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, 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! @@ -46,20 +46,21 @@ const KRYLOV_SOLVERS = Dict( :lnlq => :LnlqSolver , :craig => :CraigSolver , :craigmr => :CraigmrSolver , + :usymlqr => :UsymlqrSolver , ) "Abstract type for using Krylov solvers in-place" 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 @@ -101,14 +102,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 @@ -150,14 +151,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 @@ -193,14 +194,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 @@ -238,14 +239,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 @@ -287,14 +288,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 @@ -338,14 +339,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 @@ -383,14 +384,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 @@ -441,14 +442,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 @@ -488,14 +489,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`. """ mutable struct DqgmresSolver{T,FC,S} <: KrylovSolver{T,FC,S} @@ -541,14 +542,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`. """ mutable struct DiomSolver{T,FC,S} <: KrylovSolver{T,FC,S} @@ -592,14 +593,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 @@ -641,14 +642,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 @@ -692,14 +693,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 @@ -755,14 +756,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 @@ -826,14 +827,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 @@ -883,14 +884,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 @@ -932,14 +933,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 @@ -981,14 +982,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 @@ -1034,14 +1035,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 @@ -1089,14 +1090,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 @@ -1146,14 +1147,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 @@ -1188,14 +1189,14 @@ function CglsSolver(A, b) 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 @@ -1234,14 +1235,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 @@ -1278,14 +1279,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 @@ -1322,14 +1323,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 @@ -1370,14 +1371,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 @@ -1418,14 +1419,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 @@ -1468,14 +1469,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 @@ -1518,14 +1519,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 @@ -1568,14 +1569,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 @@ -1622,14 +1623,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`. """ mutable struct GmresSolver{T,FC,S} <: KrylovSolver{T,FC,S} @@ -1676,14 +1677,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`. """ mutable struct FgmresSolver{T,FC,S} <: KrylovSolver{T,FC,S} @@ -1730,14 +1731,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`. """ mutable struct FomSolver{T,FC,S} <: KrylovSolver{T,FC,S} @@ -1781,14 +1782,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`. """ mutable struct GpmrSolver{T,FC,S} <: KrylovSolver{T,FC,S} @@ -1845,6 +1846,57 @@ 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 + x :: S + y :: S + M⁻¹vₖ₋₁ :: S + M⁻¹vₖ :: S + N⁻¹uₖ₋₁ :: S + N⁻¹uₖ :: S + p :: S + q :: S + vₖ :: S + uₖ :: S + warm_start :: Bool + stats :: SimpleStats{T} +end + +function UsymlqrSolver(m, n, S) + FC = eltype(S) + T = real(FC) + x = S(undef, m) + y = S(undef, n) + M⁻¹vₖ₋₁ = S(undef, m) + M⁻¹vₖ = S(undef, m) + N⁻¹uₖ₋₁ = S(undef, n) + N⁻¹uₖ = S(undef, n) + p = S(undef, n) + q = S(undef, m) + vₖ = S(undef, 0) + uₖ = S(undef, 0) + stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") + solver = UsymlqrSolver{T,FC,S}(m, n, x, y, M⁻¹vₖ₋₁, M⁻¹vₖ, N⁻¹uₖ₋₁, N⁻¹uₖ, p, q, vₖ, uₖ, false, stats) + return solver +end + +function UsymlqrSolver(A, b) + m, n = size(A) + S = ktypeof(b) + UsymlqrSolver(m, n, S) +end + """ solution(solver) @@ -1932,6 +1984,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 ) ] @eval begin size(solver :: $KS) = solver.m, solver.n @@ -1959,7 +2012,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/usymlqr.jl b/src/usymlqr.jl new file mode 100644 index 000000000..29e854e98 --- /dev/null +++ b/src/usymlqr.jl @@ -0,0 +1,453 @@ +# 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 = (:(; 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 = mapreduce(extract_parameters, vcat, def_kwargs_usymlqr) + +args_usymlqr = (:A, :b, :c) +optargs_usymlqr = (:x0, :y0) +kwargs_usymlqr = (:M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) + +@eval begin + function usymlqr($(def_args_usymlqr...), $(def_optargs_usymlqr...); $(def_kwargs_usymlqr...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}} + start_time = time_ns() + solver = UsymlqrSolver(A, b) + warm_start!(solver, $(optargs_usymlqr...)) + elapsed_time = ktimer(start_time) + timemax -= elapsed_time + usymlqr!(solver, $(args_usymlqr...); $(kwargs_usymlqr...)) + solver.stats.timer += elapsed_time + return (solver.x, solver.y, solver.stats) + end + + function usymlqr($(def_args_usymlqr...); $(def_kwargs_usymlqr...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}} + start_time = time_ns() + solver = UsymlqrSolver(A, b) + elapsed_time = ktimer(start_time) + timemax -= elapsed_time + usymlqr!(solver, $(args_usymlqr...); $(kwargs_usymlqr...)) + solver.stats.timer += elapsed_time + return (solver.x, solver.y, solver.stats) + end + + 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, :vₖ, S, m) + allocate_if(!NisI, solver, :uₖ, S, n) + M⁻¹vₖ₋₁, M⁻¹vₖ, N⁻¹uₖ₋₁, N⁻¹uₖ = solver.M⁻¹vₖ₋₁, solver.M⁻¹vₖ, solver.N⁻¹uₖ₋₁, solver.N⁻¹uₖ + xₖ, yₖ, p, q = solver.x, solver.y, solver.p, solver.q + vₖ = MisI ? M⁻¹vₖ : solver.vₖ + uₖ = NisI ? N⁻¹uₖ : solver.uₖ + 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 x₀ and y₀. + warm_start && (Δx .= xₖ) + warm_start && (Δy .= yₖ) + xₖ .= zero(T) + yₖ .= zero(T) + + # Initialize preconditioned orthogonal tridiagonalization process. + M⁻¹vₖ₋₁ .= zero(T) # v₀ = 0 + N⁻¹uₖ₋₁ .= zero(T) # u₀ = 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 + + # β₁Ev₁ = b ↔ β₁v₁ = Mb + M⁻¹vₖ .= b₀ + MisI || mul!(vₖ, M, M⁻¹vₖ) + βₖ = sqrt(@kdot(m, vₖ, M⁻¹vₖ)) # β₁ = ‖v₁‖_E + if βₖ ≠ 0 + @kscal!(m, 1 / βₖ, M⁻¹vₖ) + MisI || @kscal!(m, 1 / βₖ, vₖ) + else + error("b must be nonzero") + end + + # γ₁Fu₁ = c ↔ γ₁u₁ = Nc + N⁻¹uₖ .= c₀ + NisI || mul!(uₖ, N, N⁻¹uₖ) + γₖ = sqrt(@kdot(n, uₖ, N⁻¹uₖ)) # γ₁ = ‖u₁‖_F + if γₖ ≠ 0 + @kscal!(n, 1 / γₖ, N⁻¹uₖ) + NisI || @kscal!(n, 1 / γₖ, uₖ) + else + error("c must be nonzero") + end + + (verbose > 0) && @printf(iostream, "%4s %7s %7s %7s\n", "k", "αₖ", "βₖ", "γₖ") + kdisplay(iter, verbose) && @printf(iostream, "%4d %7.1e %7.1e %7.1e\n", iter, αₖ, βₖ, γₖ) + + # Stopping criterion. + solved = rNorm ≤ ε + tired = iter ≥ itmax + status = "unknown" + ill_cond = false + user_requested_exit = false + overtimed = false + + while !(solved || tired || ill_cond || user_requested_exit || overtimed) + # Update iteration index. + iter = iter + 1 + + # Continue the orthogonal tridiagonalization process. + # AUₖ = EVₖTₖ + βₖ₊₁Evₖ₊₁(eₖ)ᵀ = EVₖ₊₁Tₖ₊₁.ₖ + # AᴴVₖ = FUₖ(Tₖ)ᴴ + γₖ₊₁Fuₖ₊₁(eₖ)ᵀ = FUₖ₊₁(Tₖ.ₖ₊₁)ᴴ + + mul!(q, A , uₖ) # Forms Evₖ₊₁ : q ← Auₖ + mul!(p, Aᴴ, vₖ) # Forms Fuₖ₊₁ : p ← Aᴴvₖ + + if iter ≥ 2 + @kaxpy!(m, -γₖ, M⁻¹vₖ₋₁, q) # q ← q - γₖ * M⁻¹vₖ₋₁ + @kaxpy!(n, -βₖ, N⁻¹uₖ₋₁, p) # p ← p - βₖ * N⁻¹uₖ₋₁ + end + + αₖ = @kdot(m, vₖ, q) # αₖ = qᴴvₖ + + @kaxpy!(m, -αₖ, M⁻¹vₖ, q) # q ← q - αₖ * M⁻¹vₖ + @kaxpy!(n, -αₖ, N⁻¹uₖ, p) # p ← p - αₖ * N⁻¹uₖ + + # Compute vₖ₊₁ and uₖ₊₁ + MisI || mul!(vₖ₊₁, M, q) # βₖ₊₁vₖ₊₁ = MAuₖ - γₖvₖ₋₁ - αₖvₖ + NisI || mul!(uₖ₊₁, N, p) # γₖ₊₁uₖ₊₁ = NAᴴvₖ - βₖuₖ₋₁ - αₖuₖ + + βₖ₊₁ = sqrt(@kdot(m, vₖ₊₁, q)) # βₖ₊₁ = ‖vₖ₊₁‖_E + γₖ₊₁ = sqrt(@kdot(n, uₖ₊₁, p)) # γₖ₊₁ = ‖uₖ₊₁‖_F + + if βₖ₊₁ ≠ 0 + @kscal!(m, one(T) / βₖ₊₁, q) + MisI || @kscal!(m, one(T) / βₖ₊₁, vₖ₊₁) + end + + if γₖ₊₁ ≠ 0 + @kscal!(n, one(T) / γₖ₊₁, p) + NisI || @kscal!(n, one(T) / γₖ₊₁, uₖ₊₁) + end + + # Continue the QR factorization of Tₖ₊₁.ₖ = Qₖ₊₁ [ Rₖ ]. + # [ Oᴴ ] + + ƛ = -cs * γₖ + ϵ = sn * γₖ + + if !solved_LS + ArNorm_qr_computed = rNorm_qr * sqrt(δbar^2 + ƛ^2) + ArNorm_qr = norm(A' * (b - A * x)) # FIXME + @debug "" ArNorm_qr_computed ArNorm_qr abs(ArNorm_qr_computed - ArNorm_qr) / ArNorm_qr + ArNorm_qr = ArNorm_qr_computed + push!(ArNorms_qr, ArNorm_qr) + + test_LS = ArNorm_qr / (Anorm * max(one(T), rNorm_qr)) + solved_lim_LS = test_LS ≤ ls_optimality_tol + solved_mach_LS = one(T) + test_LS ≤ one(T) + solved_LS = solved_mach_LS | solved_lim_LS + end + kdisplay(iter, verbose) && @printf(iostream, "%7.1e ", ArNorm_qr) + + # continue QR factorization + delta = sqrt(δbar^2 + βₖ^2) + csold = cs + snold = sn + cs = δbar/ delta + sn = βₖ / delta + + # update w (used to update x and z) + @. wold = w + @. w = cs * wbar + + if !solved_LS + # the optimality conditions of the LS problem were not triggerred + # update x and see if we have a zero residual + + ϕ = cs * ϕbar + ϕbar = sn * ϕbar + @kaxpy!(n, ϕ, w, x) + xNorm = norm(x) # FIXME + + # update least-squares residual + rNorm_qr = abs(ϕbar) + push!(rNorms_qr, rNorm_qr) + + # stopping conditions related to the least-squares problem + test_LS = rNorm_qr / (one(T) + Anorm * xNorm) + zero_resid_lim_LS = test_LS ≤ ls_zero_resid_tol + zero_resid_mach_LS = one(T) + test_LS ≤ one(T) + zero_resid_LS = zero_resid_mach_LS | zero_resid_lim_LS + solved_LS |= zero_resid_LS + + end + + # continue tridiagonalization + q = A * vₖ + @. q -= γₖ * u_prev + αₖ = @kdot(m, uₖ, q) + + # Update norm estimates + Anorm2 += αₖ * αₖ + βₖ * βₖ + γₖ * γₖ + Anorm = √Anorm2 + + # Estimate κ₂(A) based on the diagonal of L. + sigma_min = min(delta, sigma_min) + sigma_max = max(delta, sigma_max) + Acond = sigma_max / sigma_min + + # continue QR factorization of T{k+1,k} + λ = cs * ƛ + sn * αₖ + δbar= sn * ƛ - cs * αₖ + + if !solved_LN + + etaold = η + η = cs * etabar # = etak + + # compute residual of least-norm problem at y{k-1} + # TODO: use recurrence formula for LQ residual + rNorm_lq_computed = sqrt((delta * η)^2 + (ϵ * etaold)^2) + rNorm_lq = norm(A' * y - c) # FIXME + rNorm_lq = rNorm_lq_computed + push!(rNorms_lq, rNorm_lq) + + # stopping conditions related to the least-norm problem + test_LN = rNorm_lq / sqrt(cnorm^2 + Anorm2 * yNorm2) + solved_lim_LN = test_LN ≤ ln_tol + solved_mach_LN = one(T) + test_LN ≤ one(T) + solved_LN = solved_lim_LN || solved_mach_LN + + # TODO: remove this when finished + push!(tests_LN, test_LN) + + @. wbar = (vₖ - λ * w - ϵ * wold) / δbar + + if !solved_LN + + # prepare to update y and z + @. p = cs * pbar + sn * uₖ + + # update y and z + @. y += η * p + @. z -= η * w + yNorm2 += η * η + yNorm = sqrt(yNorm2) + + @. pbar = sn * pbar - cs * uₖ + etabarold = etabar + etabar = -(λ * η + ϵ * etaold) / δbar # = etabar{k+1} + + # see if CG iterate has smaller residual + # TODO: use recurrence formula for CG residual + @. yC = y + etabar * pbar + @. zC = z - etabar * wbar + yCNorm2 = yNorm2 + etabar* etabar + rNorm_cg_computed = γₖ * abs(snold * etaold - csold * etabarold) + rNorm_cg = norm(A' * yC - c) + + # if rNorm_cg < rNorm_lq + # # stopping conditions related to the least-norm problem + # test_cg = rNorm_cg / sqrt(γ₁^2 + Anorm2 * yCNorm2) + # solved_lim_LN = test_cg ≤ ln_tol + # solved_mach_LN = 1.0 + test_cg ≤ 1.0 + # solved_LN = solved_lim_LN | solved_mach_LN + # # transition_to_cg = solved_LN + # transition_to_cg = false + # end + + if transition_to_cg + # @. yC = y + etabar* pbar + # @. zC = z - etabar* wbar + end + end + end + 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) + + # Stopping conditions that apply to both problems + 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 + user_requested_exit = callback(solver) :: Bool + timer = time_ns() - start_time + overtimed = timer > timemax_ns + end + (verbose > 0) && @printf(iostream, "\n") + + # Update status + tired && (status = "maximum number of iterations exceeded") + 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") + + # Update x and y + warm_start && @kaxpy!(m, one(FC), Δx, xₖ) + warm_start && @kaxpy!(n, one(FC), Δy, yₖ) + + # 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/test/runtests.jl b/test/runtests.jl index 040d2edd1..e48b3a82b 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 b9e923c61..84d9fc371 100644 --- a/test/test_allocations.jl +++ b/test/test_allocations.jl @@ -595,6 +595,24 @@ @test inplace_trimr_bytes == 0 end + @testset "USYMLQR" begin + # USYMLQR needs: + # - _ n-vectors: ... + # - _ m-vectors: ... + storage_usymlqr(m, n) = 0 * n + 0 * 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 3e329f1b1..6f813cfb5 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 993ed8bda..f56ac25c4 100644 --- a/test/test_solvers.jl +++ b/test/test_solvers.jl @@ -48,6 +48,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 @@ -71,7 +72,7 @@ function test_solvers(FC) method ∈ (:cgls, :crls, :lslq, :lsqr, :lsmr) && @test_throws ErrorException("(solver.m, solver.n) = ($(solver.m), $(solver.n)) is inconsistent with size(A) = ($n2, $m2)") solve!(solver, Ao2, b2) 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 @@ -86,7 +87,7 @@ function test_solvers(FC) method ∈ (:cgls, :crls, :lslq, :lsqr, :lsmr) && solve!(solver, Ao, b, 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" @@ -143,7 +144,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 d688def16..c97047520 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) + :car, :minares, :fgmres, :cg_lanczos_shift, :usymlqr) @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]