Skip to content

Commit

Permalink
Clean up
Browse files Browse the repository at this point in the history
  • Loading branch information
blegat committed Sep 25, 2023
1 parent ae87405 commit 6731515
Show file tree
Hide file tree
Showing 2 changed files with 6 additions and 98 deletions.
64 changes: 6 additions & 58 deletions src/Certificate/Symmetry/block_diag.jl
Original file line number Diff line number Diff line change
@@ -1,38 +1,5 @@
import DataStructures

function _givens::T, i, j, n) where {T<:Number}
c = cos(θ)
s = sin(θ)
if T <: Complex
r = sqrt(c * conj(c) + s * conj(s))
c /= r
s /= r
end
I = collect(1:n)
J = copy(I)
V = ones(T, n)
V[i] = c
V[j] = conj(c)
push!(I, i)
push!(J, j)
push!(V, -conj(s))
push!(I, j)
push!(J, i)
push!(V, s)
return SparseArrays.sparse(I, J, V, n, n)
end

"""
_permutation_quasi_upper_triangular(S)
Given a matrix such that `S[i, j]` is zero, return a Given rotation `G` such
that `(G' * S * G)[j, i]` is zero.
"""
function _givens(A::AbstractMatrix, i, j)
n = LinearAlgebra.checksquare(A)
return _givens(atan(A[j, i] / (A[i, i] - A[j, j])), i, j, n)
end

"""
_reorder!(F::LinearAlgebra.Schur{T}) where {T}
Expand Down Expand Up @@ -130,16 +97,19 @@ function _try_integer!(A::Matrix)
end

"""
_rotate_complex(A, B)
_rotate_complex(A::AbstractMatrix{T}, B::AbstractMatrix{T}; tol = Base.rtoldefault(real(T))) where {T}
Given (quasi) upper triangular matrix `A` and `B` that have the eigenvalues in
the same order except the complex pairs which may need to be (signed) permuted,
returns an othogonal matrix `P` such that `P' * A * P = B`.
returns an othogonal matrix `P` such that `P' * A * P` and `B` have matching
low triangular part.
The upper triangular part will be dealt with by `_sign_diag`.
By (quasi), we mean that if `S` is a `Matrix{<:Real}`,
then there may be nonzero entries in `S[i+1,i]` representing
complex conjugates.
If `S` is a `Matrix{<:Complex}`, then `S` is triangular.
If `S` is a `Matrix{<:Complex}`, then `S` is upper triangular so there is
nothing to do.
"""
function _rotate_complex(A::AbstractMatrix{T}, B::AbstractMatrix{T}; tol = Base.rtoldefault(real(T))) where {T}
n = LinearAlgebra.checksquare(A)
Expand All @@ -153,30 +123,8 @@ function _rotate_complex(A::AbstractMatrix{T}, B::AbstractMatrix{T}; tol = Base.
end
pair = abs(A[i + 1, i]) > tol
if pair
#if T <: Complex
# a = A[i, i + 1]
# b = B[i, i + 1]
# else
a = (A[i + 1, i], A[i, i + 1])
b = (B[i + 1, i], B[i, i + 1])
# end
# rot = b / a
# @assert abs(rot) ≈ 1
# if T <: Complex
# V[i] = conj(rot)
# else
# θ = atan(imag(rot) / real(rot))
# c = cos(θ)
# s = sin(θ)
# V[i] = c
# V[j] = c
# push!(I, i)
# push!(J, j)
# push!(V, -s)
# push!(I, j)
# push!(J, i)
# push!(V, s)
# end
c = a[2:-1:1]
if LinearAlgebra.norm(abs.(a) .- abs.(b)) > LinearAlgebra.norm(abs.(c) .- abs.(b))
a = c
Expand Down
40 changes: 0 additions & 40 deletions test/symmetry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,46 +3,6 @@ module TestSymmetry
using LinearAlgebra, SparseArrays, Test
using SumOfSquares

function _test_givens(A)
@test A[1, 2] 0 atol = 1e-10
G = Symmetry._givens(A, 1, 2)
@test G * G' I
@test G' * G I
B = G' * A * G
@test B[2, 1] 0 atol = 1e-10
end

function _test_givens(A1, A2)
_test_givens(A1)
_test_givens(A2)
_test_givens(A1 + A2 * im)
_test_givens(A2 + A1 * im)
end

function _test_givens(T::Type)
A1 = T[
1 0
2 -1
]
A2 = T[
-2 0
-1 -3
]
A3 = T[
-2 0
1 3
]
_test_givens(A1, A2)
_test_givens(A2, A3)
_test_givens(A3, A1)
end

function test_givens()
for T in [Int, Float64]
_test_givens(T)
end
end

function test_linsolve()
x = [1, 2]
for A in [
Expand Down

0 comments on commit 6731515

Please sign in to comment.