From 8ecc2935fc22dbc42f16a0b7827c7064ea73bddc Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Mon, 20 May 2024 16:20:36 -0400 Subject: [PATCH 1/4] Add tests with LinearOperators.jl (#863) * Add tests with LinearOperators.jl * Update block_gmres for complex numbers * Update block_gmres.jl --- .buildkite/pipeline.yml | 14 ++++++++++++++ src/block_gmres.jl | 5 +++-- test/cpu/linear_operators.jl | 16 ++++++++++++++++ 3 files changed, 33 insertions(+), 2 deletions(-) create mode 100644 test/cpu/linear_operators.jl diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 120f1f736..358eb9cac 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -98,3 +98,17 @@ steps: Pkg.instantiate() include("test/cpu/component_arrays.jl")' timeout_in_minutes: 30 + + - label: "CPUs -- LinearOperators.jl" + plugins: + - JuliaCI/julia#v1: + version: "1.10" + agents: + queue: "juliaecosystem" + command: | + julia --color=yes --project -e ' + using Pkg + Pkg.add("LinearOperators") + Pkg.instantiate() + include("test/cpu/linear_operators.jl")' + timeout_in_minutes: 30 diff --git a/src/block_gmres.jl b/src/block_gmres.jl index c1e04164d..ea0b4b989 100644 --- a/src/block_gmres.jl +++ b/src/block_gmres.jl @@ -148,6 +148,7 @@ kwargs_block_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rto # Define the blocks D1 and D2 D1 = view(D, 1:p, :) D2 = view(D, p+1:2p, :) + trans = FC <: AbstractFloat ? 'T' : 'C' # Coefficients for mul! α = -one(FC) @@ -263,7 +264,7 @@ kwargs_block_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rto for i = 1 : inner_iter-1 D1 .= R[nr+i] D2 .= R[nr+i+1] - @kormqr!('L', 'T', H[i], τ[i], D) + @kormqr!('L', trans, H[i], τ[i], D) R[nr+i] .= D1 R[nr+i+1] .= D2 end @@ -276,7 +277,7 @@ kwargs_block_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rto # Update Zₖ = (Qₖ)ᴴΓE₁ = (Λ₁, ..., Λₖ, Λbarₖ₊₁) D1 .= Z[inner_iter] D2 .= zero(FC) - @kormqr!('L', 'T', H[inner_iter], τ[inner_iter], D) + @kormqr!('L', trans, H[inner_iter], τ[inner_iter], D) Z[inner_iter] .= D1 # Update residual norm estimate. diff --git a/test/cpu/linear_operators.jl b/test/cpu/linear_operators.jl new file mode 100644 index 000000000..3a58c5d20 --- /dev/null +++ b/test/cpu/linear_operators.jl @@ -0,0 +1,16 @@ +using LinearAlgebra, SparseArrays, Test +using Krylov, LinearOperators + +@testset "LinearOperators" begin + n = 50 + p = 5 + + for T in (Float64, ComplexF64) + A = rand(T, n, n) + B = rand(T, n, p) + + opA = LinearOperator(A) + x, stats = block_gmres(opA, B) + @test stats.solved + end +end From 31fe591e4758ec8d1957e12063bd4ed4f27108cc Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Mon, 20 May 2024 17:46:43 -0400 Subject: [PATCH 2/4] Support preconditioners in BiLQ and QMR (#862) * Support preconditioners in BiLQ and QMR * Fix preconditioning with ldiv=true in BiLQ and QMR * Update bilq.jl and qmr.jl --- docs/src/preconditioners.md | 5 ++- src/bilq.jl | 73 +++++++++++++++++++++++++++---------- src/krylov_solvers.jl | 12 +++++- src/qmr.jl | 70 ++++++++++++++++++++++++++--------- src/usymlq.jl | 3 +- src/usymqr.jl | 3 +- test/test_bicgstab.jl | 4 +- test/test_bilq.jl | 25 ++++++++++++- test/test_cgs.jl | 4 +- test/test_qmr.jl | 24 ++++++++++++ 10 files changed, 175 insertions(+), 48 deletions(-) diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md index 22e24ce90..38da24cc5 100644 --- a/docs/src/preconditioners.md +++ b/docs/src/preconditioners.md @@ -30,12 +30,13 @@ Krylov.jl supports both approaches thanks to the argument `ldiv` of the Krylov s ## How to use preconditioners in Krylov.jl? !!! info - - A preconditioner only need support the operation `mul!(y, P⁻¹, x)` when `ldiv=false` or `ldiv!(y, P, x)` when `ldiv=true` to be used in Krylov.jl. + - A preconditioner only needs to support the operation `mul!(y, P⁻¹, x)` when `ldiv=false` or `ldiv!(y, P, x)` when `ldiv=true` to be used in Krylov.jl. + - Additional support for `adjoint` with preconditioners is required in the methods [`BILQ`](@ref bilq) and [`QMR`](@ref qmr). - The default value of a preconditioner in Krylov.jl is the identity operator `I`. ### Square non-Hermitian linear systems -Methods concerned: [`CGS`](@ref cgs), [`BiCGSTAB`](@ref bicgstab), [`DQGMRES`](@ref dqgmres), [`GMRES`](@ref gmres), [`BLOCK-GMRES`](@ref block_gmres), [`FGMRES`](@ref fgmres), [`DIOM`](@ref diom) and [`FOM`](@ref fom). +Methods concerned: [`CGS`](@ref cgs), [`BILQ`](@ref bilq), [`QMR`](@ref qmr), [`BiCGSTAB`](@ref bicgstab), [`DQGMRES`](@ref dqgmres), [`GMRES`](@ref gmres), [`BLOCK-GMRES`](@ref block_gmres), [`FGMRES`](@ref fgmres), [`DIOM`](@ref diom) and [`FOM`](@ref fom). A Krylov method dedicated to non-Hermitian linear systems allows the three variants of preconditioning. diff --git a/src/bilq.jl b/src/bilq.jl index 2e8823e93..45a0051ef 100644 --- a/src/bilq.jl +++ b/src/bilq.jl @@ -15,8 +15,9 @@ export bilq, bilq! """ (x, stats) = bilq(A, b::AbstractVector{FC}; c::AbstractVector{FC}=b, transfer_to_bicg::Bool=true, - atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0, - timemax::Float64=Inf, verbose::Int=0, history::Bool=false, + M=I, N=I, ldiv::Bool=false, atol::T=√eps(T), + rtol::T=√eps(T), itmax::Int=0, timemax::Float64=Inf, + verbose::Int=0, history::Bool=false, callback=solver->false, iostream::IO=kstdout) `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`. @@ -30,6 +31,7 @@ Solve the square linear system Ax = b of size n using BiLQ. BiLQ is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`. The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`. When `A` is Hermitian and `b = c`, BiLQ is equivalent to SYMMLQ. +BiLQ requires support for `adjoint(M)` and `adjoint(N)` if preconditioners are provided. #### Input arguments @@ -44,6 +46,9 @@ When `A` is Hermitian and `b = c`, BiLQ is equivalent to SYMMLQ. * `c`: the second initial vector of length `n` required by the Lanczos biorthogonalization process; * `transfer_to_bicg`: transfer from the BiLQ point to the BiCG point, when it exists. The transfer is based on the residual norm; +* `M`: linear operator that models a nonsingular matrix of size `n` used for left preconditioning; +* `N`: linear operator that models a nonsingular matrix of size `n` used for right preconditioning; +* `ldiv`: define whether the preconditioners use `ldiv!` or `mul!`; * `atol`: absolute stopping tolerance based on the residual norm; * `rtol`: relative stopping tolerance based on the residual norm; * `itmax`: the maximum number of iterations. If `itmax=0`, the default number of iterations is set to `2n`; @@ -82,6 +87,9 @@ def_optargs_bilq = (:(x0::AbstractVector),) def_kwargs_bilq = (:(; c::AbstractVector{FC} = b ), :(; transfer_to_bicg::Bool = true), + :(; M = I ), + :(; N = I ), + :(; ldiv::Bool = false ), :(; atol::T = √eps(T) ), :(; rtol::T = √eps(T) ), :(; itmax::Int = 0 ), @@ -95,7 +103,7 @@ def_kwargs_bilq = mapreduce(extract_parameters, vcat, def_kwargs_bilq) args_bilq = (:A, :b) optargs_bilq = (:x0,) -kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) +kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) @eval begin function bilq($(def_args_bilq...), $(def_optargs_bilq...); $(def_kwargs_bilq...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}} @@ -131,26 +139,42 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, length(b) == m || error("Inconsistent problem size") (verbose > 0) && @printf(iostream, "BILQ: system of size %d\n", n) + # Check M = Iₙ and N = Iₙ + MisI = (M === I) + NisI = (N === I) + # Check type consistency eltype(A) == FC || @warn "eltype(A) ≠ $FC. This could lead to errors or additional allocations in operator-vector products." ktypeof(b) <: S || error("ktypeof(b) is not a subtype of $S") ktypeof(c) <: S || error("ktypeof(c) is not a subtype of $S") - # Compute the adjoint of A + # Compute the adjoint of A, M and N Aᴴ = A' + Mᴴ = M' + Nᴴ = N' # Set up workspace. + allocate_if(!MisI, solver, :t, S, n) + allocate_if(!NisI, solver, :s, S, n) uₖ₋₁, uₖ, q, vₖ₋₁, vₖ = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ p, Δx, x, d̅, stats = solver.p, solver.Δx, solver.x, solver.d̅, solver.stats warm_start = solver.warm_start rNorms = stats.residuals reset!(stats) r₀ = warm_start ? q : b + Mᴴuₖ = MisI ? uₖ : solver.t + t = MisI ? q : solver.t + Nvₖ = NisI ? vₖ : solver.s + s = NisI ? p : solver.s if warm_start mul!(r₀, A, Δx) @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₀‖. x .= zero(FC) @@ -170,10 +194,6 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, iter = 0 itmax == 0 && (itmax = 2*n) - ε = atol + rtol * bNorm - (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)) - # Initialize the Lanczos biorthogonalization process. cᴴb = @kdot(n, c, r₀) # ⟨c,r₀⟩ if cᴴb == 0 @@ -186,6 +206,10 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, return solver end + ε = atol + rtol * bNorm + (verbose > 0) && @printf(iostream, "%5s %8s %7s %5s\n", "k", "αₖ", "‖rₖ‖", "timer") + kdisplay(iter, verbose) && @printf(iostream, "%5d %8.1e %7.1e %.2fs\n", iter, cᴴb, bNorm, ktimer(start_time)) + βₖ = √(abs(cᴴb)) # β₁γ₁ = cᴴ(b - Ax₀) γₖ = cᴴb / βₖ # β₁γ₁ = cᴴ(b - Ax₀) vₖ₋₁ .= zero(FC) # v₀ = 0 @@ -214,23 +238,30 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, iter = iter + 1 # Continue the Lanczos biorthogonalization process. - # AVₖ = VₖTₖ + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ - # AᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ + # MANVₖ = VₖTₖ + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ + # NᴴAᴴMᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ - mul!(q, A , vₖ) # Forms vₖ₊₁ : q ← Avₖ - mul!(p, Aᴴ, uₖ) # Forms uₖ₊₁ : p ← Aᴴuₖ + # Forms vₖ₊₁ : q ← MANvₖ + NisI || mulorldiv!(Nvₖ, N, vₖ, ldiv) + mul!(t, A, Nvₖ) + MisI || mulorldiv!(q, M, t, ldiv) + + # Forms uₖ₊₁ : p ← NᴴAᴴMᴴuₖ + MisI || mulorldiv!(Mᴴuₖ, Mᴴ, uₖ, ldiv) + 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ₖ₋₁ - αₖ = @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 ] @@ -353,7 +384,7 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, breakdown = !solved_lq && !solved_cg && (pᴴq == 0) timer = time_ns() - start_time overtimed = timer > timemax_ns - kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %.2fs\n", iter, rNorm_lq, ktimer(start_time)) + kdisplay(iter, verbose) && @printf(iostream, "%5d %8.1e %7.1e %.2fs\n", iter, αₖ, rNorm_lq, ktimer(start_time)) end (verbose > 0) && @printf(iostream, "\n") @@ -372,6 +403,10 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, overtimed && (status = "time limit exceeded") # Update x + if !NisI + copyto!(solver.s, x) + mulorldiv!(x, N, solver.s, ldiv) + end warm_start && @kaxpy!(n, one(FC), Δx, x) solver.warm_start = false diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl index 13c98f823..d9f9fd9b1 100644 --- a/src/krylov_solvers.jl +++ b/src/krylov_solvers.jl @@ -1002,6 +1002,8 @@ mutable struct BilqSolver{T,FC,S} <: KrylovSolver{T,FC,S} Δx :: S x :: S d̅ :: S + t :: S + s :: S warm_start :: Bool stats :: SimpleStats{T} end @@ -1018,8 +1020,10 @@ function BilqSolver(m, n, S) Δx = S(undef, 0) x = S(undef, n) d̅ = S(undef, n) + t = S(undef, 0) + s = S(undef, 0) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") - solver = BilqSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, d̅, false, stats) + solver = BilqSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, d̅, t, s, false, stats) return solver end @@ -1052,6 +1056,8 @@ mutable struct QmrSolver{T,FC,S} <: KrylovSolver{T,FC,S} x :: S wₖ₋₂ :: S wₖ₋₁ :: S + t :: S + s :: S warm_start :: Bool stats :: SimpleStats{T} end @@ -1069,8 +1075,10 @@ function QmrSolver(m, n, S) x = S(undef, n) wₖ₋₂ = S(undef, n) wₖ₋₁ = S(undef, n) + t = S(undef, 0) + s = S(undef, 0) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") - solver = QmrSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, wₖ₋₂, wₖ₋₁, false, stats) + solver = QmrSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, wₖ₋₂, wₖ₋₁, t, s, false, stats) return solver end diff --git a/src/qmr.jl b/src/qmr.jl index 995392f0c..a197967a8 100644 --- a/src/qmr.jl +++ b/src/qmr.jl @@ -22,7 +22,7 @@ export qmr, qmr! """ (x, stats) = qmr(A, b::AbstractVector{FC}; - c::AbstractVector{FC}=b, atol::T=√eps(T), + c::AbstractVector{FC}=b, M=I, N=I, ldiv::Bool=false, atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0, timemax::Float64=Inf, verbose::Int=0, history::Bool=false, callback=solver->false, iostream::IO=kstdout) @@ -38,6 +38,7 @@ Solve the square linear system Ax = b of size n using QMR. QMR is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`. The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`. When `A` is Hermitian and `b = c`, QMR is equivalent to MINRES. +QMR requires support for `adjoint(M)` and `adjoint(N)` if preconditioners are provided. #### Input arguments @@ -51,6 +52,9 @@ When `A` is Hermitian and `b = c`, QMR is equivalent to MINRES. #### Keyword arguments * `c`: the second initial vector of length `n` required by the Lanczos biorthogonalization process; +* `M`: linear operator that models a nonsingular matrix of size `n` used for left preconditioning; +* `N`: linear operator that models a nonsingular matrix of size `n` used for right preconditioning; +* `ldiv`: define whether the preconditioners use `ldiv!` or `mul!`; * `atol`: absolute stopping tolerance based on the residual norm; * `rtol`: relative stopping tolerance based on the residual norm; * `itmax`: the maximum number of iterations. If `itmax=0`, the default number of iterations is set to `2n`; @@ -89,6 +93,9 @@ def_args_qmr = (:(A ), def_optargs_qmr = (:(x0::AbstractVector),) def_kwargs_qmr = (:(; c::AbstractVector{FC} = b ), + :(; M = I ), + :(; N = I ), + :(; ldiv::Bool = false ), :(; atol::T = √eps(T) ), :(; rtol::T = √eps(T) ), :(; itmax::Int = 0 ), @@ -102,7 +109,7 @@ def_kwargs_qmr = mapreduce(extract_parameters, vcat, def_kwargs_qmr) args_qmr = (:A, :b) optargs_qmr = (:x0,) -kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) +kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream) @eval begin function qmr($(def_args_qmr...), $(def_optargs_qmr...); $(def_kwargs_qmr...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}} @@ -138,26 +145,42 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, length(b) == m || error("Inconsistent problem size") (verbose > 0) && @printf(iostream, "QMR: system of size %d\n", n) + # Check M = Iₙ and N = Iₙ + MisI = (M === I) + NisI = (N === I) + # Check type consistency eltype(A) == FC || @warn "eltype(A) ≠ $FC. This could lead to errors or additional allocations in operator-vector products." ktypeof(b) <: S || error("ktypeof(b) is not a subtype of $S") ktypeof(c) <: S || error("ktypeof(c) is not a subtype of $S") - # Compute the adjoint of A + # Compute the adjoint of A, M and N Aᴴ = A' + Mᴴ = M' + Nᴴ = N' # Set up workspace. + allocate_if(!MisI, solver, :t, S, n) + allocate_if(!NisI, solver, :s, S, n) uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ, solver.p Δx, x, wₖ₋₂, wₖ₋₁, stats = solver.Δx, solver.x, solver.wₖ₋₂, solver.wₖ₋₁, solver.stats warm_start = solver.warm_start rNorms = stats.residuals reset!(stats) r₀ = warm_start ? q : b + Mᴴuₖ = MisI ? uₖ : solver.t + t = MisI ? q : solver.t + Nvₖ = NisI ? vₖ : solver.s + s = NisI ? p : solver.s if warm_start mul!(r₀, A, Δx) @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₀‖. x .= zero(FC) @@ -177,10 +200,6 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, iter = 0 itmax == 0 && (itmax = 2*n) - ε = atol + rtol * rNorm - (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)) - # Initialize the Lanczos biorthogonalization process. cᴴb = @kdot(n, c, r₀) # ⟨c,r₀⟩ if cᴴb == 0 @@ -193,6 +212,10 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, return solver end + ε = atol + rtol * rNorm + (verbose > 0) && @printf(iostream, "%5s %8s %7s %5s\n", "k", "αₖ", "‖rₖ‖", "timer") + kdisplay(iter, verbose) && @printf(iostream, "%5d %8.1e %7.1e %.2fs\n", iter, cᴴb, rNorm, ktimer(start_time)) + βₖ = √(abs(cᴴb)) # β₁γ₁ = cᴴ(b - Ax₀) γₖ = cᴴb / βₖ # β₁γ₁ = cᴴ(b - Ax₀) vₖ₋₁ .= zero(FC) # v₀ = 0 @@ -219,23 +242,30 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, iter = iter + 1 # Continue the Lanczos biorthogonalization process. - # AVₖ = VₖTₖ + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ - # AᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ + # MANVₖ = VₖTₖ + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ + # NᴴAᴴMᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ - mul!(q, A , vₖ) # Forms vₖ₊₁ : q ← Avₖ - mul!(p, Aᴴ, uₖ) # Forms uₖ₊₁ : p ← Aᴴuₖ + # Forms vₖ₊₁ : q ← MANvₖ + NisI || mulorldiv!(Nvₖ, N, vₖ, ldiv) + mul!(t, A, Nvₖ) + MisI || mulorldiv!(q, M, t, ldiv) + + # Forms uₖ₊₁ : p ← NᴴAᴴMᴴuₖ + MisI || mulorldiv!(Mᴴuₖ, Mᴴ, uₖ, ldiv) + 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ₖ₋₁ - αₖ = @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 QR factorization of Tₖ₊₁.ₖ = Qₖ [ Rₖ ]. # [ Oᵀ ] @@ -357,7 +387,7 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, breakdown = !solved && (pᴴq == 0) timer = time_ns() - start_time overtimed = timer > timemax_ns - kdisplay(iter, verbose) && @printf(iostream, "%5d %7.1e %.2fs\n", iter, rNorm, ktimer(start_time)) + kdisplay(iter, verbose) && @printf(iostream, "%5d %8.1e %7.1e %.2fs\n", iter, αₖ, rNorm, ktimer(start_time)) end (verbose > 0) && @printf(iostream, "\n") @@ -369,6 +399,10 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, overtimed && (status = "time limit exceeded") # Update x + if !NisI + copyto!(solver.s, x) + mulorldiv!(x, N, solver.s, ldiv) + end warm_start && @kaxpy!(n, one(FC), Δx, x) solver.warm_start = false diff --git a/src/usymlq.jl b/src/usymlq.jl index b80f0a622..20fe5809d 100644 --- a/src/usymlq.jl +++ b/src/usymlq.jl @@ -38,7 +38,8 @@ USYMLQ determines the least-norm solution of the consistent linear system Ax = b USYMLQ is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`. The vector `c` is only used to initialize the process and a default value can be `b` or `Aᴴb` depending on the shape of `A`. The error norm ‖x - x*‖ monotonously decreases in USYMLQ. -It's considered as a generalization of SYMMLQ. +When `A` is Hermitian and `b = c`, USYMLQ is equivalent to SYMMLQ. +USYMLQ is considered as a generalization of SYMMLQ. It can also be applied to under-determined and over-determined problems. In all cases, problems must be consistent. diff --git a/src/usymqr.jl b/src/usymqr.jl index 0aae23335..8ba7203bc 100644 --- a/src/usymqr.jl +++ b/src/usymqr.jl @@ -38,7 +38,8 @@ USYMQR solves Ax = b if it is consistent. USYMQR is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`. The vector `c` is only used to initialize the process and a default value can be `b` or `Aᴴb` depending on the shape of `A`. The residual norm ‖b - Ax‖ monotonously decreases in USYMQR. -It's considered as a generalization of MINRES. +When `A` is Hermitian and `b = c`, QMR is equivalent to MINRES. +USYMQR is considered as a generalization of MINRES. It can also be applied to under-determined and over-determined problems. USYMQR finds the minimum-norm solution if problems are inconsistent. diff --git a/test/test_bicgstab.jl b/test/test_bicgstab.jl index 6817acf3d..bde109fdd 100644 --- a/test/test_bicgstab.jl +++ b/test/test_bicgstab.jl @@ -62,7 +62,7 @@ A, b, M = square_preconditioned(FC=FC) (x, stats) = bicgstab(A, b, M=M) r = b - A * x - resid = norm(r) / norm(b) + resid = norm(M * r) / norm(M * b) @test(resid ≤ bicgstab_tol) @test(stats.solved) @@ -78,7 +78,7 @@ A, b, M, N = two_preconditioners(500, 32) (x, stats) = bicgstab(A, b, M=M, N=N) r = b - A * x - resid = norm(r) / norm(b) + resid = norm(M * r) / norm(M * b) @test(resid ≤ bicgstab_tol) @test(stats.solved) diff --git a/test/test_bilq.jl b/test/test_bilq.jl index 40b9872db..672097187 100644 --- a/test/test_bilq.jl +++ b/test/test_bilq.jl @@ -66,11 +66,34 @@ @test(resid ≤ bilq_tol) @test(stats.solved) + # Left preconditioning + A, b, M = square_preconditioned(FC=FC) + (x, stats) = bilq(A, b, M=M) + r = b - A * x + resid = norm(M * r) / norm(M * b) + @test(resid ≤ bilq_tol) + @test(stats.solved) + + # Right preconditioning + A, b, N = square_preconditioned(FC=FC) + (x, stats) = bilq(A, b, N=N) + r = b - A * x + resid = norm(r) / norm(b) + @test(resid ≤ bilq_tol) + @test(stats.solved) + + # Split preconditioning + A, b, M, N = two_preconditioners(FC=FC) + (x, stats) = bilq(A, b, M=M, N=N) + r = b - A * x + resid = norm(M * r) / norm(M * b) + @test(resid ≤ bilq_tol) + @test(stats.solved) + # Test bᴴc == 0 A, b, c = bc_breakdown(FC=FC) (x, stats) = bilq(A, b, c=c) @test stats.status == "Breakdown bᴴc = 0" - # test callback function A, b = polar_poisson(FC=FC) diff --git a/test/test_cgs.jl b/test/test_cgs.jl index 832cd76c3..e4aaca4ef 100644 --- a/test/test_cgs.jl +++ b/test/test_cgs.jl @@ -62,7 +62,7 @@ A, b, N = square_preconditioned(FC=FC) (x, stats) = cgs(A, b, N=N) r = b - A * x - resid = norm(r) / norm(b) + resid = norm(M * r) / norm(M * b) @test(resid ≤ cgs_tol) @test(stats.solved) @@ -70,7 +70,7 @@ A, b, M, N = two_preconditioners(FC=FC) (x, stats) = cgs(A, b, M=M, N=N) r = b - A * x - resid = norm(r) / norm(b) + resid = norm(M * r) / norm(M * b) @test(resid ≤ cgs_tol) @test(stats.solved) diff --git a/test/test_qmr.jl b/test/test_qmr.jl index 4a6b8c1c9..91fc77dd0 100644 --- a/test/test_qmr.jl +++ b/test/test_qmr.jl @@ -58,6 +58,30 @@ @test(resid ≤ qmr_tol) @test(stats.solved) + # Left preconditioning + A, b, M = square_preconditioned(FC=FC) + (x, stats) = qmr(A, b, M=M) + r = b - A * x + resid = norm(M * r) / norm(M * b) + @test(resid ≤ qmr_tol) + @test(stats.solved) + + # Right preconditioning + A, b, N = square_preconditioned(FC=FC) + (x, stats) = qmr(A, b, N=N) + r = b - A * x + resid = norm(r) / norm(b) + @test(resid ≤ qmr_tol) + @test(stats.solved) + + # Split preconditioning + A, b, M, N = two_preconditioners(FC=FC) + (x, stats) = qmr(A, b, M=M, N=N) + r = b - A * x + resid = norm(M * r) / norm(M * b) + @test(resid ≤ qmr_tol) + @test(stats.solved) + # Test bᴴc == 0 A, b, c = bc_breakdown(FC=FC) (x, stats) = qmr(A, b, c=c) From 5605b361552648ebc3da081854887bd4f3f094bd Mon Sep 17 00:00:00 2001 From: Alexis Montoison <35051714+amontoison@users.noreply.github.com> Date: Mon, 20 May 2024 18:48:11 -0400 Subject: [PATCH 3/4] Don't use the macro @. (#864) --- 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..5bda9802a 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, vₖ, d̅) # 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..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₁ - @. d̅ = vₖ + @kcopy!(n, vₖ, d̅) # 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..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ₖ₋₁ - @. Mv_prev = Mv # Mvₖ₋₁ ← Mvₖ + @kcopy!(n, Mv, Mv_prev) # Mvₖ₋₁ ← Mvₖ end - @. 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/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl index b523e5cc3..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ₖ₋₁ - @. Mv_prev = Mv # Mvₖ₋₁ ← Mvₖ + @kcopy!(n, Mv, Mv_prev) # Mvₖ₋₁ ← Mvₖ end - @. 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/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..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) - @. r1 = r2 - @. 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 5bc3399eb..9cb569f6a 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, vₖ, wₖ) 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..2a1ad578c 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..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) - @. Mvold = Mv + @kcopy!(m, Mv, Mvold) # Mvold ← Mv @kaxpy!(m, -α, Mv, Mv_next) - @. 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 8250e6dfc..99d714bdd 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, vₖ, gx₂ₖ₋₁) # gx₂ₖ₋₁ ← vₖ + gx₂ₖ .= -conj(δₖ) .* gx₂ₖ₋₁ + @kcopy!(n, uₖ, gy₂ₖ) # 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..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₁ - @. d̅ = uₖ + @kcopy!(n, uₖ, d̅) # 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..73d6c7114 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, uₖ, d̅) # d̅ ← vₖ 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..1d7488e36 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 a843fb3b26e35481c1f60b2913e46e3a8c3c2817 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Wed, 22 May 2024 08:31:26 -0400 Subject: [PATCH 4/4] Fix copy_triangle for Julia 1.11 --- src/block_krylov_utils.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/block_krylov_utils.jl b/src/block_krylov_utils.jl index 68ba46458..9f1a9616d 100644 --- a/src/block_krylov_utils.jl +++ b/src/block_krylov_utils.jl @@ -164,7 +164,13 @@ function copy_triangle(Q::AbstractMatrix{FC}, R::AbstractMatrix{FC}, k::Int) whe end end else - copytrito!(R, Q, 'U') + mR, nR = size(R) + mQ, nQ = size(Q) + if (mR == mQ) && (nR == nQ) + copytrito!(R, Q, 'U') + else + copytrito!(R, view(Q, 1:k, 1:k), 'U') + end end return R end