From 2558bfe937575b12a83642c272566e6c7fe8db73 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 17 Oct 2024 16:42:14 -0500 Subject: [PATCH] Restore the alignment of comments --- src/bicgstab.jl | 48 +++++++++++++++++++-------------------- src/bilq.jl | 10 ++++---- src/bilqr.jl | 14 ++++++------ src/car.jl | 6 ++--- src/cg_lanczos.jl | 18 +++++++-------- src/cg_lanczos_shift.jl | 8 +++---- src/cgls.jl | 14 ++++++------ src/cgls_lanczos_shift.jl | 24 ++++++++++---------- src/cgne.jl | 6 ++--- src/cgs.jl | 24 ++++++++++---------- src/cr.jl | 24 ++++++++++---------- src/craig.jl | 2 +- src/craigmr.jl | 8 +++---- src/crls.jl | 20 ++++++++-------- src/crmr.jl | 22 +++++++++--------- src/diom.jl | 4 ++-- src/dqgmres.jl | 4 ++-- src/fgmres.jl | 2 +- src/fom.jl | 2 +- src/gmres.jl | 6 ++--- src/lnlq.jl | 16 ++++++------- src/lsmr.jl | 8 +++---- src/minres.jl | 4 ++-- src/minres_qlp.jl | 4 ++-- src/qmr.jl | 14 ++++++------ src/symmlq.jl | 2 +- src/trilqr.jl | 22 +++++++++--------- src/usymlq.jl | 14 ++++++------ src/usymqr.jl | 30 ++++++++++++------------ 29 files changed, 190 insertions(+), 190 deletions(-) diff --git a/src/bicgstab.jl b/src/bicgstab.jl index ef05b64b1..4ebba1171 100644 --- a/src/bicgstab.jl +++ b/src/bicgstab.jl @@ -154,15 +154,15 @@ kwargs_bicgstab = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, kcopy!(n, r₀, b) # r₀ ← b end - kfill!(x, zero(FC)) # x₀ - kfill!(s, zero(FC)) # s₀ - kfill!(v, zero(FC)) # v₀ + kfill!(x, zero(FC)) # x₀ + kfill!(s, zero(FC)) # s₀ + kfill!(v, zero(FC)) # v₀ MisI || mulorldiv!(r, M, r₀, ldiv) # r₀ - kcopy!(n, p, r) # p₁ + kcopy!(n, p, r) # p₁ - α = one(FC) # α₀ - ω = one(FC) # ω₀ - ρ = one(FC) # ρ₀ + α = one(FC) # α₀ + ω = one(FC) # ω₀ + ρ = one(FC) # ρ₀ # Compute residual norm ‖r₀‖₂. rNorm = knorm(n, r) @@ -206,24 +206,24 @@ kwargs_bicgstab = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, iter = iter + 1 ρ = next_ρ - NisI || mulorldiv!(y, N, p, ldiv) # yₖ = N⁻¹pₖ - mul!(q, A, y) # qₖ = Ayₖ - mulorldiv!(v, M, q, ldiv) # vₖ = M⁻¹qₖ - α = ρ / kdot(n, c, v) # αₖ = ⟨r̅₀,rₖ₋₁⟩ / ⟨r̅₀,vₖ⟩ - kcopy!(n, s, r) # sₖ = rₖ₋₁ - kaxpy!(n, -α, v, s) # sₖ = sₖ - αₖvₖ - kaxpy!(n, α, y, x) # xₐᵤₓ = xₖ₋₁ + αₖyₖ - NisI || mulorldiv!(z, N, s, ldiv) # zₖ = N⁻¹sₖ - mul!(d, A, z) # dₖ = Azₖ - MisI || mulorldiv!(t, M, d, ldiv) # tₖ = M⁻¹dₖ + NisI || mulorldiv!(y, N, p, ldiv) # yₖ = N⁻¹pₖ + mul!(q, A, y) # qₖ = Ayₖ + mulorldiv!(v, M, q, ldiv) # vₖ = M⁻¹qₖ + α = ρ / kdot(n, c, v) # αₖ = ⟨r̅₀,rₖ₋₁⟩ / ⟨r̅₀,vₖ⟩ + kcopy!(n, s, r) # sₖ = rₖ₋₁ + kaxpy!(n, -α, v, s) # sₖ = sₖ - αₖvₖ + kaxpy!(n, α, y, x) # xₐᵤₓ = xₖ₋₁ + αₖyₖ + NisI || mulorldiv!(z, N, s, ldiv) # zₖ = N⁻¹sₖ + mul!(d, A, z) # dₖ = Azₖ + MisI || mulorldiv!(t, M, d, ldiv) # tₖ = M⁻¹dₖ ω = kdot(n, t, s) / kdot(n, t, t) # ⟨tₖ,sₖ⟩ / ⟨tₖ,tₖ⟩ - kaxpy!(n, ω, z, x) # xₖ = xₐᵤₓ + ωₖzₖ - kcopy!(n, r, s) # rₖ = sₖ - kaxpy!(n, -ω, t, r) # rₖ = rₖ - ωₖtₖ - next_ρ = kdot(n, c, r) # ρₖ₊₁ = ⟨r̅₀,rₖ⟩ - β = (next_ρ / ρ) * (α / ω) # βₖ₊₁ = (ρₖ₊₁ / ρₖ) * (αₖ / ωₖ) - kaxpy!(n, -ω, v, p) # pₐᵤₓ = pₖ - ωₖvₖ - kaxpby!(n, one(FC), r, β, p) # pₖ₊₁ = rₖ₊₁ + βₖ₊₁pₐᵤₓ + kaxpy!(n, ω, z, x) # xₖ = xₐᵤₓ + ωₖzₖ + kcopy!(n, r, s) # rₖ = sₖ + kaxpy!(n, -ω, t, r) # rₖ = rₖ - ωₖtₖ + next_ρ = kdot(n, c, r) # ρₖ₊₁ = ⟨r̅₀,rₖ⟩ + β = (next_ρ / ρ) * (α / ω) # βₖ₊₁ = (ρₖ₊₁ / ρₖ) * (αₖ / ωₖ) + kaxpy!(n, -ω, v, p) # pₐᵤₓ = pₖ - ωₖvₖ + kaxpby!(n, one(FC), r, β, p) # pₖ₊₁ = rₖ₊₁ + βₖ₊₁pₐᵤₓ # Compute residual norm ‖rₖ‖₂. rNorm = knorm(n, r) diff --git a/src/bilq.jl b/src/bilq.jl index 4d195fbcd..68fde4489 100644 --- a/src/bilq.jl +++ b/src/bilq.jl @@ -191,13 +191,13 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time βₖ = √(abs(cᴴb)) # β₁γ₁ = cᴴ(b - Ax₀) γₖ = cᴴb / βₖ # β₁γ₁ = cᴴ(b - Ax₀) - kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 - kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 + kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 + kfill!(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ₖ - kfill!(d̅, zero(FC)) # Last column of D̅ₖ = Vₖ(Qₖ)ᴴ + kfill!(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 @@ -239,8 +239,8 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time 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 / βₖ₊₁ + βₖ₊₁ = √(abs(pᴴq)) # βₖ₊₁ = √(|pᴴq|) + γₖ₊₁ = pᴴq / βₖ₊₁ # γₖ₊₁ = pᴴq / βₖ₊₁ # Update the LQ factorization of Tₖ = L̅ₖQₖ. # [ α₁ γ₂ 0 • • • 0 ] [ δ₁ 0 • • • • 0 ] diff --git a/src/bilqr.jl b/src/bilqr.jl index 2aaa895cf..35df0b98e 100644 --- a/src/bilqr.jl +++ b/src/bilqr.jl @@ -175,21 +175,21 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi # Set up workspace. βₖ = √(abs(cᴴb)) # β₁γ₁ = (c - Aᴴy₀)ᴴ(b - Ax₀) γₖ = cᴴb / βₖ # β₁γ₁ = (c - Aᴴy₀)ᴴ(b - Ax₀) - kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 - kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 + kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 + kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 vₖ .= r₀ ./ βₖ # v₁ = (b - Ax₀) / β₁ 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ₖ - kfill!(d̅, zero(FC)) # Last column of D̅ₖ = Vₖ(Qₖ)ᴴ + kfill!(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ₖ₋₁ - kfill!(wₖ₋₃, zero(FC)) # Column k-3 of Wₖ = Uₖ(Lₖ)⁻ᴴ - kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Uₖ(Lₖ)⁻ᴴ + kfill!(wₖ₋₃, zero(FC)) # Column k-3 of Wₖ = Uₖ(Lₖ)⁻ᴴ + kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Uₖ(Lₖ)⁻ᴴ τₖ = zero(T) # τₖ is used for the dual residual norm estimate # Stopping criterion. @@ -225,8 +225,8 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi 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 / βₖ₊₁ + βₖ₊₁ = √(abs(pᴴq)) # βₖ₊₁ = √(|pᴴq|) + γₖ₊₁ = pᴴq / βₖ₊₁ # γₖ₊₁ = pᴴq / βₖ₊₁ # Update the LQ factorization of Tₖ = L̅ₖQₖ. # [ α₁ γ₂ 0 • • • 0 ] [ δ₁ 0 • • • • 0 ] diff --git a/src/car.jl b/src/car.jl index 94755e0ac..37d02b03d 100644 --- a/src/car.jl +++ b/src/car.jl @@ -148,7 +148,7 @@ kwargs_car = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :ca kcopy!(n, s, q) # s ← q end - mul!(t, A, s) # t₀ = As₀ + mul!(t, A, s) # t₀ = As₀ kcopy!(n, u, t) # u₀ = Aq₀ ρ = kdotr(n, t, s) # ρ₀ = ⟨t₀ , s₀⟩ @@ -201,9 +201,9 @@ kwargs_car = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :ca solved = resid_decrease_lim || resid_decrease_mach if !solved - mul!(t, A, s) # tₖ₊₁ = A * sₖ₊₁ + mul!(t, A, s) # tₖ₊₁ = A * sₖ₊₁ ρ_next = kdotr(n, t, s) # ρₖ₊₁ = ⟨tₖ₊₁ , sₖ₊₁⟩ - β = ρ_next / ρ # βₖ = ρₖ₊₁ / ρₖ + β = ρ_next / ρ # βₖ = ρₖ₊₁ / ρₖ ρ = ρ_next kaxpby!(n, one(FC), r, β, p) # pₖ₊₁ = rₖ₊₁ + βₖ * pₖ kaxpby!(n, one(FC), s, β, q) # qₖ₊₁ = sₖ₊₁ + βₖ * qₖ diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl index 860c88765..636726c68 100644 --- a/src/cg_lanczos.jl +++ b/src/cg_lanczos.jl @@ -138,7 +138,7 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax kcopy!(n, Mv, b) # 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) @@ -185,7 +185,7 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax while ! (solved || tired || (check_curvature & indefinite) || user_requested_exit || overtimed) # Form next Lanczos vector. # βₖ₊₁Mvₖ₊₁ = Avₖ - δₖMvₖ - βₖMvₖ₋₁ - mul!(Mv_next, A, v) # Mvₖ₊₁ ← Avₖ + mul!(Mv_next, A, v) # Mvₖ₊₁ ← Avₖ δ = kdotr(n, v, Mv_next) # δₖ = vₖᴴ A vₖ # Check curvature. Exit fast if requested. @@ -194,25 +194,25 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax indefinite |= (γ ≤ 0) (check_curvature & indefinite) && continue - kaxpy!(n, -δ, Mv, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ + kaxpy!(n, -δ, Mv, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ if iter > 0 - kaxpy!(n, -β, Mv_prev, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - βₖMvₖ₋₁ - kcopy!(n, Mv_prev, Mv) # Mvₖ₋₁ ← Mvₖ + kaxpy!(n, -β, Mv_prev, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - βₖMvₖ₋₁ + kcopy!(n, Mv_prev, Mv) # Mvₖ₋₁ ← Mvₖ end kcopy!(n, Mv, Mv_next) # Mvₖ ← Mvₖ₊₁ - MisI || mulorldiv!(v, M, Mv, ldiv) # vₖ₊₁ = M⁻¹ * Mvₖ₊₁ + MisI || mulorldiv!(v, M, Mv, ldiv) # vₖ₊₁ = M⁻¹ * Mvₖ₊₁ β = 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‖₂. + Anorm2 += β_prev^2 + β^2 + δ^2 # Use ‖Tₖ₊₁‖₂ as increasing approximation of ‖A‖₂. β_prev = β # Compute next CG iterate. - kaxpy!(n, γ, p, x) # xₖ₊₁ = xₖ + γₖ * pₖ + kaxpy!(n, γ, p, x) # xₖ₊₁ = xₖ + γₖ * pₖ ω = β * γ σ = -ω * σ # σₖ₊₁ = - βₖ₊₁ * γₖ * σₖ ω = ω * ω # ωₖ = (βₖ₊₁ * γₖ)² - kaxpby!(n, σ, v, ω, p) # pₖ₊₁ = σₖ₊₁ * vₖ₊₁ + ωₖ * pₖ + kaxpby!(n, σ, v, ω, p) # pₖ₊₁ = σₖ₊₁ * vₖ₊₁ + ωₖ * pₖ rNorm = abs(σ) # ‖rₖ₊₁‖_M = |σₖ₊₁| because rₖ₊₁ = σₖ₊₁ * vₖ₊₁ and ‖vₖ₊₁‖_M = 1 history && push!(rNorms, rNorm) iter = iter + 1 diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl index 2f40d9961..0e0056f95 100644 --- a/src/cg_lanczos_shift.jl +++ b/src/cg_lanczos_shift.jl @@ -131,9 +131,9 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t for i = 1 : nshifts kfill!(x[i], zero(FC)) # x₀ end - kcopy!(n, Mv, b) # Mv₁ ← b + kcopy!(n, Mv, b) # Mv₁ ← b 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₁ kfill!(rNorms, β) if history for i = 1 : nshifts @@ -196,7 +196,7 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t while ! (solved || tired || user_requested_exit || overtimed) # Form next Lanczos vector. # βₖ₊₁Mvₖ₊₁ = Avₖ - δₖMvₖ - βₖMvₖ₋₁ - mul!(Mv_next, A, v) # Mvₖ₊₁ ← Avₖ + mul!(Mv_next, A, v) # Mvₖ₊₁ ← Avₖ δ = kdotr(n, v, Mv_next) # δₖ = vₖᴴ A vₖ kaxpy!(n, -δ, Mv, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ if iter > 0 @@ -204,7 +204,7 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t kcopy!(n, Mv_prev, Mv) # Mvₖ₋₁ ← Mvₖ end kcopy!(n, Mv, Mv_next) # Mvₖ ← Mvₖ₊₁ - MisI || mulorldiv!(v, M, Mv, ldiv) # vₖ₊₁ = M⁻¹ * Mvₖ₊₁ + MisI || mulorldiv!(v, M, Mv, ldiv) # vₖ₊₁ = M⁻¹ * Mvₖ₊₁ β = sqrt(kdotr(n, v, Mv)) # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁ kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁ MisI || kscal!(n, one(FC) / β, Mv) # Mvₖ₊₁ ← Mvₖ₊₁ / βₖ₊₁ diff --git a/src/cgls.jl b/src/cgls.jl index 29aa67b4b..76c9d7ade 100644 --- a/src/cgls.jl +++ b/src/cgls.jl @@ -147,8 +147,8 @@ kwargs_cgls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose Mq = MisI ? q : solver.Mr kfill!(x, zero(FC)) - kcopy!(m, r, b) # r ← b - bNorm = knorm(m, r) # Marginally faster than norm(b) + kcopy!(m, r, b) # r ← b + bNorm = knorm(m, r) # Marginally faster than norm(b) if bNorm == 0 stats.niter = 0 stats.solved, stats.inconsistent = true, false @@ -160,7 +160,7 @@ kwargs_cgls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose end MisI || mulorldiv!(Mr, M, r, ldiv) mul!(s, Aᴴ, Mr) - kcopy!(n, p, s) # p ← s + kcopy!(n, p, s) # p ← s γ = kdotr(n, s, s) # γ = sᴴs iter = 0 itmax == 0 && (itmax = m + n) @@ -194,12 +194,12 @@ kwargs_cgls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose on_boundary = true end - kaxpy!(n, α, p, x) # Faster than x = x + α * p - kaxpy!(m, -α, q, r) # Faster than r = r - α * q + kaxpy!(n, α, p, x) # Faster than x = x + α * p + kaxpy!(m, -α, q, r) # Faster than r = r - α * q MisI || mulorldiv!(Mr, M, r, ldiv) mul!(s, Aᴴ, Mr) - λ > 0 && kaxpy!(n, -λ, x, s) # s = A' * r - λ * x - γ_next = kdotr(n, s, s) # γ_next = sᴴs + λ > 0 && kaxpy!(n, -λ, x, s) # s = A' * r - λ * x + γ_next = kdotr(n, s, s) # γ_next = sᴴs β = γ_next / γ kaxpby!(n, one(FC), s, β, p) # p = s + βp γ = γ_next diff --git a/src/cgls_lanczos_shift.jl b/src/cgls_lanczos_shift.jl index 5e8b46801..14caba17d 100644 --- a/src/cgls_lanczos_shift.jl +++ b/src/cgls_lanczos_shift.jl @@ -137,10 +137,10 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose kfill!(x[i], zero(FC)) # x₀ end - kcopy!(m, u, b) # u ← b + kcopy!(m, u, b) # u ← b kfill!(u_prev, zero(FC)) - mul!(v, Aᴴ, u) # v₁ ← Aᴴ * b - β = sqrt(kdotr(n, v, v)) # β₁ = v₁ᵀ M v₁ + mul!(v, Aᴴ, u) # v₁ ← Aᴴ * b + β = sqrt(kdotr(n, v, v)) # β₁ = v₁ᵀ M v₁ kfill!(rNorms, β) if history for i = 1 : nshifts @@ -198,16 +198,16 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose while ! (solved || tired || user_requested_exit || overtimed) # Form next Lanczos vector. - mul!(utilde, A, v) # utildeₖ ← Avₖ - δ = kdotr(m, utilde, utilde) # δₖ = vₖᵀAᴴAvₖ - kaxpy!(m, -δ, u, utilde) # uₖ₊₁ = utildeₖ - δₖuₖ - βₖuₖ₋₁ + mul!(utilde, A, v) # utildeₖ ← Avₖ + δ = kdotr(m, utilde, utilde) # δₖ = vₖᵀAᴴAvₖ + kaxpy!(m, -δ, u, utilde) # uₖ₊₁ = utildeₖ - δₖuₖ - βₖuₖ₋₁ kaxpy!(m, -β, u_prev, utilde) - mul!(v, Aᴴ, utilde) # vₖ₊₁ = Aᴴuₖ₊₁ - β = sqrt(kdotr(n, v, v)) # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁ - kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁ - kscal!(m, one(FC) / β, utilde) # uₖ₊₁ = uₖ₊₁ / βₖ₊₁ - kcopy!(m, u_prev, u) # u_prev ← u - kcopy!(m, u, utilde) # u ← utilde + mul!(v, Aᴴ, utilde) # vₖ₊₁ = Aᴴuₖ₊₁ + β = sqrt(kdotr(n, v, v)) # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁ + kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁ + kscal!(m, one(FC) / β, utilde) # uₖ₊₁ = uₖ₊₁ / βₖ₊₁ + kcopy!(m, u_prev, u) # u_prev ← u + kcopy!(m, u, utilde) # u ← utilde MisI || (ρ = kdotr(n, v, v)) for i = 1 : nshifts diff --git a/src/cgne.jl b/src/cgne.jl index 7a775d701..357691c5f 100644 --- a/src/cgne.jl +++ b/src/cgne.jl @@ -154,7 +154,7 @@ kwargs_cgne = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor kfill!(x, zero(FC)) kcopy!(m, r, b) # r ← b NisI || mulorldiv!(z, N, r, ldiv) - rNorm = knorm(m, r) # Marginally faster than norm(r) + rNorm = knorm(m, r) # Marginally faster than norm(r) history && push!(rNorms, rNorm) if rNorm == 0 stats.niter = 0 @@ -195,8 +195,8 @@ kwargs_cgne = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor δ = kdotr(n, p, p) # Faster than dot(p, p) λ > 0 && (δ += λ * kdotr(m, s, s)) α = γ / δ - kaxpy!(n, α, p, x) # Faster than x = x + α * p - kaxpy!(m, -α, q, r) # Faster than r = r - α * q + kaxpy!(n, α, p, x) # Faster than x = x + α * p + kaxpy!(m, -α, q, r) # Faster than r = r - α * q NisI || mulorldiv!(z, N, r, ldiv) γ_next = kdotr(m, r, z) # Faster than γ_next = dot(r, z) β = γ_next / γ diff --git a/src/cgs.jl b/src/cgs.jl index 7aeb4c871..e93d46928 100644 --- a/src/cgs.jl +++ b/src/cgs.jl @@ -156,7 +156,7 @@ kwargs_cgs = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist kcopy!(n, r₀, b) # r₀ ← b end - kfill!(x, zero(FC)) # x₀ + kfill!(x, zero(FC)) # x₀ MisI || mulorldiv!(r, M, r₀, ldiv) # r₀ # Compute residual norm ‖r₀‖₂. @@ -206,22 +206,22 @@ kwargs_cgs = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist NisI || mulorldiv!(y, N, p, ldiv) # yₖ = N⁻¹pₖ mul!(t, A, y) # tₖ = Ayₖ MisI || mulorldiv!(v, M, t, ldiv) # vₖ = M⁻¹tₖ - σ = kdot(n, c, v) # σₖ = ⟨ r̅₀,M⁻¹AN⁻¹pₖ ⟩ + σ = kdot(n, c, v) # σₖ = ⟨ r̅₀,M⁻¹AN⁻¹pₖ ⟩ α = ρ / σ # αₖ = ρₖ / σₖ - kcopy!(n, q, u) # qₖ = uₖ - kaxpy!(n, -α, v, q) # qₖ = qₖ - αₖ * M⁻¹AN⁻¹pₖ - kaxpy!(n, one(FC), q, u) # uₖ₊½ = uₖ + qₖ + kcopy!(n, q, u) # qₖ = uₖ + kaxpy!(n, -α, v, q) # qₖ = qₖ - αₖ * M⁻¹AN⁻¹pₖ + kaxpy!(n, one(FC), q, u) # uₖ₊½ = uₖ + qₖ NisI || mulorldiv!(z, N, u, ldiv) # zₖ = N⁻¹uₖ₊½ - kaxpy!(n, α, z, x) # xₖ₊₁ = xₖ + αₖ * N⁻¹(uₖ + qₖ) + kaxpy!(n, α, z, x) # xₖ₊₁ = xₖ + αₖ * N⁻¹(uₖ + qₖ) mul!(s, A, z) # sₖ = Azₖ MisI || mulorldiv!(w, M, s, ldiv) # wₖ = M⁻¹sₖ - kaxpy!(n, -α, w, r) # rₖ₊₁ = rₖ - αₖ * M⁻¹AN⁻¹(uₖ + qₖ) - ρ_next = kdot(n, c, r) # ρₖ₊₁ = ⟨ r̅₀,rₖ₊₁ ⟩ + kaxpy!(n, -α, w, r) # rₖ₊₁ = rₖ - αₖ * M⁻¹AN⁻¹(uₖ + qₖ) + ρ_next = kdot(n, c, r) # ρₖ₊₁ = ⟨ r̅₀,rₖ₊₁ ⟩ β = ρ_next / ρ # βₖ = ρₖ₊₁ / ρₖ - kcopy!(n, u, r) # uₖ₊₁ = rₖ₊₁ - kaxpy!(n, β, q, u) # uₖ₊₁ = uₖ₊₁ + βₖ * qₖ - kaxpby!(n, one(FC), q, β, p) # pₐᵤₓ = qₖ + βₖ * pₖ - kaxpby!(n, one(FC), u, β, p) # pₖ₊₁ = uₖ₊₁ + βₖ * pₐᵤₓ + kcopy!(n, u, r) # uₖ₊₁ = rₖ₊₁ + kaxpy!(n, β, q, u) # uₖ₊₁ = uₖ₊₁ + βₖ * qₖ + kaxpby!(n, one(FC), q, β, p) # pₐᵤₓ = qₖ + βₖ * pₖ + kaxpby!(n, one(FC), u, β, p) # pₖ₊₁ = uₖ₊₁ + βₖ * pₐᵤₓ # Update ρ. ρ = ρ_next # ρₖ ← ρₖ₊₁ diff --git a/src/cr.jl b/src/cr.jl index d2d3fae3b..0c1ef450b 100644 --- a/src/cr.jl +++ b/src/cr.jl @@ -153,7 +153,7 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema mul!(Ar, A, r) ρ = kdotr(n, r, Ar) - rNorm = sqrt(kdotr(n, r, p)) # ‖r‖ + rNorm = sqrt(kdotr(n, r, p)) # ‖r‖ history && push!(rNorms, rNorm) # Values of ‖r‖ if ρ == 0 @@ -167,7 +167,7 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema end kcopy!(n, p, r) # p ← r kcopy!(n, q, Ar) # q ← Ar - (verbose > 0) && (m = zero(T)) # quadratic model + (verbose > 0) && (m = zero(T)) # quadratic model iter = 0 itmax == 0 && (itmax = 2 * n) @@ -221,12 +221,12 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema tr = maximum(to_boundary(n, x, r, Mq, radius; flip = false, xNorm2 = xNorm², dNorm2 = rNorm²)) (verbose > 0) && @printf(iostream, "t1 = %8.1e, t2 = %8.1e and tr = %8.1e\n", t1, t2, tr) - if abspAp ≤ γ * pNorm * knorm(n, q) # pᴴAp ≃ 0 - npcurv = true # nonpositive curvature + if abspAp ≤ γ * pNorm * knorm(n, q) # pᴴAp ≃ 0 + npcurv = true # nonpositive curvature (verbose > 0) && @printf(iostream, "pᴴAp = %8.1e ≃ 0\n", pAp) - if abspr ≤ γ * pNorm * rNorm # pᴴr ≃ 0 + if abspr ≤ γ * pNorm * rNorm # pᴴr ≃ 0 (verbose > 0) && @printf(iostream, "pᴴr = %8.1e ≃ 0, redefining p := r\n", pr) - p = r # - ∇q(x) + p = r # - ∇q(x) q = Ar # q(x + αr) = q(x) - α ‖r‖² + ½ α² rᴴAr # 1) if rᴴAr > 0, the quadratic decreases from α = 0 to α = ‖r‖² / rᴴAr @@ -245,7 +245,7 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema α = descent ? t1 : t2 ρ > 0 && (tr = min(tr, rNorm² / ρ)) Δ = -α * pr + tr * rNorm² - (tr)^2 * ρ / 2 # as pᴴAp = 0 - if Δ > 0 # direction r engenders a better decrease + if Δ > 0 # direction r engenders a better decrease (verbose > 0) && @printf(iostream, "direction r engenders a bigger decrease. q_p - q_r = %8.1e > 0\n", Δ) (verbose > 0) && @printf(iostream, "redefining p := r\n") p = r @@ -256,7 +256,7 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema end end - elseif pAp > 0 && ρ > 0 # no negative curvature + elseif pAp > 0 && ρ > 0 # no negative curvature (verbose > 0) && @printf(iostream, "positive curvatures along p and r. pᴴAp = %8.1e and rᴴAr = %8.1e\n", pAp, ρ) α = ρ / kdotr(n, q, Mq) if α ≥ t1 @@ -313,13 +313,13 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema end elseif radius == 0 - α = ρ / kdotr(n, q, Mq) # step + α = ρ / kdotr(n, q, Mq) # step end kaxpy!(n, α, p, x) xNorm = knorm(n, x) xNorm ≈ radius && (on_boundary = true) - kaxpy!(n, -α, Mq, r) # residual + kaxpy!(n, -α, Mq, r) # residual if MisI rNorm² = kdotr(n, r, r) rNorm = sqrt(rNorm²) @@ -372,9 +372,9 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema solver.warm_start = false return solver end - pr = rNorm² + β * pr - β * α * pAp # pᴴr + pr = rNorm² + β * pr - β * α * pAp # pᴴr abspr = abs(pr) - pAp = ρ + β^2 * pAp # pᴴq + pAp = ρ + β^2 * pAp # pᴴq abspAp = abs(pAp) descent = pr > 0 diff --git a/src/craig.jl b/src/craig.jl index c85111bca..de865f1c5 100644 --- a/src/craig.jl +++ b/src/craig.jl @@ -305,7 +305,7 @@ kwargs_craig = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :btol, :conlim, :at # Recur y. kaxpby!(m, one(FC), u, -θ/ρ_prev, w) # w = u - θ/ρ_prev * w - kaxpy!(m, ξ/ρ, w, y) # y = y + ξ/ρ * w + kaxpy!(m, ξ/ρ, w, y) # y = y + ξ/ρ * w Dnorm² += knorm(m, w) diff --git a/src/craigmr.jl b/src/craigmr.jl index ac36ede28..0e87622b7 100644 --- a/src/craigmr.jl +++ b/src/craigmr.jl @@ -233,9 +233,9 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver NisI || kscal!(n, one(FC)/α, Nv) # Regularization. - λₖ = λ # λ₁ = λ - cpₖ = spₖ = one(T) # Givens sines and cosines used to zero out λₖ - cdₖ = sdₖ = one(T) # Givens sines and cosines used to define λₖ₊₁ + λₖ = λ # λ₁ = λ + cpₖ = spₖ = one(T) # Givens sines and cosines used to zero out λₖ + cdₖ = sdₖ = one(T) # Givens sines and cosines used to define λₖ₊₁ λ > 0 && kcopy!(n, q, v) # Additional vector needed to update x, by definition q₀ = 0 if λ > 0 @@ -254,7 +254,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver ArNorm = α history && push!(ArNorms, ArNorm) - ɛ_c = atol + rtol * rNorm # Stopping tolerance for consistent systems. + ɛ_c = atol + rtol * rNorm # Stopping tolerance for consistent systems. ɛ_i = atol + rtol * ArNorm # Stopping tolerance for inconsistent systems. kcopy!(m, wbar, u) diff --git a/src/crls.jl b/src/crls.jl index 167b3df1f..b9c39db53 100644 --- a/src/crls.jl +++ b/src/crls.jl @@ -140,9 +140,9 @@ kwargs_crls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose MAp = MisI ? Ap : solver.Ms kfill!(x, zero(FC)) - kcopy!(m, r, b) # r ← b + kcopy!(m, r, b) # r ← b bNorm = knorm(m, r) # norm(b - A * x0) if x0 ≠ 0. - rNorm = bNorm # + λ * ‖x0‖ if x0 ≠ 0 and λ > 0. + rNorm = bNorm # + λ * ‖x0‖ if x0 ≠ 0 and λ > 0. history && push!(rNorms, rNorm) if bNorm == 0 stats.niter = 0 @@ -160,7 +160,7 @@ kwargs_crls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose kcopy!(n, p, Ar) # p ← Ar kcopy!(m, Ap, s) # Ap ← s - mul!(q, Aᴴ, Ms) # Ap + mul!(q, Aᴴ, Ms) # Ap λ > 0 && kaxpy!(n, λ, p, q) # q = q + λ * p γ = kdotr(m, s, Ms) # Faster than γ = dot(s, Ms) iter = 0 @@ -205,20 +205,20 @@ kwargs_crls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose end end - kaxpy!(n, α, p, x) # Faster than x = x + α * p - kaxpy!(n, -α, q, Ar) # Faster than Ar = Ar - α * q + kaxpy!(n, α, p, x) # Faster than x = x + α * p + kaxpy!(n, -α, q, Ar) # Faster than Ar = Ar - α * q ArNorm = knorm(n, Ar) solved = psd || on_boundary solved && continue - kaxpy!(m, -α, Ap, r) # Faster than r = r - α * Ap + kaxpy!(m, -α, Ap, r) # Faster than r = r - α * Ap mul!(s, A, Ar) MisI || mulorldiv!(Ms, M, s, ldiv) - γ_next = kdotr(m, s, Ms) # Faster than γ_next = dot(s, s) + γ_next = kdotr(m, s, Ms) # Faster than γ_next = dot(s, s) λ > 0 && (γ_next += λ * ArNorm * ArNorm) β = γ_next / γ - kaxpby!(n, one(FC), Ar, β, p) # Faster than p = Ar + β * p - kaxpby!(m, one(FC), s, β, Ap) # Faster than Ap = s + β * Ap + kaxpby!(n, one(FC), Ar, β, p) # Faster than p = Ar + β * p + kaxpby!(m, one(FC), s, β, Ap) # Faster than Ap = s + β * Ap MisI || mulorldiv!(MAp, M, Ap, ldiv) mul!(q, Aᴴ, MAp) λ > 0 && kaxpy!(n, λ, p, q) # q = q + λ * p @@ -227,7 +227,7 @@ kwargs_crls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose if λ > 0 rNorm = sqrt(kdotr(m, r, r) + λ * kdotr(n, x, x)) else - rNorm = knorm(m, r) # norm(r) + rNorm = knorm(m, r) end history && push!(rNorms, rNorm) history && push!(ArNorms, ArNorm) diff --git a/src/crmr.jl b/src/crmr.jl index b2c886fef..73dad4d1a 100644 --- a/src/crmr.jl +++ b/src/crmr.jl @@ -150,9 +150,9 @@ kwargs_crmr = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor reset!(stats) Nq = NisI ? q : solver.Nq - kfill!(x, zero(FC)) # initial estimation x = 0 + kfill!(x, zero(FC)) # initial estimation x = 0 mulorldiv!(r, N, b, ldiv) # initial residual r = N * (b - Ax) = N * b - bNorm = knorm(m, r) # norm(b - A * x0) if x0 ≠ 0. + bNorm = knorm(m, r) # norm(b - A * x0) if x0 ≠ 0. rNorm = bNorm # + λ * ‖x0‖ if x0 ≠ 0 and λ > 0. history && push!(rNorms, rNorm) if bNorm == 0 @@ -164,16 +164,16 @@ kwargs_crmr = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor return solver end λ > 0 && kcopy!(m, s, r) # s ← r - mul!(Aᴴr, Aᴴ, r) # - λ * x0 if x0 ≠ 0. - kcopy!(n, p, Aᴴr) # p ← Aᴴr - γ = kdotr(n, Aᴴr, Aᴴr) # Faster than γ = dot(Aᴴr, Aᴴr) + mul!(Aᴴr, Aᴴ, r) # - λ * x0 if x0 ≠ 0. + kcopy!(n, p, Aᴴr) # p ← Aᴴr + γ = kdotr(n, Aᴴr, Aᴴr) # Faster than γ = dot(Aᴴr, Aᴴr) λ > 0 && (γ += λ * rNorm * rNorm) iter = 0 itmax == 0 && (itmax = m + n) ArNorm = sqrt(γ) history && push!(ArNorms, ArNorm) - ɛ_c = atol + rtol * rNorm # Stopping tolerance for consistent systems. + ɛ_c = atol + rtol * rNorm # Stopping tolerance for consistent systems. ɛ_i = atol + rtol * ArNorm # Stopping tolerance for inconsistent systems. (verbose > 0) && @printf(iostream, "%5s %8s %8s %5s\n", "k", "‖Aᴴr‖", "‖r‖", "timer") kdisplay(iter, verbose) && @printf(iostream, "%5d %8.2e %8.2e %.2fs\n", iter, ArNorm, rNorm, ktimer(start_time)) @@ -189,10 +189,10 @@ kwargs_crmr = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor mul!(q, A, p) λ > 0 && kaxpy!(m, λ, s, q) # q = q + λ * s NisI || mulorldiv!(Nq, N, q, ldiv) - α = γ / kdotr(m, q, Nq) # Compute qᴴ * N * q - kaxpy!(n, α, p, x) # Faster than x = x + α * p - kaxpy!(m, -α, Nq, r) # Faster than r = r - α * Nq - rNorm = knorm(m, r) # norm(r) + α = γ / kdotr(m, q, Nq) # Compute qᴴ * N * q + kaxpy!(n, α, p, x) # Faster than x = x + α * p + kaxpy!(m, -α, Nq, r) # Faster than r = r - α * Nq + rNorm = knorm(m, r) # norm(r) mul!(Aᴴr, Aᴴ, r) γ_next = kdotr(n, Aᴴr, Aᴴr) # Faster than γ_next = dot(Aᴴr, Aᴴr) λ > 0 && (γ_next += λ * rNorm * rNorm) @@ -200,7 +200,7 @@ kwargs_crmr = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor kaxpby!(n, one(FC), Aᴴr, β, p) # Faster than p = Aᴴr + β * p if λ > 0 - kaxpby!(m, one(FC), r, β, s) # s = r + β * s + kaxpby!(m, one(FC), r, β, s) # s = r + β * s end γ = γ_next diff --git a/src/diom.jl b/src/diom.jl index 7d166be72..2c3d2331f 100644 --- a/src/diom.jl +++ b/src/diom.jl @@ -149,7 +149,7 @@ kwargs_diom = (:M, :N, :ldiv, :reorthogonalization, :atol, :rtol, :itmax, :timem kcopy!(n, t, b) # t ← b end MisI || mulorldiv!(r₀, M, t, ldiv) # M(b - Ax₀) - rNorm = knorm(n, r₀) # β = ‖r₀‖₂ + rNorm = knorm(n, r₀) # β = ‖r₀‖₂ history && push!(rNorms, rNorm) if rNorm == 0 stats.niter = 0 @@ -225,7 +225,7 @@ kwargs_diom = (:M, :N, :ldiv, :reorthogonalization, :atol, :rtol, :itmax, :timem end # Compute hₖ₊₁.ₖ and vₖ₊₁. - Haux = knorm(n, w) # hₖ₊₁.ₖ = ‖vₖ₊₁‖₂ + Haux = knorm(n, w) # hₖ₊₁.ₖ = ‖vₖ₊₁‖₂ if Haux ≠ 0 # hₖ₊₁.ₖ = 0 ⇒ "lucky breakdown" V[next_pos] .= w ./ Haux # vₖ₊₁ = w / hₖ₊₁.ₖ end diff --git a/src/dqgmres.jl b/src/dqgmres.jl index a4d0161a1..06b78d5e9 100644 --- a/src/dqgmres.jl +++ b/src/dqgmres.jl @@ -149,7 +149,7 @@ kwargs_dqgmres = (:M, :N, :ldiv, :reorthogonalization, :atol, :rtol, :itmax, :ti kcopy!(n, t, b) # t ← b end MisI || mulorldiv!(r₀, M, t, ldiv) # M(b - Ax₀) - rNorm = knorm(n, r₀) # β = ‖r₀‖₂ + rNorm = knorm(n, r₀) # β = ‖r₀‖₂ history && push!(rNorms, rNorm) if rNorm == 0 stats.niter = 0 @@ -227,7 +227,7 @@ kwargs_dqgmres = (:M, :N, :ldiv, :reorthogonalization, :atol, :rtol, :itmax, :ti end # Compute hₖ₊₁.ₖ and vₖ₊₁. - Haux = knorm(n, w) # hₖ₊₁.ₖ = ‖vₖ₊₁‖₂ + Haux = knorm(n, w) # hₖ₊₁.ₖ = ‖vₖ₊₁‖₂ if Haux ≠ 0 # hₖ₊₁.ₖ = 0 ⇒ "lucky breakdown" V[next_pos] .= w ./ Haux # vₖ₊₁ = w / hₖ₊₁.ₖ end diff --git a/src/fgmres.jl b/src/fgmres.jl index 85c5bc2b9..6b5a3a23b 100644 --- a/src/fgmres.jl +++ b/src/fgmres.jl @@ -155,7 +155,7 @@ kwargs_fgmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :i kcopy!(n, w, b) # w ← b end MisI || mulorldiv!(r₀, M, w, ldiv) # r₀ = M(b - Ax₀) - β = knorm(n, r₀) # β = ‖r₀‖₂ + β = knorm(n, r₀) # β = ‖r₀‖₂ rNorm = β history && push!(rNorms, β) diff --git a/src/fom.jl b/src/fom.jl index dc6943cda..05450862d 100644 --- a/src/fom.jl +++ b/src/fom.jl @@ -150,7 +150,7 @@ kwargs_fom = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itma kcopy!(n, w, b) # w ← b end MisI || mulorldiv!(r₀, M, w, ldiv) # r₀ = M(b - Ax₀) - β = knorm(n, r₀) # β = ‖r₀‖₂ + β = knorm(n, r₀) # β = ‖r₀‖₂ rNorm = β history && push!(rNorms, β) diff --git a/src/gmres.jl b/src/gmres.jl index 5a3a29d35..ce07382e8 100644 --- a/src/gmres.jl +++ b/src/gmres.jl @@ -150,7 +150,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it kcopy!(n, w, b) # w ← b end MisI || mulorldiv!(r₀, M, w, ldiv) # r₀ = M(b - Ax₀) - β = knorm(n, r₀) # β = ‖r₀‖₂ + β = knorm(n, r₀) # β = ‖r₀‖₂ rNorm = β history && push!(rNorms, β) @@ -241,8 +241,8 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it mul!(w, A, p) # w ← ANvₖ MisI || mulorldiv!(q, M, w, ldiv) # q ← MANvₖ for i = 1 : inner_iter - R[nr+i] = kdot(n, V[i], q) # hᵢₖ = (vᵢ)ᴴq - kaxpy!(n, -R[nr+i], V[i], q) # q ← q - hᵢₖvᵢ + R[nr+i] = kdot(n, V[i], q) # hᵢₖ = (vᵢ)ᴴq + kaxpy!(n, -R[nr+i], V[i], q) # q ← q - hᵢₖvᵢ end # Reorthogonalization of the Krylov basis. diff --git a/src/lnlq.jl b/src/lnlq.jl index 6123f232b..58a3cf611 100644 --- a/src/lnlq.jl +++ b/src/lnlq.jl @@ -228,7 +228,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly # β₁Mu₁ = b. kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) # u₁ = M⁻¹ * Mu₁ - βₖ = sqrt(kdotr(m, u, Mu)) # β₁ = ‖u₁‖_M + βₖ = sqrt(kdotr(m, u, Mu)) # β₁ = ‖u₁‖_M if βₖ ≠ 0 kscal!(m, one(FC) / βₖ, u) MisI || kscal!(m, one(FC) / βₖ, Mu) @@ -238,22 +238,22 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly mul!(Aᴴu, Aᴴ, u) kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) # v₁ = N⁻¹ * Nv₁ - αₖ = sqrt(kdotr(n, v, Nv)) # α₁ = ‖v₁‖_N + αₖ = sqrt(kdotr(n, v, Nv)) # α₁ = ‖v₁‖_N if αₖ ≠ 0 kscal!(n, one(FC) / αₖ, v) NisI || kscal!(n, one(FC) / αₖ, Nv) end - kcopy!(m, w̄, u) # Direction w̄₁ + kcopy!(m, w̄, u) # Direction w̄₁ cₖ = zero(T) # Givens cosines used for the LQ factorization of (Lₖ)ᴴ sₖ = zero(FC) # Givens sines used for the LQ factorization of (Lₖ)ᴴ ζₖ₋₁ = zero(FC) # ζₖ₋₁ and ζbarₖ are the last components of z̅ₖ ηₖ = zero(FC) # Coefficient of M̅ₖ # Variable used for the regularization. - λₖ = λ # λ₁ = λ - cpₖ = spₖ = one(T) # Givens sines and cosines used to zero out λₖ - cdₖ = sdₖ = one(FC) # Givens sines and cosines used to define λₖ₊₁ + λₖ = λ # λ₁ = λ + cpₖ = spₖ = one(T) # Givens sines and cosines used to zero out λₖ + cdₖ = sdₖ = one(FC) # Givens sines and cosines used to define λₖ₊₁ λ > 0 && kcopy!(n, q, v) # Additional vector needed to update x, by definition q₀ = 0 # Initialize the regularization. @@ -346,7 +346,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly mul!(Av, A, v) kaxpby!(m, one(FC), Av, -αₖ, Mu) MisI || mulorldiv!(u, M, Mu, ldiv) # uₖ₊₁ = M⁻¹ * Muₖ₊₁ - βₖ₊₁ = sqrt(kdotr(m, u, Mu)) # βₖ₊₁ = ‖uₖ₊₁‖_M + βₖ₊₁ = sqrt(kdotr(m, u, Mu)) # βₖ₊₁ = ‖uₖ₊₁‖_M if βₖ₊₁ ≠ 0 kscal!(m, one(FC) / βₖ₊₁, u) MisI || kscal!(m, one(FC) / βₖ₊₁, Mu) @@ -356,7 +356,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly mul!(Aᴴu, Aᴴ, u) kaxpby!(n, one(FC), Aᴴu, -βₖ₊₁, Nv) NisI || mulorldiv!(v, N, Nv, ldiv) # vₖ₊₁ = N⁻¹ * Nvₖ₊₁ - αₖ₊₁ = sqrt(kdotr(n, v, Nv)) # αₖ₊₁ = ‖vₖ₊₁‖_N + αₖ₊₁ = sqrt(kdotr(n, v, Nv)) # αₖ₊₁ = ‖vₖ₊₁‖_N if αₖ₊₁ ≠ 0 kscal!(n, one(FC) / αₖ₊₁, v) NisI || kscal!(n, one(FC) / αₖ₊₁, Nv) diff --git a/src/lsmr.jl b/src/lsmr.jl index 352d81ee7..14bb646c7 100644 --- a/src/lsmr.jl +++ b/src/lsmr.jl @@ -332,8 +332,8 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, iter ≥ window && (err_lbnd = knorm(window, err_vec)) # Update h, hbar and x. - δ = θbar * ρ / (ρold * ρbarold) # δₖ = θbarₖ * ρₖ / (ρₖ₋₁ * ρbarₖ₋₁) - kaxpby!(n, one(FC), h, -δ, hbar) # ĥₖ = hₖ - δₖ * ĥₖ₋₁ + δ = θbar * ρ / (ρold * ρbarold) # δₖ = θbarₖ * ρₖ / (ρₖ₋₁ * ρbarₖ₋₁) + kaxpby!(n, one(FC), h, -δ, hbar) # ĥₖ = hₖ - δₖ * ĥₖ₋₁ # if a trust-region constraint is given, compute step to the boundary # the step ϕ/ρ is not necessarily positive @@ -345,8 +345,8 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, σ = σ > 0 ? min(σ, tmax) : max(σ, tmin) end - kaxpy!(n, σ, hbar, x) # xₖ = xₖ₋₁ + σₖ * ĥₖ - kaxpby!(n, one(FC), v, -θnew / ρ, h) # hₖ₊₁ = vₖ₊₁ - (θₖ₊₁/ρₖ) * hₖ + kaxpy!(n, σ, hbar, x) # xₖ = xₖ₋₁ + σₖ * ĥₖ + kaxpby!(n, one(FC), v, -θnew / ρ, h) # hₖ₊₁ = vₖ₊₁ - (θₖ₊₁/ρₖ) * hₖ # Estimate ‖r‖. βacute = chat * βdd diff --git a/src/minres.jl b/src/minres.jl index 8fd711b28..2b5f6ce51 100644 --- a/src/minres.jl +++ b/src/minres.jl @@ -237,9 +237,9 @@ kwargs_minres = (:M, :ldiv, :λ, :atol, :rtol, :etol, :conlim, :itmax, :timemax, # Generate next Lanczos vector. mul!(y, A, v) - λ ≠ 0 && kaxpy!(n, λ, v, y) # (y = y + λ * v) + λ ≠ 0 && kaxpy!(n, λ, v, y) # (y = y + λ * v) kscal!(n, one(FC) / β, y) - iter ≥ 2 && kaxpy!(n, -β / oldβ, r1, y) # (y = y - β / oldβ * r1) + iter ≥ 2 && kaxpy!(n, -β / oldβ, r1, y) # (y = y - β / oldβ * r1) α = kdotr(n, v, y) / β kaxpy!(n, -α / β, r2, y) # y = y - α / β * r2 diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl index 78532c284..085f7d728 100644 --- a/src/minres_qlp.jl +++ b/src/minres_qlp.jl @@ -216,13 +216,13 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve # M(A + λI)Vₖ = Vₖ₊₁Tₖ₊₁.ₖ # βₖ₊₁vₖ₊₁ = M(A + λI)vₖ - αₖvₖ - βₖvₖ₋₁ - mul!(p, A, vₖ) # p ← Avₖ + mul!(p, A, vₖ) # p ← Avₖ if λ ≠ 0 kaxpy!(n, λ, vₖ, p) # p ← p + λvₖ end if iter ≥ 2 - kaxpy!(n, -βₖ, M⁻¹vₖ₋₁, p) # p ← p - βₖ * M⁻¹vₖ₋₁ + kaxpy!(n, -βₖ, M⁻¹vₖ₋₁, p) # p ← p - βₖ * M⁻¹vₖ₋₁ end αₖ = kdotr(n, vₖ, p) # αₖ = ⟨vₖ,p⟩ diff --git a/src/qmr.jl b/src/qmr.jl index 22ed8f29f..a05c052c6 100644 --- a/src/qmr.jl +++ b/src/qmr.jl @@ -197,16 +197,16 @@ kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist βₖ = √(abs(cᴴb)) # β₁γ₁ = cᴴ(b - Ax₀) γₖ = cᴴb / βₖ # β₁γ₁ = cᴴ(b - Ax₀) - kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 - kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 + kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 + kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 vₖ .= r₀ ./ βₖ # v₁ = (b - Ax₀) / β₁ uₖ .= c ./ conj(γₖ) # u₁ = c / γ̄₁ cₖ₋₂ = cₖ₋₁ = cₖ = zero(T) # Givens cosines used for the QR factorization of Tₖ₊₁.ₖ sₖ₋₂ = sₖ₋₁ = sₖ = zero(FC) # Givens sines used for the QR factorization of Tₖ₊₁.ₖ - kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Vₖ(Rₖ)⁻¹ - kfill!(wₖ₋₁, zero(FC)) # Column k-1 of Wₖ = Vₖ(Rₖ)⁻¹ + kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Vₖ(Rₖ)⁻¹ + kfill!(wₖ₋₁, zero(FC)) # Column k-1 of Wₖ = Vₖ(Rₖ)⁻¹ ζbarₖ = βₖ # ζbarₖ is the last component of z̅ₖ = (Qₖ)ᴴβ₁e₁ - τₖ = kdotr(n, vₖ, vₖ) # τₖ is used for the residual norm estimate + τₖ = kdotr(n, vₖ, vₖ) # τₖ is used for the residual norm estimate # Stopping criterion. solved = rNorm ≤ ε @@ -243,8 +243,8 @@ kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist 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 / βₖ₊₁ + βₖ₊₁ = √(abs(pᴴq)) # βₖ₊₁ = √(|pᴴq|) + γₖ₊₁ = pᴴq / βₖ₊₁ # γₖ₊₁ = pᴴq / βₖ₊₁ # Update the QR factorization of Tₖ₊₁.ₖ = Qₖ [ Rₖ ]. # [ Oᵀ ] diff --git a/src/symmlq.jl b/src/symmlq.jl index 31b3a53d5..d19559c42 100644 --- a/src/symmlq.jl +++ b/src/symmlq.jl @@ -238,7 +238,7 @@ kwargs_symmlq = (:M, :ldiv, :transfer_to_cg, :λ, :λest, :atol, :rtol, :etol, : cw = ρbar / ρ sw = β / ρ - history && push!(errors, abs(β₁/λest)) + history && push!(errors, abs(β₁ / λest)) if γbar ≠ 0 history && push!(errorscg, sqrt(errors[1]^2 - ζbar^2)) else diff --git a/src/trilqr.jl b/src/trilqr.jl index 90a8fbb33..b609cd914 100644 --- a/src/trilqr.jl +++ b/src/trilqr.jl @@ -160,10 +160,10 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %7.1e %.2fs\n", iter, bNorm, cNorm, ktimer(start_time)) # Set up workspace. - βₖ = knorm(m, r₀) # β₁ = ‖r₀‖ = ‖v₁‖ - γₖ = knorm(n, s₀) # γ₁ = ‖s₀‖ = ‖u₁‖ - kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 - kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 + βₖ = knorm(m, r₀) # β₁ = ‖r₀‖ = ‖v₁‖ + γₖ = knorm(n, s₀) # γ₁ = ‖s₀‖ = ‖u₁‖ + kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 + kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 vₖ .= r₀ ./ βₖ # v₁ = (b - Ax₀) / β₁ uₖ .= s₀ ./ γₖ # u₁ = (c - Aᴴy₀) / γ₁ cₖ₋₁ = cₖ = -one(T) # Givens cosines used for the LQ factorization of Tₖ @@ -174,8 +174,8 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, δ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₁ ϵₖ₋₃ = λₖ₋₂ = zero(FC) # Components of Lₖ₋₁ - kfill!(wₖ₋₃, zero(FC)) # Column k-3 of Wₖ = Vₖ(Lₖ)⁻ᴴ - kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Vₖ(Lₖ)⁻ᴴ + kfill!(wₖ₋₃, zero(FC)) # Column k-3 of Wₖ = Vₖ(Lₖ)⁻ᴴ + kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Vₖ(Lₖ)⁻ᴴ # Stopping criterion. inconsistent = false @@ -204,13 +204,13 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, kaxpy!(m, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁ kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - βₖ * uₖ₋₁ - αₖ = kdot(m, vₖ, q) # αₖ = ⟨vₖ,q⟩ + αₖ = kdot(m, vₖ, q) # αₖ = ⟨vₖ,q⟩ - kaxpy!(m, - αₖ , vₖ, q) # q ← q - αₖ * vₖ - kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ + kaxpy!(m, - αₖ , vₖ, q) # q ← q - αₖ * vₖ + kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ - βₖ₊₁ = knorm(m, q) # βₖ₊₁ = ‖q‖ - γₖ₊₁ = knorm(n, p) # γₖ₊₁ = ‖p‖ + βₖ₊₁ = knorm(m, q) # βₖ₊₁ = ‖q‖ + γₖ₊₁ = knorm(n, p) # γₖ₊₁ = ‖p‖ # Update the LQ factorization of Tₖ = L̅ₖQₖ. # [ α₁ γ₂ 0 • • • 0 ] [ δ₁ 0 • • • • 0 ] diff --git a/src/usymlq.jl b/src/usymlq.jl index b16b72145..db7569e4a 100644 --- a/src/usymlq.jl +++ b/src/usymlq.jl @@ -166,15 +166,15 @@ kwargs_usymlq = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, (verbose > 0) && @printf(iostream, "%5s %7s %5s\n", "k", "‖rₖ‖", "timer") kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %.2fs\n", iter, bNorm, ktimer(start_time)) - βₖ = knorm(m, r₀) # β₁ = ‖v₁‖ = ‖r₀‖ - γₖ = knorm(n, c) # γ₁ = ‖u₁‖ = ‖c‖ - kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 - kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 + βₖ = knorm(m, r₀) # β₁ = ‖v₁‖ = ‖r₀‖ + γₖ = knorm(n, c) # γ₁ = ‖u₁‖ = ‖c‖ + kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 + kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 vₖ .= r₀ ./ βₖ # v₁ = (b - Ax₀) / β₁ uₖ .= c ./ γₖ # 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ₖ - kfill!(d̅, zero(FC)) # Last column of D̅ₖ = Uₖ(Qₖ)ᴴ + kfill!(d̅, zero(FC)) # Last column of D̅ₖ = Uₖ(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 @@ -203,8 +203,8 @@ kwargs_usymlq = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, αₖ = kdot(m, vₖ, q) # αₖ = ⟨vₖ,q⟩ - kaxpy!(m, - αₖ , vₖ, q) # q ← q - αₖ * vₖ - kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ + kaxpy!(m, - αₖ , vₖ, q) # q ← q - αₖ * vₖ + kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ βₖ₊₁ = knorm(m, q) # βₖ₊₁ = ‖q‖ γₖ₊₁ = knorm(n, p) # γₖ₊₁ = ‖p‖ diff --git a/src/usymqr.jl b/src/usymqr.jl index f366c63e7..812f30c81 100644 --- a/src/usymqr.jl +++ b/src/usymqr.jl @@ -170,16 +170,16 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, (verbose > 0) && @printf(iostream, "%5s %7s %8s %5s\n", "k", "‖rₖ‖", "‖Aᴴrₖ₋₁‖", "timer") kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %8s %.2fs\n", iter, rNorm, " ✗ ✗ ✗ ✗", ktimer(start_time)) - βₖ = knorm(m, r₀) # β₁ = ‖v₁‖ = ‖r₀‖ - γₖ = knorm(n, c) # γ₁ = ‖u₁‖ = ‖c‖ - kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 - kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 + βₖ = knorm(m, r₀) # β₁ = ‖v₁‖ = ‖r₀‖ + γₖ = knorm(n, c) # γ₁ = ‖u₁‖ = ‖c‖ + kfill!(vₖ₋₁, zero(FC)) # v₀ = 0 + kfill!(uₖ₋₁, zero(FC)) # u₀ = 0 vₖ .= r₀ ./ βₖ # v₁ = (b - Ax₀) / β₁ uₖ .= c ./ γₖ # u₁ = c / γ₁ cₖ₋₂ = cₖ₋₁ = cₖ = one(T) # Givens cosines used for the QR factorization of Tₖ₊₁.ₖ sₖ₋₂ = sₖ₋₁ = sₖ = zero(FC) # Givens sines used for the QR factorization of Tₖ₊₁.ₖ - kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Uₖ(Rₖ)⁻¹ - kfill!(wₖ₋₁, zero(FC)) # Column k-1 of Wₖ = Uₖ(Rₖ)⁻¹ + kfill!(wₖ₋₂, zero(FC)) # Column k-2 of Wₖ = Uₖ(Rₖ)⁻¹ + kfill!(wₖ₋₁, zero(FC)) # Column k-1 of Wₖ = Uₖ(Rₖ)⁻¹ ζbarₖ = βₖ # ζbarₖ is the last component of z̅ₖ = (Qₖ)ᴴβ₁e₁ # Stopping criterion. @@ -201,16 +201,16 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, mul!(q, A , uₖ) # Forms vₖ₊₁ : q ← Auₖ mul!(p, Aᴴ, vₖ) # Forms uₖ₊₁ : p ← Aᴴvₖ - kaxpy!(m, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁ - kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - βₖ * uₖ₋₁ + kaxpy!(m, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁ + kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - βₖ * uₖ₋₁ - αₖ = kdot(m, vₖ, q) # αₖ = ⟨vₖ,q⟩ + αₖ = kdot(m, vₖ, q) # αₖ = ⟨vₖ,q⟩ - kaxpy!(m, - αₖ , vₖ, q) # q ← q - αₖ * vₖ - kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ + kaxpy!(m, - αₖ , vₖ, q) # q ← q - αₖ * vₖ + kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ - βₖ₊₁ = knorm(m, q) # βₖ₊₁ = ‖q‖ - γₖ₊₁ = knorm(n, p) # γₖ₊₁ = ‖p‖ + βₖ₊₁ = knorm(m, q) # βₖ₊₁ = ‖q‖ + γₖ₊₁ = knorm(n, p) # γₖ₊₁ = ‖p‖ # Update the QR factorization of Tₖ₊₁.ₖ = Qₖ [ Rₖ ]. # [ Oᵀ ] @@ -298,10 +298,10 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ if βₖ₊₁ ≠ zero(T) - vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q + vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q end if γₖ₊₁ ≠ zero(T) - uₖ .= p ./ γₖ₊₁ # γₖ₊₁uₖ₊₁ = p + uₖ .= p ./ γₖ₊₁ # γₖ₊₁uₖ₊₁ = p end # Update directions for x.