Skip to content

[USYMLQR] Try to support elliptic norms #1068

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
amontoison wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

[USYMLQR] Try to support elliptic norms #1068

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 11 additions & 1 deletion docs/src/examples/usymlqr.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
```@example usymlqr
using LinearAlgebra, Printf, SparseArrays
using Krylov
using Krylov, LinearOperators

# Identity matrix.
eye(n::Int) = sparse(1.0 * I, n, n)
Expand Down Expand Up @@ -39,4 +39,14 @@ d = [0*b; c]
r = d - K * [x; y]
resid = norm(r)
@printf("USYMLQR: Relative residual: %8.1e\n", resid)

# [D A] [x] = [b]
# [Aᴴ 0] [y] [c]
opH = BlockDiagonalOperator(inv(D), eye(n))
(x, y, stats) = usymlqr(A, b, c, M=inv(D))
K = [D A; A' zeros(n,n)]
d = [b; c]
r = d - K * [x; y]
resid = sqrt(dot(r, opH * r))
@printf("USYMLQR: Relative residual: %8.1e\n", resid)
```
6 changes: 4 additions & 2 deletions docs/src/preconditioners.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -133,7 +135,7 @@ Methods concerned: [`CGNE`](@ref cgne), [`CRMR`](@ref crmr), [`LNLQ`](@ref lnlq)
| Arguments | `C` and `E` with `ldiv=false` | `C` and `E` with `ldiv=true` | `D` and `F` with `ldiv=false` | `D` and `F` with `ldiv=true` |

!!! note
Our implementations of [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr), [`USYMLQR`](@ref usymlqr) and [`TriLQR`](@ref trilqr) don't support preconditioning.
Our implementations of [`BiLQ`](@ref bilq), [`QMR`](@ref qmr), [`BiLQR`](@ref bilqr), [`USYMLQ`](@ref usymlq), [`USYMQR`](@ref usymqr) and [`TriLQR`](@ref trilqr) don't support preconditioning.

## Packages that provide preconditioners

Expand Down
34 changes: 20 additions & 14 deletions src/krylov_workspaces.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3239,17 +3239,19 @@ mutable struct UsymlqrWorkspace{T,FC,Sm,Sn} <: _KrylovWorkspace{T,FC,Sm,Sn}
x :: Sm
y :: Sn
z :: Sn
vₖ₋₁ :: Sm
vₖ :: Sm
uₖ₋₁ :: Sn
uₖ :: Sn
M⁻¹vₖ₋₁ :: Sm
M⁻¹vₖ :: Sm
N⁻¹uₖ₋₁ :: Sn
N⁻¹uₖ :: Sn
p :: Sn
q :: Sm
d̅ :: Sm
wₖ₋₂ :: Sn
wₖ₋₁ :: Sn
Δx :: Sm
Δy :: Sn
vₖ :: Sm
uₖ :: Sn
warm_start :: Bool
stats :: SimpleStats{T}
end
Expand All @@ -3264,19 +3266,21 @@ function UsymlqrWorkspace(kc::KrylovConstructor{Sm,Sn}) where {Sm,Sn}
x = similar(kc.vm)
y = similar(kc.vn)
z = similar(kc.vn)
vₖ₋₁ = similar(kc.vm)
vₖ = similar(kc.vm)
uₖ₋₁ = similar(kc.vn)
uₖ = similar(kc.vn)
M⁻¹vₖ₋₁ = similar(kc.vm)
M⁻¹vₖ = similar(kc.vm)
N⁻¹uₖ₋₁ = similar(kc.vn)
N⁻¹uₖ = similar(kc.vn)
p = similar(kc.vn)
q = similar(kc.vm)
d̅ = similar(kc.vm)
wₖ₋₂ = similar(kc.vn)
wₖ₋₁ = similar(kc.vn)
Δx = similar(kc.vm_empty)
Δy = similar(kc.vn_empty)
vₖ = similar(kc.vm_empty)
uₖ = similar(kc.vn_empty)
stats = SimpleStats(0, false, false, false, 0, T[], T[], T[], 0.0, 0.0, "unknown")
workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, vₖ₋₁, vₖ, uₖ₋₁, uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, false, stats)
workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, M⁻¹vₖ₋₁, M⁻¹vₖ, N⁻¹uₖ₋₁, N⁻¹uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, vₖ, uₖ, false, stats)
workspace.stats.allocation_timer = start_allocation_time |> ktimer
return workspace
end
Expand All @@ -3289,21 +3293,23 @@ function UsymlqrWorkspace(m::Integer, n::Integer, Sm::Type, Sn::Type)
x = Sm(undef, m)
y = Sn(undef, n)
z = Sn(undef, n)
vₖ₋₁ = Sm(undef, m)
vₖ = Sm(undef, m)
uₖ₋₁ = Sn(undef, n)
uₖ = Sn(undef, n)
M⁻¹vₖ₋₁ = Sm(undef, m)
M⁻¹vₖ = Sm(undef, m)
N⁻¹uₖ₋₁ = Sn(undef, n)
N⁻¹uₖ = Sn(undef, n)
p = Sn(undef, n)
q = Sm(undef, m)
d̅ = Sm(undef, m)
wₖ₋₂ = Sn(undef, n)
wₖ₋₁ = Sn(undef, n)
Δx = Sm(undef, 0)
Δy = Sn(undef, 0)
vₖ = Sm(undef, 0)
uₖ = Sn(undef, 0)
Sm = isconcretetype(Sm) ? Sm : typeof(x)
Sn = isconcretetype(Sn) ? Sn : typeof(y)
stats = SimpleStats(0, false, false, false, 0, T[], T[], T[], 0.0, 0.0, "unknown")
workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, vₖ₋₁, vₖ, uₖ₋₁, uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, false, stats)
workspace = UsymlqrWorkspace{T,FC,Sm,Sn}(m, n, r, x, y, z, M⁻¹vₖ₋₁, M⁻¹vₖ, N⁻¹uₖ₋₁, N⁻¹uₖ, p, q, d̅, wₖ₋₂, wₖ₋₁, Δx, Δy, vₖ, uₖ, false, stats)
workspace.stats.allocation_timer = start_allocation_time |> ktimer
return workspace
end
Expand Down
Loading
Loading