diff --git a/src/algorithms/lyapunov.jl b/src/algorithms/lyapunov.jl index 0d149206..2a95839d 100644 --- a/src/algorithms/lyapunov.jl +++ b/src/algorithms/lyapunov.jl @@ -27,64 +27,69 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, end # timeit_debug @timeit_debug timer "Solve" begin - X, solved, i, reached_tol = solve_lyapunov_equation(A, C, + X, i, reached_tol = solve_lyapunov_equation(A, C, Val(lyapunov_algorithm), - tol = tol, + # tol = tol, timer = timer) if verbose - println("Lyapunov equation - converged to tol $tol: $solved; iterations: $i; reached tol: $reached_tol; algorithm: $lyapunov_algorithm") + println("Lyapunov equation - converged to tol $tol: $(reached_tol < tol); iterations: $i; reached tol: $reached_tol; algorithm: $lyapunov_algorithm") end - if !solved + if reached_tol > tol if (reached_tol < sqrt(tol) || A isa AbstractSparseMatrix) && lyapunov_algorithm ≠ :bicgstab C = collect(C) - X, solved, i, reached_tol = solve_lyapunov_equation(A, C, + X, i, reached_tol = solve_lyapunov_equation(A, C, Val(:bicgstab), - tol = tol, + # tol = tol, timer = timer) if verbose - println("Lyapunov equation - converged to tol $tol: $solved; iterations: $i; reached tol: $reached_tol; algorithm: gmres") + println("Lyapunov equation - converged to tol $tol: $(reached_tol < tol); iterations: $i; reached tol: $reached_tol; algorithm: gmres") end else A = collect(A) C = collect(C) - X, solved, i, reached_tol = solve_lyapunov_equation(A, C, + X, i, reached_tol = solve_lyapunov_equation(A, C, Val(:lyapunov), - tol = tol, + # tol = tol, timer = timer) if verbose - println("Lyapunov equation - converged to tol $tol: $solved; iterations: $i; reached tol: $reached_tol; algorithm: lyapunov") + println("Lyapunov equation - converged to tol $tol: $(reached_tol < tol); iterations: $i; reached tol: $reached_tol; algorithm: lyapunov") end end end end # timeit_debug end # timeit_debug + + if (reached_tol > tol) println("Lyapunov failed: $reached_tol") end - - return X, solved + return X, reached_tol < tol end function rrule(::typeof(solve_lyapunov_equation), A::AbstractMatrix{Float64}, C::AbstractMatrix{Float64}; lyapunov_algorithm::Symbol = :doubling, - tol::AbstractFloat = 1e-9, + tol::AbstractFloat = 1e-12, verbose::Bool = false, timer::TimerOutput = TimerOutput()) - P, solved = solve_lyapunov_equation(A, C, lyapunov_algorithm = lyapunov_algorithm, tol = tol, verbose = verbose) + P, solved = solve_lyapunov_equation(A, C, lyapunov_algorithm = lyapunov_algorithm, + # tol = tol, + verbose = verbose) # pullback # https://arxiv.org/abs/2011.11430 function solve_lyapunov_equation_pullback(∂P) - ∂C, solved = solve_lyapunov_equation(A', ∂P[1], lyapunov_algorithm = lyapunov_algorithm, tol = tol, verbose = verbose) + ∂C, solved = solve_lyapunov_equation(A', ∂P[1], lyapunov_algorithm = lyapunov_algorithm, + # tol = tol, + verbose = verbose) ∂A = ∂C * A * P' + ∂C' * A * P @@ -99,14 +104,16 @@ end function solve_lyapunov_equation( A::AbstractMatrix{ℱ.Dual{Z,S,N}}, C::AbstractMatrix{ℱ.Dual{Z,S,N}}; lyapunov_algorithm::Symbol = :doubling, - tol::AbstractFloat = 1e-9, + tol::AbstractFloat = 1e-12, verbose::Bool = false, timer::TimerOutput = TimerOutput()) where {Z,S,N} # unpack: AoS -> SoA Â = ℱ.value.(A) Ĉ = ℱ.value.(C) - P̂, solved = solve_lyapunov_equation(Â, Ĉ, lyapunov_algorithm = lyapunov_algorithm, tol = tol, verbose = verbose) + P̂, solved = solve_lyapunov_equation(Â, Ĉ, lyapunov_algorithm = lyapunov_algorithm, + # tol = tol, + verbose = verbose) Ã = copy(Â) C̃ = copy(Ĉ) @@ -122,7 +129,9 @@ function solve_lyapunov_equation( A::AbstractMatrix{ℱ.Dual{Z,S,N}}, if ℒ.norm(X) < eps() continue end - P, solved = solve_lyapunov_equation(Â, X, lyapunov_algorithm = lyapunov_algorithm, tol = tol, verbose = verbose) + P, solved = solve_lyapunov_equation(Â, X, lyapunov_algorithm = lyapunov_algorithm, + # tol = tol, + verbose = verbose) P̃[:,i] = vec(P) end @@ -153,11 +162,11 @@ function solve_lyapunov_equation(A::Union{ℒ.Adjoint{Float64,Matrix{Float64}},D reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) - if reached_tol > tol - println("Lyapunov: lyapunov $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: lyapunov $reached_tol") + # end - return 𝐂, reached_tol < tol, 0, reached_tol # return info on convergence + return 𝐂, 0, reached_tol # return info on convergence end @@ -165,7 +174,7 @@ end function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; - tol::Float64 = 1e-12, + tol::Float64 = 1e-14, timer::TimerOutput = TimerOutput()) 𝐂 = copy(C) 𝐀 = copy(A) @@ -201,18 +210,18 @@ function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) - if reached_tol > tol - println("Lyapunov: doubling $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end function solve_lyapunov_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, C::AbstractSparseMatrix{Float64}, ::Val{:doubling}; - tol::Float64 = 1e-12, + tol::Float64 = 1e-14, timer::TimerOutput = TimerOutput()) 𝐂 = copy(C) 𝐀 = copy(A) @@ -251,18 +260,18 @@ function solve_lyapunov_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) - if reached_tol > tol - println("Lyapunov: doubling $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:doubling}; - tol::Float64 = 1e-12, + tol::Float64 = 1e-14, timer::TimerOutput = TimerOutput()) 𝐂 = copy(C) 𝐀 = copy(A) @@ -311,11 +320,11 @@ function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) - if reached_tol > tol - println("Lyapunov: doubling $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -324,7 +333,7 @@ end function solve_lyapunov_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:doubling}; - tol::Float64 = 1e-12, + tol::Float64 = 1e-14, timer::TimerOutput = TimerOutput()) 𝐂 = copy(C) 𝐂¹ = copy(C) @@ -368,11 +377,11 @@ function solve_lyapunov_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) - if reached_tol > tol - println("Lyapunov: doubling $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -381,7 +390,7 @@ end function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:bicgstab}; - tol::Float64 = 1e-8, + tol::Float64 = 1e-14, timer::TimerOutput = TimerOutput()) tmp̄ = similar(C) 𝐗 = similar(C) @@ -395,7 +404,7 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, lyapunov = LinearOperators.LinearOperator(Float64, length(C), length(C), true, true, lyapunov!) - 𝐂, info = Krylov.bicgstab(lyapunov, [vec(C);], rtol = tol / 100, atol = tol / 100) + 𝐂, info = Krylov.bicgstab(lyapunov, [vec(C);], rtol = tol, atol = tol) copyto!(𝐗, 𝐂) @@ -410,18 +419,18 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐗 * A' + C - 𝐗) / ℒ.norm(𝐗) - if reached_tol > tol - println("Lyapunov: bicgstab $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: bicgstab $reached_tol") + # end - return 𝐗, reached_tol < tol, info.niter, reached_tol + return 𝐗, info.niter, reached_tol end function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:gmres}; - tol::Float64 = 1e-8, + tol::Float64 = 1e-14, timer::TimerOutput = TimerOutput()) tmp̄ = similar(C) 𝐗 = similar(C) @@ -437,7 +446,7 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, lyapunov = LinearOperators.LinearOperator(Float64, length(C), length(C), true, true, lyapunov!) - 𝐂, info = Krylov.gmres(lyapunov, [vec(C);], rtol = tol / 100, atol = tol / 100) + 𝐂, info = Krylov.gmres(lyapunov, [vec(C);], rtol = tol, atol = tol) copyto!(𝐗, 𝐂) @@ -452,18 +461,18 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐗 * A' + C - 𝐗) / ℒ.norm(𝐗) - if reached_tol > tol - println("Lyapunov: gmres $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: gmres $reached_tol") + # end - return 𝐗, reached_tol < tol, info.niter, reached_tol + return 𝐗, info.niter, reached_tol end function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:iterative}; - tol::AbstractFloat = 1e-12, + tol::AbstractFloat = 1e-14, timer::TimerOutput = TimerOutput()) 𝐂 = copy(C) 𝐂¹ = copy(C) @@ -502,18 +511,18 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) - if reached_tol > tol - println("Lyapunov: iterative $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: iterative $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end function solve_lyapunov_equation(A::AbstractMatrix{Float64}, C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}}, ::Val{:speedmapping}; - tol::AbstractFloat = 1e-12, + tol::AbstractFloat = 1e-14, timer::TimerOutput = TimerOutput()) 𝐂A = similar(C) @@ -528,9 +537,9 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) - if reached_tol > tol - println("Lyapunov: speedmapping $reached_tol") - end + # if reached_tol > tol + # println("Lyapunov: speedmapping $reached_tol") + # end - return 𝐂, reached_tol < tol, soll.maps, reached_tol + return 𝐂, soll.maps, reached_tol end diff --git a/src/algorithms/quadratic_matrix_equation.jl b/src/algorithms/quadratic_matrix_equation.jl index 5577fc03..ee820c94 100644 --- a/src/algorithms/quadratic_matrix_equation.jl +++ b/src/algorithms/quadratic_matrix_equation.jl @@ -13,7 +13,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, initial_guess::AbstractMatrix{R} = zeros(0,0), quadratic_matrix_equation_solver::Symbol = :doubling, timer::TimerOutput = TimerOutput(), - tol::AbstractFloat = 1e-13, # 1e-14 is too tight + tol::AbstractFloat = 1e-10, # 1e-14 is too tight verbose::Bool = false) where R <: Real if length(initial_guess) > 0 X = initial_guess @@ -35,41 +35,41 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, end end - sol, solved, iterations, reached_tol = solve_quadratic_matrix_equation(A, B, C, + sol, iterations, reached_tol = solve_quadratic_matrix_equation(A, B, C, Val(quadratic_matrix_equation_solver), T; initial_guess = initial_guess, - tol = tol, + # tol = tol, timer = timer, verbose = verbose) - if verbose println("Quadratic matrix equation solver: $quadratic_matrix_equation_solver - converged: $solved in $iterations iterations to tolerance: $reached_tol") end + if verbose println("Quadratic matrix equation solver: $quadratic_matrix_equation_solver - converged: $(reached_tol < tol) in $iterations iterations to tolerance: $reached_tol") end - if !solved + if reached_tol > tol if quadratic_matrix_equation_solver ≠ :schur # try schur if previous one didn't solve it - sol, solved, iterations, reached_tol = solve_quadratic_matrix_equation(A, B, C, + sol, iterations, reached_tol = solve_quadratic_matrix_equation(A, B, C, Val(:schur), T; initial_guess = initial_guess, - tol = tol, + # tol = tol, timer = timer, verbose = verbose) - if verbose println("Quadratic matrix equation solver: schur - converged: $solved in $iterations iterations to tolerance: $reached_tol") end + if verbose println("Quadratic matrix equation solver: schur - converged: $(reached_tol < tol) in $iterations iterations to tolerance: $reached_tol") end else quadratic_matrix_equation_solver ≠ :doubling - sol, solved, iterations, reached_tol = solve_quadratic_matrix_equation(A, B, C, + sol, iterations, reached_tol = solve_quadratic_matrix_equation(A, B, C, Val(:doubling), T; initial_guess = initial_guess, - tol = tol, + # tol = tol, timer = timer, verbose = verbose) - if verbose println("Quadratic matrix equation solver: doubling - converged: $solved in $iterations iterations to tolerance: $reached_tol") end + if verbose println("Quadratic matrix equation solver: doubling - converged: $(reached_tol < tol) in $iterations iterations to tolerance: $reached_tol") end end end - return sol, solved + return sol, reached_tol < tol end function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, @@ -193,11 +193,11 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, reached_tol = ℒ.norm(AXX) / AXXnorm - if reached_tol > tol - println("QME: schur $reached_tol") - end + # if reached_tol > tol + # println("QME: schur $reached_tol") + # end - return X, reached_tol < tol, iter, reached_tol # schur can fail + return X, iter, reached_tol # schur can fail end @@ -383,13 +383,11 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, reached_tol = ℒ.norm(AXX) / AXXnorm - converged = reached_tol < tol - - if reached_tol > tol - println("QME: doubling $reached_tol") - end + # if reached_tol > tol + # println("QME: doubling $reached_tol") + # end - return X, converged, iter, reached_tol + return X, iter, reached_tol end @@ -400,7 +398,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, ::Val{:linear_time_iteration}, T::timings; initial_guess::AbstractMatrix{R} = zeros(0,0), - tol::AbstractFloat = 1e-8, # lower tol not possible for NAWM (and probably other models this size) + tol::AbstractFloat = 1e-14, # lower tol not possible for NAWM (and probably other models this size) timer::TimerOutput = TimerOutput(), verbose::Bool = false, max_iter::Int = 5000) where R <: Real @@ -438,15 +436,13 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, reached_tol = ℒ.norm(AXX) / AXXnorm - converged = reached_tol < tol - end # timeit_debug - if reached_tol > tol - println("QME: linear time iteration $reached_tol") - end + # if reached_tol > tol + # println("QME: linear time iteration $reached_tol") + # end - return X, converged, sol.maps, reached_tol + return X, sol.maps, reached_tol end @@ -457,7 +453,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, ::Val{:quadratic_iteration}, T::timings; initial_guess::AbstractMatrix{R} = zeros(0,0), - tol::AbstractFloat = 1e-8, # lower tol not possible for NAWM (and probably other models this size) + tol::AbstractFloat = 1e-14, # lower tol not possible for NAWM (and probably other models this size) timer::TimerOutput = TimerOutput(), verbose::Bool = false, max_iter::Int = 50000) where R <: Real @@ -496,13 +492,11 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, reached_tol = ℒ.norm(AXX) / AXXnorm - converged = reached_tol < tol - - if reached_tol > tol - println("QME: quadratic iteration $reached_tol") - end + # if reached_tol > tol + # println("QME: quadratic iteration $reached_tol") + # end - return X, converged, sol.maps, reached_tol + return X, sol.maps, reached_tol end @@ -512,7 +506,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{ℱ.Dual{Z,S,N}}, C::AbstractMatrix{ℱ.Dual{Z,S,N}}, T::timings; initial_guess::AbstractMatrix{<:Real} = zeros(0,0), - tol::AbstractFloat = 1e-12, + tol::AbstractFloat = 1e-10, quadratic_matrix_equation_solver::Symbol = :doubling, timer::TimerOutput = TimerOutput(), verbose::Bool = false) where {Z,S,N} @@ -524,6 +518,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{ℱ.Dual{Z,S,N}}, X, solved = solve_quadratic_matrix_equation(Â, B̂, Ĉ, Val(quadratic_matrix_equation_solver), T; + tol = tol, initial_guess = initial_guess, timer = timer, verbose = verbose) diff --git a/src/algorithms/sylvester.jl b/src/algorithms/sylvester.jl index ee428149..59b3e9b0 100644 --- a/src/algorithms/sylvester.jl +++ b/src/algorithms/sylvester.jl @@ -15,7 +15,7 @@ function solve_sylvester_equation(A::M, C::O; initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), sylvester_algorithm::Symbol = :doubling, - tol::AbstractFloat = 1e-14, + tol::AbstractFloat = 1e-10, timer::TimerOutput = TimerOutput(), verbose::Bool = false) where {M <: AbstractMatrix{Float64}, N <: AbstractMatrix{Float64}, O <: AbstractMatrix{Float64}} @timeit_debug timer "Choose matrix formats" begin @@ -56,27 +56,27 @@ function solve_sylvester_equation(A::M, end # timeit_debug @timeit_debug timer "Solve sylvester equation" begin - x, solved, i, reached_tol = solve_sylvester_equation(a, b, c, Val(sylvester_algorithm), + x, i, reached_tol = solve_sylvester_equation(a, b, c, Val(sylvester_algorithm), initial_guess = initial_guess, - tol = tol, + # tol = tol, verbose = verbose, timer = timer) if verbose && i != 0 - println("Sylvester equation - converged to tol $tol: $solved; iterations: $i; reached tol: $reached_tol; algorithm: $sylvester_algorithm") + println("Sylvester equation - converged to tol $tol: $(reached_tol < tol); iterations: $i; reached tol: $reached_tol; algorithm: $sylvester_algorithm") end - if !solved && sylvester_algorithm ≠ :sylvester && length(B) < 5e7 # try sylvester if previous one didn't solve it + if !(reached_tol < tol) && sylvester_algorithm ≠ :sylvester && length(B) < 5e7 # try sylvester if previous one didn't solve it aa = collect(A) bb = collect(B) cc = collect(C) - x, solved, i, reached_tol = solve_sylvester_equation(aa, bb, cc, + x, i, reached_tol = solve_sylvester_equation(aa, bb, cc, Val(:sylvester), initial_guess = zeros(0,0), - tol = tol, + # tol = tol, verbose = verbose, timer = timer) @@ -85,15 +85,15 @@ function solve_sylvester_equation(A::M, end end - if !solved && reached_tol < sqrt(tol) + if !(reached_tol < tol) && reached_tol < sqrt(tol) aa = collect(A) cc = collect(C) - X, solved, i, Reached_tol = solve_sylvester_equation(aa, b, cc, + X, i, Reached_tol = solve_sylvester_equation(aa, b, cc, Val(:dqgmres), initial_guess = x, - tol = tol, + # tol = tol, verbose = verbose, timer = timer) if Reached_tol < reached_tol @@ -102,36 +102,36 @@ function solve_sylvester_equation(A::M, end if verbose# && i != 0 - println("Sylvester equation - converged to tol $tol: $solved; iterations: $i; reached tol: $Reached_tol; algorithm: dqgmres (refinement of previous solution)") + println("Sylvester equation - converged to tol $tol: $(reached_tol < tol); iterations: $i; reached tol: $Reached_tol; algorithm: dqgmres (refinement of previous solution)") end end - if !solved && sylvester_algorithm ≠ :gmres + if !(reached_tol < tol) && sylvester_algorithm ≠ :gmres aa = collect(A) cc = collect(C) - x, solved, i, reached_tol = solve_sylvester_equation(aa, b, cc, + x, i, reached_tol = solve_sylvester_equation(aa, b, cc, Val(:gmres), initial_guess = zeros(0,0), - tol = tol, + # tol = tol, verbose = verbose, timer = timer) if verbose# && i != 0 - println("Sylvester equation - converged to tol $tol: $solved; iterations: $i; reached tol: $reached_tol; algorithm: gmres") + println("Sylvester equation - converged to tol $tol: $(reached_tol < tol); iterations: $i; reached tol: $reached_tol; algorithm: gmres") end end - if !solved && reached_tol < sqrt(tol) + if !(reached_tol < tol) && reached_tol < sqrt(tol) aa = collect(A) cc = collect(C) - X, solved, i, Reached_tol = solve_sylvester_equation(aa, b, cc, + X, i, Reached_tol = solve_sylvester_equation(aa, b, cc, Val(:dqgmres), initial_guess = x, - tol = tol, + # tol = tol, verbose = verbose, timer = timer) if Reached_tol < reached_tol @@ -203,7 +203,9 @@ function solve_sylvester_equation(A::M, X = choose_matrix_format(x)# |> sparse - return X, solved + if (reached_tol > tol) println("Sylvester failed: $reached_tol") end + + return X, reached_tol < tol end @@ -292,7 +294,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 # guess_provided = true @@ -343,11 +345,11 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -359,7 +361,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 # guess_provided = true @@ -426,11 +428,11 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -443,7 +445,7 @@ function solve_sylvester_equation( A::Matrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 @timeit_debug timer "Doubling solve" begin @timeit_debug timer "Setup buffers" begin @@ -531,11 +533,11 @@ function solve_sylvester_equation( A::Matrix{Float64}, end # timeit_debug end # timeit_debug - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -546,7 +548,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 On Smith-type iterative algorithms for the Stein matrix equation # guess_provided = true @@ -611,11 +613,11 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -627,7 +629,7 @@ function solve_sylvester_equation( A::Matrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 # guess_provided = true @@ -693,11 +695,11 @@ function solve_sylvester_equation( A::Matrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -709,7 +711,7 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 # guess_provided = true @@ -775,11 +777,11 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -790,7 +792,7 @@ function solve_sylvester_equation( A::Matrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 # guess_provided = true @@ -856,11 +858,11 @@ function solve_sylvester_equation( A::Matrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -871,7 +873,7 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-12) + tol::Float64 = 1e-14) # see doi:10.1016/j.aml.2009.01.012 @timeit_debug timer "Setup buffers" begin @@ -954,11 +956,11 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} end # timeit_debug - if reached_tol > tol - println("Sylvester: doubling $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: doubling $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -995,11 +997,11 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: sylvester $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: sylvester $reached_tol") + # end - return 𝐂, reached_tol < tol, -1, reached_tol # return info on convergence + return 𝐂, -1, reached_tol # return info on convergence end @@ -1010,7 +1012,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-8) + tol::Float64 = 1e-14) @timeit_debug timer "Preallocate matrices" begin # guess_provided = true @@ -1115,11 +1117,11 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, 𝐗 = choose_matrix_format(𝐗, density_threshold = 1.0) end - if reached_tol > tol - println("Sylvester: bicgstab $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: bicgstab $reached_tol") + # end - return 𝐗, reached_tol < tol, info.niter, reached_tol + return 𝐗, info.niter, reached_tol end @@ -1130,7 +1132,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-8) + tol::Float64 = 1e-14) @timeit_debug timer "Preallocate matrices" begin # guess_provided = true @@ -1221,11 +1223,11 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, end # timeit_debug - if reached_tol > tol - println("Sylvester: dqgmres $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: dqgmres $reached_tol") + # end - return 𝐗, reached_tol < tol, info.niter, reached_tol + return 𝐗, info.niter, reached_tol end @@ -1236,7 +1238,7 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::Float64 = 1e-8) + tol::Float64 = 1e-14) @timeit_debug timer "Preallocate matrices" begin # guess_provided = true @@ -1327,11 +1329,11 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, end # timeit_debug - if reached_tol > tol - println("Sylvester: gmres $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: gmres $reached_tol") + # end - return 𝐗, reached_tol < tol, info.niter, reached_tol + return 𝐗, info.niter, reached_tol end @@ -1342,7 +1344,7 @@ function solve_sylvester_equation(A::AbstractMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::AbstractFloat = 1e-12) + tol::AbstractFloat = 1e-14) # guess_provided = true if length(initial_guess) == 0 @@ -1394,11 +1396,11 @@ function solve_sylvester_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: iterative $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: iterative $reached_tol") + # end - return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence + return 𝐂, iters, reached_tol # return info on convergence end @@ -1409,7 +1411,7 @@ function solve_sylvester_equation(A::AbstractMatrix{Float64}, initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), timer::TimerOutput = TimerOutput(), verbose::Bool = false, - tol::AbstractFloat = 1e-12) + tol::AbstractFloat = 1e-14) # guess_provided = true if length(initial_guess) == 0 @@ -1439,9 +1441,9 @@ function solve_sylvester_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * B + C - 𝐂) / max(ℒ.norm(𝐂), ℒ.norm(C)) - if reached_tol > tol - println("Sylvester: speedmapping $reached_tol") - end + # if reached_tol > tol + # println("Sylvester: speedmapping $reached_tol") + # end - return 𝐂, reached_tol < tol, soll.maps, reached_tol + return 𝐂, soll.maps, reached_tol end diff --git a/src/filter/find_shocks.jl b/src/filter/find_shocks.jl index e5c1a1ff..353d63dd 100644 --- a/src/filter/find_shocks.jl +++ b/src/filter/find_shocks.jl @@ -128,7 +128,7 @@ function find_shocks(::Val{:LagrangeNewton}, # println(ℒ.norm(Δxλ)) # println(ℒ.norm(Δxλ) / ℒ.norm(xλ)) if !(ℒ.norm(x̂) / max(norm1,norm2) < tol && ℒ.norm(Δxλ) / ℒ.norm(xλ) < sqrt(tol)) - println("Norm 1: $(ℒ.norm(x̂) / max(norm1,norm2)); Norm 2: $(ℒ.norm(Δxλ) / ℒ.norm(xλ))") + println("Find shocks failed. Norm 1: $(ℒ.norm(x̂) / max(norm1,norm2)); Norm 2: $(ℒ.norm(Δxλ) / ℒ.norm(xλ))") end return x, ℒ.norm(x̂) / max(norm1,norm2) < tol && ℒ.norm(Δxλ) / ℒ.norm(xλ) < sqrt(tol) end @@ -352,9 +352,9 @@ function find_shocks(::Val{:LagrangeNewton}, # println(ℒ.norm(x̂) / max(norm1,norm2) < tol && ℒ.norm(Δxλ) / ℒ.norm(xλ) < tol) if !(ℒ.norm(x̂) / max(norm1,norm2) < tol && ℒ.norm(Δxλ) / ℒ.norm(xλ) < sqrt(tol)) - println("Norm 1: $(ℒ.norm(x̂) / max(norm1,norm2)); Norm 2: $(ℒ.norm(Δxλ) / ℒ.norm(xλ))") + println("Find shocks failed. Norm 1: $(ℒ.norm(x̂) / max(norm1,norm2)); Norm 2: $(ℒ.norm(Δxλ) / ℒ.norm(xλ))") end - + return x, ℒ.norm(x̂) / max(norm1,norm2) < tol && ℒ.norm(Δxλ) / ℒ.norm(xλ) < sqrt(tol) end