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

Remove old macros #907

Merged
merged 5 commits into from
Oct 17, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
62 changes: 31 additions & 31 deletions src/bicgstab.jl
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) # ρ₀
Ξ± = 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 @@ -206,27 +206,27 @@ 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β‚–
Ο‰ = @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ₐᡀₓ
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ₐᡀₓ

# 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
50 changes: 25 additions & 25 deletions src/bilq.jl
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,17 +230,17 @@ 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⟩
Ξ²β‚–β‚Šβ‚ = √(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 @@ -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
Loading