Skip to content

Commit

Permalink
[code] Use Aᴴ instead of Aᵀ
Browse files Browse the repository at this point in the history
  • Loading branch information
amontoison committed Sep 7, 2022
1 parent 604c960 commit 73e95a7
Show file tree
Hide file tree
Showing 30 changed files with 339 additions and 339 deletions.
6 changes: 3 additions & 3 deletions src/bicgstab.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,10 @@ export bicgstab, bicgstab!
Solve the square linear system Ax = b using the BICGSTAB method.
BICGSTAB requires two initial vectors `b` and `c`.
The relation `bᵀc ≠ 0` must be satisfied and by default `c = b`.
The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
The Biconjugate Gradient Stabilized method is a variant of BiCG, like CGS,
but using different updates for the Aᵀ-sequence in order to obtain smoother
but using different updates for the Aᴴ-sequence in order to obtain smoother
convergence than CGS.
If BICGSTAB stagnates, we recommend DQGMRES and BiLQ as alternative methods for unsymmetric square systems.
Expand Down Expand Up @@ -157,7 +157,7 @@ function bicgstab!(solver :: BicgstabSolver{T,FC,S}, A, b :: AbstractVector{FC};
if next_ρ == 0
stats.niter = 0
stats.solved, stats.inconsistent = false, false
stats.status = "Breakdown bᵀc = 0"
stats.status = "Breakdown bᴴc = 0"
solver.warm_start = false
return solver
end
Expand Down
36 changes: 18 additions & 18 deletions src/bilq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ export bilq, bilq!
Solve the square linear system Ax = b using the BiLQ method.
BiLQ is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`.
The relation `bᵀc ≠ 0` must be satisfied and by default `c = b`.
The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
When `A` is symmetric and `b = c`, BiLQ is equivalent to SYMMLQ.
An option gives the possibility of transferring to the BiCG point,
Expand Down Expand Up @@ -90,7 +90,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
ktypeof(c) == S || error("ktypeof(c) ≠ $S")

# Compute the adjoint of A
Aᵀ = A'
Aᴴ = A'

# Set up workspace.
uₖ₋₁, uₖ, q, vₖ₋₁, vₖ = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ
Expand Down Expand Up @@ -127,25 +127,25 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
kdisplay(iter, verbose) && @printf("%5d %7.1e\n", iter, bNorm)

# Initialize the Lanczos biorthogonalization process.
cᵗb = @kdot(n, c, r₀) # ⟨c,r₀⟩
if cᵗb == 0
cᴴb = @kdot(n, c, r₀) # ⟨c,r₀⟩
if cᴴb == 0
stats.niter = 0
stats.solved = false
stats.inconsistent = false
stats.status = "Breakdown bᵀc = 0"
stats.status = "Breakdown bᴴc = 0"
solver.warm_start = false
return solver
end

βₖ = (abs(cᵗb)) # β₁γ₁ = cᵀ(b - Ax₀)
γₖ = cᵗb / βₖ # β₁γ₁ = cᵀ(b - Ax₀)
βₖ = (abs(cᴴb)) # β₁γ₁ = cᴴ(b - Ax₀)
γₖ = cᴴb / βₖ # β₁γ₁ = cᴴ(b - Ax₀)
vₖ₋₁ .= zero(FC) # v₀ = 0
uₖ₋₁ .= zero(FC) # u₀ = 0
vₖ .= r₀ ./ βₖ # v₁ = (b - Ax₀) / β₁
uₖ .= c ./ conj(γₖ) # u₁ = c / γ̄₁
cₖ₋₁ = cₖ = -one(T) # Givens cosines used for the LQ factorization of Tₖ
sₖ₋₁ = sₖ = zero(FC) # Givens sines used for the LQ factorization of Tₖ
d̅ .= zero(FC) # Last column of D̅ₖ = Vₖ(Qₖ)
d̅ .= zero(FC) # Last column of D̅ₖ = Vₖ(Qₖ)
ζₖ₋₁ = ζbarₖ = zero(FC) # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
ζₖ₋₂ = ηₖ = zero(FC) # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
δbarₖ₋₁ = δbarₖ = zero(FC) # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations
Expand All @@ -165,10 +165,10 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab

# Continue the Lanczos biorthogonalization process.
# AVₖ = VₖTₖ + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
# AᵀUₖ = Uₖ(Tₖ) + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)
# AᴴUₖ = Uₖ(Tₖ) + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)

mul!(q, A , vₖ) # Forms vₖ₊₁ : q ← Avₖ
mul!(p, Aᵀ, uₖ) # Forms uₖ₊₁ : p ← Aᵀuₖ
mul!(p, Aᴴ, uₖ) # Forms uₖ₊₁ : p ← Aᴴuₖ

@kaxpy!(n, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁
@kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - β̄ₖ * uₖ₋₁
Expand All @@ -178,9 +178,9 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
@kaxpy!(n, - αₖ , vₖ, q) # q ← q - αₖ * vₖ
@kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ

pᵗq = @kdot(n, p, q) # pᵗq = ⟨p,q⟩
βₖ₊₁ = (abs(pᵗq)) # βₖ₊₁ = √(|pᵗq|)
γₖ₊₁ = pᵗq / βₖ₊₁ # γₖ₊₁ = pᵗq / βₖ₊₁
pᴴq = @kdot(n, p, q) # pᴴq = ⟨p,q⟩
βₖ₊₁ = (abs(pᴴq)) # βₖ₊₁ = √(|pᴴq|)
γₖ₊₁ = pᴴq / βₖ₊₁ # γₖ₊₁ = pᴴq / βₖ₊₁

# Update the LQ factorization of Tₖ = L̅ₖQₖ.
# [ α₁ γ₂ 0 • • • 0 ] [ δ₁ 0 • • • • 0 ]
Expand Down Expand Up @@ -235,7 +235,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
ηₖ = -ϵₖ₋₂ * ζₖ₋₂ - λₖ₋₁ * ζₖ₋₁
end

# Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ).
# Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ).
# [d̅ₖ₋₁ vₖ] [cₖ s̄ₖ] = [dₖ₋₁ d̅ₖ] ⟷ dₖ₋₁ = cₖ * d̅ₖ₋₁ + sₖ * vₖ
# [sₖ -cₖ] ⟷ d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ
if iter 2
Expand All @@ -258,13 +258,13 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
@. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ
@. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ

if pᵗq 0
if pᴴq 0
@. vₖ = q / βₖ₊₁ # βₖ₊₁vₖ₊₁ = q
@. uₖ = p / conj(γₖ₊₁) # γ̄ₖ₊₁uₖ₊₁ = p
end

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

# Compute BiLQ residual norm
Expand All @@ -274,7 +274,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
else
μₖ = βₖ * (sₖ₋₁ * ζₖ₋₂ - cₖ₋₁ * cₖ * ζₖ₋₁) + αₖ * sₖ * ζₖ₋₁
ωₖ = βₖ₊₁ * sₖ * ζₖ₋₁
θₖ = conj(μₖ) * ωₖ * vₖᵀvₖ₊₁
θₖ = conj(μₖ) * ωₖ * vₖᴴvₖ₊₁
rNorm_lq = sqrt(abs2(μₖ) * norm_vₖ^2 + abs2(ωₖ) * norm_vₖ₊₁^2 + 2 * real(θₖ))
end
history && push!(rNorms, rNorm_lq)
Expand All @@ -300,7 +300,7 @@ function bilq!(solver :: BilqSolver{T,FC,S}, A, b :: AbstractVector{FC}; c :: Ab
solved_lq = rNorm_lq ε
solved_cg = transfer_to_bicg && (abs(δbarₖ) > eps(T)) && (rNorm_cg ε)
tired = iter itmax
breakdown = !solved_lq && !solved_cg && (pᵗq == 0)
breakdown = !solved_lq && !solved_cg && (pᴴq == 0)
kdisplay(iter, verbose) && @printf("%5d %7.1e\n", iter, rNorm_lq)
end
(verbose > 0) && @printf("\n")
Expand Down
54 changes: 27 additions & 27 deletions src/bilqr.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# An implementation of BILQR for the solution of square
# consistent linear adjoint systems Ax = b and Aᵀy = c.
# consistent linear adjoint systems Ax = b and Aᴴy = c.
#
# This method is described in
#
Expand All @@ -24,11 +24,11 @@ export bilqr, bilqr!
Combine BiLQ and QMR to solve adjoint systems.
[0 A] [y] = [b]
[Aᵀ 0] [x] [c]
[Aᴴ 0] [x] [c]
The relation `bᵀc ≠ 0` must be satisfied.
The relation `bᴴc ≠ 0` must be satisfied.
BiLQ is used for solving primal system `Ax = b`.
QMR is used for solving dual system `Aᵀy = c`.
QMR is used for solving dual system `Aᴴy = c`.
An option gives the possibility of transferring from the BiLQ point to the
BiCG point, when it exists. The transfer is based on the residual norm.
Expand Down Expand Up @@ -94,7 +94,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
ktypeof(c) == S || error("ktypeof(c) ≠ $S")

# Compute the adjoint of A
Aᵀ = A'
Aᴴ = A'

# Set up workspace.
uₖ₋₁, uₖ, q, vₖ₋₁, vₖ = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ
Expand All @@ -109,15 +109,15 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
if warm_start
mul!(r₀, A, Δx)
@kaxpby!(n, one(FC), b, -one(FC), r₀)
mul!(s₀, Aᵀ, Δy)
mul!(s₀, Aᴴ, Δy)
@kaxpby!(n, one(FC), c, -one(FC), s₀)
end

# Initial solution x₀ and residual norm ‖r₀‖ = ‖b - Ax₀‖.
x .= zero(FC) # x₀
bNorm = @knrm2(n, r₀) # rNorm = ‖r₀‖

# Initial solution t₀ and residual norm ‖s₀‖ = ‖c - Aᵀy₀‖.
# Initial solution t₀ and residual norm ‖s₀‖ = ‖c - Aᴴy₀‖.
t .= zero(FC) # t₀
cNorm = @knrm2(n, s₀) # sNorm = ‖s₀‖

Expand All @@ -132,34 +132,34 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
kdisplay(iter, verbose) && @printf("%5d %7.1e %7.1e\n", iter, bNorm, cNorm)

# Initialize the Lanczos biorthogonalization process.
cᵗb = @kdot(n, s₀, r₀) # ⟨s₀,r₀⟩ = ⟨c - Aᵀy₀,b - Ax₀⟩
if cᵗb == 0
cᴴb = @kdot(n, s₀, r₀) # ⟨s₀,r₀⟩ = ⟨c - Aᴴy₀,b - Ax₀⟩
if cᴴb == 0
stats.niter = 0
stats.solved_primal = false
stats.solved_dual = false
stats.status = "Breakdown bᵀc = 0"
stats.status = "Breakdown bᴴc = 0"
solver.warm_start = false
return solver
end

# Set up workspace.
βₖ = (abs(cᵗb)) # β₁γ₁ = (c - Aᵀy₀)ᵀ(b - Ax₀)
γₖ = cᵗb / βₖ # β₁γ₁ = (c - Aᵀy₀)ᵀ(b - Ax₀)
βₖ = (abs(cᴴb)) # β₁γ₁ = (c - Aᴴy₀)ᴴ(b - Ax₀)
γₖ = cᴴb / βₖ # β₁γ₁ = (c - Aᴴy₀)ᴴ(b - Ax₀)
vₖ₋₁ .= zero(FC) # v₀ = 0
uₖ₋₁ .= zero(FC) # u₀ = 0
vₖ .= r₀ ./ βₖ # v₁ = (b - Ax₀) / β₁
uₖ .= s₀ ./ conj(γₖ) # u₁ = (c - Aᵀy₀) / γ̄₁
uₖ .= s₀ ./ conj(γₖ) # u₁ = (c - Aᴴy₀) / γ̄₁
cₖ₋₁ = cₖ = -one(T) # Givens cosines used for the LQ factorization of Tₖ
sₖ₋₁ = sₖ = zero(FC) # Givens sines used for the LQ factorization of Tₖ
d̅ .= zero(FC) # Last column of D̅ₖ = Vₖ(Qₖ)
d̅ .= zero(FC) # Last column of D̅ₖ = Vₖ(Qₖ)
ζₖ₋₁ = ζbarₖ = zero(FC) # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ = (L̅ₖ)⁻¹β₁e₁
ζₖ₋₂ = ηₖ = zero(FC) # ζₖ₋₂ and ηₖ are used to update ζₖ₋₁ and ζbarₖ
δbarₖ₋₁ = δbarₖ = zero(FC) # Coefficients of Lₖ₋₁ and L̅ₖ modified over the course of two iterations
ψbarₖ₋₁ = ψₖ₋₁ = zero(FC) # ψₖ₋₁ and ψbarₖ are the last components of h̅ₖ = Qₖγ̄₁e₁
norm_vₖ = bNorm / βₖ # ‖vₖ‖ is used for residual norm estimates
ϵₖ₋₃ = λₖ₋₂ = zero(FC) # Components of Lₖ₋₁
wₖ₋₃ .= zero(FC) # Column k-3 of Wₖ = Uₖ(Lₖ)⁻
wₖ₋₂ .= zero(FC) # Column k-2 of Wₖ = Uₖ(Lₖ)⁻
wₖ₋₃ .= zero(FC) # Column k-3 of Wₖ = Uₖ(Lₖ)⁻
wₖ₋₂ .= zero(FC) # Column k-2 of Wₖ = Uₖ(Lₖ)⁻
τₖ = zero(T) # τₖ is used for the dual residual norm estimate

# Stopping criterion.
Expand All @@ -180,10 +180,10 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::

# Continue the Lanczos biorthogonalization process.
# AVₖ = VₖTₖ + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
# AᵀUₖ = Uₖ(Tₖ) + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)
# AᴴUₖ = Uₖ(Tₖ) + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)

mul!(q, A , vₖ) # Forms vₖ₊₁ : q ← Avₖ
mul!(p, Aᵀ, uₖ) # Forms uₖ₊₁ : p ← Aᵀuₖ
mul!(p, Aᴴ, uₖ) # Forms uₖ₊₁ : p ← Aᴴuₖ

@kaxpy!(n, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁
@kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - β̄ₖ * uₖ₋₁
Expand All @@ -193,9 +193,9 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
@kaxpy!(n, - αₖ , vₖ, q) # q ← q - αₖ * vₖ
@kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ

pᵗq = @kdot(n, p, q) # pᵗq = ⟨p,q⟩
βₖ₊₁ = (abs(pᵗq)) # βₖ₊₁ = √(|pᵗq|)
γₖ₊₁ = pᵗq / βₖ₊₁ # γₖ₊₁ = pᵗq / βₖ₊₁
pᴴq = @kdot(n, p, q) # pᴴq = ⟨p,q⟩
βₖ₊₁ = (abs(pᴴq)) # βₖ₊₁ = √(|pᴴq|)
γₖ₊₁ = pᴴq / βₖ₊₁ # γₖ₊₁ = pᴴq / βₖ₊₁

# Update the LQ factorization of Tₖ = L̅ₖQₖ.
# [ α₁ γ₂ 0 • • • 0 ] [ δ₁ 0 • • • • 0 ]
Expand Down Expand Up @@ -251,7 +251,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
ηₖ = -ϵₖ₋₂ * ζₖ₋₂ - λₖ₋₁ * ζₖ₋₁
end

# Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ).
# Relations for the directions dₖ₋₁ and d̅ₖ, the last two columns of D̅ₖ = Vₖ(Qₖ).
# [d̅ₖ₋₁ vₖ] [cₖ s̄ₖ] = [dₖ₋₁ d̅ₖ] ⟷ dₖ₋₁ = cₖ * d̅ₖ₋₁ + sₖ * vₖ
# [sₖ -cₖ] ⟷ d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ
if iter 2
Expand All @@ -271,7 +271,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
end

# Compute ⟨vₖ,vₖ₊₁⟩ and ‖vₖ₊₁‖
vₖᵀvₖ₊₁ = @kdot(n, vₖ, q) / βₖ₊₁
vₖᴴvₖ₊₁ = @kdot(n, vₖ, q) / βₖ₊₁
norm_vₖ₊₁ = @knrm2(n, q) / βₖ₊₁

# Compute BiLQ residual norm
Expand All @@ -281,7 +281,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
else
μₖ = βₖ * (sₖ₋₁ * ζₖ₋₂ - cₖ₋₁ * cₖ * ζₖ₋₁) + αₖ * sₖ * ζₖ₋₁
ωₖ = βₖ₊₁ * sₖ * ζₖ₋₁
θₖ = conj(μₖ) * ωₖ * vₖᵀvₖ₊₁
θₖ = conj(μₖ) * ωₖ * vₖᴴvₖ₊₁
rNorm_lq = sqrt(abs2(μₖ) * norm_vₖ^2 + abs2(ωₖ) * norm_vₖ₊₁^2 + 2 * real(θₖ))
end
history && push!(rNorms, rNorm_lq)
Expand Down Expand Up @@ -318,7 +318,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
ψbarₖ = sₖ * ψbarₖ₋₁
end

# Compute the direction wₖ₋₁, the last column of Wₖ₋₁ = (Uₖ₋₁)(Lₖ₋₁)⁻ ⟷ (L̄ₖ₋₁)(Wₖ₋₁)ᵀ = (Uₖ₋₁)ᵀ.
# Compute the direction wₖ₋₁, the last column of Wₖ₋₁ = (Uₖ₋₁)(Lₖ₋₁)⁻ ⟷ (L̄ₖ₋₁)(Wₖ₋₁)ᵀ = (Uₖ₋₁)ᵀ.
# w₁ = u₁ / δ̄₁
if iter == 2
wₖ₋₁ = wₖ₋₂
Expand Down Expand Up @@ -372,7 +372,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::
@. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ
@. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ

if pᵗq zero(FC)
if pᴴq zero(FC)
@. vₖ = q / βₖ₊₁ # βₖ₊₁vₖ₊₁ = q
@. uₖ = p / conj(γₖ₊₁) # γ̄ₖ₊₁uₖ₊₁ = p
end
Expand All @@ -392,7 +392,7 @@ function bilqr!(solver :: BilqrSolver{T,FC,S}, A, b :: AbstractVector{FC}, c ::

user_requested_exit = callback(solver) :: Bool
tired = iter itmax
breakdown = !solved_lq && !solved_cg && (pᵗq == 0)
breakdown = !solved_lq && !solved_cg && (pᴴq == 0)

kdisplay(iter, verbose) && solved_primal && !solved_dual && @printf("%5d %7s %7.1e\n", iter, "", sNorm)
kdisplay(iter, verbose) && !solved_primal && solved_dual && @printf("%5d %7.1e %7s\n", iter, rNorm_lq, "")
Expand Down
8 changes: 4 additions & 4 deletions src/cg_lanczos.jl
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ function cg_lanczos!(solver :: CgLanczosSolver{T,FC,S}, A, b :: AbstractVector{F
Mv .= b
end
MisI || mulorldiv!(v, M, Mv, ldiv) # v₁ = M⁻¹r₀
β = sqrt(@kdotr(n, v, Mv)) # β₁ = v₁ M v₁
β = sqrt(@kdotr(n, v, Mv)) # β₁ = v₁ M v₁
σ = β
rNorm = σ
history && push!(rNorms, rNorm)
Expand Down Expand Up @@ -157,10 +157,10 @@ function cg_lanczos!(solver :: CgLanczosSolver{T,FC,S}, A, b :: AbstractVector{F
# Form next Lanczos vector.
# βₖ₊₁Mvₖ₊₁ = Avₖ - δₖMvₖ - βₖMvₖ₋₁
mul!(Mv_next, A, v) # Mvₖ₊₁ ← Avₖ
δ = @kdotr(n, v, Mv_next) # δₖ = vₖᵀ A vₖ
δ = @kdotr(n, v, Mv_next) # δₖ = vₖᴴ A vₖ

# Check curvature. Exit fast if requested.
# It is possible to show that σₖ² (δₖ - ωₖ₋₁ / γₖ₋₁) = pₖᵀ A pₖ.
# It is possible to show that σₖ² (δₖ - ωₖ₋₁ / γₖ₋₁) = pₖᴴ A pₖ.
γ = one(T) /- ω / γ) # γₖ = 1 / (δₖ - ωₖ₋₁ / γₖ₋₁)
indefinite |= 0)
(check_curvature & indefinite) && continue
Expand All @@ -172,7 +172,7 @@ function cg_lanczos!(solver :: CgLanczosSolver{T,FC,S}, A, b :: AbstractVector{F
end
@. Mv = Mv_next # Mvₖ ← Mvₖ₊₁
MisI || mulorldiv!(v, M, Mv, ldiv) # vₖ₊₁ = M⁻¹ * Mvₖ₊₁
β = sqrt(@kdotr(n, v, Mv)) # βₖ₊₁ = vₖ₊₁ M vₖ₊₁
β = sqrt(@kdotr(n, v, Mv)) # βₖ₊₁ = vₖ₊₁ M vₖ₊₁
@kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁
MisI || @kscal!(n, one(FC) / β, Mv) # Mvₖ₊₁ ← Mvₖ₊₁ / βₖ₊₁
Anorm2 += β_prev^2 + β^2 + δ^2 # Use ‖Tₖ₊₁‖₂ as increasing approximation of ‖A‖₂.
Expand Down
Loading

0 comments on commit 73e95a7

Please sign in to comment.