Skip to content

Commit

Permalink
Remove old macros
Browse files Browse the repository at this point in the history
  • Loading branch information
@amontoison
amontoison committed Oct 17, 2024
1 parent 40039bf commit c0b4539
Show file tree
Hide file tree
Showing 46 changed files with 1,080 additions and 1,127 deletions.
42 changes: 21 additions & 21 deletions src/bicgstab.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,23 +149,23 @@ kwargs_bicgstab = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose,

if warm_start
mul!(r₀, A, Δx)
@kaxpby!(n, one(FC), b, -one(FC), r₀)
kaxpby!(n, one(FC), b, -one(FC), r₀)
else
@kcopy!(n, r₀, b) # r₀ ← b
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) # ρ₀

# Compute residual norm ‖r₀‖₂.
rNorm = @knrm2(n, r)
rNorm = knorm(n, r)
history && push!(rNorms, rNorm)
if rNorm == 0
stats.niter = 0
Expand All @@ -183,7 +183,7 @@ kwargs_bicgstab = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose,
(verbose > 0) && @printf(iostream, "%5s %7s %8s %8s %5s\n", "k", "‖rₖ‖", "|αₖ|", "|ωₖ|", "timer")
kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %8.1e %8.1e %.2fs\n", iter, rNorm, abs(α), abs(ω), ktimer(start_time))

next_ρ = @kdot(n, c, r) # ρ₁ = ⟨r̅₀,r₀⟩
next_ρ = kdot(n, c, r) # ρ₁ = ⟨r̅₀,r₀⟩
if next_ρ == 0
stats.niter = 0
stats.solved, stats.inconsistent = false, false
Expand All @@ -209,24 +209,24 @@ kwargs_bicgstab = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose,
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ₖ
α = ρ / 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ₖ⟩
ω = 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, -ω, v, p) # pₐᵤₓ = pₖ - ωₖvₖ
kaxpby!(n, one(FC), r, β, p) # pₖ₊₁ = rₖ₊₁ + βₖ₊₁pₐᵤₓ

# Compute residual norm ‖rₖ‖₂.
rNorm = @knrm2(n, r)
rNorm = knorm(n, r)
history && push!(rNorms, rNorm)

# Stopping conditions that do not depend on user input.
Expand All @@ -253,7 +253,7 @@ kwargs_bicgstab = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose,
overtimed && (status = "time limit exceeded")

# Update x
warm_start && @kaxpy!(n, one(FC), Δx, x)
warm_start && kaxpy!(n, one(FC), Δx, x)
solver.warm_start = false

# Update stats
Expand Down
46 changes: 23 additions & 23 deletions src/bilq.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,16 +148,16 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time

if warm_start
mul!(r₀, A, Δx)
@kaxpby!(n, one(FC), b, -one(FC), r₀)
kaxpby!(n, one(FC), b, -one(FC), r₀)
end
if !MisI
mulorldiv!(solver.t, M, r₀, ldiv)
r₀ = solver.t
end

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

history && push!(rNorms, bNorm)
if bNorm == 0
Expand All @@ -174,7 +174,7 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time
itmax == 0 && (itmax = 2*n)

# Initialize the Lanczos biorthogonalization process.
cᴴb = @kdot(n, c, r₀) # ⟨c,r₀⟩
cᴴb = kdot(n, c, r₀) # ⟨c,r₀⟩
if cᴴb == 0
stats.niter = 0
stats.solved = false
Expand All @@ -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
Expand Down Expand Up @@ -230,15 +230,15 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time
mul!(s, Aᴴ, Mᴴuₖ)
NisI || mulorldiv!(p, Nᴴ, s, ldiv)

@kaxpy!(n, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁
@kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - β̄ₖ * uₖ₋₁
kaxpy!(n, -γₖ, vₖ₋₁, q) # q ← q - γₖ * vₖ₋₁
kaxpy!(n, -βₖ, uₖ₋₁, p) # p ← p - β̄ₖ * uₖ₋₁

αₖ = @kdot(n, uₖ, q) # αₖ = ⟨uₖ,q⟩
αₖ = kdot(n, uₖ, q) # αₖ = ⟨uₖ,q⟩

@kaxpy!(n, - αₖ , vₖ, q) # q ← q - αₖ * vₖ
@kaxpy!(n, -conj(αₖ), uₖ, p) # p ← p - ᾱₖ * uₖ
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⟩
pᴴq = kdot(n, p, q) # pᴴq = ⟨p,q⟩
βₖ₊₁ = √(abs(pᴴq)) # βₖ₊₁ = √(|pᴴq|)
γₖ₊₁ = pᴴq / βₖ₊₁ # γₖ₊₁ = pᴴq / βₖ₊₁

Expand Down Expand Up @@ -301,31 +301,31 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time
if iter ≥ 2
# Compute solution xₖ.
# (xᴸ)ₖ₋₁ ← (xᴸ)ₖ₋₂ + ζₖ₋₁ * dₖ₋₁
@kaxpy!(n, ζₖ₋₁ * cₖ, d̅, x)
@kaxpy!(n, ζₖ₋₁ * sₖ, vₖ, x)
kaxpy!(n, ζₖ₋₁ * cₖ, d̅, x)
kaxpy!(n, ζₖ₋₁ * sₖ, vₖ, x)
end

# Compute d̅ₖ.
if iter == 1
# d̅₁ = v₁
@kcopy!(n, d̅, vₖ) # d̅ ← vₖ
kcopy!(n, d̅, vₖ) # d̅ ← vₖ
else
# d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ
@kaxpby!(n, -cₖ, vₖ, conj(sₖ), d̅)
kaxpby!(n, -cₖ, vₖ, conj(sₖ), d̅)
end

# Compute vₖ₊₁ and uₖ₊₁.
@kcopy!(n, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ
@kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ
kcopy!(n, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ
kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ

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ₖ)
norm_vₖ₊₁ = @knrm2(n, vₖ)
vₖᴴvₖ₊₁ = kdot(n, vₖ₋₁, vₖ)
norm_vₖ₊₁ = knorm(n, vₖ)

# Compute BiLQ residual norm
# ‖rₖ‖ = √(|μₖ|²‖vₖ‖² + |ωₖ|²‖vₖ₊₁‖² + μ̄ₖωₖ⟨vₖ,vₖ₊₁⟩ + μₖω̄ₖ⟨vₖ₊₁,vₖ⟩)
Expand Down Expand Up @@ -370,7 +370,7 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time
# Compute BICG point
# (xᶜ)ₖ ← (xᴸ)ₖ₋₁ + ζbarₖ * d̅ₖ
if solved_cg
@kaxpy!(n, ζbarₖ, d̅, x)
kaxpy!(n, ζbarₖ, d̅, x)
end

# Termination status
Expand All @@ -386,7 +386,7 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time
copyto!(solver.s, x)
mulorldiv!(x, N, solver.s, ldiv)
end
warm_start && @kaxpy!(n, one(FC), Δx, x)
warm_start && kaxpy!(n, one(FC), Δx, x)
solver.warm_start = false

# Update stats
Expand Down
Loading

0 comments on commit c0b4539

Please sign in to comment.