diff --git a/src/endgame_tracker.jl b/src/endgame_tracker.jl index 78051253..fb92cf1b 100644 --- a/src/endgame_tracker.jl +++ b/src/endgame_tracker.jl @@ -40,6 +40,7 @@ These parameters control the behaviour during the endgame. * `val_finite_tol = 1e-3`: Tolerance on the valuation which has to be satisfied before the endgame is started. * `sing_cond = 1e14`: value for the condition number above which a solution is considered singular. * `sing_accuracy = 1e-12`: value for the accuracy number above which a solution is considered singular. +* `scaling_threshold = -30.0`: Row scaling of matrices is only applied to rows with `e < scaling_threshold`, where is the norm of the row is estimated to be `2^e`. See [`skeel_row_scaling`](@skeel_row_scaling) for details. * `refine_steps = 3`: number of steps for refining solutions at the end. """ Base.@kwdef mutable struct EndgameOptions @@ -63,6 +64,7 @@ Base.@kwdef mutable struct EndgameOptions # singular solutions parameters sing_cond::Float64 = 1e14 sing_accuracy::Float64 = 1e-12 + scaling_threshold::Float64 = -30.0 # refinement parameters refine_steps::Int = 3 @@ -386,7 +388,7 @@ function step!(endgame_tracker::EndgameTracker, debug::Bool = false) update!(state.val, tracker.predictor, t) debug && println(state.val) if check_finite!(state, options) - switch_to_singular!(state, tracker; debug = debug) + switch_to_singular!(state, tracker, options; debug = debug) return state.code elseif check_at_infinity!(state, tracker, options; debug = debug) return state.code @@ -442,6 +444,7 @@ function check_at_infinity!(state, tracker, options; debug::Bool = false) state.row_scaling, workspace(tracker.state.jacobian), state.col_scaling, + options.scaling_threshold, ) end κ = LA.cond(tracker.state.jacobian, state.row_scaling, state.col_scaling) @@ -494,7 +497,7 @@ function check_at_infinity!(state, tracker, options; debug::Bool = false) false end -function switch_to_singular!(state, tracker; debug::Bool) +function switch_to_singular!(state, tracker, options; debug::Bool) state.singular_endgame = true t = real(tracker.state.t) state.singular_start = t @@ -506,6 +509,7 @@ function switch_to_singular!(state, tracker; debug::Bool) state.row_scaling, workspace(tracker.state.jacobian), state.col_scaling, + options.scaling_threshold, ) end add_sample!(state, tracker, t) @@ -703,6 +707,7 @@ function tracking_stopped!(endgame_tracker::EndgameTracker) state.row_scaling, workspace(tracker.state.jacobian), state.col_scaling, + options.scaling_threshold, ) state.cond = LA.cond(tracker, state.solution, 0.0, state.row_scaling, state.col_scaling) diff --git a/src/linear_algebra.jl b/src/linear_algebra.jl index 081658cb..7bf235ab 100644 --- a/src/linear_algebra.jl +++ b/src/linear_algebra.jl @@ -418,19 +418,21 @@ end """ skeel_row_scaling!(W::MatrixWorkspace, c) - skeel_row_scaling!(d, A, c) + skeel_row_scaling!(d, A, c; + scaling_threshold = -30.0) Compute optimal scaling factors `d` for the matrix `A` following Skeel [^S79] if `c` is approximately of the order of the solution of the linear system of interest. -The scaling factors are rounded to powers of the base radix 2. +The scaling factors are rounded to powers `e` of the base radix 2. Row scaling is only applied to rows with `e - m ≥ scaling_threshold` (this is an inequality of norms at log-scale) to prevent zero rows from being scaled, where `2^m` is approximately the maximum norm of a row of `A`. [^S79]: Skeel, Robert D. "Scaling for numerical stability in Gaussian elimination." Journal of the ACM (JACM) 26.3 (1979): 494-526. """ function skeel_row_scaling!( d::AbstractVector{<:Real}, A::AbstractMatrix{<:Complex}, - c::AbstractVector{<:Real}, + c::AbstractVector{<:Real}; + scaling_threshold::Float64 = -30.0, ) n = length(c) @inbounds d .= zero(eltype(d)) @@ -440,8 +442,16 @@ function skeel_row_scaling!( d[i] += fast_abs(A[i, j]) * cj end end + + m = maximum(d) + s = scaling_threshold + m @inbounds for i = 1:n - d[i] = exp2(-last(frexp(d[i]))) + e = last(frexp(d[i])) + if e < s + d[i] = 1.0 + else + d[i] = exp2(-e) + end end d @@ -456,10 +466,11 @@ function row_scaling!( d::AbstractVector{<:Real}, WS::MatrixWorkspace, c::AbstractVector{<:Real}, + scaling_threshold::Float64, ) m, n = size(WS) if m == n - skeel_row_scaling!(d, WS.A, c) + skeel_row_scaling!(d, WS.A, c; scaling_threshold = scaling_threshold) else d .= 1.0 end @@ -754,6 +765,7 @@ function LA.cond( rmax = max(rmax, rᵢ) rmin = min(rmin, rᵢ) end + #@show (rmax, rmin) rmax / rmin else inverse_inf_norm_est(WS, d_l, d_r) * inf_norm(WS, d_l, d_r) diff --git a/src/numerical_irreducible_decomposition.jl b/src/numerical_irreducible_decomposition.jl index 246701d0..43cf85a6 100644 --- a/src/numerical_irreducible_decomposition.jl +++ b/src/numerical_irreducible_decomposition.jl @@ -449,7 +449,7 @@ function regeneration!( endgame_options = EndgameOptions(; max_endgame_steps = 100, max_endgame_extended_steps = 100, - sing_accuracy = 1e-10, + sing_cond = 1e12, ), threading::Bool = true, seed = rand(UInt32), diff --git a/test/nid_test.jl b/test/nid_test.jl index ee20d9a0..8c4a43c1 100644 --- a/test/nid_test.jl +++ b/test/nid_test.jl @@ -28,6 +28,9 @@ N = nid(F; seed = nothing) @test isnothing(seed(N)) + # bad seed + N = nid(F; seed = 0xc770fa47) + # progress N = nid(F; show_progress = false) @test isa(N, NumericalIrreducibleDecomposition)