From c45a492f4375e8e5648e58c2fcd2fa13b9e6ebbd Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Tue, 27 Sep 2022 20:45:49 -0400 Subject: [PATCH] Move callback_utils.jl file --- src/Krylov.jl | 2 - src/callback_utils.jl | 50 -------------- test/callback_utils.jl | 152 +++++++++++++++++++++++++++++++++++++++++ test/test_utils.jl | 108 +---------------------------- 4 files changed, 153 insertions(+), 159 deletions(-) delete mode 100644 src/callback_utils.jl create mode 100644 test/callback_utils.jl diff --git a/src/Krylov.jl b/src/Krylov.jl index 7c480896f..e3903e124 100644 --- a/src/Krylov.jl +++ b/src/Krylov.jl @@ -50,6 +50,4 @@ include("lnlq.jl") include("craig.jl") include("craigmr.jl") -include("callback_utils.jl") - end diff --git a/src/callback_utils.jl b/src/callback_utils.jl deleted file mode 100644 index eac362e5d..000000000 --- a/src/callback_utils.jl +++ /dev/null @@ -1,50 +0,0 @@ -export StorageGetxRestartedGmres - -export get_x_restarted_gmres! - -mutable struct StorageGetxRestartedGmres{S} - x::S - y::S - p::S -end -StorageGetxRestartedGmres(solver::GmresSolver; N = I) = - StorageGetxRestartedGmres(similar(solver.x), similar(solver.z), (N === I) ? similar(solver.p) : similar(solver.x)) - -function get_x_restarted_gmres!(solver::GmresSolver{T,FC,S}, A, - stor::StorageGetxRestartedGmres{S}, N) where {T,FC,S} - NisI = (N === I) - x2, y2, p2 = stor.x, stor.y, stor.p - n = size(A, 2) - # Compute yₖ by solving Rₖyₖ = zₖ with backward substitution. - nr = sum(1:solver.inner_iter) - y = solver.z # yᵢ = zᵢ - y2 .= y - R = solver.R - V = solver.V - x2 .= solver.Δx - for i = solver.inner_iter : -1 : 1 - pos = nr + i - solver.inner_iter # position of rᵢ.ₖ - for j = solver.inner_iter : -1 : i+1 - y2[i] = y2[i] - R[pos] * y2[j] # yᵢ ← yᵢ - rᵢⱼyⱼ - pos = pos - j + 1 # position of rᵢ.ⱼ₋₁ - end - # Rₖ can be singular if the system is inconsistent - if abs(R[pos]) ≤ eps(T)^(3/4) - y2[i] = zero(FC) - inconsistent = true - else - y2[i] = y2[i] / R[pos] # yᵢ ← yᵢ / rᵢᵢ - end - end - - # Form xₖ = N⁻¹Vₖyₖ - for i = 1 : solver.inner_iter - @kaxpy!(n, y2[i], V[i], x2) - end - if !NisI - p2 .= solver.p - p2 .= x2 - mul!(x2, N, p2) - end - x2 .+= solver.x -end diff --git a/test/callback_utils.jl b/test/callback_utils.jl new file mode 100644 index 000000000..c5993c2a3 --- /dev/null +++ b/test/callback_utils.jl @@ -0,0 +1,152 @@ +mutable struct StorageGetxRestartedGmres{S} + x::S + y::S + p::S +end +StorageGetxRestartedGmres(solver::GmresSolver; N = I) = + StorageGetxRestartedGmres(similar(solver.x), similar(solver.z), (N === I) ? similar(solver.p) : similar(solver.x)) + +function get_x_restarted_gmres!(solver::GmresSolver{T,FC,S}, A, + stor::StorageGetxRestartedGmres{S}, N) where {T,FC,S} + NisI = (N === I) + x2, y2, p2 = stor.x, stor.y, stor.p + n = size(A, 2) + # Compute yₖ by solving Rₖyₖ = zₖ with backward substitution. + nr = sum(1:solver.inner_iter) + y = solver.z # yᵢ = zᵢ + y2 .= y + R = solver.R + V = solver.V + x2 .= solver.Δx + for i = solver.inner_iter : -1 : 1 + pos = nr + i - solver.inner_iter # position of rᵢ.ₖ + for j = solver.inner_iter : -1 : i+1 + y2[i] = y2[i] - R[pos] * y2[j] # yᵢ ← yᵢ - rᵢⱼyⱼ + pos = pos - j + 1 # position of rᵢ.ⱼ₋₁ + end + # Rₖ can be singular if the system is inconsistent + if abs(R[pos]) ≤ eps(T)^(3/4) + y2[i] = zero(FC) + inconsistent = true + else + y2[i] = y2[i] / R[pos] # yᵢ ← yᵢ / rᵢᵢ + end + end + + # Form xₖ = N⁻¹Vₖyₖ + for i = 1 : solver.inner_iter + Krylov.@kaxpy!(n, y2[i], V[i], x2) + end + if !NisI + p2 .= solver.p + p2 .= x2 + mul!(x2, N, p2) + end + x2 .+= solver.x +end + +mutable struct TestCallbackN2{T, S, M} + A::M + b::S + storage_vec::S + tol::T +end +TestCallbackN2(A, b; tol = 0.1) = TestCallbackN2(A, b, similar(b), tol) + +function (cb_n2::TestCallbackN2)(solver) + mul!(cb_n2.storage_vec, cb_n2.A, solver.x) + cb_n2.storage_vec .-= cb_n2.b + return norm(cb_n2.storage_vec) ≤ cb_n2.tol +end + +mutable struct TestCallbackN2Adjoint{T, S, M} + A::M + b::S + c::S + storage_vec1::S + storage_vec2::S + tol::T +end +TestCallbackN2Adjoint(A, b, c; tol = 0.1) = TestCallbackN2Adjoint(A, b, c, similar(b), similar(c), tol) + +function (cb_n2::TestCallbackN2Adjoint)(solver) + mul!(cb_n2.storage_vec1, cb_n2.A, solver.x) + cb_n2.storage_vec1 .-= cb_n2.b + mul!(cb_n2.storage_vec2, cb_n2.A', solver.y) + cb_n2.storage_vec2 .-= cb_n2.c + return (norm(cb_n2.storage_vec1) ≤ cb_n2.tol && norm(cb_n2.storage_vec2) ≤ cb_n2.tol) +end + +mutable struct TestCallbackN2Shifts{T, S, M} + A::M + b::S + shifts::Vector{T} + tol::T +end +TestCallbackN2Shifts(A, b, shifts; tol = 0.1) = TestCallbackN2Shifts(A, b, shifts, tol) + +function (cb_n2::TestCallbackN2Shifts)(solver) + r = residuals(cb_n2.A, cb_n2.b, cb_n2.shifts, solver.x) + return all(map(norm, r) .≤ cb_n2.tol) +end + +mutable struct TestCallbackN2LS{T, S, M} + A::M + b::S + λ::T + storage_vec1::S + storage_vec2::S + tol::T +end +TestCallbackN2LS(A, b, λ; tol = 0.1) = TestCallbackN2LS(A, b, λ, similar(b), similar(b, size(A, 2)), tol) + +function (cb_n2::TestCallbackN2LS)(solver) + mul!(cb_n2.storage_vec1, cb_n2.A, solver.x) + cb_n2.storage_vec1 .-= cb_n2.b + mul!(cb_n2.storage_vec2, cb_n2.A', cb_n2.storage_vec1) + cb_n2.storage_vec2 .+= cb_n2.λ .* solver.x + return norm(cb_n2.storage_vec2) ≤ cb_n2.tol +end + +mutable struct TestCallbackN2LN{T, S, M} + A::M + b::S + λ::T + storage_vec::S + tol::T +end +TestCallbackN2LN(A, b, λ; tol = 0.1) = TestCallbackN2LN(A, b, λ, similar(b), tol) + +function (cb_n2::TestCallbackN2LN)(solver) + mul!(cb_n2.storage_vec, cb_n2.A, solver.x) + cb_n2.storage_vec .-= cb_n2.b + cb_n2.λ != 0 && (cb_n2.storage_vec .+= sqrt(cb_n2.λ) .* solver.s) + return norm(cb_n2.storage_vec) ≤ cb_n2.tol +end + +mutable struct TestCallbackN2SaddlePts{T, S, M} + A::M + b::S + c::S + storage_vec1::S + storage_vec2::S + tol::T +end +TestCallbackN2SaddlePts(A, b, c; tol = 0.1) = + TestCallbackN2SaddlePts(A, b, c, similar(b), similar(c), tol) + +function (cb_n2::TestCallbackN2SaddlePts)(solver) + mul!(cb_n2.storage_vec1, cb_n2.A, solver.y) + cb_n2.storage_vec1 .+= solver.x .- cb_n2.b + mul!(cb_n2.storage_vec2, cb_n2.A', solver.x) + cb_n2.storage_vec2 .-= solver.y .+ cb_n2.c + return (norm(cb_n2.storage_vec1) ≤ cb_n2.tol && norm(cb_n2.storage_vec2) ≤ cb_n2.tol) +end + +function restarted_gmres_callback_n2(solver::GmresSolver, A, b, stor, N, storage_vec, tol) + get_x_restarted_gmres!(solver, A, stor, N) + x = stor.x + mul!(storage_vec, A, x) + storage_vec .-= b + return (norm(storage_vec) ≤ tol) +end diff --git a/test/test_utils.jl b/test/test_utils.jl index fbfe2e4e0..0ac2e1538 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -1,6 +1,7 @@ include("get_div_grad.jl") include("gen_lsq.jl") include("check_min_norm.jl") +include("callback_utils.jl") # Symmetric and positive definite systems. function symmetric_definite(n :: Int=10; FC=Float64) @@ -363,110 +364,3 @@ function check_reset(stats :: KS) where KS <: Krylov.KrylovStats end end end - -# Test callback -mutable struct TestCallbackN2{T, S, M} - A::M - b::S - storage_vec::S - tol::T -end -TestCallbackN2(A, b; tol = 0.1) = TestCallbackN2(A, b, similar(b), tol) - -function (cb_n2::TestCallbackN2)(solver) - mul!(cb_n2.storage_vec, cb_n2.A, solver.x) - cb_n2.storage_vec .-= cb_n2.b - return norm(cb_n2.storage_vec) ≤ cb_n2.tol -end - -mutable struct TestCallbackN2Adjoint{T, S, M} - A::M - b::S - c::S - storage_vec1::S - storage_vec2::S - tol::T -end -TestCallbackN2Adjoint(A, b, c; tol = 0.1) = TestCallbackN2Adjoint(A, b, c, similar(b), similar(c), tol) - -function (cb_n2::TestCallbackN2Adjoint)(solver) - mul!(cb_n2.storage_vec1, cb_n2.A, solver.x) - cb_n2.storage_vec1 .-= cb_n2.b - mul!(cb_n2.storage_vec2, cb_n2.A', solver.y) - cb_n2.storage_vec2 .-= cb_n2.c - return (norm(cb_n2.storage_vec1) ≤ cb_n2.tol && norm(cb_n2.storage_vec2) ≤ cb_n2.tol) -end - -mutable struct TestCallbackN2Shifts{T, S, M} - A::M - b::S - shifts::Vector{T} - tol::T -end -TestCallbackN2Shifts(A, b, shifts; tol = 0.1) = TestCallbackN2Shifts(A, b, shifts, tol) - -function (cb_n2::TestCallbackN2Shifts)(solver) - r = residuals(cb_n2.A, cb_n2.b, cb_n2.shifts, solver.x) - return all(map(norm, r) .≤ cb_n2.tol) -end - -mutable struct TestCallbackN2LS{T, S, M} - A::M - b::S - λ::T - storage_vec1::S - storage_vec2::S - tol::T -end -TestCallbackN2LS(A, b, λ; tol = 0.1) = TestCallbackN2LS(A, b, λ, similar(b), similar(b, size(A, 2)), tol) - -function (cb_n2::TestCallbackN2LS)(solver) - mul!(cb_n2.storage_vec1, cb_n2.A, solver.x) - cb_n2.storage_vec1 .-= cb_n2.b - mul!(cb_n2.storage_vec2, cb_n2.A', cb_n2.storage_vec1) - cb_n2.storage_vec2 .+= cb_n2.λ .* solver.x - return norm(cb_n2.storage_vec2) ≤ cb_n2.tol -end - -mutable struct TestCallbackN2LN{T, S, M} - A::M - b::S - λ::T - storage_vec::S - tol::T -end -TestCallbackN2LN(A, b, λ; tol = 0.1) = TestCallbackN2LN(A, b, λ, similar(b), tol) - -function (cb_n2::TestCallbackN2LN)(solver) - mul!(cb_n2.storage_vec, cb_n2.A, solver.x) - cb_n2.storage_vec .-= cb_n2.b - cb_n2.λ != 0 && (cb_n2.storage_vec .+= sqrt(cb_n2.λ) .* solver.s) - return norm(cb_n2.storage_vec) ≤ cb_n2.tol -end - -mutable struct TestCallbackN2SaddlePts{T, S, M} - A::M - b::S - c::S - storage_vec1::S - storage_vec2::S - tol::T -end -TestCallbackN2SaddlePts(A, b, c; tol = 0.1) = - TestCallbackN2SaddlePts(A, b, c, similar(b), similar(c), tol) - -function (cb_n2::TestCallbackN2SaddlePts)(solver) - mul!(cb_n2.storage_vec1, cb_n2.A, solver.y) - cb_n2.storage_vec1 .+= solver.x .- cb_n2.b - mul!(cb_n2.storage_vec2, cb_n2.A', solver.x) - cb_n2.storage_vec2 .-= solver.y .+ cb_n2.c - return (norm(cb_n2.storage_vec1) ≤ cb_n2.tol && norm(cb_n2.storage_vec2) ≤ cb_n2.tol) -end - -function restarted_gmres_callback_n2(solver::GmresSolver, A, b, stor, N, storage_vec, tol) - get_x_restarted_gmres!(solver, A, stor, N) - x = stor.x - mul!(storage_vec, A, x) - storage_vec .-= b - return (norm(storage_vec) ≤ tol) -end