From ffdd08c839ebe4f8b3f259f4979d5ce534b6b6a7 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 20 May 2024 17:26:27 -0400 Subject: [PATCH 1/4] Don't use the macro @. --- src/bilq.jl | 10 +++++----- src/bilqr.jl | 16 ++++++++-------- src/cg_lanczos.jl | 4 ++-- src/cg_lanczos_shift.jl | 4 ++-- src/fgmres.jl | 4 ++-- src/fom.jl | 4 ++-- src/gmres.jl | 4 ++-- src/gpmr.jl | 8 ++++---- src/minres.jl | 4 ++-- src/minres_qlp.jl | 4 ++-- src/qmr.jl | 14 +++++++------- src/symmlq.jl | 4 ++-- src/tricg.jl | 18 +++++++++--------- src/trilqr.jl | 16 ++++++++-------- src/trimr.jl | 22 +++++++++++----------- src/usymlq.jl | 10 +++++----- src/usymqr.jl | 14 +++++++------- 17 files changed, 80 insertions(+), 80 deletions(-) diff --git a/src/bilq.jl b/src/bilq.jl index 45a0051ef..6153c59ab 100644 --- a/src/bilq.jl +++ b/src/bilq.jl @@ -329,19 +329,19 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time # Compute d̅ₖ. if iter == 1 # d̅₁ = v₁ - @. d̅ = vₖ + @kcopy!(n, d̅, vₖ) else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ @kaxpby!(n, -cₖ, vₖ, conj(sₖ), d̅) end # Compute vₖ₊₁ and uₖ₊₁. - @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ - @. 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 + vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q + uₖ .= p ./ conj(γₖ₊₁) # γ̄ₖ₊₁uₖ₊₁ = p end # Compute ⟨vₖ,vₖ₊₁⟩ and ‖vₖ₊₁‖ diff --git a/src/bilqr.jl b/src/bilqr.jl index 486ccceec..01bd7b8fb 100644 --- a/src/bilqr.jl +++ b/src/bilqr.jl @@ -316,7 +316,7 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi # Compute d̅ₖ. if iter == 1 # d̅₁ = v₁ - @. d̅ = vₖ + @kcopy!(n, d̅, vₖ) else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ @kaxpby!(n, -cₖ, vₖ, conj(sₖ), d̅) @@ -375,14 +375,14 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi if iter == 2 wₖ₋₁ = wₖ₋₂ @kaxpy!(n, one(FC), uₖ₋₁, wₖ₋₁) - @. wₖ₋₁ = uₖ₋₁ / conj(δₖ₋₁) + wₖ₋₁ .= uₖ₋₁ ./ conj(δₖ₋₁) end # w₂ = (u₂ - λ̄₁w₁) / δ̄₂ if iter == 3 wₖ₋₁ = wₖ₋₃ @kaxpy!(n, one(FC), uₖ₋₁, wₖ₋₁) @kaxpy!(n, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁) - @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁) + wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁) end # wₖ₋₁ = (uₖ₋₁ - λ̄ₖ₋₂wₖ₋₂ - ϵ̄ₖ₋₃wₖ₋₃) / δ̄ₖ₋₁ if iter ≥ 4 @@ -390,7 +390,7 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi wₖ₋₁ = wₖ₋₃ @kaxpy!(n, one(FC), uₖ₋₁, wₖ₋₁) @kaxpy!(n, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁) - @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁) + wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁) end if iter ≥ 3 @@ -421,12 +421,12 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi end # Compute vₖ₊₁ and uₖ₊₁. - @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ - @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ + @kcopy!(n, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ + @kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ if pᴴq ≠ zero(FC) - @. vₖ = q / βₖ₊₁ # βₖ₊₁vₖ₊₁ = q - @. uₖ = p / conj(γₖ₊₁) # γ̄ₖ₊₁uₖ₊₁ = p + vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q + uₖ .= p ./ conj(γₖ₊₁) # γ̄ₖ₊₁uₖ₊₁ = p end # Update ϵₖ₋₃, λₖ₋₂, δbarₖ₋₁, cₖ₋₁, sₖ₋₁, γₖ and βₖ. diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl index 2c5d72a64..892a37527 100644 --- a/src/cg_lanczos.jl +++ b/src/cg_lanczos.jl @@ -218,9 +218,9 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax @kaxpy!(n, -δ, Mv, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ if iter > 0 @kaxpy!(n, -β, Mv_prev, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - βₖMvₖ₋₁ - @. Mv_prev = Mv # Mvₖ₋₁ ← Mvₖ + @kcopy!(n, Mv_prev, Mv) # Mvₖ₋₁ ← Mvₖ end - @. Mv = Mv_next # Mvₖ ← Mvₖ₊₁ + @kcopy!(n, Mv, Mv_next) # Mvₖ ← 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ₖ₊₁ / βₖ₊₁ diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl index b523e5cc3..b38accf2c 100644 --- a/src/cg_lanczos_shift.jl +++ b/src/cg_lanczos_shift.jl @@ -211,9 +211,9 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t @kaxpy!(n, -δ, Mv, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ if iter > 0 @kaxpy!(n, -β, Mv_prev, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - βₖMvₖ₋₁ - @. Mv_prev = Mv # Mvₖ₋₁ ← Mvₖ + kcopy!(n, Mv_prev, Mv) # Mvₖ₋₁ ← Mvₖ end - @. Mv = Mv_next # Mvₖ ← Mvₖ₊₁ + kcopy!(n, Mv, Mv_next) # Mvₖ ← 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ₖ₊₁ / βₖ₊₁ diff --git a/src/fgmres.jl b/src/fgmres.jl index 1a68aac6c..0972f9274 100644 --- a/src/fgmres.jl +++ b/src/fgmres.jl @@ -241,7 +241,7 @@ kwargs_fgmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :i # Initial ζ₁ and V₁ β = @knrm2(n, r₀) z[1] = β - @. V[1] = r₀ / rNorm + V[1] .= r₀ ./ rNorm npass = npass + 1 solver.inner_iter = 0 @@ -332,7 +332,7 @@ kwargs_fgmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :i push!(V, S(undef, n)) push!(z, zero(FC)) end - @. V[inner_iter+1] = q / Hbis # hₖ₊₁.ₖvₖ₊₁ = q + V[inner_iter+1] .= q ./ Hbis # hₖ₊₁.ₖvₖ₊₁ = q z[inner_iter+1] = ζₖ₊₁ end end diff --git a/src/fom.jl b/src/fom.jl index 351fb246f..c27485f2f 100644 --- a/src/fom.jl +++ b/src/fom.jl @@ -233,7 +233,7 @@ kwargs_fom = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itma # Initial ζ₁ and V₁ β = @knrm2(n, r₀) z[1] = β - @. V[1] = r₀ / rNorm + V[1] .= r₀ ./ rNorm npass = npass + 1 inner_iter = 0 @@ -314,7 +314,7 @@ kwargs_fom = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itma if !restart && (inner_iter ≥ mem) push!(V, S(undef, n)) end - @. V[inner_iter+1] = q / Hbis # hₖ₊₁.ₖvₖ₊₁ = q + V[inner_iter+1] .= q ./ Hbis # hₖ₊₁.ₖvₖ₊₁ = q end end diff --git a/src/gmres.jl b/src/gmres.jl index 7ee6e2341..a172be6d9 100644 --- a/src/gmres.jl +++ b/src/gmres.jl @@ -235,7 +235,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it # Initial ζ₁ and V₁ β = @knrm2(n, r₀) z[1] = β - @. V[1] = r₀ / rNorm + V[1] .= r₀ ./ rNorm npass = npass + 1 solver.inner_iter = 0 @@ -325,7 +325,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it push!(V, S(undef, n)) push!(z, zero(FC)) end - @. V[inner_iter+1] = q / Hbis # hₖ₊₁.ₖvₖ₊₁ = q + V[inner_iter+1] .= q ./ Hbis # hₖ₊₁.ₖvₖ₊₁ = q z[inner_iter+1] = ζₖ₊₁ end end diff --git a/src/gpmr.jl b/src/gpmr.jl index 1049c3b50..b47d9d32e 100644 --- a/src/gpmr.jl +++ b/src/gpmr.jl @@ -261,12 +261,12 @@ kwargs_gpmr = (:C, :D, :E, :F, :ldiv, :gsp, :λ, :μ, :reorthogonalization, :ato # βv₁ = Cb β = @knrm2(m, b₀) β ≠ 0 || error("b must be nonzero") - @. V[1] = b₀ / β + V[1] .= b₀ ./ β # γu₁ = Dc γ = @knrm2(n, c₀) γ ≠ 0 || error("c must be nonzero") - @. U[1] = c₀ / γ + U[1] .= c₀ ./ γ # Compute ‖r₀‖² = γ² + β² rNorm = sqrt(γ^2 + β^2) @@ -479,7 +479,7 @@ kwargs_gpmr = (:C, :D, :E, :F, :ldiv, :gsp, :λ, :μ, :reorthogonalization, :ato # hₖ₊₁.ₖ ≠ 0 if Haux > btol - @. V[k+1] = q / Haux # hₖ₊₁.ₖvₖ₊₁ = q + V[k+1] .= q ./ Haux # hₖ₊₁.ₖvₖ₊₁ = q else # Breakdown -- hₖ₊₁.ₖ = ‖q‖₂ = 0 and Auₖ ∈ Span{v₁, ..., vₖ} V[k+1] .= zero(FC) # vₖ₊₁ = 0 such that vₖ₊₁ ⊥ Span{v₁, ..., vₖ} @@ -487,7 +487,7 @@ kwargs_gpmr = (:C, :D, :E, :F, :ldiv, :gsp, :λ, :μ, :reorthogonalization, :ato # fₖ₊₁.ₖ ≠ 0 if Faux > btol - @. U[k+1] = p / Faux # fₖ₊₁.ₖuₖ₊₁ = p + U[k+1] .= p ./ Faux # fₖ₊₁.ₖuₖ₊₁ = p else # Breakdown -- fₖ₊₁.ₖ = ‖p‖₂ = 0 and Bvₖ ∈ Span{u₁, ..., uₖ} U[k+1] .= zero(FC) # uₖ₊₁ = 0 such that uₖ₊₁ ⊥ Span{u₁, ..., uₖ} diff --git a/src/minres.jl b/src/minres.jl index d6ca85d8c..30f248ec3 100644 --- a/src/minres.jl +++ b/src/minres.jl @@ -276,8 +276,8 @@ kwargs_minres = (:M, :ldiv, :λ, :atol, :rtol, :etol, :conlim, :itmax, :timemax, end @kaxpy!(n, one(FC) / β, v, w) - @. r1 = r2 - @. r2 = y + @kcopy!(n, r1, r2) + @kcopy!(n, r2, y) MisI || mulorldiv!(v, M, r2, ldiv) oldβ = β β = @kdotr(n, r2, v) diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl index 5bc3399eb..48f189c99 100644 --- a/src/minres_qlp.jl +++ b/src/minres_qlp.jl @@ -369,12 +369,12 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve # Compute directions wₖ₋₂, ẘₖ₋₁ and w̄ₖ, last columns of Wₖ = Vₖ(Pₖ)ᴴ if iter == 1 # w̅₁ = v₁ - @. wₖ = vₖ + @kcopy!(n, wₖ, vₖ) elseif iter == 2 # [w̅ₖ₋₁ vₖ] [cpₖ spₖ] = [ẘₖ₋₁ w̅ₖ] ⟷ ẘₖ₋₁ = cpₖ * w̅ₖ₋₁ + spₖ * vₖ # [spₖ -cpₖ] ⟷ w̅ₖ = spₖ * w̅ₖ₋₁ - cpₖ * vₖ @kswap(wₖ₋₁, wₖ) - @. wₖ = spₖ * wₖ₋₁ - cpₖ * vₖ + wₖ .= spₖ .* wₖ₋₁ .- cpₖ .* vₖ @kaxpby!(n, spₖ, vₖ, cpₖ, wₖ₋₁) else # [ẘₖ₋₂ w̄ₖ₋₁ vₖ] [cpₖ 0 spₖ] [1 0 0 ] = [wₖ₋₂ ẘₖ₋₁ w̄ₖ] ⟷ wₖ₋₂ = cpₖ * ẘₖ₋₂ + spₖ * vₖ diff --git a/src/qmr.jl b/src/qmr.jl index a197967a8..7ce36b82f 100644 --- a/src/qmr.jl +++ b/src/qmr.jl @@ -326,14 +326,14 @@ kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist if iter == 1 wₖ = wₖ₋₁ @kaxpy!(n, one(FC), vₖ, wₖ) - @. wₖ = wₖ / δₖ + wₖ .= wₖ ./ δₖ end # w₂ = (v₂ - λ₁w₁) / δ₂ if iter == 2 wₖ = wₖ₋₂ @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) @kaxpy!(n, one(FC), vₖ, wₖ) - @. wₖ = wₖ / δₖ + wₖ .= wₖ ./ δₖ end # wₖ = (vₖ - λₖ₋₁wₖ₋₁ - ϵₖ₋₂wₖ₋₂) / δₖ if iter ≥ 3 @@ -341,7 +341,7 @@ kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist wₖ = wₖ₋₂ @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) @kaxpy!(n, one(FC), vₖ, wₖ) - @. wₖ = wₖ / δₖ + wₖ .= wₖ ./ δₖ end # Compute solution xₖ. @@ -349,12 +349,12 @@ kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist @kaxpy!(n, ζₖ, wₖ, x) # Compute vₖ₊₁ and uₖ₊₁. - @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ - @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ + @kcopy!(n, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ + @kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ if pᴴq ≠ zero(FC) - @. vₖ = q / βₖ₊₁ # βₖ₊₁vₖ₊₁ = q - @. uₖ = p / conj(γₖ₊₁) # γ̄ₖ₊₁uₖ₊₁ = p + vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q + uₖ .= p ./ conj(γₖ₊₁) # γ̄ₖ₊₁uₖ₊₁ = p end # Compute τₖ₊₁ = τₖ + ‖vₖ₊₁‖² diff --git a/src/symmlq.jl b/src/symmlq.jl index 604698525..3e1a57a7a 100644 --- a/src/symmlq.jl +++ b/src/symmlq.jl @@ -303,9 +303,9 @@ kwargs_symmlq = (:M, :ldiv, :transfer_to_cg, :λ, :λest, :atol, :rtol, :etol, : mul!(Mv_next, A, v) α = @kdotr(m, v, Mv_next) + λ @kaxpy!(m, -oldβ, Mvold, Mv_next) - @. Mvold = Mv + @kcopy!(m, Mvold, Mv) @kaxpy!(m, -α, Mv, Mv_next) - @. Mv = Mv_next + @kcopy!(m, Mv, Mv_next) MisI || mulorldiv!(v, M, Mv, ldiv) β = @kdotr(m, v, Mv) β < 0 && error("Preconditioner is not positive definite") diff --git a/src/tricg.jl b/src/tricg.jl index 8250e6dfc..e516fcd3e 100644 --- a/src/tricg.jl +++ b/src/tricg.jl @@ -369,22 +369,22 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax if iter == 1 # [ 1 0 ] [ gx₁ gy₁ ] = [ v₁ 0 ] # [ δ̄₁ 1 ] [ gx₂ gy₂ ] [ 0 u₁ ] - @. gx₂ₖ₋₁ = vₖ - @. gx₂ₖ = - conj(δₖ) * gx₂ₖ₋₁ - @. gy₂ₖ = uₖ + @kcopy!(m, gx₂ₖ₋₁, vₖ) + gx₂ₖ .= -conj(δₖ) .* gx₂ₖ₋₁ + @kcopy!(n, gy₂ₖ, uₖ) else # [ 0 σ̄ₖ 1 0 ] [ gx₂ₖ₋₃ gy₂ₖ₋₃ ] = [ vₖ 0 ] # [ η̄ₖ λ̄ₖ δ̄ₖ 1 ] [ gx₂ₖ₋₂ gy₂ₖ₋₂ ] [ 0 uₖ ] # [ gx₂ₖ₋₁ gy₂ₖ₋₁ ] # [ gx₂ₖ gy₂ₖ ] - @. gx₂ₖ₋₁ = conj(ηₖ) * gx₂ₖ₋₁ + conj(λₖ) * gx₂ₖ - @. gy₂ₖ₋₁ = conj(ηₖ) * gy₂ₖ₋₁ + conj(λₖ) * gy₂ₖ + gx₂ₖ₋₁ .= conj(ηₖ) .* gx₂ₖ₋₁ .+ conj(λₖ) .* gx₂ₖ + gy₂ₖ₋₁ .= conj(ηₖ) .* gy₂ₖ₋₁ .+ conj(λₖ) .* gy₂ₖ - @. gx₂ₖ = vₖ - conj(σₖ) * gx₂ₖ - @. gy₂ₖ = - conj(σₖ) * gy₂ₖ + gx₂ₖ .= vₖ .- conj(σₖ) .* gx₂ₖ + gy₂ₖ .= .- conj(σₖ) .* gy₂ₖ - @. gx₂ₖ₋₁ = - gx₂ₖ₋₁ - conj(δₖ) * gx₂ₖ - @. gy₂ₖ₋₁ = uₖ - gy₂ₖ₋₁ - conj(δₖ) * gy₂ₖ + gx₂ₖ₋₁ .= .- gx₂ₖ₋₁ .- conj(δₖ) .* gx₂ₖ + gy₂ₖ₋₁ .= uₖ .- gy₂ₖ₋₁ .- conj(δₖ) .* gy₂ₖ # g₂ₖ₋₃ == g₂ₖ and g₂ₖ₋₂ == g₂ₖ₋₁ @kswap(gx₂ₖ₋₁, gx₂ₖ) diff --git a/src/trilqr.jl b/src/trilqr.jl index 2b584c216..6e1f22833 100644 --- a/src/trilqr.jl +++ b/src/trilqr.jl @@ -300,7 +300,7 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, # Compute d̅ₖ. if iter == 1 # d̅₁ = u₁ - @. d̅ = uₖ + @kcopy!(n, d̅, uₖ) else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅) @@ -351,14 +351,14 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, if iter == 2 wₖ₋₁ = wₖ₋₂ @kaxpy!(m, one(FC), vₖ₋₁, wₖ₋₁) - @. wₖ₋₁ = vₖ₋₁ / conj(δₖ₋₁) + wₖ₋₁ .= vₖ₋₁ ./ conj(δₖ₋₁) end # w₂ = (v₂ - λ̄₁w₁) / δ̄₂ if iter == 3 wₖ₋₁ = wₖ₋₃ @kaxpy!(m, one(FC), vₖ₋₁, wₖ₋₁) @kaxpy!(m, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁) - @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁) + wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁) end # wₖ₋₁ = (vₖ₋₁ - λ̄ₖ₋₂wₖ₋₂ - ϵ̄ₖ₋₃wₖ₋₃) / δ̄ₖ₋₁ if iter ≥ 4 @@ -366,7 +366,7 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, wₖ₋₁ = wₖ₋₃ @kaxpy!(m, one(FC), vₖ₋₁, wₖ₋₁) @kaxpy!(m, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁) - @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁) + wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁) end if iter ≥ 3 @@ -399,14 +399,14 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, end # Compute uₖ₊₁ and uₖ₊₁. - @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ - @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ + @kcopy!(m, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ + @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 ϵₖ₋₃, λₖ₋₂, δbarₖ₋₁, cₖ₋₁, sₖ₋₁, γₖ and βₖ. diff --git a/src/trimr.jl b/src/trimr.jl index ae61b785a..780265e4a 100644 --- a/src/trimr.jl +++ b/src/trimr.jl @@ -447,9 +447,9 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : if iter == 1 # [ δ₁ 0 ] [ gx₁ gy₁ ] = [ v₁ 0 ] # [ σ₁ δ₂ ] [ gx₂ gy₂ ] [ 0 u₁ ] - @. gx₂ₖ₋₁ = vₖ / δ₂ₖ₋₁ - @. gx₂ₖ = - σ₂ₖ₋₁ / δ₂ₖ * gx₂ₖ₋₁ - @. gy₂ₖ = uₖ / δ₂ₖ + gx₂ₖ₋₁ .= vₖ ./ δ₂ₖ₋₁ + gx₂ₖ .= -(σ₂ₖ₋₁ / δ₂ₖ) .* gx₂ₖ₋₁ + gy₂ₖ .= uₖ ./ δ₂ₖ elseif iter == 2 # [ η₁ σ₂ δ₃ 0 ] [ gx₁ gy₁ ] = [ v₂ 0 ] # [ λ₁ η₂ σ₃ δ₄ ] [ gx₂ gy₂ ] [ 0 u₂ ] @@ -458,23 +458,23 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : @kswap(gx₂ₖ₋₃, gx₂ₖ₋₁) @kswap(gx₂ₖ₋₂, gx₂ₖ) @kswap(gy₂ₖ₋₂, gy₂ₖ) - @. gx₂ₖ₋₁ = (vₖ - η₂ₖ₋₃ * gx₂ₖ₋₃ - σ₂ₖ₋₂ * gx₂ₖ₋₂ ) / δ₂ₖ₋₁ - @. gx₂ₖ = ( - λ₂ₖ₋₃ * gx₂ₖ₋₃ - η₂ₖ₋₂ * gx₂ₖ₋₂ - σ₂ₖ₋₁ * gx₂ₖ₋₁) / δ₂ₖ - @. gy₂ₖ₋₁ = ( - η₂ₖ₋₃ * gy₂ₖ₋₃ - σ₂ₖ₋₂ * gy₂ₖ₋₂ ) / δ₂ₖ₋₁ - @. gy₂ₖ = (uₖ - λ₂ₖ₋₃ * gy₂ₖ₋₃ - η₂ₖ₋₂ * gy₂ₖ₋₂ - σ₂ₖ₋₁ * gy₂ₖ₋₁) / δ₂ₖ + gx₂ₖ₋₁ .= (vₖ .- η₂ₖ₋₃ .* gx₂ₖ₋₃ .- σ₂ₖ₋₂ .* gx₂ₖ₋₂ ) ./ δ₂ₖ₋₁ + gx₂ₖ .= ( .- λ₂ₖ₋₃ .* gx₂ₖ₋₃ .- η₂ₖ₋₂ .* gx₂ₖ₋₂ .- σ₂ₖ₋₁ .* gx₂ₖ₋₁) ./ δ₂ₖ + gy₂ₖ₋₁ .= ( .- η₂ₖ₋₃ .* gy₂ₖ₋₃ .- σ₂ₖ₋₂ .* gy₂ₖ₋₂ ) ./ δ₂ₖ₋₁ + gy₂ₖ .= (uₖ .- λ₂ₖ₋₃ .* gy₂ₖ₋₃ .- η₂ₖ₋₂ .* gy₂ₖ₋₂ .- σ₂ₖ₋₁ .* gy₂ₖ₋₁) ./ δ₂ₖ else # μ₂ₖ₋₅ * gx₂ₖ₋₅ + λ₂ₖ₋₄ * gx₂ₖ₋₄ + η₂ₖ₋₃ * gx₂ₖ₋₃ + σ₂ₖ₋₂ * gx₂ₖ₋₂ + δ₂ₖ₋₁ * gx₂ₖ₋₁ = vₖ # μ₂ₖ₋₄ * gx₂ₖ₋₄ + λ₂ₖ₋₃ * gx₂ₖ₋₃ + η₂ₖ₋₂ * gx₂ₖ₋₂ + σ₂ₖ₋₁ * gx₂ₖ₋₁ + δ₂ₖ * gx₂ₖ = 0 g₂ₖ₋₁ = g₂ₖ₋₅ = gx₂ₖ₋₃; g₂ₖ = g₂ₖ₋₄ = gx₂ₖ₋₂; g₂ₖ₋₃ = gx₂ₖ₋₁; g₂ₖ₋₂ = gx₂ₖ - @. g₂ₖ₋₁ = (vₖ - μ₂ₖ₋₅ * g₂ₖ₋₅ - λ₂ₖ₋₄ * g₂ₖ₋₄ - η₂ₖ₋₃ * g₂ₖ₋₃ - σ₂ₖ₋₂ * g₂ₖ₋₂ ) / δ₂ₖ₋₁ - @. g₂ₖ = ( - μ₂ₖ₋₄ * g₂ₖ₋₄ - λ₂ₖ₋₃ * g₂ₖ₋₃ - η₂ₖ₋₂ * g₂ₖ₋₂ - σ₂ₖ₋₁ * g₂ₖ₋₁) / δ₂ₖ + g₂ₖ₋₁ .= (vₖ .- μ₂ₖ₋₅ .* g₂ₖ₋₅ .- λ₂ₖ₋₄ .* g₂ₖ₋₄ .- η₂ₖ₋₃ .* g₂ₖ₋₃ .- σ₂ₖ₋₂ .* g₂ₖ₋₂ ) ./ δ₂ₖ₋₁ + g₂ₖ .= ( .- μ₂ₖ₋₄ .* g₂ₖ₋₄ .- λ₂ₖ₋₃ .* g₂ₖ₋₃ .- η₂ₖ₋₂ .* g₂ₖ₋₂ .- σ₂ₖ₋₁ .* g₂ₖ₋₁) ./ δ₂ₖ @kswap(gx₂ₖ₋₃, gx₂ₖ₋₁) @kswap(gx₂ₖ₋₂, gx₂ₖ) # μ₂ₖ₋₅ * gy₂ₖ₋₅ + λ₂ₖ₋₄ * gy₂ₖ₋₄ + η₂ₖ₋₃ * gy₂ₖ₋₃ + σ₂ₖ₋₂ * gy₂ₖ₋₂ + δ₂ₖ₋₁ * gy₂ₖ₋₁ = 0 # μ₂ₖ₋₄ * gy₂ₖ₋₄ + λ₂ₖ₋₃ * gy₂ₖ₋₃ + η₂ₖ₋₂ * gy₂ₖ₋₂ + σ₂ₖ₋₁ * gy₂ₖ₋₁ + δ₂ₖ * gy₂ₖ = uₖ g₂ₖ₋₁ = g₂ₖ₋₅ = gy₂ₖ₋₃; g₂ₖ = g₂ₖ₋₄ = gy₂ₖ₋₂; g₂ₖ₋₃ = gy₂ₖ₋₁; g₂ₖ₋₂ = gy₂ₖ - @. g₂ₖ₋₁ = ( - μ₂ₖ₋₅ * g₂ₖ₋₅ - λ₂ₖ₋₄ * g₂ₖ₋₄ - η₂ₖ₋₃ * g₂ₖ₋₃ - σ₂ₖ₋₂ * g₂ₖ₋₂ ) / δ₂ₖ₋₁ - @. g₂ₖ = (uₖ - μ₂ₖ₋₄ * g₂ₖ₋₄ - λ₂ₖ₋₃ * g₂ₖ₋₃ - η₂ₖ₋₂ * g₂ₖ₋₂ - σ₂ₖ₋₁ * g₂ₖ₋₁) / δ₂ₖ + g₂ₖ₋₁ .= ( .- μ₂ₖ₋₅ .* g₂ₖ₋₅ .- λ₂ₖ₋₄ .* g₂ₖ₋₄ .- η₂ₖ₋₃ .* g₂ₖ₋₃ .- σ₂ₖ₋₂ .* g₂ₖ₋₂ ) ./ δ₂ₖ₋₁ + g₂ₖ .= (uₖ .- μ₂ₖ₋₄ .* g₂ₖ₋₄ .- λ₂ₖ₋₃ .* g₂ₖ₋₃ .- η₂ₖ₋₂ .* g₂ₖ₋₂ .- σ₂ₖ₋₁ .* g₂ₖ₋₁) ./ δ₂ₖ @kswap(gy₂ₖ₋₃, gy₂ₖ₋₁) @kswap(gy₂ₖ₋₂, gy₂ₖ) end diff --git a/src/usymlq.jl b/src/usymlq.jl index 20fe5809d..f15138e6b 100644 --- a/src/usymlq.jl +++ b/src/usymlq.jl @@ -296,21 +296,21 @@ kwargs_usymlq = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, # Compute d̅ₖ. if iter == 1 # d̅₁ = u₁ - @. d̅ = uₖ + @kcopy!(n, d̅, uₖ) else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅) end # Compute uₖ₊₁ and uₖ₊₁. - @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ - @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ + @kcopy!(m, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ + @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 # Compute USYMLQ residual norm diff --git a/src/usymqr.jl b/src/usymqr.jl index 8ba7203bc..adaff466c 100644 --- a/src/usymqr.jl +++ b/src/usymqr.jl @@ -279,14 +279,14 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, if iter == 1 wₖ = wₖ₋₁ @kaxpy!(n, one(FC), uₖ, wₖ) - @. wₖ = wₖ / δₖ + wₖ .= wₖ ./ δₖ end # w₂ = (u₂ - λ₁w₁) / δ₂ if iter == 2 wₖ = wₖ₋₂ @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) @kaxpy!(n, one(FC), uₖ, wₖ) - @. wₖ = wₖ / δₖ + wₖ .= wₖ ./ δₖ end # wₖ = (uₖ - λₖ₋₁wₖ₋₁ - ϵₖ₋₂wₖ₋₂) / δₖ if iter ≥ 3 @@ -294,7 +294,7 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, wₖ = wₖ₋₂ @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ) @kaxpy!(n, one(FC), uₖ, wₖ) - @. wₖ = wₖ / δₖ + wₖ .= wₖ ./ δₖ end # Compute solution xₖ. @@ -310,14 +310,14 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, history && push!(AᴴrNorms, AᴴrNorm) # Compute uₖ₊₁ and uₖ₊₁. - @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ - @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ + @kcopy!(m, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ + @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. From 9dc980f8420773f53ca28999b27628556dde0c61 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 20 May 2024 18:05:56 -0400 Subject: [PATCH 2/4] Use the correct order of arguments for kcopy! --- src/bilq.jl | 6 +++--- src/bilqr.jl | 6 +++--- src/cg_lanczos.jl | 4 ++-- src/minres.jl | 4 ++-- src/minres_qlp.jl | 2 +- src/qmr.jl | 4 ++-- src/symmlq.jl | 4 ++-- src/tricg.jl | 4 ++-- src/trilqr.jl | 6 +++--- src/usymlq.jl | 6 +++--- src/usymqr.jl | 4 ++-- 11 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/bilq.jl b/src/bilq.jl index 6153c59ab..5bda9802a 100644 --- a/src/bilq.jl +++ b/src/bilq.jl @@ -329,15 +329,15 @@ kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :time # Compute d̅ₖ. if iter == 1 # d̅₁ = v₁ - @kcopy!(n, d̅, vₖ) + @kcopy!(n, vₖ, d̅) # d̅ ← vₖ else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ @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 diff --git a/src/bilqr.jl b/src/bilqr.jl index 01bd7b8fb..1abc60bce 100644 --- a/src/bilqr.jl +++ b/src/bilqr.jl @@ -316,7 +316,7 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi # Compute d̅ₖ. if iter == 1 # d̅₁ = v₁ - @kcopy!(n, d̅, vₖ) + @kcopy!(n, vₖ, d̅) # d̅ ← vₖ else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ @kaxpby!(n, -cₖ, vₖ, conj(sₖ), d̅) @@ -421,8 +421,8 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi 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 ≠ zero(FC) vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl index 892a37527..778a73d6f 100644 --- a/src/cg_lanczos.jl +++ b/src/cg_lanczos.jl @@ -218,9 +218,9 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax @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ₖ + @kcopy!(n, Mv, Mv_prev) # Mvₖ₋₁ ← Mvₖ end - @kcopy!(n, Mv, Mv_next) # Mvₖ ← Mvₖ₊₁ + @kcopy!(n, Mv_next, Mv) # Mvₖ ← 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ₖ₊₁ / βₖ₊₁ diff --git a/src/minres.jl b/src/minres.jl index 30f248ec3..2e0dd114b 100644 --- a/src/minres.jl +++ b/src/minres.jl @@ -276,8 +276,8 @@ kwargs_minres = (:M, :ldiv, :λ, :atol, :rtol, :etol, :conlim, :itmax, :timemax, end @kaxpy!(n, one(FC) / β, v, w) - @kcopy!(n, r1, r2) - @kcopy!(n, r2, y) + @kcopy!(n, r2, r1) # r1 ← r2 + @kcopy!(n, y, r2) # r2 ← y MisI || mulorldiv!(v, M, r2, ldiv) oldβ = β β = @kdotr(n, r2, v) diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl index 48f189c99..9cb569f6a 100644 --- a/src/minres_qlp.jl +++ b/src/minres_qlp.jl @@ -369,7 +369,7 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve # Compute directions wₖ₋₂, ẘₖ₋₁ and w̄ₖ, last columns of Wₖ = Vₖ(Pₖ)ᴴ if iter == 1 # w̅₁ = v₁ - @kcopy!(n, wₖ, vₖ) + @kcopy!(n, vₖ, wₖ) elseif iter == 2 # [w̅ₖ₋₁ vₖ] [cpₖ spₖ] = [ẘₖ₋₁ w̅ₖ] ⟷ ẘₖ₋₁ = cpₖ * w̅ₖ₋₁ + spₖ * vₖ # [spₖ -cpₖ] ⟷ w̅ₖ = spₖ * w̅ₖ₋₁ - cpₖ * vₖ diff --git a/src/qmr.jl b/src/qmr.jl index 7ce36b82f..2a1ad578c 100644 --- a/src/qmr.jl +++ b/src/qmr.jl @@ -349,8 +349,8 @@ kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist @kaxpy!(n, ζₖ, wₖ, x) # 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 ≠ zero(FC) vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q diff --git a/src/symmlq.jl b/src/symmlq.jl index 3e1a57a7a..bfc866944 100644 --- a/src/symmlq.jl +++ b/src/symmlq.jl @@ -303,9 +303,9 @@ kwargs_symmlq = (:M, :ldiv, :transfer_to_cg, :λ, :λest, :atol, :rtol, :etol, : mul!(Mv_next, A, v) α = @kdotr(m, v, Mv_next) + λ @kaxpy!(m, -oldβ, Mvold, Mv_next) - @kcopy!(m, Mvold, Mv) + @kcopy!(m, Mv, Mvold) # Mvold ← Mv @kaxpy!(m, -α, Mv, Mv_next) - @kcopy!(m, Mv, Mv_next) + @kcopy!(m, Mv_next, Mv) # Mv ← Mv_next MisI || mulorldiv!(v, M, Mv, ldiv) β = @kdotr(m, v, Mv) β < 0 && error("Preconditioner is not positive definite") diff --git a/src/tricg.jl b/src/tricg.jl index e516fcd3e..99d714bdd 100644 --- a/src/tricg.jl +++ b/src/tricg.jl @@ -369,9 +369,9 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax if iter == 1 # [ 1 0 ] [ gx₁ gy₁ ] = [ v₁ 0 ] # [ δ̄₁ 1 ] [ gx₂ gy₂ ] [ 0 u₁ ] - @kcopy!(m, gx₂ₖ₋₁, vₖ) + @kcopy!(m, vₖ, gx₂ₖ₋₁) # gx₂ₖ₋₁ ← vₖ gx₂ₖ .= -conj(δₖ) .* gx₂ₖ₋₁ - @kcopy!(n, gy₂ₖ, uₖ) + @kcopy!(n, uₖ, gy₂ₖ) # gy₂ₖ ← uₖ else # [ 0 σ̄ₖ 1 0 ] [ gx₂ₖ₋₃ gy₂ₖ₋₃ ] = [ vₖ 0 ] # [ η̄ₖ λ̄ₖ δ̄ₖ 1 ] [ gx₂ₖ₋₂ gy₂ₖ₋₂ ] [ 0 uₖ ] diff --git a/src/trilqr.jl b/src/trilqr.jl index 6e1f22833..042db2604 100644 --- a/src/trilqr.jl +++ b/src/trilqr.jl @@ -300,7 +300,7 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, # Compute d̅ₖ. if iter == 1 # d̅₁ = u₁ - @kcopy!(n, d̅, uₖ) + @kcopy!(n, uₖ, d̅) # d̅ ← uₖ else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅) @@ -399,8 +399,8 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, end # Compute uₖ₊₁ and uₖ₊₁. - @kcopy!(m, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ - @kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ + @kcopy!(m, vₖ, vₖ₋₁) # vₖ₋₁ ← vₖ + @kcopy!(n, uₖ, uₖ₋₁) # uₖ₋₁ ← uₖ if βₖ₊₁ ≠ zero(T) vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q diff --git a/src/usymlq.jl b/src/usymlq.jl index f15138e6b..73d6c7114 100644 --- a/src/usymlq.jl +++ b/src/usymlq.jl @@ -296,15 +296,15 @@ kwargs_usymlq = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose, # Compute d̅ₖ. if iter == 1 # d̅₁ = u₁ - @kcopy!(n, d̅, uₖ) + @kcopy!(n, uₖ, d̅) # d̅ ← vₖ else # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅) end # Compute uₖ₊₁ and uₖ₊₁. - @kcopy!(m, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ - @kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ + @kcopy!(m, vₖ, vₖ₋₁) # vₖ₋₁ ← vₖ + @kcopy!(n, uₖ, uₖ₋₁) # uₖ₋₁ ← uₖ if βₖ₊₁ ≠ zero(T) vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q diff --git a/src/usymqr.jl b/src/usymqr.jl index adaff466c..1d7488e36 100644 --- a/src/usymqr.jl +++ b/src/usymqr.jl @@ -310,8 +310,8 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, history && push!(AᴴrNorms, AᴴrNorm) # Compute uₖ₊₁ and uₖ₊₁. - @kcopy!(m, vₖ₋₁, vₖ) # vₖ₋₁ ← vₖ - @kcopy!(n, uₖ₋₁, uₖ) # uₖ₋₁ ← uₖ + @kcopy!(m, vₖ, vₖ₋₁) # vₖ₋₁ ← vₖ + @kcopy!(n, uₖ, uₖ₋₁) # uₖ₋₁ ← uₖ if βₖ₊₁ ≠ zero(T) vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q From ffea212b78ceb68ab77a4bfe6f91b34fa2672eeb Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 20 May 2024 18:18:18 -0400 Subject: [PATCH 3/4] Fix a typo in cg_lanczos_shift.jl --- src/cg_lanczos_shift.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl index b38accf2c..909a32c78 100644 --- a/src/cg_lanczos_shift.jl +++ b/src/cg_lanczos_shift.jl @@ -211,9 +211,9 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t @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ₖ + @kcopy!(n, Mv_prev, Mv) # Mvₖ₋₁ ← Mvₖ end - kcopy!(n, Mv, Mv_next) # Mvₖ ← Mvₖ₊₁ + @kcopy!(n, Mv, Mv_next) # Mvₖ ← 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ₖ₊₁ / βₖ₊₁ From 865ee18ffc26294fd66c322263149a37cfa0277c Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Mon, 20 May 2024 18:27:34 -0400 Subject: [PATCH 4/4] fix a typo --- src/cg_lanczos_shift.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl index 909a32c78..e47c882e6 100644 --- a/src/cg_lanczos_shift.jl +++ b/src/cg_lanczos_shift.jl @@ -211,9 +211,9 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t @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ₖ + @kcopy!(n, Mv, Mv_prev) # Mvₖ₋₁ ← Mvₖ end - @kcopy!(n, Mv, Mv_next) # Mvₖ ← Mvₖ₊₁ + @kcopy!(n, Mv_next, Mv) # Mvₖ ← 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ₖ₊₁ / βₖ₊₁