diff --git a/SL(3,Z).jl b/SL(3,Z).jl index d756640..3e4352f 100644 --- a/SL(3,Z).jl +++ b/SL(3,Z).jl @@ -2,11 +2,11 @@ using JuMP import SCS: SCSSolver import Mosek: MosekSolver -push!(LOAD_PATH, "./") +workers_processes = addprocs() +@everywhere push!(LOAD_PATH, "./") using GroupAlgebras -include("property(T).jl") - +@everywhere include("property(T).jl") function E(i::Int, j::Int, N::Int=3) @assert i≠j @@ -39,33 +39,41 @@ const TOL=10.0^-7 # κ, A = solve_for_property_T(S₁, solver, verbose=VERBOSE) -product_matrix = readdlm("SL3Z.product_matrix", Int) -L = readdlm("SL3Z.delta.coefficients")[:, 1] -Δ = GroupAlgebraElement(L, product_matrix) +const product_matrix = readdlm("SL3Z.product_matrix", Int) +const L = readdlm("SL3Z.delta.coefficients")[:, 1] +const Δ = GroupAlgebraElement(L, product_matrix) -A = readdlm("SL3Z.SDPmatrixA.Mosek") -κ = readdlm("SL3Z.kappa.Mosek")[1] +const A = readdlm("SL3Z.SDPmatrixA.Mosek") +const κ = readdlm("SL3Z.kappa.Mosek")[1] @assert isapprox(eigvals(A), abs(eigvals(A)), atol=TOL) @assert A == Symmetric(A) const A_sqrt = real(sqrtm(A)) -SOS_fp_diff, SOS_fp_L₁_distance = check_solution(κ, A_sqrt, Δ) +const SOS_fp_diff, SOS_fp_L₁_distance = check_solution(κ, A_sqrt, Δ) @show SOS_fp_L₁_distance @show GroupAlgebras.ɛ(SOS_fp_diff) -κ_rational = rationalize(BigInt, κ;) -A_sqrt_rational = rationalize(BigInt, A_sqrt) -Δ_rational = rationalize(BigInt, Δ) +const κ_rational = rationalize(BigInt, κ, tol=TOL) +const A_sqrt_rational = rationalize(BigInt, A_sqrt, tol=TOL) +const Δ_rational = rationalize(BigInt, Δ, tol=TOL) -SOS_rational_diff, SOS_rat_L₁_distance = check_solution(κ_rational, A_sqrt_rational, Δ_rational) +const SOS_rational_diff, SOS_rat_L₁_distance = check_solution(κ_rational, A_sqrt_rational, Δ_rational) @assert isa(SOS_rat_L₁_distance, Rational{BigInt}) @show float(SOS_rat_L₁_distance) @show float(GroupAlgebras.ɛ(SOS_rational_diff)) -A_sqrt_augmented = correct_to_augmentation_ideal(A_sqrt_rational) +const A_sqrt_augmented = correct_to_augmentation_ideal(A_sqrt_rational) -SOS_rational_diff_aug, SOS_rat_L₁_distance_aug = check_solution(κ_rational, A_sqrt_augmented, Δ_rational) +const SOS_rational_aug_diff, SOS_aug_rat_L₁_distance = check_solution(κ_rational, A_sqrt_augmented, Δ_rational) + +@assert isa(SOS_aug_rat_L₁_distance, Rational{BigInt}) +@assert GroupAlgebras.ɛ(SOS_rational_aug_diff) == 0//1 + +@show float(SOS_aug_rat_L₁_distance) +@show float(κ_rational - 2^3*SOS_aug_rat_L₁_distance) + +rmprocs(workers_processes) diff --git a/property(T).jl b/property(T).jl index 38f6ecf..d43b275 100644 --- a/property(T).jl +++ b/property(T).jl @@ -1,5 +1,6 @@ using JuMP import Base: rationalize +using GroupAlgebras function products{T<:Real}(S1::Array{Array{T,2},1}, S2::Array{Array{T,2},1}) result = [0*similar(S1[1])] @@ -132,16 +133,19 @@ function EOI{T<:Number}(Δ::GroupAlgebraElement{T}, κ::T) return Δ*Δ - κ*Δ end -function resulting_SOS{T<:Number}(sqrt_matrix::Array{T,2}, - elt::GroupAlgebraElement{T}) - result = zeros(elt.coefficients) +@everywhere function square(vector, elt) zzz = zeros(elt.coefficients) + zzz[1:length(vector)] = vector +# new_base_elt = GroupAlgebraElement(zzz, elt.product_matrix) +# return (new_base_elt*new_base_elt).coefficients + return GroupAlgebras.algebra_multiplication(zzz, zzz, elt.product_matrix) +end + +function compute_SOS{T<:Number}(sqrt_matrix::Array{T,2}, + elt::GroupAlgebraElement{T}) L = size(sqrt_matrix,2) - for i in 1:L - info("$i of $L") - zzz[1:L] = view(sqrt_matrix, :,i) - new_base = GroupAlgebraElement(zzz, elt.product_matrix) - result += (new_base*new_base).coefficients + result = @parallel (+) for i in 1:L + square(sqrt_matrix[:,i], elt) end return GroupAlgebraElement{T}(result, elt.product_matrix) end @@ -161,14 +165,13 @@ function check_solution{T<:Number}(κ::T, sqrt_matrix::Array{T,2}, Δ::GroupAlgebraElement{T}) eoi = EOI(Δ, κ) - result = resulting_SOS(sqrt_matrix, Δ) + result = compute_SOS(sqrt_matrix, Δ) L₁_dist = norm(result - eoi,1) return eoi - result, L₁_dist end function rationalize{T<:Integer, S<:Real}(::Type{T}, X::AbstractArray{S}; tol::Real=eps(eltype(X))) - r(x) = rationalize(T, x, tol=tol) return r.(X) end;