Add an implementation of USYMLQR

Author: amontoison

docs/make.jl

@@ -41,6 +41,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",

docs/src/api.md

@@ -31,6 +31,7 @@ DqgmresWorkspace
 GmresWorkspace
 UsymlqWorkspace
 UsymqrWorkspace
+UsymlqrWorkspace
 TricgWorkspace
 TrimrWorkspace
 TrilqrWorkspace

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 = -3*b
+
+# [D   A] [x] = [b]
+# [Aᴴ  0] [y]   [c]
+opH = BlockDiagonalOperator(inv(D), eye(n))
+(x, y, stats) = usymlqr(A, b, c, M=inv(D))
+K = [D A; A' zeros(n,n)]
+B = [b; c]
+r = B - K * [x; y]
+resid = sqrt(dot(r, opH * r))
+@printf("USYMLQR: Relative residual: %8.1e\n", resid)
+```

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

docs/src/preconditioners.md

@@ -111,14 +111,16 @@ 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},
 ```
 | Preconditioners | $E^{-1}$              | $E$                  | $F^{-1}$              | $F$                  |
 |:---------------:|:---------------------:|:--------------------:|:---------------------:|:--------------------:|
 | Arguments       | `M` with `ldiv=false` | `M` with `ldiv=true` | `N` with `ldiv=false` | `N` with `ldiv=true` |
 
+In the special case of [USYMLQR](@ref usymlqr), $\tau = 1$ and $\nu = 0$.
+
 !!! warning
     The preconditioners `M` and `N` must be hermitian and positive definite.
 

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}$.

docs/src/solvers/sp_sqd.md

@@ -15,3 +15,10 @@ tricg!
 trimr
 trimr!
 ```
+
+## USYMLQR
+
+```@docs
+usymlqr!
+usymlqr
+```

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
 

src/Krylov.jl

@@ -40,6 +40,7 @@ include("gpmr.jl")
 
 include("usymlq.jl")
 include("usymqr.jl")
+include("usymlqr.jl")
 include("tricg.jl")
 include("trimr.jl")
 include("trilqr.jl")

src/interface.jl

@@ -97,6 +97,7 @@ for (workspace, krylov, args, def_args, optargs, def_optargs, kwargs, def_kwargs
   (:FgmresWorkspace   , :fgmres    , args_fgmres    , def_args_fgmres    , optargs_fgmres    , def_optargs_fgmres    , kwargs_fgmres    , def_kwargs_fgmres    )
   (:FomWorkspace      , :fom       , args_fom       , def_args_fom       , optargs_fom       , def_optargs_fom       , kwargs_fom       , def_kwargs_fom       )
   (:GpmrWorkspace     , :gpmr      , args_gpmr      , def_args_gpmr      , optargs_gpmr      , def_optargs_gpmr      , kwargs_gpmr      , def_kwargs_gpmr      )
+  (:UsymlqrWorkspace  , :usymlqr   , args_usymlqr   , def_args_usymlqr   , optargs_usymlqr   , def_optargs_usymlqr   , kwargs_usymlqr   , def_kwargs_usymlqr   )
   (:CgLanczosShiftWorkspace  , :cg_lanczos_shift  , args_cg_lanczos_shift  , def_args_cg_lanczos_shift  , (), (), kwargs_cg_lanczos_shift  , def_kwargs_cg_lanczos_shift  )
   (:CglsLanczosShiftWorkspace, :cgls_lanczos_shift, args_cgls_lanczos_shift, def_args_cgls_lanczos_shift, (), (), kwargs_cgls_lanczos_shift, def_kwargs_cgls_lanczos_shift)
 ]