From 4dafdf1a287cb3b750fbce21283345e9303ff54f Mon Sep 17 00:00:00 2001 From: thorek1 Date: Sun, 20 Oct 2024 22:15:39 +0200 Subject: [PATCH] handle sufficient init guess outside algos (qme, sylvester) --- src/algorithms/quadratic_matrix_equation.jl | 18 +++-- src/algorithms/sylvester.jl | 88 +++++---------------- 2 files changed, 31 insertions(+), 75 deletions(-) diff --git a/src/algorithms/quadratic_matrix_equation.jl b/src/algorithms/quadratic_matrix_equation.jl index 721d9f8d..969b70e3 100644 --- a/src/algorithms/quadratic_matrix_equation.jl +++ b/src/algorithms/quadratic_matrix_equation.jl @@ -13,8 +13,18 @@ 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-12, + tol::AbstractFloat = 1e-14, verbose::Bool = false) where R <: Real + if length(initial_guess) > 0 + guess_ϵ = ℒ.norm(A * initial_guess ^ 2 + B * initial_guess + C) / ℒ.norm(A * initial_guess ^ 2) + + if guess_ϵ < tol # 1e-12 is too large eps is too small + if verbose println("Quadratic matrix equation solver previous solution has tolerance: $reached_tol") end + + return initial_guess, true + end + end + sol, solved, iterations, reached_tol = solve_quadratic_matrix_equation(A, B, C, Val(quadratic_matrix_equation_solver), T; @@ -197,12 +207,6 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, initial_guess = zero(A) end - guess_ϵ = ℒ.norm(A * initial_guess ^ 2 + B * initial_guess + C) / ℒ.norm(A * initial_guess ^ 2) - - if guess_ϵ < tol # 1e-12 is too large eps is too small - return initial_guess, true, 0, guess_ϵ - end - E = copy(C) F = copy(A) diff --git a/src/algorithms/sylvester.jl b/src/algorithms/sylvester.jl index 16b53500..ea4ac178 100644 --- a/src/algorithms/sylvester.jl +++ b/src/algorithms/sylvester.jl @@ -17,6 +17,23 @@ function solve_sylvester_equation(A::M, tol::AbstractFloat = 1e-14, timer::TimerOutput = TimerOutput(), verbose::Bool = false) where {M <: AbstractMatrix{Float64}, N <: AbstractMatrix{Float64}, O <: AbstractMatrix{Float64}} + @timeit_debug timer "Check if guess solves it already" begin + + if length(initial_guess) > 0 + 𝐂 = A * initial_guess * B + C - initial_guess + + reached_tol = ℒ.norm(𝐂) / ℒ.norm(initial_guess) + + if reached_tol < tol + if verbose println("Sylvester equation - previous solution achieves relative tol of $reached_tol") end + + # X = choose_matrix_format(initial_guess) + + return X, true + end + end + + end # timeit_debug @timeit_debug timer "Choose matrix formats" begin if sylvester_algorithm == :sylvester @@ -205,11 +222,6 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) max_iter = 500 @@ -275,11 +287,6 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) 𝐂¹ = similar(𝐂) 𝐂B = copy(C) @@ -363,11 +370,6 @@ function solve_sylvester_equation( A::Matrix{Float64}, # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) 𝐂¹ = similar(𝐂) 𝐂B = similar(C) @@ -464,11 +466,6 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) 𝐂¹ = similar(𝐂) 𝐂B = copy(C) @@ -546,14 +543,9 @@ function solve_sylvester_equation( A::Matrix{Float64}, # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) 𝐂¹ = similar(𝐂) - 𝐂B = copy(C) + # 𝐂B = copy(C) max_iter = 500 @@ -629,11 +621,6 @@ function solve_sylvester_equation( A::AbstractSparseMatrix{Float64}, # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) 𝐂¹ = similar(𝐂) # 𝐂B = copy(C) @@ -711,14 +698,9 @@ function solve_sylvester_equation( A::Matrix{Float64}, # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) 𝐂¹ = similar(𝐂) - 𝐂B = copy(C) + # 𝐂B = copy(C) max_iter = 500 @@ -795,11 +777,6 @@ function solve_sylvester_equation( A::Union{ℒ.Adjoint{Float64,Matrix{Float64} # 𝐂 = length(init) == 0 ? copy(C) : copy(init) 𝐂 = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - # ℒ.rmul!(𝐂, -1) 𝐂¹ = similar(𝐂) 𝐂B = copy(C) @@ -886,11 +863,6 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, 𝐂¹ = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂¹) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂¹) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - 𝐂 = try MatrixEquations.sylvd(-A, B, 𝐂¹) catch @@ -930,11 +902,6 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, 𝐂¹ = A * initial_guess * B + C - initial_guess #copy(C) - if ℒ.norm(𝐂¹) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂¹) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - tmp̄ = zero(C) 𝐗 = zero(C) end # timeit_debug @@ -1072,11 +1039,6 @@ function solve_sylvester_equation(A::DenseMatrix{Float64}, 𝐂¹ = A * initial_guess * B + C - initial_guess # 𝐂¹ = copy(C) - if ℒ.norm(𝐂¹) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂¹) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - tmp̄ = similar(C) 𝐗 = similar(C) end # timeit_debug @@ -1172,12 +1134,7 @@ function solve_sylvester_equation(A::AbstractMatrix{Float64}, 𝐂 = A * initial_guess * B + C - initial_guess 𝐂⁰ = copy(𝐂) # 𝐂 = copy(C) - - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - + 𝐂¹ = similar(C) 𝐂B = similar(C) @@ -1240,11 +1197,6 @@ function solve_sylvester_equation(A::AbstractMatrix{Float64}, 𝐂 = A * initial_guess * B + C - initial_guess # 𝐂 = copy(C) - if ℒ.norm(𝐂) / ℒ.norm(initial_guess) < tol - if verbose println("Previous solution of sylvester equation achieves relative tol of $(ℒ.norm(𝐂) / ℒ.norm(initial_guess))") end - return initial_guess, true, 0, 0.0 - end - if !(C isa DenseMatrix) C = collect(C) end