Skip to content

Commit

Permalink
handle sufficient init guess outside algos (qme, sylvester)
Browse files Browse the repository at this point in the history
  • Loading branch information
@thorek1
thorek1 committed Oct 20, 2024
1 parent 646254d commit 4dafdf1
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 75 deletions.
18 changes: 11 additions & 7 deletions src/algorithms/quadratic_matrix_equation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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)

Expand Down
88 changes: 20 additions & 68 deletions src/algorithms/sylvester.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 4dafdf1

Please sign in to comment.