Add an implementation of cgls-shift (#853)

Author: tmigot

docs/make.jl

@@ -46,6 +46,7 @@ makedocs(
                            "CRAIG" => "examples/craig.md",
                            "CRAIGMR" => "examples/craigmr.md",
                            "CGLS" => "examples/cgls.md",
+                           "CGLS-LANCZOS-SHIFT" => "examples/cgls_lanczos_shift.md",
                            "CRLS" => "examples/crls.md",
                            "LSQR" => "examples/lsqr.md",
                            "LSMR" => "examples/lsmr.md"],

docs/src/api.md

@@ -41,6 +41,7 @@ BilqSolver
 QmrSolver
 BilqrSolver
 CglsSolver
+CglsLanczosShiftSolver
 CrlsSolver
 CgneSolver
 CrmrSolver

docs/src/examples/cgls_lanczos_shift.md

@@ -0,0 +1,33 @@
+```@example cgls_lanczos_shift
+using MatrixMarket, SuiteSparseMatrixCollection
+using Krylov, LinearOperators
+using LinearAlgebra, Printf
+
+function residuals(A, b, shifts, x)
+  nshifts = length(shifts)
+  r = [ A' * (A * x[i] - b) + shifts[i] * x[i] for i = 1 : nshifts ]
+  return r
+end
+ssmc = ssmc_db(verbose=false)
+matrix = ssmc_matrices(ssmc, "HB", "well1033")
+path = fetch_ssmc(matrix, format="MM")
+
+A = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1]).mtx"))
+b = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1])_b.mtx"))[:]
+(m, n) = size(A)
+@printf("System size: %d rows and %d columns\n", m, n)
+
+# Define regularization parameters.
+shifts = [1.0, 2.0, 3.0, 4.0]
+
+(x, stats) = cgls_lanczos_shift(A, b, shifts)
+show(stats)
+r = residuals(A, b, shifts, x)
+
+resids = map(norm, r) / norm(b)
+@printf("CGLS: Relative residuals with shifts:\n")
+for resid in resids
+  @printf(" %8.1e", resid)
+end
+@printf("\n")
+```

docs/src/inplace.md

@@ -16,7 +16,7 @@ Given an operator `A` and a right-hand side `b`, you can create a `KrylovSolver`
 For example, use `S = Vector{Float64}` if you want to solve linear systems in double precision on the CPU and `S = CuVector{Float32}` if you want to solve linear systems in single precision on an Nvidia GPU.
 
 !!! note
-    `DiomSolver`, `FomSolver`, `DqgmresSolver`, `GmresSolver`, `BlockGmresSolver`, `FgmresSolver`, `GpmrSolver` and `CgLanczosShiftSolver` require an additional argument (`memory` or `nshifts`).
+    `DiomSolver`, `FomSolver`, `DqgmresSolver`, `GmresSolver`, `BlockGmresSolver`, `FgmresSolver`, `GpmrSolver`, `CgLanczosShiftSolver` and `CglsLanczosShiftSolver` require an additional argument (`memory` or `nshifts`).
 
 The workspace is always the first argument of the in-place methods:
 

docs/src/matrix_free.md

@@ -41,6 +41,7 @@ Some methods only require `A * v` products, whereas other ones also require `A'
 | 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                             |
+| CG-LANCZOS-SHIFT                                | CGLS-LANCZOS-SHIFT                       |
 
 !!! info
     GPMR is the only method that requires `A * v` and `B * w` products.

docs/src/preconditioners.md

@@ -61,7 +61,7 @@ However, there is no need to specify $L$ and one may specify $P_c = LL^H$ or its
 
 ### Linear least-squares problems
 
-Methods concerned: [`CGLS`](@ref cgls), [`CRLS`](@ref crls), [`LSLQ`](@ref lslq), [`LSQR`](@ref lsqr) and [`LSMR`](@ref lsmr).
+Methods concerned: [`CGLS`](@ref cgls), [`CGLS-LANCZOS-SHIFT`](@ref cgls_lanczos_shift), [`CRLS`](@ref crls), [`LSLQ`](@ref lslq), [`LSQR`](@ref lsqr) and [`LSMR`](@ref lsmr).
 
 | Formulation           | Without preconditioning              | With preconditioning                        |
 |:---------------------:|:------------------------------------:|:-------------------------------------------:|

docs/src/processes.md

@@ -91,7 +91,7 @@ Note that depending on how we normalize the vectors that compose $V_k$, $T_{k+1,
 
 The function [`hermitian_lanczos`](@ref hermitian_lanczos(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)) returns $V_{k+1}$, $\beta_1$ and $T_{k+1,k}$.
 
-Related methods: [`SYMMLQ`](@ref symmlq), [`CG`](@ref cg), [`CR`](@ref cr), [`CAR`](@ref car), [`MINRES`](@ref minres), [`MINRES-QLP`](@ref minres_qlp), [`MINARES`](@ref minares), [`CGLS`](@ref cgls), [`CRLS`](@ref crls), [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`CG-LANCZOS`](@ref cg_lanczos) and [`CG-LANCZOS-SHIFT`](@ref cg_lanczos_shift).
+Related methods: [`SYMMLQ`](@ref symmlq), [`CG`](@ref cg), [`CR`](@ref cr), [`CAR`](@ref car), [`MINRES`](@ref minres), [`MINRES-QLP`](@ref minres_qlp), [`MINARES`](@ref minares), [`CGLS`](@ref cgls), [`CGLS-LANCZOS-SHIFT`](@ref cgls_lanczos_shift), [`CRLS`](@ref crls), [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`CG-LANCZOS`](@ref cg_lanczos) and [`CG-LANCZOS-SHIFT`](@ref cg_lanczos_shift).
 
 ```@docs
 hermitian_lanczos(::Any, ::AbstractVector{FC}, ::Int) where FC <: (Union{Complex{T}, T} where T <: AbstractFloat)

docs/src/solvers/ls.md

@@ -9,6 +9,13 @@ cgls
 cgls!
 ```
 
+## CGLS-LANCZOS-SHIFT
+
+```@docs
+cgls_lanczos_shift
+cgls_lanczos_shift!
+```
+
 ## CRLS
 
 ```@docs

docs/src/storage.md

@@ -39,7 +39,7 @@ This section provides the storage requirements of all Krylov methods available i
 
 We denote by $m$ and $n$ the number of rows and columns of the linear problem.
 The memory parameter of DIOM, FOM, DQGMRES, GMRES, FGMRES and GPMR is $k$.
-The numbers of shifts of CG-LANCZOS-SHIFT is $p$.
+The numbers of shifts of CG-LANCZOS-SHIFT and CGLS-LANCZOS-SHIFT is $p$.
 
 ## Theoretical storage requirements
 
@@ -81,9 +81,9 @@ Each table summarizes the storage requirements of Krylov methods recommended to
 
 #### Least-squares problems
 
-| Methods | [`USYMQR`](@ref usymqr) | [`CGLS`](@ref cgls) | [`CRLS`](@ref crls) | [`LSLQ`](@ref lslq) | [`LSQR`](@ref lsqr) | [`LSMR`](@ref lsmr) |
-|:-------:|:-----------------------:|:-------------------:|:-------------------:|:-------------------:|:-------------------:|:-------------------:|
-| Storage | $6n + 3m$               | $3n + 2m$           | $4n + 3m$           | $4n + 2m$           | $4n + 2m$           | $5n + 2m$           |
+| Methods | [`USYMQR`](@ref usymqr) | [`CGLS`](@ref cgls) | [`CG-LANCZOS-SHIFT`](@ref cg_lanczos_shift) | [`CRLS`](@ref crls) | [`LSLQ`](@ref lslq) | [`LSQR`](@ref lsqr) | [`LSMR`](@ref lsmr) |
+|:-------:|:-----------------------:|:-------------------:|:-------------------------------------------:|:-------------------:|:-------------------:|:-------------------:|
+| Storage | $6n + 3m$               | $3n + 2m$           | $3n + 2m + 5p + 2np$                        | $4n + 3m$           | $4n + 2m$           | $4n + 2m$           | $5n + 2m$           |
 
 #### Adjoint systems
 

src/Krylov.jl

@@ -20,7 +20,6 @@ include("car.jl")
 
 include("symmlq.jl")
 include("cg_lanczos.jl")
-include("cg_lanczos_shift.jl")
 include("minres.jl")
 include("minres_qlp.jl")
 include("minares.jl")
@@ -60,5 +59,8 @@ include("lnlq.jl")
 include("craig.jl")
 include("craigmr.jl")
 
+include("cg_lanczos_shift.jl")
+include("cgls_lanczos_shift.jl")
+
 include("krylov_solve.jl")
 end