From b80e32f3a8c23d6f86da73db0ea77864a38b31cd Mon Sep 17 00:00:00 2001
From: kalmar <kalmar@amu.edu.pl>
Date: Sat, 14 Jan 2017 15:24:16 +0100
Subject: [PATCH] Parallel compute_SOS

---
 SL(3,Z).jl     | 38 +++++++++++++++++++++++---------------
 property(T).jl | 23 +++++++++++++----------
 2 files changed, 36 insertions(+), 25 deletions(-)

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;