diff --git a/src/algorithms/lyapunov.jl b/src/algorithms/lyapunov.jl index 21aba007..0d149206 100644 --- a/src/algorithms/lyapunov.jl +++ b/src/algorithms/lyapunov.jl @@ -153,6 +153,10 @@ 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 + return 𝐂, reached_tol < tol, 0, reached_tol # return info on convergence end @@ -197,6 +201,10 @@ function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) + if reached_tol > tol + println("Lyapunov: doubling $reached_tol") + end + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -243,6 +251,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -299,6 +311,10 @@ function solve_lyapunov_equation( A::AbstractSparseMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) + if reached_tol > tol + println("Lyapunov: doubling $reached_tol") + end + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -352,6 +368,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -390,6 +410,10 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐗 * A' + C - 𝐗) / ℒ.norm(𝐗) + if reached_tol > tol + println("Lyapunov: bicgstab $reached_tol") + end + return 𝐗, reached_tol < tol, info.niter, reached_tol end @@ -428,6 +452,10 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐗 * A' + C - 𝐗) / ℒ.norm(𝐗) + if reached_tol > tol + println("Lyapunov: gmres $reached_tol") + end + return 𝐗, reached_tol < tol, info.niter, reached_tol end @@ -474,6 +502,10 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64}, reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂) + if reached_tol > tol + println("Lyapunov: iterative $reached_tol") + end + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -496,5 +528,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 + return 𝐂, reached_tol < tol, soll.maps, reached_tol end diff --git a/src/algorithms/quadratic_matrix_equation.jl b/src/algorithms/quadratic_matrix_equation.jl index 4e38f5e3..5577fc03 100644 --- a/src/algorithms/quadratic_matrix_equation.jl +++ b/src/algorithms/quadratic_matrix_equation.jl @@ -193,6 +193,10 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, reached_tol = ℒ.norm(AXX) / AXXnorm + if reached_tol > tol + println("QME: schur $reached_tol") + end + return X, reached_tol < tol, iter, reached_tol # schur can fail end @@ -381,6 +385,10 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, converged = reached_tol < tol + if reached_tol > tol + println("QME: doubling $reached_tol") + end + return X, converged, iter, reached_tol end @@ -434,6 +442,10 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, end # timeit_debug + if reached_tol > tol + println("QME: linear time iteration $reached_tol") + end + return X, converged, sol.maps, reached_tol end @@ -486,6 +498,10 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, converged = reached_tol < tol + if reached_tol > tol + println("QME: quadratic iteration $reached_tol") + end + return X, converged, sol.maps, reached_tol end diff --git a/src/algorithms/sylvester.jl b/src/algorithms/sylvester.jl index 2546d222..ee428149 100644 --- a/src/algorithms/sylvester.jl +++ b/src/algorithms/sylvester.jl @@ -343,6 +343,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -422,6 +426,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -523,6 +531,10 @@ function solve_sylvester_equation( A::Matrix{Float64}, end # timeit_debug end # timeit_debug + if reached_tol > tol + println("Sylvester: doubling $reached_tol") + end + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -599,6 +611,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -677,6 +693,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -755,6 +775,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -832,6 +856,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -926,6 +954,10 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} end # timeit_debug + if reached_tol > tol + println("Sylvester: doubling $reached_tol") + end + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -963,6 +995,10 @@ 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 + return 𝐂, reached_tol < tol, -1, reached_tol # return info on convergence end @@ -1079,6 +1115,10 @@ 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 + return 𝐗, reached_tol < tol, info.niter, reached_tol end @@ -1181,6 +1221,10 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, end # timeit_debug + if reached_tol > tol + println("Sylvester: dqgmres $reached_tol") + end + return 𝐗, reached_tol < tol, info.niter, reached_tol end @@ -1283,6 +1327,10 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, end # timeit_debug + if reached_tol > tol + println("Sylvester: gmres $reached_tol") + end + return 𝐗, reached_tol < tol, info.niter, reached_tol end @@ -1346,6 +1394,10 @@ 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 + return 𝐂, reached_tol < tol, iters, reached_tol # return info on convergence end @@ -1387,5 +1439,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 + return 𝐂, reached_tol < tol, soll.maps, reached_tol end diff --git a/src/filter/find_shocks.jl b/src/filter/find_shocks.jl index f2f36e80..e5c1a1ff 100644 --- a/src/filter/find_shocks.jl +++ b/src/filter/find_shocks.jl @@ -127,6 +127,9 @@ function find_shocks(::Val{:LagrangeNewton}, # println("Norm: $(ℒ.norm(x̂) / max(norm1,norm2))") # 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λ))") + end return x, ℒ.norm(x̂) / max(norm1,norm2) < tol && ℒ.norm(Δxλ) / ℒ.norm(xλ) < sqrt(tol) end @@ -347,6 +350,11 @@ function find_shocks(::Val{:LagrangeNewton}, # println("Norm: $(ℒ.norm(x̂) / max(norm1,norm2))") # println(ℒ.norm(Δxλ)) # 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λ))") + end + return x, ℒ.norm(x̂) / max(norm1,norm2) < tol && ℒ.norm(Δxλ) / ℒ.norm(xλ) < sqrt(tol) end