From 63c080661f685c7d13a783ba58bcefd9db2b949d Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Thu, 17 Oct 2024 17:26:29 -0500 Subject: [PATCH] Add a function knorm_elliptic --- src/car.jl | 4 ++-- src/cg_lanczos.jl | 4 ++-- src/cg_lanczos_shift.jl | 4 ++-- src/cgls_lanczos_shift.jl | 8 ++++---- src/cr.jl | 4 ++-- src/craig.jl | 6 +++--- src/craigmr.jl | 8 ++++---- src/krylov_utils.jl | 2 ++ src/lnlq.jl | 8 ++++---- src/lslq.jl | 8 ++++---- src/lsmr.jl | 8 ++++---- src/lsqr.jl | 6 +++--- src/minres_qlp.jl | 4 ++-- src/tricg.jl | 8 ++++---- src/trimr.jl | 8 ++++---- 15 files changed, 46 insertions(+), 44 deletions(-) diff --git a/src/car.jl b/src/car.jl index 37d02b03d..5d2ed68b6 100644 --- a/src/car.jl +++ b/src/car.jl @@ -157,7 +157,7 @@ kwargs_car = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :ca history && push!(rNorms, rNorm) # Compute ‖Ar₀‖ - ArNorm = MisI ? knorm(n, s) : sqrt(kdotr(n, r, u)) + ArNorm = MisI ? knorm(n, s) : knorm_elliptic(n, r, u) history && push!(ArNorms, ArNorm) if rNorm == 0 @@ -210,7 +210,7 @@ kwargs_car = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :ca kaxpby!(n, one(FC), t, β, u) # uₖ₊₁ = tₖ₊₁ + βₖ * uₖ # Compute ‖Arₖ‖ - ArNorm = MisI ? knorm(n, s) : sqrt(kdotr(n, r, u)) + ArNorm = MisI ? knorm(n, s) : knorm_elliptic(n, r, u) history && push!(ArNorms, ArNorm) end diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl index 636726c68..469492b25 100644 --- a/src/cg_lanczos.jl +++ b/src/cg_lanczos.jl @@ -138,7 +138,7 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax kcopy!(n, Mv, b) # Mv ← b end MisI || mulorldiv!(v, M, Mv, ldiv) # v₁ = M⁻¹r₀ - β = sqrt(kdotr(n, v, Mv)) # β₁ = v₁ᴴ M v₁ + β = knorm_elliptic(n, v, Mv) # β₁ = v₁ᴴ M v₁ σ = β rNorm = σ history && push!(rNorms, rNorm) @@ -201,7 +201,7 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax end kcopy!(n, Mv, Mv_next) # Mvₖ ← Mvₖ₊₁ MisI || mulorldiv!(v, M, Mv, ldiv) # vₖ₊₁ = M⁻¹ * Mvₖ₊₁ - β = sqrt(kdotr(n, v, Mv)) # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁ + β = knorm_elliptic(n, v, Mv) # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁ kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁ MisI || kscal!(n, one(FC) / β, Mv) # Mvₖ₊₁ ← Mvₖ₊₁ / βₖ₊₁ Anorm2 += β_prev^2 + β^2 + δ^2 # Use ‖Tₖ₊₁‖₂ as increasing approximation of ‖A‖₂. diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl index 0e0056f95..9e278ce4e 100644 --- a/src/cg_lanczos_shift.jl +++ b/src/cg_lanczos_shift.jl @@ -133,7 +133,7 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t end kcopy!(n, Mv, b) # Mv₁ ← b MisI || mulorldiv!(v, M, Mv, ldiv) # v₁ = M⁻¹ * Mv₁ - β = sqrt(kdotr(n, v, Mv)) # β₁ = v₁ᴴ M v₁ + β = knorm_elliptic(n, v, Mv) # β₁ = v₁ᴴ M v₁ kfill!(rNorms, β) if history for i = 1 : nshifts @@ -205,7 +205,7 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t end kcopy!(n, Mv, Mv_next) # Mvₖ ← Mvₖ₊₁ MisI || mulorldiv!(v, M, Mv, ldiv) # vₖ₊₁ = M⁻¹ * Mvₖ₊₁ - β = sqrt(kdotr(n, v, Mv)) # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁ + β = knorm_elliptic(n, v, Mv) # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁ kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁ MisI || kscal!(n, one(FC) / β, Mv) # Mvₖ₊₁ ← Mvₖ₊₁ / βₖ₊₁ diff --git a/src/cgls_lanczos_shift.jl b/src/cgls_lanczos_shift.jl index 14caba17d..02ced4abb 100644 --- a/src/cgls_lanczos_shift.jl +++ b/src/cgls_lanczos_shift.jl @@ -137,10 +137,10 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose kfill!(x[i], zero(FC)) # x₀ end - kcopy!(m, u, b) # 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₁ + mul!(v, Aᴴ, u) # v₁ ← Aᴴ * b + β = knorm_elliptic(n, v, v) # β₁ = v₁ᵀ M v₁ kfill!(rNorms, β) if history for i = 1 : nshifts @@ -203,7 +203,7 @@ kwargs_cgls_lanczos_shift = (:M, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose kaxpy!(m, -δ, u, utilde) # uₖ₊₁ = utildeₖ - δₖuₖ - βₖuₖ₋₁ kaxpy!(m, -β, u_prev, utilde) mul!(v, Aᴴ, utilde) # vₖ₊₁ = Aᴴuₖ₊₁ - β = sqrt(kdotr(n, v, v)) # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁ + β = knorm_elliptic(n, v, v) # βₖ₊₁ = vₖ₊₁ᵀ M vₖ₊₁ kscal!(n, one(FC) / β, v) # vₖ₊₁ ← vₖ₊₁ / βₖ₊₁ kscal!(m, one(FC) / β, utilde) # uₖ₊₁ = uₖ₊₁ / βₖ₊₁ kcopy!(m, u_prev, u) # u_prev ← u diff --git a/src/cr.jl b/src/cr.jl index 0c1ef450b..0fb897141 100644 --- a/src/cr.jl +++ b/src/cr.jl @@ -153,8 +153,8 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema mul!(Ar, A, r) ρ = kdotr(n, r, Ar) - rNorm = sqrt(kdotr(n, r, p)) # ‖r‖ - history && push!(rNorms, rNorm) # Values of ‖r‖ + rNorm = knorm_elliptic(n, r, p) # ‖r‖ + history && push!(rNorms, rNorm) # Values of ‖r‖ if ρ == 0 stats.niter = 0 diff --git a/src/craig.jl b/src/craig.jl index de865f1c5..b331fcfa4 100644 --- a/src/craig.jl +++ b/src/craig.jl @@ -204,7 +204,7 @@ kwargs_craig = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :btol, :conlim, :at kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) - β₁ = sqrt(kdotr(m, u, Mu)) + β₁ = knorm_elliptic(m, u, Mu) rNorm = β₁ history && push!(rNorms, rNorm) if β₁ == 0 @@ -270,7 +270,7 @@ kwargs_craig = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :btol, :conlim, :at mul!(Aᴴu, Aᴴ, u) kaxpby!(n, one(FC), Aᴴu, -β, Nv) NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) + α = knorm_elliptic(n, v, Nv) if α == 0 inconsistent = true continue @@ -313,7 +313,7 @@ kwargs_craig = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :btol, :conlim, :at mul!(Av, A, v) kaxpby!(m, one(FC), Av, -α, Mu) MisI || mulorldiv!(u, M, Mu, ldiv) - β = sqrt(kdotr(m, u, Mu)) + β = knorm_elliptic(m, u, Mu) if β ≠ 0 kscal!(m, one(FC) / β, u) MisI || kscal!(m, one(FC) / β, Mu) diff --git a/src/craigmr.jl b/src/craigmr.jl index 0e87622b7..7628e55a5 100644 --- a/src/craigmr.jl +++ b/src/craigmr.jl @@ -191,7 +191,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver kfill!(y, zero(FC)) kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) - β = sqrt(kdotr(m, u, Mu)) + β = knorm_elliptic(m, u, Mu) if β == 0 stats.niter = 0 stats.solved, stats.inconsistent = true, false @@ -210,7 +210,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver mul!(Aᴴu, Aᴴ, u) kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) + α = knorm_elliptic(n, v, Nv) Anorm² = α * α iter = 0 @@ -277,7 +277,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver mul!(Av, A, v) kaxpby!(m, one(FC), Av, -α, Mu) MisI || mulorldiv!(u, M, Mu, ldiv) - β = sqrt(kdotr(m, u, Mu)) + β = knorm_elliptic(m, u, Mu) if β ≠ 0 kscal!(m, one(FC)/β, u) MisI || kscal!(m, one(FC)/β, Mu) @@ -339,7 +339,7 @@ kwargs_craigmr = (:M, :N, :ldiv, :sqd, :λ, :atol, :rtol, :itmax, :timemax, :ver mul!(Aᴴu, Aᴴ, u) kaxpby!(n, one(FC), Aᴴu, -β, Nv) NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) + α = knorm_elliptic(n, v, Nv) Anorm² = Anorm² + α * α # = ‖Lₖ‖ ArNorm = α * β * abs(ζ/ρ) history && push!(ArNorms, ArNorm) diff --git a/src/krylov_utils.jl b/src/krylov_utils.jl index 14c0e430c..d88c5226e 100644 --- a/src/krylov_utils.jl +++ b/src/krylov_utils.jl @@ -315,6 +315,8 @@ kdotr(n :: Integer, x :: AbstractVector{Complex{T}}, y :: AbstractVector{Complex knorm(n :: Integer, x :: Vector{T}) where T <: BLAS.BlasFloat = BLAS.nrm2(n, x, 1) knorm(n :: Integer, x :: AbstractVector{T}) where T <: FloatOrComplex = norm(x) +knorm_elliptic(n :: Integer, x :: AbstractVector{T}, y :: AbstractVector{T}) where T <: FloatOrComplex = (x === y) ? knorm(n, x) : sqrt(kdotr(n, x, y)) + kscal!(n :: Integer, s :: T, x :: Vector{T}) where T <: BLAS.BlasFloat = BLAS.scal!(n, s, x, 1) kscal!(n :: Integer, s :: T, x :: AbstractVector{T}) where T <: FloatOrComplex = (x .*= s) kscal!(n :: Integer, s :: T, x :: AbstractVector{Complex{T}}) where T <: AbstractFloat = kscal!(n, Complex{T}(s), x) diff --git a/src/lnlq.jl b/src/lnlq.jl index 58a3cf611..a721be25a 100644 --- a/src/lnlq.jl +++ b/src/lnlq.jl @@ -228,7 +228,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly # β₁Mu₁ = b. kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) # u₁ = M⁻¹ * Mu₁ - βₖ = sqrt(kdotr(m, u, Mu)) # β₁ = ‖u₁‖_M + βₖ = knorm_elliptic(m, u, Mu) # β₁ = ‖u₁‖_M if βₖ ≠ 0 kscal!(m, one(FC) / βₖ, u) MisI || kscal!(m, one(FC) / βₖ, Mu) @@ -238,7 +238,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly mul!(Aᴴu, 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 + αₖ = knorm_elliptic(n, v, Nv) # α₁ = ‖v₁‖_N if αₖ ≠ 0 kscal!(n, one(FC) / αₖ, v) NisI || kscal!(n, one(FC) / αₖ, Nv) @@ -346,7 +346,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly mul!(Av, A, v) kaxpby!(m, one(FC), Av, -αₖ, Mu) MisI || mulorldiv!(u, M, Mu, ldiv) # uₖ₊₁ = M⁻¹ * Muₖ₊₁ - βₖ₊₁ = sqrt(kdotr(m, u, Mu)) # βₖ₊₁ = ‖uₖ₊₁‖_M + βₖ₊₁ = knorm_elliptic(m, u, Mu) # βₖ₊₁ = ‖uₖ₊₁‖_M if βₖ₊₁ ≠ 0 kscal!(m, one(FC) / βₖ₊₁, u) MisI || kscal!(m, one(FC) / βₖ₊₁, Mu) @@ -356,7 +356,7 @@ kwargs_lnlq = (:M, :N, :ldiv, :transfer_to_craig, :sqd, :λ, :σ, :utolx, :utoly mul!(Aᴴu, Aᴴ, u) kaxpby!(n, one(FC), Aᴴu, -βₖ₊₁, Nv) NisI || mulorldiv!(v, N, Nv, ldiv) # vₖ₊₁ = N⁻¹ * Nvₖ₊₁ - αₖ₊₁ = sqrt(kdotr(n, v, Nv)) # αₖ₊₁ = ‖vₖ₊₁‖_N + αₖ₊₁ = knorm_elliptic(n, v, Nv) # αₖ₊₁ = ‖vₖ₊₁‖_N if αₖ₊₁ ≠ 0 kscal!(n, one(FC) / αₖ₊₁, v) NisI || kscal!(n, one(FC) / αₖ₊₁, Nv) diff --git a/src/lslq.jl b/src/lslq.jl index c6b5bd42a..0aac95338 100644 --- a/src/lslq.jl +++ b/src/lslq.jl @@ -228,7 +228,7 @@ kwargs_lslq = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :σ, :etol, :utol, : # β₁ M u₁ = b. kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) - β₁ = sqrt(kdotr(m, u, Mu)) + β₁ = knorm_elliptic(m, u, Mu) if β₁ == 0 stats.niter = 0 stats.solved, stats.inconsistent = true, false @@ -246,7 +246,7 @@ kwargs_lslq = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :σ, :etol, :utol, : mul!(Aᴴu, Aᴴ, u) kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) # = α₁ + α = knorm_elliptic(n, v, Nv) # = α₁ # Aᴴb = 0 so x = 0 is a minimum least-squares solution if α == 0 @@ -326,7 +326,7 @@ kwargs_lslq = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :σ, :etol, :utol, : mul!(Av, A, v) kaxpby!(m, one(FC), Av, -α, Mu) MisI || mulorldiv!(u, M, Mu, ldiv) - β = sqrt(kdotr(m, u, Mu)) + β = knorm_elliptic(m, u, Mu) if β ≠ 0 kscal!(m, one(FC)/β, u) MisI || kscal!(m, one(FC)/β, Mu) @@ -335,7 +335,7 @@ kwargs_lslq = (:M, :N, :ldiv, :transfer_to_lsqr, :sqd, :λ, :σ, :etol, :utol, : mul!(Aᴴu, Aᴴ, u) kaxpby!(n, one(FC), Aᴴu, -β, Nv) NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) + α = knorm_elliptic(n, v, Nv) if α ≠ 0 kscal!(n, one(FC)/α, v) NisI || kscal!(n, one(FC)/α, Nv) diff --git a/src/lsmr.jl b/src/lsmr.jl index 14bb646c7..64b91e647 100644 --- a/src/lsmr.jl +++ b/src/lsmr.jl @@ -202,7 +202,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, # β₁ M u₁ = b. kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) - β₁ = sqrt(kdotr(m, u, Mu)) + β₁ = knorm_elliptic(m, u, Mu) if β₁ == 0 stats.niter = 0 stats.solved, stats.inconsistent = true, false @@ -219,7 +219,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, mul!(Aᴴu, Aᴴ, u) kcopy!(n, Nv, Aᴴu) # Nv ← Aᴴu NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) + α = knorm_elliptic(n, v, Nv) ζbar = α * β αbar = α @@ -295,7 +295,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, mul!(Av, A, v) kaxpby!(m, one(FC), Av, -α, Mu) MisI || mulorldiv!(u, M, Mu, ldiv) - β = sqrt(kdotr(m, u, Mu)) + β = knorm_elliptic(m, u, Mu) if β ≠ 0 kscal!(m, one(FC)/β, u) MisI || kscal!(m, one(FC)/β, Mu) @@ -304,7 +304,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, mul!(Aᴴu, Aᴴ, u) kaxpby!(n, one(FC), Aᴴu, -β, Nv) NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) + α = knorm_elliptic(n, v, Nv) if α ≠ 0 kscal!(n, one(FC)/α, v) NisI || kscal!(n, one(FC)/α, Nv) diff --git a/src/lsqr.jl b/src/lsqr.jl index e1493dcce..d980434e1 100644 --- a/src/lsqr.jl +++ b/src/lsqr.jl @@ -199,7 +199,7 @@ kwargs_lsqr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, # β₁ M u₁ = b. kcopy!(m, Mu, b) # Mu ← b MisI || mulorldiv!(u, M, Mu, ldiv) - β₁ = sqrt(kdotr(m, u, Mu)) + β₁ = knorm_elliptic(m, u, Mu) if β₁ == 0 stats.niter = 0 stats.solved, stats.inconsistent = true, false @@ -283,7 +283,7 @@ kwargs_lsqr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, mul!(Av, A, v) kaxpby!(m, one(FC), Av, -α, Mu) MisI || mulorldiv!(u, M, Mu, ldiv) - β = sqrt(kdotr(m, u, Mu)) + β = knorm_elliptic(m, u, Mu) if β ≠ 0 kscal!(m, one(FC)/β, u) MisI || kscal!(m, one(FC)/β, Mu) @@ -294,7 +294,7 @@ kwargs_lsqr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim, mul!(Aᴴu, Aᴴ, u) kaxpby!(n, one(FC), Aᴴu, -β, Nv) NisI || mulorldiv!(v, N, Nv, ldiv) - α = sqrt(kdotr(n, v, Nv)) + α = knorm_elliptic(n, v, Nv) if α ≠ 0 kscal!(n, one(FC)/α, v) NisI || kscal!(n, one(FC)/α, Nv) diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl index 085f7d728..abcfbbd59 100644 --- a/src/minres_qlp.jl +++ b/src/minres_qlp.jl @@ -151,7 +151,7 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve # β₁v₁ = Mb MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv) - βₖ = sqrt(kdotr(n, vₖ, M⁻¹vₖ)) + βₖ = knorm_elliptic(n, vₖ, M⁻¹vₖ) if βₖ ≠ 0 kscal!(n, one(FC) / βₖ, M⁻¹vₖ) MisI || kscal!(n, one(FC) / βₖ, vₖ) @@ -231,7 +231,7 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve MisI || mulorldiv!(vₖ₊₁, M, p, ldiv) # βₖ₊₁vₖ₊₁ = MAvₖ - γₖvₖ₋₁ - αₖvₖ - βₖ₊₁ = sqrt(kdotr(m, vₖ₊₁, p)) + βₖ₊₁ = knorm_elliptic(m, vₖ₊₁, p) # βₖ₊₁.ₖ ≠ 0 if βₖ₊₁ > btol diff --git a/src/tricg.jl b/src/tricg.jl index 66a9ebba0..cd2fd3c45 100644 --- a/src/tricg.jl +++ b/src/tricg.jl @@ -214,7 +214,7 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax # β₁Ev₁ = b ↔ β₁v₁ = Mb kcopy!(m, M⁻¹vₖ, b₀) # M⁻¹vₖ ← b₀ MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv) - βₖ = sqrt(kdotr(m, vₖ, M⁻¹vₖ)) # β₁ = ‖v₁‖_E + βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ) # β₁ = ‖v₁‖_E if βₖ ≠ 0 kscal!(m, one(FC) / βₖ, M⁻¹vₖ) MisI || kscal!(m, one(FC) / βₖ, vₖ) @@ -225,7 +225,7 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax # γ₁Fu₁ = c ↔ γ₁u₁ = Nc kcopy!(n, N⁻¹uₖ, c₀) # M⁻¹uₖ ← c₀ NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv) - γₖ = sqrt(kdotr(n, uₖ, N⁻¹uₖ)) # γ₁ = ‖u₁‖_F + γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ)) # γ₁ = ‖u₁‖_F if γₖ ≠ 0 kscal!(n, one(FC) / γₖ, N⁻¹uₖ) NisI || kscal!(n, one(FC) / γₖ, uₖ) @@ -382,8 +382,8 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax MisI || mulorldiv!(vₖ₊₁, M, q, ldiv) # βₖ₊₁vₖ₊₁ = MAuₖ - γₖvₖ₋₁ - αₖvₖ NisI || mulorldiv!(uₖ₊₁, N, p, ldiv) # γₖ₊₁uₖ₊₁ = NAᴴvₖ - βₖuₖ₋₁ - ᾱₖuₖ - βₖ₊₁ = sqrt(kdotr(m, vₖ₊₁, q)) # βₖ₊₁ = ‖vₖ₊₁‖_E - γₖ₊₁ = sqrt(kdotr(n, uₖ₊₁, p)) # γₖ₊₁ = ‖uₖ₊₁‖_F + βₖ₊₁ = knorm_elliptic(m, vₖ₊₁, q) # βₖ₊₁ = ‖vₖ₊₁‖_E + γₖ₊₁ = knorm_elliptic(n, uₖ₊₁, p) # γₖ₊₁ = ‖uₖ₊₁‖_F # βₖ₊₁ ≠ 0 if βₖ₊₁ > btol diff --git a/src/trimr.jl b/src/trimr.jl index 5b0124e8b..5093c7beb 100644 --- a/src/trimr.jl +++ b/src/trimr.jl @@ -220,7 +220,7 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : # β₁Ev₁ = b ↔ β₁v₁ = Mb kcopy!(m, M⁻¹vₖ, b₀) # M⁻¹vₖ ← b₀ MisI || mulorldiv!(vₖ, M, M⁻¹vₖ, ldiv) - βₖ = sqrt(kdotr(m, vₖ, M⁻¹vₖ)) # β₁ = ‖v₁‖_E + βₖ = knorm_elliptic(m, vₖ, M⁻¹vₖ) # β₁ = ‖v₁‖_E if βₖ ≠ 0 kscal!(m, one(FC) / βₖ, M⁻¹vₖ) MisI || kscal!(m, one(FC) / βₖ, vₖ) @@ -231,7 +231,7 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : # γ₁Fu₁ = c ↔ γ₁u₁ = Nc kcopy!(n, N⁻¹uₖ, c₀) # N⁻¹uₖ ← c₀ NisI || mulorldiv!(uₖ, N, N⁻¹uₖ, ldiv) - γₖ = sqrt(kdotr(n, uₖ, N⁻¹uₖ)) # γ₁ = ‖u₁‖_F + γₖ = knorm_elliptic(n, uₖ, N⁻¹uₖ) # γ₁ = ‖u₁‖_F if γₖ ≠ 0 kscal!(n, one(FC) / γₖ, N⁻¹uₖ) NisI || kscal!(n, one(FC) / γₖ, uₖ) @@ -302,8 +302,8 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, : MisI || mulorldiv!(vₖ₊₁, M, q, ldiv) # βₖ₊₁vₖ₊₁ = MAuₖ - γₖvₖ₋₁ - αₖvₖ NisI || mulorldiv!(uₖ₊₁, N, p, ldiv) # γₖ₊₁uₖ₊₁ = NAᴴvₖ - βₖuₖ₋₁ - ᾱₖuₖ - βₖ₊₁ = sqrt(kdotr(m, vₖ₊₁, q)) # βₖ₊₁ = ‖vₖ₊₁‖_E - γₖ₊₁ = sqrt(kdotr(n, uₖ₊₁, p)) # γₖ₊₁ = ‖uₖ₊₁‖_F + βₖ₊₁ = knorm_elliptic(m, vₖ₊₁, q) # βₖ₊₁ = ‖vₖ₊₁‖_E + γₖ₊₁ = knorm_elliptic(n, uₖ₊₁, p) # γₖ₊₁ = ‖uₖ₊₁‖_F # βₖ₊₁ ≠ 0 if βₖ₊₁ > btol