diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl index 2a2eef1d5..905eda930 100644 --- a/src/krylov_solvers.jl +++ b/src/krylov_solvers.jl @@ -55,9 +55,26 @@ const KRYLOV_SOLVERS = Dict( KrylovConstructor(vm; vm_empty=vm) KrylovConstructor(vm, vn; vm_empty=vm, vn_empty=vn) -A type containing vectors used to determine the ones that need to be allocated in the Krylov solvers using `Krylov.ksimilar`. -By default, `Krylov.ksimilar(v)` simply calls `similar(v)`. -If a different behavior is required, the function `Krylov.ksimilar` can be specialized with additional methods for your specific type. +Krylov methods require a workspace containing vectors of length `m` and `n` to solve linear problems of size `m × n`. +The `KrylovConstructor` facilitates the allocation of these vectors using `similar`. + +For square problems (`m == n`), use the first constructor with a single vector `vm`. +For rectangular problems (`m ≠ n`), use the second constructor with `vm` and `vn`. + +#### Input arguments + +* `vm`: a vector of length `m`; +* `vn`: a vector of length `n`. + +#### Keyword arguments + +- `vm_empty`: an empty vector that may be replaced with a vector of length `m`; +- `vn_empty`: an empty vector that may be replaced with a vector of length `n`. + +#### Note + +Empty vectors `vm_empty` and `vn_empty` reduce storage requirements when features such as warm-start or preconditioners are unused. +These empty vectors will be replaced within a [`KrylovSolver`](@ref) only if required, such as when preconditioners are provided. """ struct KrylovConstructor{S} vm::S @@ -110,14 +127,14 @@ function MinresSolver(kc::KrylovConstructor; window :: Int=5) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - r1 = ksimilar(kc.vn) - r2 = ksimilar(kc.vn) - w1 = ksimilar(kc.vn) - w2 = ksimilar(kc.vn) - y = ksimilar(kc.vn) - v = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + r1 = similar(kc.vn) + r2 = similar(kc.vn) + w1 = similar(kc.vn) + w2 = similar(kc.vn) + y = similar(kc.vn) + v = similar(kc.vn_empty) err_vec = zeros(T, window) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = MinresSolver{T,FC,S}(m, n, Δx, x, r1, r2, w1, w2, y, v, err_vec, false, stats) @@ -180,15 +197,15 @@ function MinaresSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - vₖ = ksimilar(kc.vn) - vₖ₊₁ = ksimilar(kc.vn) - x = ksimilar(kc.vn) - wₖ₋₂ = ksimilar(kc.vn) - wₖ₋₁ = ksimilar(kc.vn) - dₖ₋₂ = ksimilar(kc.vn) - dₖ₋₁ = ksimilar(kc.vn) - q = ksimilar(kc.vn) + Δx = similar(kc.vn_empty) + vₖ = similar(kc.vn) + vₖ₊₁ = similar(kc.vn) + x = similar(kc.vn) + wₖ₋₂ = similar(kc.vn) + wₖ₋₁ = similar(kc.vn) + dₖ₋₂ = similar(kc.vn) + dₖ₋₁ = similar(kc.vn) + q = similar(kc.vn) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = MinaresSolver{T,FC,S}(m, n, Δx, vₖ, vₖ₊₁, x, wₖ₋₂, wₖ₋₁, dₖ₋₂, dₖ₋₁, q, false, stats) return solver @@ -247,12 +264,12 @@ function CgSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - r = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Ap = ksimilar(kc.vn) - z = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + r = similar(kc.vn) + p = similar(kc.vn) + Ap = similar(kc.vn) + z = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CgSolver{T,FC,S}(m, n, Δx, x, r, p, Ap, z, false, stats) return solver @@ -309,13 +326,13 @@ function CrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - r = ksimilar(kc.vn) - p = ksimilar(kc.vn) - q = ksimilar(kc.vn) - Ar = ksimilar(kc.vn) - Mq = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + r = similar(kc.vn) + p = similar(kc.vn) + q = similar(kc.vn) + Ar = similar(kc.vn) + Mq = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CrSolver{T,FC,S}(m, n, Δx, x, r, p, q, Ar, Mq, false, stats) return solver @@ -375,15 +392,15 @@ function CarSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - r = ksimilar(kc.vn) - p = ksimilar(kc.vn) - s = ksimilar(kc.vn) - q = ksimilar(kc.vn) - t = ksimilar(kc.vn) - u = ksimilar(kc.vn) - Mu = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + r = similar(kc.vn) + p = similar(kc.vn) + s = similar(kc.vn) + q = similar(kc.vn) + t = similar(kc.vn) + u = similar(kc.vn) + Mu = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CarSolver{T,FC,S}(m, n, Δx, x, r, p, s, q, t, u, Mu, false, stats) return solver @@ -446,13 +463,13 @@ function SymmlqSolver(kc::KrylovConstructor; window :: Int=5) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - Mvold = ksimilar(kc.vn) - Mv = ksimilar(kc.vn) - Mv_next = ksimilar(kc.vn) - w̅ = ksimilar(kc.vn) - v = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + Mvold = similar(kc.vn) + Mv = similar(kc.vn) + Mv_next = similar(kc.vn) + w̅ = similar(kc.vn) + v = similar(kc.vn_empty) clist = zeros(T, window) zlist = zeros(T, window) sprod = ones(T, window) @@ -516,13 +533,13 @@ function CgLanczosSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - Mv = ksimilar(kc.vn) - Mv_prev = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Mv_next = ksimilar(kc.vn) - v = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + Mv = similar(kc.vn) + Mv_prev = similar(kc.vn) + p = similar(kc.vn) + Mv_next = similar(kc.vn) + v = similar(kc.vn_empty) stats = LanczosStats(0, false, T[], false, T(NaN), T(NaN), 0.0, "unknown") solver = CgLanczosSolver{T,FC,S}(m, n, Δx, x, Mv, Mv_prev, p, Mv_next, v, false, stats) return solver @@ -586,12 +603,12 @@ function CgLanczosShiftSolver(kc::KrylovConstructor, nshifts) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Mv = ksimilar(kc.vn) - Mv_prev = ksimilar(kc.vn) - Mv_next = ksimilar(kc.vn) - v = ksimilar(kc.vn_empty) - x = S[ksimilar(kc.vn) for i = 1 : nshifts] - p = S[ksimilar(kc.vn) for i = 1 : nshifts] + Mv = similar(kc.vn) + Mv_prev = similar(kc.vn) + Mv_next = similar(kc.vn) + v = similar(kc.vn_empty) + x = S[similar(kc.vn) for i = 1 : nshifts] + p = S[similar(kc.vn) for i = 1 : nshifts] σ = Vector{T}(undef, nshifts) δhat = Vector{T}(undef, nshifts) ω = Vector{T}(undef, nshifts) @@ -665,14 +682,14 @@ function MinresQlpSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - wₖ₋₁ = ksimilar(kc.vn) - wₖ = ksimilar(kc.vn) - M⁻¹vₖ₋₁ = ksimilar(kc.vn) - M⁻¹vₖ = ksimilar(kc.vn) - x = ksimilar(kc.vn) - p = ksimilar(kc.vn) - vₖ = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + wₖ₋₁ = similar(kc.vn) + wₖ = similar(kc.vn) + M⁻¹vₖ₋₁ = similar(kc.vn) + M⁻¹vₖ = similar(kc.vn) + x = similar(kc.vn) + p = similar(kc.vn) + vₖ = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = MinresQlpSolver{T,FC,S}(m, n, Δx, wₖ₋₁, wₖ, M⁻¹vₖ₋₁, M⁻¹vₖ, x, p, vₖ, false, stats) return solver @@ -737,13 +754,13 @@ function DqgmresSolver(kc::KrylovConstructor, memory = 20) m = length(kc.vm) n = length(kc.vn) memory = min(m, memory) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - t = ksimilar(kc.vn) - z = ksimilar(kc.vn_empty) - w = ksimilar(kc.vn_empty) - P = S[ksimilar(kc.vn) for i = 1 : memory] - V = S[ksimilar(kc.vn) for i = 1 : memory] + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + t = similar(kc.vn) + z = similar(kc.vn_empty) + w = similar(kc.vn_empty) + P = S[similar(kc.vn) for i = 1 : memory] + V = S[similar(kc.vn) for i = 1 : memory] c = Vector{T}(undef, memory) s = Vector{FC}(undef, memory) H = Vector{FC}(undef, memory+1) @@ -813,13 +830,13 @@ function DiomSolver(kc::KrylovConstructor, memory = 20) m = length(kc.vm) n = length(kc.vn) memory = min(m, memory) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - t = ksimilar(kc.vn) - z = ksimilar(kc.vn_empty) - w = ksimilar(kc.vn_empty) - P = S[ksimilar(kc.vn) for i = 1 : memory-1] - V = S[ksimilar(kc.vn) for i = 1 : memory] + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + t = similar(kc.vn) + z = similar(kc.vn_empty) + w = similar(kc.vn_empty) + P = S[similar(kc.vn) for i = 1 : memory-1] + V = S[similar(kc.vn) for i = 1 : memory] L = Vector{FC}(undef, memory-1) H = Vector{FC}(undef, memory) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") @@ -884,15 +901,15 @@ function UsymlqSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - uₖ₋₁ = ksimilar(kc.vn) - uₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - d̅ = ksimilar(kc.vn) - vₖ₋₁ = ksimilar(kc.vm) - vₖ = ksimilar(kc.vm) - q = ksimilar(kc.vm) + uₖ₋₁ = similar(kc.vn) + uₖ = similar(kc.vn) + p = similar(kc.vn) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + d̅ = similar(kc.vn) + vₖ₋₁ = similar(kc.vm) + vₖ = similar(kc.vm) + q = similar(kc.vm) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = UsymlqSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, p, Δx, x, d̅, vₖ₋₁, vₖ, q, false, stats) return solver @@ -955,16 +972,16 @@ function UsymqrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - vₖ₋₁ = ksimilar(kc.vm) - vₖ = ksimilar(kc.vm) - q = ksimilar(kc.vm) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - wₖ₋₂ = ksimilar(kc.vn) - wₖ₋₁ = ksimilar(kc.vn) - uₖ₋₁ = ksimilar(kc.vn) - uₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) + vₖ₋₁ = similar(kc.vm) + vₖ = similar(kc.vm) + q = similar(kc.vm) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + wₖ₋₂ = similar(kc.vn) + wₖ₋₁ = similar(kc.vn) + uₖ₋₁ = similar(kc.vn) + uₖ = similar(kc.vn) + p = similar(kc.vn) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = UsymqrSolver{T,FC,S}(m, n, vₖ₋₁, vₖ, q, Δx, x, wₖ₋₂, wₖ₋₁, uₖ₋₁, uₖ, p, false, stats) return solver @@ -1034,22 +1051,22 @@ function TricgSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - y = ksimilar(kc.vn) - N⁻¹uₖ₋₁ = ksimilar(kc.vn) - N⁻¹uₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) - gy₂ₖ₋₁ = ksimilar(kc.vn) - gy₂ₖ = ksimilar(kc.vn) - x = ksimilar(kc.vm) - M⁻¹vₖ₋₁ = ksimilar(kc.vm) - M⁻¹vₖ = ksimilar(kc.vm) - q = ksimilar(kc.vm) - gx₂ₖ₋₁ = ksimilar(kc.vm) - gx₂ₖ = ksimilar(kc.vm) - Δx = ksimilar(kc.vm_empty) - Δy = ksimilar(kc.vn_empty) - uₖ = ksimilar(kc.vn_empty) - vₖ = ksimilar(kc.vm_empty) + y = similar(kc.vn) + N⁻¹uₖ₋₁ = similar(kc.vn) + N⁻¹uₖ = similar(kc.vn) + p = similar(kc.vn) + gy₂ₖ₋₁ = similar(kc.vn) + gy₂ₖ = similar(kc.vn) + x = similar(kc.vm) + M⁻¹vₖ₋₁ = similar(kc.vm) + M⁻¹vₖ = similar(kc.vm) + q = similar(kc.vm) + gx₂ₖ₋₁ = similar(kc.vm) + gx₂ₖ = similar(kc.vm) + Δx = similar(kc.vm_empty) + Δy = similar(kc.vn_empty) + uₖ = similar(kc.vn_empty) + vₖ = similar(kc.vm_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = TricgSolver{T,FC,S}(m, n, y, N⁻¹uₖ₋₁, N⁻¹uₖ, p, gy₂ₖ₋₁, gy₂ₖ, x, M⁻¹vₖ₋₁, M⁻¹vₖ, q, gx₂ₖ₋₁, gx₂ₖ, Δx, Δy, uₖ, vₖ, false, stats) return solver @@ -1129,26 +1146,26 @@ function TrimrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - y = ksimilar(kc.vn) - N⁻¹uₖ₋₁ = ksimilar(kc.vn) - N⁻¹uₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) - gy₂ₖ₋₃ = ksimilar(kc.vn) - gy₂ₖ₋₂ = ksimilar(kc.vn) - gy₂ₖ₋₁ = ksimilar(kc.vn) - gy₂ₖ = ksimilar(kc.vn) - x = ksimilar(kc.vm) - M⁻¹vₖ₋₁ = ksimilar(kc.vm) - M⁻¹vₖ = ksimilar(kc.vm) - q = ksimilar(kc.vm) - gx₂ₖ₋₃ = ksimilar(kc.vm) - gx₂ₖ₋₂ = ksimilar(kc.vm) - gx₂ₖ₋₁ = ksimilar(kc.vm) - gx₂ₖ = ksimilar(kc.vm) - Δx = ksimilar(kc.vm_empty) - Δy = ksimilar(kc.vn_empty) - uₖ = ksimilar(kc.vn_empty) - vₖ = ksimilar(kc.vm_empty) + y = similar(kc.vn) + N⁻¹uₖ₋₁ = similar(kc.vn) + N⁻¹uₖ = similar(kc.vn) + p = similar(kc.vn) + gy₂ₖ₋₃ = similar(kc.vn) + gy₂ₖ₋₂ = similar(kc.vn) + gy₂ₖ₋₁ = similar(kc.vn) + gy₂ₖ = similar(kc.vn) + x = similar(kc.vm) + M⁻¹vₖ₋₁ = similar(kc.vm) + M⁻¹vₖ = similar(kc.vm) + q = similar(kc.vm) + gx₂ₖ₋₃ = similar(kc.vm) + gx₂ₖ₋₂ = similar(kc.vm) + gx₂ₖ₋₁ = similar(kc.vm) + gx₂ₖ = similar(kc.vm) + Δx = similar(kc.vm_empty) + Δy = similar(kc.vn_empty) + uₖ = similar(kc.vn_empty) + vₖ = similar(kc.vm_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = TrimrSolver{T,FC,S}(m, n, y, N⁻¹uₖ₋₁, N⁻¹uₖ, p, gy₂ₖ₋₃, gy₂ₖ₋₂, gy₂ₖ₋₁, gy₂ₖ, x, M⁻¹vₖ₋₁, M⁻¹vₖ, q, gx₂ₖ₋₃, gx₂ₖ₋₂, gx₂ₖ₋₁, gx₂ₖ, Δx, Δy, uₖ, vₖ, false, stats) return solver @@ -1225,19 +1242,19 @@ function TrilqrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - uₖ₋₁ = ksimilar(kc.vn) - uₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) - d̅ = ksimilar(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - vₖ₋₁ = ksimilar(kc.vm) - vₖ = ksimilar(kc.vm) - q = ksimilar(kc.vm) - Δy = ksimilar(kc.vm_empty) - y = ksimilar(kc.vm) - wₖ₋₃ = ksimilar(kc.vm) - wₖ₋₂ = ksimilar(kc.vm) + uₖ₋₁ = similar(kc.vn) + uₖ = similar(kc.vn) + p = similar(kc.vn) + d̅ = similar(kc.vn) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + vₖ₋₁ = similar(kc.vm) + vₖ = similar(kc.vm) + q = similar(kc.vm) + Δy = similar(kc.vm_empty) + y = similar(kc.vm) + wₖ₋₃ = similar(kc.vm) + wₖ₋₂ = similar(kc.vm) stats = AdjointStats(0, false, false, T[], T[], 0.0, "unknown") solver = TrilqrSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, p, d̅, Δx, x, vₖ₋₁, vₖ, q, Δy, y, wₖ₋₃, wₖ₋₂, false, stats) return solver @@ -1303,15 +1320,15 @@ function CgsSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - r = ksimilar(kc.vn) - u = ksimilar(kc.vn) - p = ksimilar(kc.vn) - q = ksimilar(kc.vn) - ts = ksimilar(kc.vn) - yz = ksimilar(kc.vn_empty) - vw = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + r = similar(kc.vn) + u = similar(kc.vn) + p = similar(kc.vn) + q = similar(kc.vn) + ts = similar(kc.vn) + yz = similar(kc.vn_empty) + vw = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CgsSolver{T,FC,S}(m, n, Δx, x, r, u, p, q, ts, yz, vw, false, stats) return solver @@ -1373,15 +1390,15 @@ function BicgstabSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - r = ksimilar(kc.vn) - p = ksimilar(kc.vn) - v = ksimilar(kc.vn) - s = ksimilar(kc.vn) - qd = ksimilar(kc.vn) - yz = ksimilar(kc.vn_empty) - t = ksimilar(kc.vn_empty) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + r = similar(kc.vn) + p = similar(kc.vn) + v = similar(kc.vn) + s = similar(kc.vn) + qd = similar(kc.vn) + yz = similar(kc.vn_empty) + t = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = BicgstabSolver{T,FC,S}(m, n, Δx, x, r, p, v, s, qd, yz, t, false, stats) return solver @@ -1445,17 +1462,17 @@ function BilqSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - uₖ₋₁ = ksimilar(kc.vn) - uₖ = ksimilar(kc.vn) - q = ksimilar(kc.vn) - vₖ₋₁ = ksimilar(kc.vn) - vₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - d̅ = ksimilar(kc.vn) - t = ksimilar(kc.vn_empty) - s = ksimilar(kc.vn_empty) + uₖ₋₁ = similar(kc.vn) + uₖ = similar(kc.vn) + q = similar(kc.vn) + vₖ₋₁ = similar(kc.vn) + vₖ = similar(kc.vn) + p = similar(kc.vn) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + d̅ = similar(kc.vn) + t = similar(kc.vn_empty) + s = similar(kc.vn_empty) 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̅, t, s, false, stats) return solver @@ -1522,18 +1539,18 @@ function QmrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - uₖ₋₁ = ksimilar(kc.vn) - uₖ = ksimilar(kc.vn) - q = ksimilar(kc.vn) - vₖ₋₁ = ksimilar(kc.vn) - vₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - wₖ₋₂ = ksimilar(kc.vn) - wₖ₋₁ = ksimilar(kc.vn) - t = ksimilar(kc.vn_empty) - s = ksimilar(kc.vn_empty) + uₖ₋₁ = similar(kc.vn) + uₖ = similar(kc.vn) + q = similar(kc.vn) + vₖ₋₁ = similar(kc.vn) + vₖ = similar(kc.vn) + p = similar(kc.vn) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + wₖ₋₂ = similar(kc.vn) + wₖ₋₁ = similar(kc.vn) + t = similar(kc.vn_empty) + s = similar(kc.vn_empty) 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ₖ₋₁, t, s, false, stats) return solver @@ -1602,19 +1619,19 @@ function BilqrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - uₖ₋₁ = ksimilar(kc.vn) - uₖ = ksimilar(kc.vn) - q = ksimilar(kc.vn) - vₖ₋₁ = ksimilar(kc.vn) - vₖ = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - Δy = ksimilar(kc.vn_empty) - y = ksimilar(kc.vn) - d̅ = ksimilar(kc.vn) - wₖ₋₃ = ksimilar(kc.vn) - wₖ₋₂ = ksimilar(kc.vn) + uₖ₋₁ = similar(kc.vn) + uₖ = similar(kc.vn) + q = similar(kc.vn) + vₖ₋₁ = similar(kc.vn) + vₖ = similar(kc.vn) + p = similar(kc.vn) + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + Δy = similar(kc.vn_empty) + y = similar(kc.vn) + d̅ = similar(kc.vn) + wₖ₋₃ = similar(kc.vn) + wₖ₋₂ = similar(kc.vn) stats = AdjointStats(0, false, false, T[], T[], 0.0, "unknown") solver = BilqrSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, Δy, y, d̅, wₖ₋₃, wₖ₋₂, false, stats) return solver @@ -1676,12 +1693,12 @@ function CglsSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - p = ksimilar(kc.vn) - s = ksimilar(kc.vn) - r = ksimilar(kc.vm) - q = ksimilar(kc.vm) - Mr = ksimilar(kc.vm_empty) + x = similar(kc.vn) + p = similar(kc.vn) + s = similar(kc.vn) + r = similar(kc.vm) + q = similar(kc.vm) + Mr = similar(kc.vm_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CglsSolver{T,FC,S}(m, n, x, p, s, r, q, Mr, stats) return solver @@ -1745,13 +1762,13 @@ function CglsLanczosShiftSolver(kc::KrylovConstructor, nshifts) T = real(FC) m = length(kc.vm) n = length(kc.vn) - Mv = ksimilar(kc.vn) - u_prev = ksimilar(kc.vm) - u_next = ksimilar(kc.vm) - u = ksimilar(kc.vm) - v = ksimilar(kc.vn_empty) - x = S[ksimilar(kc.vn) for i = 1 : nshifts] - p = S[ksimilar(kc.vn) for i = 1 : nshifts] + Mv = similar(kc.vn) + u_prev = similar(kc.vm) + u_next = similar(kc.vm) + u = similar(kc.vm) + v = similar(kc.vn_empty) + x = S[similar(kc.vn) for i = 1 : nshifts] + p = S[similar(kc.vn) for i = 1 : nshifts] σ = Vector{T}(undef, nshifts) δhat = Vector{T}(undef, nshifts) ω = Vector{T}(undef, nshifts) @@ -1825,14 +1842,14 @@ function CrlsSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Ar = ksimilar(kc.vn) - q = ksimilar(kc.vn) - r = ksimilar(kc.vm) - Ap = ksimilar(kc.vm) - s = ksimilar(kc.vm) - Ms = ksimilar(kc.vm_empty) + x = similar(kc.vn) + p = similar(kc.vn) + Ar = similar(kc.vn) + q = similar(kc.vn) + r = similar(kc.vm) + Ap = similar(kc.vm) + s = similar(kc.vm) + Ms = similar(kc.vm_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CrlsSolver{T,FC,S}(m, n, x, p, Ar, q, r, Ap, s, Ms, stats) return solver @@ -1890,13 +1907,13 @@ function CgneSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Aᴴz = ksimilar(kc.vn) - r = ksimilar(kc.vm) - q = ksimilar(kc.vm) - s = ksimilar(kc.vm_empty) - z = ksimilar(kc.vm_empty) + x = similar(kc.vn) + p = similar(kc.vn) + Aᴴz = similar(kc.vn) + r = similar(kc.vm) + q = similar(kc.vm) + s = similar(kc.vm_empty) + z = similar(kc.vm_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CgneSolver{T,FC,S}(m, n, x, p, Aᴴz, r, q, s, z, stats) return solver @@ -1953,13 +1970,13 @@ function CrmrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - p = ksimilar(kc.vn) - Aᴴr = ksimilar(kc.vn) - r = ksimilar(kc.vm) - q = ksimilar(kc.vm) - Nq = ksimilar(kc.vm_empty) - s = ksimilar(kc.vm_empty) + x = similar(kc.vn) + p = similar(kc.vn) + Aᴴr = similar(kc.vn) + r = similar(kc.vm) + q = similar(kc.vm) + Nq = similar(kc.vm_empty) + s = similar(kc.vm_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CrmrSolver{T,FC,S}(m, n, x, p, Aᴴr, r, q, Nq, s, stats) return solver @@ -2018,14 +2035,14 @@ function LslqSolver(kc::KrylovConstructor; window :: Int=5) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - Nv = ksimilar(kc.vn) - Aᴴu = ksimilar(kc.vn) - w̄ = ksimilar(kc.vn) - Mu = ksimilar(kc.vm) - Av = ksimilar(kc.vm) - u = ksimilar(kc.vm_empty) - v = ksimilar(kc.vn_empty) + x = similar(kc.vn) + Nv = similar(kc.vn) + Aᴴu = similar(kc.vn) + w̄ = similar(kc.vn) + Mu = similar(kc.vm) + Av = similar(kc.vm) + u = similar(kc.vm_empty) + v = similar(kc.vn_empty) err_vec = zeros(T, window) stats = LSLQStats(0, false, false, T[], T[], T[], false, T[], T[], 0.0, "unknown") solver = LslqSolver{T,FC,S}(m, n, x, Nv, Aᴴu, w̄, Mu, Av, u, v, err_vec, stats) @@ -2087,14 +2104,14 @@ function LsqrSolver(kc::KrylovConstructor; window :: Int=5) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - Nv = ksimilar(kc.vn) - Aᴴu = ksimilar(kc.vn) - w = ksimilar(kc.vn) - Mu = ksimilar(kc.vm) - Av = ksimilar(kc.vm) - u = ksimilar(kc.vm_empty) - v = ksimilar(kc.vn_empty) + x = similar(kc.vn) + Nv = similar(kc.vn) + Aᴴu = similar(kc.vn) + w = similar(kc.vn) + Mu = similar(kc.vm) + Av = similar(kc.vm) + u = similar(kc.vm_empty) + v = similar(kc.vn_empty) err_vec = zeros(T, window) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = LsqrSolver{T,FC,S}(m, n, x, Nv, Aᴴu, w, Mu, Av, u, v, err_vec, stats) @@ -2157,15 +2174,15 @@ function LsmrSolver(kc::KrylovConstructor; window :: Int=5) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - Nv = ksimilar(kc.vn) - Aᴴu = ksimilar(kc.vn) - h = ksimilar(kc.vn) - hbar = ksimilar(kc.vn) - Mu = ksimilar(kc.vm) - Av = ksimilar(kc.vm) - u = ksimilar(kc.vm_empty) - v = ksimilar(kc.vn_empty) + x = similar(kc.vn) + Nv = similar(kc.vn) + Aᴴu = similar(kc.vn) + h = similar(kc.vn) + hbar = similar(kc.vn) + Mu = similar(kc.vm) + Av = similar(kc.vm) + u = similar(kc.vm_empty) + v = similar(kc.vn_empty) err_vec = zeros(T, window) stats = LsmrStats(0, false, false, T[], T[], zero(T), zero(T), zero(T), zero(T), zero(T), 0.0, "unknown") solver = LsmrSolver{T,FC,S}(m, n, x, Nv, Aᴴu, h, hbar, Mu, Av, u, v, err_vec, stats) @@ -2229,16 +2246,16 @@ function LnlqSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - Nv = ksimilar(kc.vn) - Aᴴu = ksimilar(kc.vn) - y = ksimilar(kc.vm) - w̄ = ksimilar(kc.vm) - Mu = ksimilar(kc.vm) - Av = ksimilar(kc.vm) - u = ksimilar(kc.vm_empty) - v = ksimilar(kc.vn_empty) - q = ksimilar(kc.vn_empty) + x = similar(kc.vn) + Nv = similar(kc.vn) + Aᴴu = similar(kc.vn) + y = similar(kc.vm) + w̄ = similar(kc.vm) + Mu = similar(kc.vm) + Av = similar(kc.vm) + u = similar(kc.vm_empty) + v = similar(kc.vn_empty) + q = similar(kc.vn_empty) stats = LNLQStats(0, false, T[], false, T[], T[], 0.0, "unknown") solver = LnlqSolver{T,FC,S}(m, n, x, Nv, Aᴴu, y, w̄, Mu, Av, u, v, q, stats) return solver @@ -2301,16 +2318,16 @@ function CraigSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - Nv = ksimilar(kc.vn) - Aᴴu = ksimilar(kc.vn) - y = ksimilar(kc.vm) - w = ksimilar(kc.vm) - Mu = ksimilar(kc.vm) - Av = ksimilar(kc.vm) - u = ksimilar(kc.vm_empty) - v = ksimilar(kc.vn_empty) - w2 = ksimilar(kc.vn_empty) + x = similar(kc.vn) + Nv = similar(kc.vn) + Aᴴu = similar(kc.vn) + y = similar(kc.vm) + w = similar(kc.vm) + Mu = similar(kc.vm) + Av = similar(kc.vm) + u = similar(kc.vm_empty) + v = similar(kc.vn_empty) + w2 = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CraigSolver{T,FC,S}(m, n, x, Nv, Aᴴu, y, w, Mu, Av, u, v, w2, stats) return solver @@ -2375,18 +2392,18 @@ function CraigmrSolver(kc::KrylovConstructor) T = real(FC) m = length(kc.vm) n = length(kc.vn) - x = ksimilar(kc.vn) - Nv = ksimilar(kc.vn) - Aᴴu = ksimilar(kc.vn) - d = ksimilar(kc.vn) - y = ksimilar(kc.vm) - Mu = ksimilar(kc.vm) - w = ksimilar(kc.vm) - wbar = ksimilar(kc.vm) - Av = ksimilar(kc.vm) - u = ksimilar(kc.vm_empty) - v = ksimilar(kc.vn_empty) - q = ksimilar(kc.vn_empty) + x = similar(kc.vn) + Nv = similar(kc.vn) + Aᴴu = similar(kc.vn) + d = similar(kc.vn) + y = similar(kc.vm) + Mu = similar(kc.vm) + w = similar(kc.vm) + wbar = similar(kc.vm) + Av = similar(kc.vm) + u = similar(kc.vm_empty) + v = similar(kc.vn_empty) + q = similar(kc.vn_empty) stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown") solver = CraigmrSolver{T,FC,S}(m, n, x, Nv, Aᴴu, d, y, Mu, w, wbar, Av, u, v, q, stats) return solver @@ -2456,12 +2473,12 @@ function GmresSolver(kc::KrylovConstructor, memory = 20) m = length(kc.vm) n = length(kc.vn) memory = min(m, memory) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - w = ksimilar(kc.vn) - p = ksimilar(kc.vn_empty) - q = ksimilar(kc.vn_empty) - V = S[ksimilar(kc.vn) for i = 1 : memory] + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + w = similar(kc.vn) + p = similar(kc.vn_empty) + q = similar(kc.vn_empty) + V = S[similar(kc.vn) for i = 1 : memory] c = Vector{T}(undef, memory) s = Vector{FC}(undef, memory) z = Vector{FC}(undef, memory) @@ -2534,12 +2551,12 @@ function FgmresSolver(kc::KrylovConstructor, memory = 20) m = length(kc.vm) n = length(kc.vn) memory = min(m, memory) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - w = ksimilar(kc.vn) - q = ksimilar(kc.vn_empty) - V = S[ksimilar(kc.vn) for i = 1 : memory] - Z = S[ksimilar(kc.vn) for i = 1 : memory] + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + w = similar(kc.vn) + q = similar(kc.vn_empty) + V = S[similar(kc.vn) for i = 1 : memory] + Z = S[similar(kc.vn) for i = 1 : memory] c = Vector{T}(undef, memory) s = Vector{FC}(undef, memory) z = Vector{FC}(undef, memory) @@ -2610,12 +2627,12 @@ function FomSolver(kc::KrylovConstructor, memory = 20) m = length(kc.vm) n = length(kc.vn) memory = min(m, memory) - Δx = ksimilar(kc.vn_empty) - x = ksimilar(kc.vn) - w = ksimilar(kc.vn) - p = ksimilar(kc.vn_empty) - q = ksimilar(kc.vn_empty) - V = S[ksimilar(kc.vn) for i = 1 : memory] + Δx = similar(kc.vn_empty) + x = similar(kc.vn) + w = similar(kc.vn) + p = similar(kc.vn_empty) + q = similar(kc.vn_empty) + V = S[similar(kc.vn) for i = 1 : memory] l = Vector{FC}(undef, memory) z = Vector{FC}(undef, memory) U = Vector{FC}(undef, div(memory * (memory+1), 2)) @@ -2691,18 +2708,18 @@ function GpmrSolver(kc::KrylovConstructor, memory = 20) m = length(kc.vm) n = length(kc.vn) memory = min(n + m, memory) - wA = ksimilar(kc.vn_empty) - wB = ksimilar(kc.vm_empty) - dA = ksimilar(kc.vm) - dB = ksimilar(kc.vn) - Δx = ksimilar(kc.vm_empty) - Δy = ksimilar(kc.vn_empty) - x = ksimilar(kc.vm) - y = ksimilar(kc.vn) - q = ksimilar(kc.vm_empty) - p = ksimilar(kc.vn_empty) - V = S[ksimilar(kc.vm) for i = 1 : memory] - U = S[ksimilar(kc.vn) for i = 1 : memory] + wA = similar(kc.vn_empty) + wB = similar(kc.vm_empty) + dA = similar(kc.vm) + dB = similar(kc.vn) + Δx = similar(kc.vm_empty) + Δy = similar(kc.vn_empty) + x = similar(kc.vm) + y = similar(kc.vn) + q = similar(kc.vm_empty) + p = similar(kc.vn_empty) + V = S[similar(kc.vm) for i = 1 : memory] + U = S[similar(kc.vn) for i = 1 : memory] gs = Vector{FC}(undef, 4 * memory) gc = Vector{T}(undef, 4 * memory) zt = Vector{FC}(undef, 2 * memory) diff --git a/src/krylov_utils.jl b/src/krylov_utils.jl index 7f995d860..dee289ca7 100644 --- a/src/krylov_utils.jl +++ b/src/krylov_utils.jl @@ -296,11 +296,9 @@ Create a vector of storage type `S` of length `n` only composed of one. """ kones(S, n) = fill!(S(undef, n), one(eltype(S))) -kisempty(v) = isempty(v) -ksimilar(v) = similar(v) -allocate_if(bool, solver, v, S, u) = bool && kisempty(solver.:($v)::S) && (solver.:($v)::S = ksimilar(u)) -# allocate_if(bool, solver, v, S, n::Int) = bool && kisempty(solver.:($v)::S) && (solver.:($v)::S = S(undef, n)) -allocate_if(bool, solver, v, S, m::Int, n::Int) = bool && kisempty(solver.:($v)::S) && (solver.:($v)::S = S(undef, m, n)) +allocate_if(bool, solver, v, S, u) = bool && isempty(solver.:($v)::S) && (solver.:($v)::S = similar(u)) +# allocate_if(bool, solver, v, S, n::Int) = bool && isempty(solver.:($v)::S) && (solver.:($v)::S = S(undef, n)) +allocate_if(bool, solver, v, S, m::Int, n::Int) = bool && isempty(solver.:($v)::S) && (solver.:($v)::S = S(undef, m, n)) kdisplay(iter, verbose) = (verbose > 0) && (mod(iter, verbose) == 0)