From 9f50cff8ace3f3a9afd378b1705ef67f0274af73 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Wed, 16 Oct 2024 13:55:21 -0500 Subject: [PATCH 1/5] Use the macro kcopy! instead of broadcast --- src/bicgstab.jl | 4 +- src/car.jl | 12 ++-- src/cg.jl | 6 +- src/cg_lanczos.jl | 6 +- src/cg_lanczos_shift.jl | 8 +-- src/cgls.jl | 4 +- src/cgls_lanczos_shift.jl | 18 ++--- src/cgne.jl | 4 +- src/cgs.jl | 8 +-- src/cr.jl | 8 +-- src/craig.jl | 2 +- src/craigmr.jl | 14 ++-- src/crls.jl | 6 +- src/crmr.jl | 4 +- src/diom.jl | 2 +- src/dqgmres.jl | 2 +- src/fgmres.jl | 2 +- src/fom.jl | 4 +- src/gmres.jl | 4 +- src/gpmr.jl | 4 +- src/krylov_solve.jl | 68 +++++++++--------- src/krylov_solvers.jl | 140 +++++++++++++++++++------------------- src/lnlq.jl | 8 +-- src/lslq.jl | 6 +- src/lsmr.jl | 6 +- src/lsqr.jl | 6 +- src/minares.jl | 2 +- src/minres.jl | 4 +- src/minres_qlp.jl | 8 +-- src/symmlq.jl | 4 +- src/tricg.jl | 12 ++-- src/trimr.jl | 16 ++--- 32 files changed, 201 insertions(+), 201 deletions(-) diff --git a/src/bicgstab.jl b/src/bicgstab.jl index b3fa1d370..54cc02ea9 100644 --- a/src/bicgstab.jl +++ b/src/bicgstab.jl @@ -151,14 +151,14 @@ kwargs_bicgstab = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, mul!(r₀, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), r₀) else - r₀ .= b + @kcopy!(n, r₀, b) # r₀ ← b end @kfill!(x, zero(FC)) # x₀ @kfill!(s, zero(FC)) # s₀ @kfill!(v, zero(FC)) # v₀ MisI || mulorldiv!(r, M, r₀, ldiv) # r₀ - p .= r # p₁ + @kcopy!(n, p, r) # p₁ α = one(FC) # α₀ ω = one(FC) # ω₀ diff --git a/src/car.jl b/src/car.jl index 030885947..464a6458c 100644 --- a/src/car.jl +++ b/src/car.jl @@ -127,29 +127,29 @@ kwargs_car = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :ca mul!(r, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), r) else - r .= b + @kcopy!(n, r, b) # r ← b end # p₀ = r₀ = M(b - Ax₀) if MisI - p .= r + @kcopy!(n, p, r) # p ← r else mulorldiv!(p, M, r, ldiv) - r .= p + @kcopy!(n, r, p) # r ← p end mul!(s, A, r) # s₀ = Ar₀ # q₀ = MAp₀ and s₀ = MAr₀ if MisI - q .= s + @kcopy!(n, q, s) # q ← s else mulorldiv!(q, M, s, ldiv) - s .= q + @kcopy!(n, s, q) # s ← q end mul!(t, A, s) # t₀ = As₀ - u .= t # u₀ = Aq₀ + @kcopy!(n, u, t) # u₀ = Aq₀ ρ = @kdotr(n, t, s) # ρ₀ = ⟨t₀ , s₀⟩ # Compute ‖r₀‖ diff --git a/src/cg.jl b/src/cg.jl index 8edb709fb..769bb6c3c 100644 --- a/src/cg.jl +++ b/src/cg.jl @@ -138,10 +138,10 @@ kwargs_cg = (:M, :ldiv, :radius, :linesearch, :atol, :rtol, :itmax, :timemax, :v mul!(r, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), r) else - r .= b + @kcopy!(n, r, b) # r ← b end MisI || mulorldiv!(z, M, r, ldiv) - p .= z + @kcopy!(n, p, z) # p ← z γ = @kdotr(n, r, z) rNorm = sqrt(γ) history && push!(rNorms, rNorm) @@ -182,7 +182,7 @@ kwargs_cg = (:M, :ldiv, :radius, :linesearch, :atol, :rtol, :itmax, :timemax, :v inconsistent = !linesearch end if linesearch - iter == 0 && (x .= b) + iter == 0 && @kcopy!(n, x, b) # x ← b solved = true end end diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl index c08164c7c..0582561b9 100644 --- a/src/cg_lanczos.jl +++ b/src/cg_lanczos.jl @@ -135,7 +135,7 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax mul!(Mv, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), Mv) else - Mv .= b + @kcopy!(n, Mv, b) # Mv ← b end MisI || mulorldiv!(v, M, Mv, ldiv) # v₁ = M⁻¹r₀ β = sqrt(@kdotr(n, v, Mv)) # β₁ = v₁ᴴ M v₁ @@ -152,13 +152,13 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax solver.warm_start = false return solver end - p .= v + @kcopy!(n, p, v) # p ← v # Initialize Lanczos process. # β₁Mv₁ = b @kscal!(n, one(FC) / β, v) # v₁ ← v₁ / β₁ MisI || @kscal!(n, one(FC) / β, Mv) # Mv₁ ← Mv₁ / β₁ - Mv_prev .= Mv + @kcopy!(n, Mv_prev, Mv) # Mv_prev ← Mv iter = 0 itmax == 0 && (itmax = 2 * n) diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl index d4e146287..97adc860f 100644 --- a/src/cg_lanczos_shift.jl +++ b/src/cg_lanczos_shift.jl @@ -131,10 +131,10 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t for i = 1 : nshifts @kfill!(x[i], zero(FC)) # x₀ end - 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₁ - rNorms .= β + @kfill!(rNorms, β) if history for i = 1 : nshifts push!(rNorms_history[i], rNorms[i]) @@ -155,14 +155,14 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t # Initialize each p to v. for i = 1 : nshifts - p[i] .= v + @kcopy!(n, p[i], v) # pᵢ ← v end # Initialize Lanczos process. # β₁Mv₁ = b @kscal!(n, one(FC) / β, v) # v₁ ← v₁ / β₁ MisI || @kscal!(n, one(FC) / β, Mv) # Mv₁ ← Mv₁ / β₁ - Mv_prev .= Mv + @kcopy!(n, Mv_prev, Mv) # Mv_prev ← Mv # Initialize some constants used in recursions below. ρ = one(T) diff --git a/src/cgls.jl b/src/cgls.jl index bfbd9d43d..551268317 100644 --- a/src/cgls.jl +++ b/src/cgls.jl @@ -147,7 +147,7 @@ kwargs_cgls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose Mq = MisI ? q : solver.Mr @kfill!(x, zero(FC)) - r .= b + @copy!(m, r, b) # r ← b bNorm = @knrm2(m, r) # Marginally faster than norm(b) if bNorm == 0 stats.niter = 0 @@ -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) - p .= s + @kcopy!(n, p, s) # p ← s γ = @kdotr(n, s, s) # γ = sᴴs iter = 0 itmax == 0 && (itmax = m + n) diff --git a/src/cgls_lanczos_shift.jl b/src/cgls_lanczos_shift.jl index 407d6fb2a..c03a0bc90 100644 --- a/src/cgls_lanczos_shift.jl +++ b/src/cgls_lanczos_shift.jl @@ -137,11 +137,11 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose @kfill!(x[i], zero(FC)) # x₀ end - 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₁ - rNorms .= β + mul!(v, Aᴴ, u) # v₁ ← Aᴴ * b + β = sqrt(@kdotr(n, v, v)) # β₁ = v₁ᵀ M v₁ + @kcopy!(nshifts, rNorms, β) # rNorms ← β if history for i = 1 : nshifts push!(rNorms_history[i], rNorms[i]) @@ -158,17 +158,17 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose # Initialize each p to v. for i = 1 : nshifts - p[i] .= v + @kcopy!(n, p[i], v) # pᵢ ← v end # Initialize Lanczos process. # β₁v₁ = b - @kscal!(n, one(FC) / β, v) # v₁ ← v₁ / β₁ + @kscal!(n, one(FC) / β, v) # v₁ ← v₁ / β₁ @kscal!(m, one(FC) / β, u) # Initialize some constants used in recursions below. ρ = one(T) - σ .= β + @kcopy!(nshifts, σ, β) # σ ← β @kfill!(δhat, zero(T)) @kfill!(ω, zero(T)) @kfill!(γ, one(T)) @@ -206,8 +206,8 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose β = sqrt(@kdotr(n, v, v)) # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁ @kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁ @kscal!(m, one(FC) / β, utilde) # uₖ₊₁ = uₖ₊₁ / βₖ₊₁ - u_prev .= u - u .= utilde + @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 307f53095..906f51c01 100644 --- a/src/cgne.jl +++ b/src/cgne.jl @@ -152,7 +152,7 @@ kwargs_cgne = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor z = NisI ? r : solver.z @kfill!(x, zero(FC)) - r .= b + @kcopy!(m, r, b) # r ← b NisI || mulorldiv!(z, N, r, ldiv) rNorm = @knrm2(m, r) # Marginally faster than norm(r) history && push!(rNorms, rNorm) @@ -163,7 +163,7 @@ kwargs_cgne = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor stats.status = "x = 0 is a zero-residual solution" return solver end - λ > 0 && (s .= r) + λ > 0 && @kcopy!(m, s, r) # s ← r mul!(p, Aᴴ, z) # Use ‖p‖ to detect inconsistent system. diff --git a/src/cgs.jl b/src/cgs.jl index 38216fd88..692aaa010 100644 --- a/src/cgs.jl +++ b/src/cgs.jl @@ -153,7 +153,7 @@ kwargs_cgs = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist mul!(r₀, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), r₀) else - r₀ .= b + @kcopy!(n, r₀, b) # r₀ ← b end @kfill!(x, zero(FC)) # x₀ @@ -189,9 +189,9 @@ kwargs_cgs = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist (verbose > 0) && @printf(iostream, "%5s %7s %5s\n", "k", "‖rₖ‖", "timer") kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %.2fs\n", iter, rNorm, ktimer(start_time)) - u .= r # u₀ - p .= r # p₀ - @kfill!(q, zero(FC)) # q₋₁ + @kcopy!(n, u, r) # u₀ + @kcopy!(n, p, r) # p₀ + @kfill!(q, zero(FC)) # q₋₁ # Stopping criterion. solved = rNorm ≤ ε diff --git a/src/cr.jl b/src/cr.jl index d8981049b..14dc46f44 100644 --- a/src/cr.jl +++ b/src/cr.jl @@ -146,9 +146,9 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema mul!(p, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), p) else - p .= b + @kcopy!(n, p, b) # p ← b end - MisI && (r .= p) + MisI && @kcopy!(n, r, p) # r ← p MisI || mulorldiv!(r, M, p, ldiv) mul!(Ar, A, r) ρ = @kdotr(n, r, Ar) @@ -165,8 +165,8 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema solver.warm_start = false return solver end - p .= r - q .= Ar + @kcopy!(n, p, r) # p ← r + @kcopy!(n, q, Ar) # q ← Ar (verbose > 0) && (m = zero(T)) # quadratic model iter = 0 diff --git a/src/craig.jl b/src/craig.jl index 305771cd0..1311ad35e 100644 --- a/src/craig.jl +++ b/src/craig.jl @@ -202,7 +202,7 @@ kwargs_craig = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :btol, :conlim, :at @kfill!(x, zero(FC)) @kfill!(y, zero(FC)) - Mu .= b + @kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) β₁ = sqrt(@kdotr(m, u, Mu)) rNorm = β₁ diff --git a/src/craigmr.jl b/src/craigmr.jl index ff9f3e322..341ec4420 100644 --- a/src/craigmr.jl +++ b/src/craigmr.jl @@ -189,7 +189,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver # Compute y such that AAᴴy = b. Then recover x = Aᴴy. @kfill!(x, zero(FC)) @kfill!(y, zero(FC)) - Mu .= b + @kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) β = sqrt(@kdotr(m, u, Mu)) if β == 0 @@ -208,7 +208,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver MisI || @kscal!(m, one(FC)/β, Mu) # α₁Nv₁ = Aᴴu₁. mul!(Aᴴu, Aᴴ, u) - Nv .= Aᴴu + @kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) α = sqrt(@kdotr(n, v, Nv)) Anorm² = α * α @@ -233,10 +233,10 @@ 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 λₖ₊₁ - λ > 0 && (q .= v) # Additional vector needed to update x, by definition q₀ = 0 + λₖ = λ # λ₁ = λ + 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 (cpₖ, spₖ, αhat) = sym_givens(α, λₖ) @@ -257,7 +257,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver ɛ_c = atol + rtol * rNorm # Stopping tolerance for consistent systems. ɛ_i = atol + rtol * ArNorm # Stopping tolerance for inconsistent systems. - wbar .= u + @kcopy!(m, wbar, u) @kscal!(m, one(FC)/αhat, wbar) @kfill!(w, zero(FC)) @kfill!(d, zero(FC)) diff --git a/src/crls.jl b/src/crls.jl index 7dffa152e..2e4f4f177 100644 --- a/src/crls.jl +++ b/src/crls.jl @@ -140,7 +140,7 @@ kwargs_crls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose MAp = MisI ? Ap : solver.Ms @kfill!(x, zero(FC)) - r .= b + @kcopy!(m, r, b) # r ← b bNorm = @knrm2(m, r) # norm(b - A * x0) if x0 ≠ 0. rNorm = bNorm # + λ * ‖x0‖ if x0 ≠ 0 and λ > 0. history && push!(rNorms, rNorm) @@ -158,8 +158,8 @@ kwargs_crls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose mul!(s, A, Ar) MisI || mulorldiv!(Ms, M, s, ldiv) - p .= Ar - Ap .= s + @kcopy!(n, p, Ar) # p ← Ar + @kcopy!(m, Ap, s) # Ap ← s mul!(q, Aᴴ, Ms) # Ap λ > 0 && @kaxpy!(n, λ, p, q) # q = q + λ * p γ = @kdotr(m, s, Ms) # Faster than γ = dot(s, Ms) diff --git a/src/crmr.jl b/src/crmr.jl index b63608a33..204f595f5 100644 --- a/src/crmr.jl +++ b/src/crmr.jl @@ -163,9 +163,9 @@ kwargs_crmr = (:N, :ldiv, :λ, :atol, :rtol, :itmax, :timemax, :verbose, :histor history && push!(ArNorms, zero(T)) return solver end - λ > 0 && (s .= r) + λ > 0 && @kcopy!(m, s, r) # s ← r mul!(Aᴴr, Aᴴ, r) # - λ * x0 if x0 ≠ 0. - p .= Aᴴr + @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 diff --git a/src/diom.jl b/src/diom.jl index d28b076dc..e491f64e5 100644 --- a/src/diom.jl +++ b/src/diom.jl @@ -146,7 +146,7 @@ kwargs_diom = (:M, :N, :ldiv, :reorthogonalization, :atol, :rtol, :itmax, :timem mul!(t, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), t) else - t .= b + @kcopy!(n, t, b) # t ← b end MisI || mulorldiv!(r₀, M, t, ldiv) # M(b - Ax₀) rNorm = @knrm2(n, r₀) # β = ‖r₀‖₂ diff --git a/src/dqgmres.jl b/src/dqgmres.jl index 99ff7901f..53fa3ac76 100644 --- a/src/dqgmres.jl +++ b/src/dqgmres.jl @@ -146,7 +146,7 @@ kwargs_dqgmres = (:M, :N, :ldiv, :reorthogonalization, :atol, :rtol, :itmax, :ti mul!(t, A, Δx) @kaxpby!(n, one(FC), b, -one(FC), t) else - t .= b + @kcopy!(n, t, b) # t ← b end MisI || mulorldiv!(r₀, M, t, ldiv) # M(b - Ax₀) rNorm = @knrm2(n, r₀) # β = ‖r₀‖₂ diff --git a/src/fgmres.jl b/src/fgmres.jl index 8e3cea789..4482897a9 100644 --- a/src/fgmres.jl +++ b/src/fgmres.jl @@ -152,7 +152,7 @@ kwargs_fgmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :i @kaxpby!(n, one(FC), b, -one(FC), w) restart && @kaxpy!(n, one(FC), Δx, x) else - w .= b + @kcopy!(n, w, b) # w ← b end MisI || mulorldiv!(r₀, M, w, ldiv) # r₀ = M(b - Ax₀) β = @knrm2(n, r₀) # β = ‖r₀‖₂ diff --git a/src/fom.jl b/src/fom.jl index 3d9773b6c..ee6874a09 100644 --- a/src/fom.jl +++ b/src/fom.jl @@ -147,7 +147,7 @@ kwargs_fom = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itma @kaxpby!(n, one(FC), b, -one(FC), w) restart && @kaxpy!(n, one(FC), Δx, x) else - w .= b + @kcopy!(n, w, b) # w ← b end MisI || mulorldiv!(r₀, M, w, ldiv) # r₀ = M(b - Ax₀) β = @knrm2(n, r₀) # β = ‖r₀‖₂ @@ -314,7 +314,7 @@ kwargs_fom = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itma @kaxpy!(n, y[i], V[i], xr) end if !NisI - solver.p .= xr + @kcopy!(n, solver.p, xr) # p ← xr mulorldiv!(xr, N, solver.p, ldiv) end restart && @kaxpy!(n, one(FC), xr, x) diff --git a/src/gmres.jl b/src/gmres.jl index 8c7dd7e7a..4c8d14d45 100644 --- a/src/gmres.jl +++ b/src/gmres.jl @@ -147,7 +147,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it @kaxpby!(n, one(FC), b, -one(FC), w) restart && @kaxpy!(n, one(FC), Δx, x) else - w .= b + @kcopy!(n, w, b) # w ← b end MisI || mulorldiv!(r₀, M, w, ldiv) # r₀ = M(b - Ax₀) β = @knrm2(n, r₀) # β = ‖r₀‖₂ @@ -331,7 +331,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it @kaxpy!(n, y[i], V[i], xr) end if !NisI - solver.p .= xr + @kcopy!(n, solver.p, xr) # p ← xr mulorldiv!(xr, N, solver.p, ldiv) end restart && @kaxpy!(n, one(FC), xr, x) diff --git a/src/gpmr.jl b/src/gpmr.jl index 0c3f9166f..f0db40825 100644 --- a/src/gpmr.jl +++ b/src/gpmr.jl @@ -500,11 +500,11 @@ kwargs_gpmr = (:C, :D, :E, :F, :ldiv, :gsp, :λ, :μ, :reorthogonalization, :ato @kaxpy!(n, zt[2i] , U[i], y) # xₖ = ζ₂u₁ + ζ₄u₂ + ••• + ζ₂ₖuₖ end if !EisI - wB .= x + @kcopy!(m, wB, x) # wB ← x mulorldiv!(x, E, wB, ldiv) end if !FisI - wA .= y + @kcopy!(n, wA, y) # wA ← y mulorldiv!(y, F, wA, ldiv) end warm_start && @kaxpy!(m, one(FC), Δx, x) diff --git a/src/krylov_solve.jl b/src/krylov_solve.jl index bfc14715c..8a96bef19 100644 --- a/src/krylov_solve.jl +++ b/src/krylov_solve.jl @@ -7,40 +7,40 @@ function solve! end # Krylov methods for (KS, fun, fun2, args, def_args, optargs, def_optargs, kwargs, def_kwargs) in [ - (:LsmrSolver , :lsmr! , :lsmr , args_lsmr , def_args_lsmr , () , () , kwargs_lsmr , def_kwargs_lsmr ) - (:CgsSolver , :cgs! , :cgs , args_cgs , def_args_cgs , optargs_cgs , def_optargs_cgs , kwargs_cgs , def_kwargs_cgs ) - (:UsymlqSolver , :usymlq! , :usymlq , args_usymlq , def_args_usymlq , optargs_usymlq , def_optargs_usymlq , kwargs_usymlq , def_kwargs_usymlq ) - (:LnlqSolver , :lnlq! , :lnlq , args_lnlq , def_args_lnlq , () , () , kwargs_lnlq , def_kwargs_lnlq ) - (:BicgstabSolver , :bicgstab! , :bicgstab , args_bicgstab , def_args_bicgstab , optargs_bicgstab , def_optargs_bicgstab , kwargs_bicgstab , def_kwargs_bicgstab ) - (:CrlsSolver , :crls! , :crls , args_crls , def_args_crls , () , () , kwargs_crls , def_kwargs_crls ) - (:LsqrSolver , :lsqr! , :lsqr , args_lsqr , def_args_lsqr , () , () , kwargs_lsqr , def_kwargs_lsqr ) - (:MinresSolver , :minres! , :minres , args_minres , def_args_minres , optargs_minres , def_optargs_minres , kwargs_minres , def_kwargs_minres ) - (:MinaresSolver , :minares! , :minares , args_minares , def_args_minares , optargs_minares , def_optargs_minares , kwargs_minares , def_kwargs_minares ) - (:CgneSolver , :cgne! , :cgne , args_cgne , def_args_cgne , () , () , kwargs_cgne , def_kwargs_cgne ) - (:DqgmresSolver , :dqgmres! , :dqgmres , args_dqgmres , def_args_dqgmres , optargs_dqgmres , def_optargs_dqgmres , kwargs_dqgmres , def_kwargs_dqgmres ) - (:SymmlqSolver , :symmlq! , :symmlq , args_symmlq , def_args_symmlq , optargs_symmlq , def_optargs_symmlq , kwargs_symmlq , def_kwargs_symmlq ) - (:TrimrSolver , :trimr! , :trimr , args_trimr , def_args_trimr , optargs_trimr , def_optargs_trimr , kwargs_trimr , def_kwargs_trimr ) - (:UsymqrSolver , :usymqr! , :usymqr , args_usymqr , def_args_usymqr , optargs_usymqr , def_optargs_usymqr , kwargs_usymqr , def_kwargs_usymqr ) - (:BilqrSolver , :bilqr! , :bilqr , args_bilqr , def_args_bilqr , optargs_bilqr , def_optargs_bilqr , kwargs_bilqr , def_kwargs_bilqr ) - (:CrSolver , :cr! , :cr , args_cr , def_args_cr , optargs_cr , def_optargs_cr , kwargs_cr , def_kwargs_cr ) - (:CarSolver , :car! , :car , args_car , def_args_car , optargs_car , def_optargs_car , kwargs_car , def_kwargs_car ) - (:CraigmrSolver , :craigmr! , :craigmr , args_craigmr , def_args_craigmr , () , () , kwargs_craigmr , def_kwargs_craigmr ) - (:TricgSolver , :tricg! , :tricg , args_tricg , def_args_tricg , optargs_tricg , def_optargs_tricg , kwargs_tricg , def_kwargs_tricg ) - (:CraigSolver , :craig! , :craig , args_craig , def_args_craig , () , () , kwargs_craig , def_kwargs_craig ) - (:DiomSolver , :diom! , :diom , args_diom , def_args_diom , optargs_diom , def_optargs_diom , kwargs_diom , def_kwargs_diom ) - (:LslqSolver , :lslq! , :lslq , args_lslq , def_args_lslq , () , () , kwargs_lslq , def_kwargs_lslq ) - (:TrilqrSolver , :trilqr! , :trilqr , args_trilqr , def_args_trilqr , optargs_trilqr , def_optargs_trilqr , kwargs_trilqr , def_kwargs_trilqr ) - (:CrmrSolver , :crmr! , :crmr , args_crmr , def_args_crmr , () , () , kwargs_crmr , def_kwargs_crmr ) - (:CgSolver , :cg! , :cg , args_cg , def_args_cg , optargs_cg , def_optargs_cg , kwargs_cg , def_kwargs_cg ) - (:CglsSolver , :cgls! , :cgls , args_cgls , def_args_cgls , () , () , kwargs_cgls , def_kwargs_cgls ) - (:CgLanczosSolver, :cg_lanczos! , :cg_lanczos , args_cg_lanczos , def_args_cg_lanczos , optargs_cg_lanczos, def_optargs_cg_lanczos, kwargs_cg_lanczos , def_kwargs_cg_lanczos ) - (:BilqSolver , :bilq! , :bilq , args_bilq , def_args_bilq , optargs_bilq , def_optargs_bilq , kwargs_bilq , def_kwargs_bilq ) - (:MinresQlpSolver, :minres_qlp! , :minres_qlp , args_minres_qlp , def_args_minres_qlp , optargs_minres_qlp, def_optargs_minres_qlp, kwargs_minres_qlp , def_kwargs_minres_qlp ) - (:QmrSolver , :qmr! , :qmr , args_qmr , def_args_qmr , optargs_qmr , def_optargs_qmr , kwargs_qmr , def_kwargs_qmr ) - (:GmresSolver , :gmres! , :gmres , args_gmres , def_args_gmres , optargs_gmres , def_optargs_gmres , kwargs_gmres , def_kwargs_gmres ) - (:FgmresSolver , :fgmres! , :fgmres , args_fgmres , def_args_fgmres , optargs_fgmres , def_optargs_fgmres , kwargs_fgmres , def_kwargs_fgmres ) - (:FomSolver , :fom! , :fom , args_fom , def_args_fom , optargs_fom , def_optargs_fom , kwargs_fom , def_kwargs_fom ) - (:GpmrSolver , :gpmr! , :gpmr , args_gpmr , def_args_gpmr , optargs_gpmr , def_optargs_gpmr , kwargs_gpmr , def_kwargs_gpmr ) + (:LsmrSolver , :lsmr! , :lsmr , args_lsmr , def_args_lsmr , () , () , kwargs_lsmr , def_kwargs_lsmr ) + (:CgsSolver , :cgs! , :cgs , args_cgs , def_args_cgs , optargs_cgs , def_optargs_cgs , kwargs_cgs , def_kwargs_cgs ) + (:UsymlqSolver , :usymlq! , :usymlq , args_usymlq , def_args_usymlq , optargs_usymlq , def_optargs_usymlq , kwargs_usymlq , def_kwargs_usymlq ) + (:LnlqSolver , :lnlq! , :lnlq , args_lnlq , def_args_lnlq , () , () , kwargs_lnlq , def_kwargs_lnlq ) + (:BicgstabSolver , :bicgstab! , :bicgstab , args_bicgstab , def_args_bicgstab , optargs_bicgstab , def_optargs_bicgstab , kwargs_bicgstab , def_kwargs_bicgstab ) + (:CrlsSolver , :crls! , :crls , args_crls , def_args_crls , () , () , kwargs_crls , def_kwargs_crls ) + (:LsqrSolver , :lsqr! , :lsqr , args_lsqr , def_args_lsqr , () , () , kwargs_lsqr , def_kwargs_lsqr ) + (:MinresSolver , :minres! , :minres , args_minres , def_args_minres , optargs_minres , def_optargs_minres , kwargs_minres , def_kwargs_minres ) + (:MinaresSolver , :minares! , :minares , args_minares , def_args_minares , optargs_minares , def_optargs_minares , kwargs_minares , def_kwargs_minares ) + (:CgneSolver , :cgne! , :cgne , args_cgne , def_args_cgne , () , () , kwargs_cgne , def_kwargs_cgne ) + (:DqgmresSolver , :dqgmres! , :dqgmres , args_dqgmres , def_args_dqgmres , optargs_dqgmres , def_optargs_dqgmres , kwargs_dqgmres , def_kwargs_dqgmres ) + (:SymmlqSolver , :symmlq! , :symmlq , args_symmlq , def_args_symmlq , optargs_symmlq , def_optargs_symmlq , kwargs_symmlq , def_kwargs_symmlq ) + (:TrimrSolver , :trimr! , :trimr , args_trimr , def_args_trimr , optargs_trimr , def_optargs_trimr , kwargs_trimr , def_kwargs_trimr ) + (:UsymqrSolver , :usymqr! , :usymqr , args_usymqr , def_args_usymqr , optargs_usymqr , def_optargs_usymqr , kwargs_usymqr , def_kwargs_usymqr ) + (:BilqrSolver , :bilqr! , :bilqr , args_bilqr , def_args_bilqr , optargs_bilqr , def_optargs_bilqr , kwargs_bilqr , def_kwargs_bilqr ) + (:CrSolver , :cr! , :cr , args_cr , def_args_cr , optargs_cr , def_optargs_cr , kwargs_cr , def_kwargs_cr ) + (:CarSolver , :car! , :car , args_car , def_args_car , optargs_car , def_optargs_car , kwargs_car , def_kwargs_car ) + (:CraigmrSolver , :craigmr! , :craigmr , args_craigmr , def_args_craigmr , () , () , kwargs_craigmr , def_kwargs_craigmr ) + (:TricgSolver , :tricg! , :tricg , args_tricg , def_args_tricg , optargs_tricg , def_optargs_tricg , kwargs_tricg , def_kwargs_tricg ) + (:CraigSolver , :craig! , :craig , args_craig , def_args_craig , () , () , kwargs_craig , def_kwargs_craig ) + (:DiomSolver , :diom! , :diom , args_diom , def_args_diom , optargs_diom , def_optargs_diom , kwargs_diom , def_kwargs_diom ) + (:LslqSolver , :lslq! , :lslq , args_lslq , def_args_lslq , () , () , kwargs_lslq , def_kwargs_lslq ) + (:TrilqrSolver , :trilqr! , :trilqr , args_trilqr , def_args_trilqr , optargs_trilqr , def_optargs_trilqr , kwargs_trilqr , def_kwargs_trilqr ) + (:CrmrSolver , :crmr! , :crmr , args_crmr , def_args_crmr , () , () , kwargs_crmr , def_kwargs_crmr ) + (:CgSolver , :cg! , :cg , args_cg , def_args_cg , optargs_cg , def_optargs_cg , kwargs_cg , def_kwargs_cg ) + (:CglsSolver , :cgls! , :cgls , args_cgls , def_args_cgls , () , () , kwargs_cgls , def_kwargs_cgls ) + (:CgLanczosSolver, :cg_lanczos! , :cg_lanczos , args_cg_lanczos , def_args_cg_lanczos , optargs_cg_lanczos, def_optargs_cg_lanczos, kwargs_cg_lanczos , def_kwargs_cg_lanczos) + (:BilqSolver , :bilq! , :bilq , args_bilq , def_args_bilq , optargs_bilq , def_optargs_bilq , kwargs_bilq , def_kwargs_bilq ) + (:MinresQlpSolver, :minres_qlp! , :minres_qlp , args_minres_qlp , def_args_minres_qlp , optargs_minres_qlp, def_optargs_minres_qlp, kwargs_minres_qlp , def_kwargs_minres_qlp) + (:QmrSolver , :qmr! , :qmr , args_qmr , def_args_qmr , optargs_qmr , def_optargs_qmr , kwargs_qmr , def_kwargs_qmr ) + (:GmresSolver , :gmres! , :gmres , args_gmres , def_args_gmres , optargs_gmres , def_optargs_gmres , kwargs_gmres , def_kwargs_gmres ) + (:FgmresSolver , :fgmres! , :fgmres , args_fgmres , def_args_fgmres , optargs_fgmres , def_optargs_fgmres , kwargs_fgmres , def_kwargs_fgmres ) + (:FomSolver , :fom! , :fom , args_fom , def_args_fom , optargs_fom , def_optargs_fom , kwargs_fom , def_kwargs_fom ) + (:GpmrSolver , :gpmr! , :gpmr , args_gpmr , def_args_gpmr , optargs_gpmr , def_optargs_gpmr , kwargs_gpmr , def_kwargs_gpmr ) (:CgLanczosShiftSolver , :cg_lanczos_shift! , :cg_lanczos_shift , args_cg_lanczos_shift , def_args_cg_lanczos_shift , (), (), kwargs_cg_lanczos_shift , def_kwargs_cg_lanczos_shift ) (:CglsLanczosShiftSolver, :cgls_lanczos_shift!, :cgls_lanczos_shift, args_cgls_lanczos_shift, def_args_cgls_lanczos_shift, (), (), kwargs_cgls_lanczos_shift, def_kwargs_cgls_lanczos_shift) ] diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl index 9a64253bd..51f3eaeb7 100644 --- a/src/krylov_solvers.jl +++ b/src/krylov_solvers.jl @@ -11,42 +11,42 @@ niterations, Aprod, Atprod, Bprod, warm_start! import Base.size, Base.sizeof, Base.format_bytes const KRYLOV_SOLVERS = Dict( - :cg => :CgSolver , - :cr => :CrSolver , - :car => :CarSolver , - :symmlq => :SymmlqSolver , - :cg_lanczos => :CgLanczosSolver , - :cg_lanczos_shift => :CgLanczosShiftSolver, - :minares => :MinaresSolver , - :minres => :MinresSolver , - :minres_qlp => :MinresQlpSolver , - :diom => :DiomSolver , - :fom => :FomSolver , - :dqgmres => :DqgmresSolver , - :gmres => :GmresSolver , - :fgmres => :FgmresSolver , - :gpmr => :GpmrSolver , - :usymlq => :UsymlqSolver , - :usymqr => :UsymqrSolver , - :tricg => :TricgSolver , - :trimr => :TrimrSolver , - :trilqr => :TrilqrSolver , - :cgs => :CgsSolver , - :bicgstab => :BicgstabSolver , - :bilq => :BilqSolver , - :qmr => :QmrSolver , - :bilqr => :BilqrSolver , - :cgls => :CglsSolver , + :cg => :CgSolver , + :cr => :CrSolver , + :car => :CarSolver , + :symmlq => :SymmlqSolver , + :cg_lanczos => :CgLanczosSolver, + :minares => :MinaresSolver , + :minres => :MinresSolver , + :minres_qlp => :MinresQlpSolver, + :diom => :DiomSolver , + :fom => :FomSolver , + :dqgmres => :DqgmresSolver , + :gmres => :GmresSolver , + :fgmres => :FgmresSolver , + :gpmr => :GpmrSolver , + :usymlq => :UsymlqSolver , + :usymqr => :UsymqrSolver , + :tricg => :TricgSolver , + :trimr => :TrimrSolver , + :trilqr => :TrilqrSolver , + :cgs => :CgsSolver , + :bicgstab => :BicgstabSolver , + :bilq => :BilqSolver , + :qmr => :QmrSolver , + :bilqr => :BilqrSolver , + :cgls => :CglsSolver , + :crls => :CrlsSolver , + :cgne => :CgneSolver , + :crmr => :CrmrSolver , + :lslq => :LslqSolver , + :lsqr => :LsqrSolver , + :lsmr => :LsmrSolver , + :lnlq => :LnlqSolver , + :craig => :CraigSolver , + :craigmr => :CraigmrSolver , + :cg_lanczos_shift => :CgLanczosShiftSolver , :cgls_lanczos_shift => :CglsLanczosShiftSolver, - :crls => :CrlsSolver , - :cgne => :CgneSolver , - :crmr => :CrmrSolver , - :lslq => :LslqSolver , - :lsqr => :LsqrSolver , - :lsmr => :LsmrSolver , - :lnlq => :LnlqSolver , - :craig => :CraigSolver , - :craigmr => :CraigmrSolver , ) "Abstract type for using Krylov solvers in-place" @@ -1973,42 +1973,42 @@ For example, instead of `x, stats = cg(A, b)`, you can use: function results end for (KS, fun, nsol, nA, nAt, warm_start) in [ - (:CarSolver , :car! , 1, 1, 0, true ) - (:LsmrSolver , :lsmr! , 1, 1, 1, false) - (:CgsSolver , :cgs! , 1, 2, 0, true ) - (:UsymlqSolver , :usymlq! , 1, 1, 1, true ) - (:LnlqSolver , :lnlq! , 2, 1, 1, false) - (:BicgstabSolver , :bicgstab! , 1, 2, 0, true ) - (:CrlsSolver , :crls! , 1, 1, 1, false) - (:LsqrSolver , :lsqr! , 1, 1, 1, false) - (:MinresSolver , :minres! , 1, 1, 0, true ) - (:MinaresSolver , :minares! , 1, 1, 0, true ) - (:CgneSolver , :cgne! , 1, 1, 1, false) - (:DqgmresSolver , :dqgmres! , 1, 1, 0, true ) - (:SymmlqSolver , :symmlq! , 1, 1, 0, true ) - (:TrimrSolver , :trimr! , 2, 1, 1, true ) - (:UsymqrSolver , :usymqr! , 1, 1, 1, true ) - (:BilqrSolver , :bilqr! , 2, 1, 1, true ) - (:CrSolver , :cr! , 1, 1, 0, true ) - (:CraigmrSolver , :craigmr! , 2, 1, 1, false) - (:TricgSolver , :tricg! , 2, 1, 1, true ) - (:CraigSolver , :craig! , 2, 1, 1, false) - (:DiomSolver , :diom! , 1, 1, 0, true ) - (:LslqSolver , :lslq! , 1, 1, 1, false) - (:TrilqrSolver , :trilqr! , 2, 1, 1, true ) - (:CrmrSolver , :crmr! , 1, 1, 1, false) - (:CgSolver , :cg! , 1, 1, 0, true ) - (:CgLanczosShiftSolver, :cg_lanczos_shift!, 1, 1, 0, false) - (:CglsSolver , :cgls! , 1, 1, 1, false) + (:CarSolver , :car! , 1, 1, 0, true ) + (:LsmrSolver , :lsmr! , 1, 1, 1, false) + (:CgsSolver , :cgs! , 1, 2, 0, true ) + (:UsymlqSolver , :usymlq! , 1, 1, 1, true ) + (:LnlqSolver , :lnlq! , 2, 1, 1, false) + (:BicgstabSolver , :bicgstab! , 1, 2, 0, true ) + (:CrlsSolver , :crls! , 1, 1, 1, false) + (:LsqrSolver , :lsqr! , 1, 1, 1, false) + (:MinresSolver , :minres! , 1, 1, 0, true ) + (:MinaresSolver , :minares! , 1, 1, 0, true ) + (:CgneSolver , :cgne! , 1, 1, 1, false) + (:DqgmresSolver , :dqgmres! , 1, 1, 0, true ) + (:SymmlqSolver , :symmlq! , 1, 1, 0, true ) + (:TrimrSolver , :trimr! , 2, 1, 1, true ) + (:UsymqrSolver , :usymqr! , 1, 1, 1, true ) + (:BilqrSolver , :bilqr! , 2, 1, 1, true ) + (:CrSolver , :cr! , 1, 1, 0, true ) + (:CraigmrSolver , :craigmr! , 2, 1, 1, false) + (:TricgSolver , :tricg! , 2, 1, 1, true ) + (:CraigSolver , :craig! , 2, 1, 1, false) + (:DiomSolver , :diom! , 1, 1, 0, true ) + (:LslqSolver , :lslq! , 1, 1, 1, false) + (:TrilqrSolver , :trilqr! , 2, 1, 1, true ) + (:CrmrSolver , :crmr! , 1, 1, 1, false) + (:CgSolver , :cg! , 1, 1, 0, true ) + (:CglsSolver , :cgls! , 1, 1, 1, false) + (:CgLanczosSolver, :cg_lanczos!, 1, 1, 0, true ) + (:BilqSolver , :bilq! , 1, 1, 1, true ) + (:MinresQlpSolver, :minres_qlp!, 1, 1, 0, true ) + (:QmrSolver , :qmr! , 1, 1, 1, true ) + (:GmresSolver , :gmres! , 1, 1, 0, true ) + (:FgmresSolver , :fgmres! , 1, 1, 0, true ) + (:FomSolver , :fom! , 1, 1, 0, true ) + (:GpmrSolver , :gpmr! , 2, 1, 0, true ) + (:CgLanczosShiftSolver , :cg_lanczos_shift! , 1, 1, 0, false) (:CglsLanczosShiftSolver, :cgls_lanczos_shift!, 1, 1, 1, false) - (:CgLanczosSolver , :cg_lanczos! , 1, 1, 0, true ) - (:BilqSolver , :bilq! , 1, 1, 1, true ) - (:MinresQlpSolver , :minres_qlp! , 1, 1, 0, true ) - (:QmrSolver , :qmr! , 1, 1, 1, true ) - (:GmresSolver , :gmres! , 1, 1, 0, true ) - (:FgmresSolver , :fgmres! , 1, 1, 0, true ) - (:FomSolver , :fom! , 1, 1, 0, true ) - (:GpmrSolver , :gpmr! , 2, 1, 0, true ) ] @eval begin size(solver :: $KS) = solver.m, solver.n diff --git a/src/lnlq.jl b/src/lnlq.jl index b03b1e68a..1ca359b54 100644 --- a/src/lnlq.jl +++ b/src/lnlq.jl @@ -226,7 +226,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly # Initialize generalized Golub-Kahan bidiagonalization. # β₁Mu₁ = b. - Mu .= b + @kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) # u₁ = M⁻¹ * Mu₁ βₖ = sqrt(@kdotr(m, u, Mu)) # β₁ = ‖u₁‖_M if βₖ ≠ 0 @@ -236,7 +236,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly # α₁Nv₁ = Aᴴu₁. mul!(Aᴴu, Aᴴ, u) - Nv .= 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 if αₖ ≠ 0 @@ -244,7 +244,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly NisI || @kscal!(n, one(FC) / αₖ, Nv) end - 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̅ₖ @@ -254,7 +254,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly λₖ = λ # λ₁ = λ cpₖ = spₖ = one(T) # Givens sines and cosines used to zero out λₖ cdₖ = sdₖ = one(FC) # Givens sines and cosines used to define λₖ₊₁ - λ > 0 && (q .= v) # Additional vector needed to update x, by definition q₀ = 0 + λ > 0 && @kcopy!(n, q, v) # Additional vector needed to update x, by definition q₀ = 0 # Initialize the regularization. if λ > 0 diff --git a/src/lslq.jl b/src/lslq.jl index 01c98c9cc..76ba2b052 100644 --- a/src/lslq.jl +++ b/src/lslq.jl @@ -226,7 +226,7 @@ kwargs_lslq = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :σ, :etol, :utol, : # Initialize Golub-Kahan process. # β₁ M u₁ = b. - Mu .= b + @kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) β₁ = sqrt(@kdotr(m, u, Mu)) if β₁ == 0 @@ -244,7 +244,7 @@ kwargs_lslq = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :σ, :etol, :utol, : @kscal!(m, one(FC)/β₁, u) MisI || @kscal!(m, one(FC)/β₁, Mu) mul!(Aᴴu, Aᴴ, u) - Nv .= Aᴴu + @kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) α = sqrt(@kdotr(n, v, Nv)) # = α₁ @@ -275,7 +275,7 @@ kwargs_lslq = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :σ, :etol, :utol, : xcgNorm = zero(T) xcgNorm² = zero(T) - w̄ .= v # w̄₁ = v₁ + @kcopy!(n, w̄, v) # w̄₁ = v₁ err_lbnd = zero(T) window = length(err_vec) diff --git a/src/lsmr.jl b/src/lsmr.jl index 6b068c87a..36db8b5af 100644 --- a/src/lsmr.jl +++ b/src/lsmr.jl @@ -200,7 +200,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, # Initialize Golub-Kahan process. # β₁ M u₁ = b. - Mu .= b + @kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) β₁ = sqrt(@kdotr(m, u, Mu)) if β₁ == 0 @@ -217,7 +217,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, @kscal!(m, one(FC)/β₁, u) MisI || @kscal!(m, one(FC)/β₁, Mu) mul!(Aᴴu, Aᴴ, u) - Nv .= Aᴴu + @kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) α = sqrt(@kdotr(n, v, Nv)) @@ -274,7 +274,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, @kscal!(n, one(FC)/α, v) NisI || @kscal!(n, one(FC)/α, Nv) - h .= v + @kcopy!(n, h, v) # h ← v @kfill!(hbar, zero(FC)) status = "unknown" diff --git a/src/lsqr.jl b/src/lsqr.jl index c5cf13ac0..120174105 100644 --- a/src/lsqr.jl +++ b/src/lsqr.jl @@ -197,7 +197,7 @@ kwargs_lsqr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, # Initialize Golub-Kahan process. # β₁ M u₁ = b. - Mu .= b + @kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) β₁ = sqrt(@kdotr(m, u, Mu)) if β₁ == 0 @@ -214,7 +214,7 @@ kwargs_lsqr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, @kscal!(m, one(FC)/β₁, u) MisI || @kscal!(m, one(FC)/β₁, Mu) mul!(Aᴴu, Aᴴ, u) - Nv .= Aᴴu + @kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) Anorm² = @kdotr(n, v, Nv) Anorm = sqrt(Anorm²) @@ -255,7 +255,7 @@ kwargs_lsqr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, end @kscal!(n, one(FC)/α, v) NisI || @kscal!(n, one(FC)/α, Nv) - w .= v + @kcopy!(n, w, v) # w ← v # Initialize other constants. ϕbar = β₁ diff --git a/src/minares.jl b/src/minares.jl index c3f74ab25..1a185f46f 100644 --- a/src/minares.jl +++ b/src/minares.jl @@ -140,7 +140,7 @@ kwargs_minares = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :verbo (λ ≠ 0) && @kaxpy!(n, λ, Δx, vₖ) @kaxpby!(n, one(FC), b, -one(FC), vₖ) else - vₖ .= b # r₀ = b + @kcopy!(n, vₖ, b) # r₀ = b end βₖ = @knrm2(n, vₖ) # β₁ = ‖v₁‖ if βₖ ≠ 0 diff --git a/src/minres.jl b/src/minres.jl index 3fd079893..472a8cf2d 100644 --- a/src/minres.jl +++ b/src/minres.jl @@ -166,12 +166,12 @@ kwargs_minres = (:M, :ldiv, :λ, :atol, :rtol, :etol, :conlim, :itmax, :timemax, (λ ≠ 0) && @kaxpy!(n, λ, Δx, r1) @kaxpby!(n, one(FC), b, -one(FC), r1) else - r1 .= b + @kcopy!(n, r1, b) # r1 ← b end # Initialize Lanczos process. # β₁ M v₁ = b. - r2 .= r1 + @kcopy!(n, r2, r1) # r2 ← r1 MisI || mulorldiv!(v, M, r1, ldiv) β₁ = @kdotr(m, r1, v) β₁ < 0 && error("Preconditioner is not positive definite") diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl index 3bc46e6aa..b664965c1 100644 --- a/src/minres_qlp.jl +++ b/src/minres_qlp.jl @@ -146,7 +146,7 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve (λ ≠ 0) && @kaxpy!(n, λ, Δx, M⁻¹vₖ) @kaxpby!(n, one(FC), b, -one(FC), M⁻¹vₖ) else - M⁻¹vₖ .= b + @kcopy!(n, M⁻¹vₖ, b) # M⁻¹vₖ ← b end # β₁v₁ = Mb @@ -373,9 +373,9 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve end # Update vₖ, M⁻¹vₖ₋₁, M⁻¹vₖ - MisI || (vₖ .= vₖ₊₁) - M⁻¹vₖ₋₁ .= M⁻¹vₖ - M⁻¹vₖ .= p + MisI || @kcopy!(n, vₖ, vₖ₊₁) # vₖ ← vₖ₊₁ + @kcopy!(n, M⁻¹vₖ₋₁, M⁻¹vₖ) # M⁻¹vₖ₋₁ ← M⁻¹vₖ + @kcopy!(n, M⁻¹vₖ, p) # M⁻¹vₖ ← p # Update ‖rₖ‖ estimate # ‖ rₖ ‖ = |ζbarₖ₊₁| diff --git a/src/symmlq.jl b/src/symmlq.jl index 94ad3655a..de3cf34f9 100644 --- a/src/symmlq.jl +++ b/src/symmlq.jl @@ -154,7 +154,7 @@ kwargs_symmlq = (:M, :ldiv, :transfer_to_cg, :λ, :λest, :atol, :rtol, :etol, : (λ ≠ 0) && @kaxpy!(n, λ, Δx, Mvold) @kaxpby!(n, one(FC), b, -one(FC), Mvold) else - Mvold .= b + @kcopy!(n, Mvold, b) # Mvold ← b end # Initialize Lanczos process. @@ -178,7 +178,7 @@ kwargs_symmlq = (:M, :ldiv, :transfer_to_cg, :λ, :λest, :atol, :rtol, :etol, : @kscal!(m, one(FC) / β, vold) MisI || @kscal!(m, one(FC) / β, Mvold) - w̅ .= vold + @kcopy!(n, w̅, vold) # w̅ ← vold mul!(Mv, A, vold) α = @kdotr(m, vold, Mv) + λ diff --git a/src/tricg.jl b/src/tricg.jl index 1bb4e5eea..a728cc2bf 100644 --- a/src/tricg.jl +++ b/src/tricg.jl @@ -212,7 +212,7 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax end # β₁Ev₁ = b ↔ β₁v₁ = Mb - M⁻¹vₖ .= b₀ + @kcopy!(m, M⁻¹vₖ, b₀) # M⁻¹vₖ ← b₀ MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv) βₖ = sqrt(@kdotr(m, vₖ, M⁻¹vₖ)) # β₁ = ‖v₁‖_E if βₖ ≠ 0 @@ -223,7 +223,7 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax end # γ₁Fu₁ = c ↔ γ₁u₁ = Nc - N⁻¹uₖ .= c₀ + @kcopy!(n, N⁻¹uₖ, c₀) # M⁻¹uₖ ← c₀ NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv) γₖ = sqrt(@kdotr(n, uₖ, N⁻¹uₖ)) # γ₁ = ‖u₁‖_F if γₖ ≠ 0 @@ -285,8 +285,8 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax @kaxpy!(n, -conj(αₖ), N⁻¹uₖ, p) # p ← p - ᾱₖ * N⁻¹uₖ # Update M⁻¹vₖ₋₁ and N⁻¹uₖ₋₁ - M⁻¹vₖ₋₁ .= M⁻¹vₖ - N⁻¹uₖ₋₁ .= N⁻¹uₖ + @kcopy!(m, M⁻¹vₖ₋₁, M⁻¹vₖ) # M⁻¹vₖ₋₁ ← M⁻¹vₖ + @kcopy!(n, N⁻¹uₖ₋₁, N⁻¹uₖ) # N⁻¹uₖ₋₁ ← N⁻¹uₖ # Notations : Wₖ = [w₁ ••• wₖ] = [v₁ 0 ••• vₖ 0 ] # [0 u₁ ••• 0 uₖ] @@ -398,8 +398,8 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax end # Update M⁻¹vₖ and N⁻¹uₖ - M⁻¹vₖ .= q - N⁻¹uₖ .= p + @kcopy!(m, M⁻¹vₖ, q) # M⁻¹vₖ ← q + @kcopy!(n, N⁻¹uₖ, p) # N⁻¹uₖ ← p # Compute ‖rₖ‖² = |γₖ₊₁ζ₂ₖ₋₁|² + |βₖ₊₁ζ₂ₖ|² ζ₂ₖ₋₁ = π₂ₖ₋₁ - conj(δₖ) * π₂ₖ diff --git a/src/trimr.jl b/src/trimr.jl index 56965d5cb..8177c8fa2 100644 --- a/src/trimr.jl +++ b/src/trimr.jl @@ -218,7 +218,7 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : end # β₁Ev₁ = b ↔ β₁v₁ = Mb - M⁻¹vₖ .= b₀ + @kcopy!(m, M⁻¹vₖ, b₀) # M⁻¹vₖ ← b₀ MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv) βₖ = sqrt(@kdotr(m, vₖ, M⁻¹vₖ)) # β₁ = ‖v₁‖_E if βₖ ≠ 0 @@ -229,7 +229,7 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : end # γ₁Fu₁ = c ↔ γ₁u₁ = Nc - N⁻¹uₖ .= c₀ + @kcopy!(n, N⁻¹uₖ, c₀) # N⁻¹uₖ ← c₀ NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv) γₖ = sqrt(@kdotr(n, uₖ, N⁻¹uₖ)) # γ₁ = ‖u₁‖_F if γₖ ≠ 0 @@ -484,16 +484,16 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : history && push!(rNorms, rNorm) # Update vₖ and uₖ - MisI || (vₖ .= vₖ₊₁) - NisI || (uₖ .= uₖ₊₁) + MisI || @kcopy!(m, vₖ, vₖ₊₁) # vₖ ← vₖ₊₁ + NisI || @kcopy!(n, uₖ, uₖ₊₁) # uₖ ← uₖ₊₁ # Update M⁻¹vₖ₋₁ and N⁻¹uₖ₋₁ - M⁻¹vₖ₋₁ .= M⁻¹vₖ - N⁻¹uₖ₋₁ .= N⁻¹uₖ + @kcopy!(m, M⁻¹vₖ₋₁, M⁻¹vₖ) # M⁻¹vₖ₋₁ ← M⁻¹vₖ + @kcopy!(n, N⁻¹uₖ₋₁, N⁻¹uₖ) # N⁻¹uₖ₋₁ ← N⁻¹uₖ # Update M⁻¹vₖ and N⁻¹uₖ - M⁻¹vₖ .= q - N⁻¹uₖ .= p + @kcopy!(m, M⁻¹vₖ, q) # M⁻¹vₖ ← q + @kcopy!(n, N⁻¹uₖ, p) # N⁻¹uₖ ← p # Update cosines and sines old_s₁ₖ = s₁ₖ From 9caf674f9a0afb906e877b75aa3535882269dcf5 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Wed, 16 Oct 2024 13:58:07 -0500 Subject: [PATCH 2/5] Fix a typo in cgls.jl --- src/cgls.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cgls.jl b/src/cgls.jl index 551268317..ad8992d17 100644 --- a/src/cgls.jl +++ b/src/cgls.jl @@ -147,7 +147,7 @@ kwargs_cgls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose Mq = MisI ? q : solver.Mr @kfill!(x, zero(FC)) - @copy!(m, r, b) # r ← b + @kcopy!(m, r, b) # r ← b bNorm = @knrm2(m, r) # Marginally faster than norm(b) if bNorm == 0 stats.niter = 0 From a65c74e40c4e8fd2f056edee057c1d7274f5f478 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Wed, 16 Oct 2024 14:04:26 -0500 Subject: [PATCH 3/5] Fix a typo in cgls_lanczos_shift.jl --- src/cgls_lanczos_shift.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cgls_lanczos_shift.jl b/src/cgls_lanczos_shift.jl index c03a0bc90..6f24dc3cd 100644 --- a/src/cgls_lanczos_shift.jl +++ b/src/cgls_lanczos_shift.jl @@ -141,7 +141,7 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose @kfill!(u_prev, zero(FC)) mul!(v, Aᴴ, u) # v₁ ← Aᴴ * b β = sqrt(@kdotr(n, v, v)) # β₁ = v₁ᵀ M v₁ - @kcopy!(nshifts, rNorms, β) # rNorms ← β + @kfill!(rNorms, β) if history for i = 1 : nshifts push!(rNorms_history[i], rNorms[i]) From e55ddd592a6ba0733cf60f4f792b430c0603343b Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Wed, 16 Oct 2024 14:07:26 -0500 Subject: [PATCH 4/5] Update src/cgs.jl --- src/cgs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cgs.jl b/src/cgs.jl index 692aaa010..fbe03ceeb 100644 --- a/src/cgs.jl +++ b/src/cgs.jl @@ -190,7 +190,7 @@ kwargs_cgs = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :hist kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %.2fs\n", iter, rNorm, ktimer(start_time)) @kcopy!(n, u, r) # u₀ - @kcopy!(n, p, r) # p₀ + @kcopy!(n, p, r) # p₀ @kfill!(q, zero(FC)) # q₋₁ # Stopping criterion. From 2d8d344bde63c6b928326173fe808d136ad904a7 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Wed, 16 Oct 2024 14:19:28 -0500 Subject: [PATCH 5/5] Fix an error in cgls_lanczos_shift.jl --- src/cgls_lanczos_shift.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cgls_lanczos_shift.jl b/src/cgls_lanczos_shift.jl index 6f24dc3cd..f3e751b74 100644 --- a/src/cgls_lanczos_shift.jl +++ b/src/cgls_lanczos_shift.jl @@ -168,7 +168,7 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose # Initialize some constants used in recursions below. ρ = one(T) - @kcopy!(nshifts, σ, β) # σ ← β + @kfill!(σ, β) @kfill!(δhat, zero(T)) @kfill!(ω, zero(T)) @kfill!(γ, one(T))