Skip to content
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.

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

BiLQR & TriLQR #159

Merged
merged 13 commits into from
Dec 18, 2019
1 change: 1 addition & 0 deletions docs/src/api.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ Krylov.KrylovStats
Krylov.SimpleStats
Krylov.LanczosStats
Krylov.SymmlqStats
Krylov.AdjointStats
```

## Utilities
Expand Down
1 change: 1 addition & 0 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ Given a linear operator ``A`` and a right-hand side ``b``, solve ``Ax = b``, whi
Krylov can be installed and tested through the Julia package manager:

```julia
julia> import Pkg
julia> Pkg.add("Krylov")
julia> Pkg.test("Krylov")
```
Expand Down
2 changes: 2 additions & 0 deletions docs/src/solvers.md
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,11 @@ diom
dqgmres
usymlq
usymqr
trilqr
bilq
cgs
qmr
bilqr
cgls
crls
cgne
Expand Down
21 changes: 21 additions & 0 deletions src/Krylov.jl
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,15 @@ mutable struct SymmlqStats{T} <: KrylovStats{T}
status :: String
end

"Type for statistics returned by adjoint systems solvers"
mutable struct AdjointStats{T} <: KrylovStats{T}
solved_primal :: Bool
solved_dual :: Bool
residuals_primal :: Vector{T}
residuals_dual :: Vector{T}
status :: String
end

import Base.show

function show(io :: IO, stats :: SimpleStats)
Expand Down Expand Up @@ -72,6 +81,16 @@ function show(io :: IO, stats :: SymmlqStats)
print(io, s)
end

function show(io :: IO, stats :: AdjointStats)
s = "\nAdjoint stats\n"
s *= @sprintf(" solved primal: %s\n", stats.solved_primal)
s *= @sprintf(" solved dual: %s\n", stats.solved_dual)
s *= @sprintf(" residuals primal: %s\n", vec2str(stats.residuals_primal))
s *= @sprintf(" residuals dual: %s\n", vec2str(stats.residuals_dual))
s *= @sprintf(" status: %s\n", stats.status)
print(io, s)
end

include("krylov_aux.jl")
include("krylov_utils.jl")

Expand All @@ -88,10 +107,12 @@ include("diom.jl")

include("usymlq.jl")
include("usymqr.jl")
include("trilqr.jl")

include("bilq.jl")
include("cgs.jl")
include("qmr.jl")
include("bilqr.jl")

include("cgls.jl")
include("crls.jl")
Expand Down
6 changes: 3 additions & 3 deletions src/bilq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ function bilq(A :: AbstractLinearOperator{T}, b :: AbstractVector{T}; c :: Abstr
d̅ = zeros(T, n) # Last column of D̅ₖ = Vₖ(Qₖ)ᵀ
ζₖ₋₁ = ζbarₖ = zero(T) # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
ζₖ₋₂ = ηₖ = zero(T) # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
δbarₖ₋₁ = δbarₖ = zero(T) # Coefficients of Lₖ₋₁ and L̅ₖ modified during two iterations
δbarₖ₋₁ = δbarₖ = zero(T) # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations
norm_vₖ = bNorm / βₖ # ‖vₖ‖ is used for residual norm estimates

# Stopping criterion.
Expand Down Expand Up @@ -178,7 +178,7 @@ function bilq(A :: AbstractLinearOperator{T}, b :: AbstractVector{T}; c :: Abstr
end

# Compute ⟨vₖ,vₖ₊₁⟩ and ‖vₖ₊₁‖
dot_vₖ_vₖ₊₁ = @kdot(n, vₖ₋₁, vₖ)
vₖᵀvₖ₊₁ = @kdot(n, vₖ₋₁, vₖ)
norm_vₖ₊₁ = @knrm2(n, vₖ)

# Compute BiLQ residual norm
Expand All @@ -188,7 +188,7 @@ function bilq(A :: AbstractLinearOperator{T}, b :: AbstractVector{T}; c :: Abstr
else
μₖ = βₖ * (sₖ₋₁ * ζₖ₋₂ - cₖ₋₁ * cₖ * ζₖ₋₁) + αₖ * sₖ * ζₖ₋₁
ωₖ = βₖ₊₁ * sₖ * ζₖ₋₁
rNorm_lq = sqrt(μₖ^2 * norm_vₖ^2 + ωₖ^2 * norm_vₖ₊₁^2 + 2 * μₖ * ωₖ * dot_vₖ_vₖ₊₁)
rNorm_lq = sqrt(μₖ^2 * norm_vₖ^2 + ωₖ^2 * norm_vₖ₊₁^2 + 2 * μₖ * ωₖ * vₖᵀvₖ₊₁)
end
push!(rNorms, rNorm_lq)

Expand Down
Loading