From 4e135d81d436bd4cee0f277057b3d2a0417fe3d1 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Tue, 29 Jun 2021 16:44:18 +0200 Subject: [PATCH 1/9] Adding support for BigFLoat cones --- Project.toml | 3 ++- src/COSMO.jl | 2 +- src/convexset.jl | 6 ++++++ src/interface.jl | 2 +- 4 files changed, 10 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 1735d92e..82c0d905 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.8.1" AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" COSMOAccelerators = "bbd8fffe-5ad0-4d78-a55e-85575421b4ac" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" @@ -29,8 +30,8 @@ DataStructures = "^0.17.0, ^0.18.0" IterTools = "^1" MathOptInterface = "~0.9.16" QDLDL = "~0.1.4" -Requires = "^1" Reexport = "0.2, ^1" +Requires = "^1" UnsafeArrays = "0.3, 1" julia = "^1" diff --git a/src/COSMO.jl b/src/COSMO.jl index e8cf5acc..88a2559a 100644 --- a/src/COSMO.jl +++ b/src/COSMO.jl @@ -1,7 +1,7 @@ __precompile__() module COSMO -using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg, DataStructures, Requires, Printf, IterTools +using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg, DataStructures, Requires, Printf, IterTools, GenericLinearAlgebra using Reexport using COSMOAccelerators diff --git a/src/convexset.jl b/src/convexset.jl index 29a31498..026e5ab5 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -203,6 +203,12 @@ function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T} rank_k_update!(X, ws) end +function _project!(X::AbstractMatrix{BigFloat}, ws::PsdBlasWorkspace{BigFloat}) + w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) + D = LinearAlgebra.Diagonal(max.(w,0)) + X .= Z'D'Z +end + function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where {T} n = size(X, 1) @. X = zero(T) diff --git a/src/interface.jl b/src/interface.jl index 302dbe44..93e36072 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -388,7 +388,7 @@ function type_checks(constraints::Vector{COSMO.Constraint{T}}) where {T <: Abstr return nothing end type_checks(convex_set::AbstractConvexSet) = nothing -type_checks(convex_set::Union{PsdCone{BigFloat}, PsdConeTriangle{BigFloat}}) = throw(ArgumentError("COSMO currently does not support the combination of PSD constraints and BigFloat.")) +# type_checks(convex_set::Union{PsdCone{BigFloat}, PsdConeTriangle{BigFloat}}) = throw(ArgumentError("COSMO currently does not support the combination of PSD constraints and BigFloat.")) From 86070814b7413e751eef312888b1eac6517a7f26 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Tue, 29 Jun 2021 16:56:26 +0200 Subject: [PATCH 2/9] removing completely the typecheck --- src/interface.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/interface.jl b/src/interface.jl index 93e36072..f61a8c5e 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -388,8 +388,6 @@ function type_checks(constraints::Vector{COSMO.Constraint{T}}) where {T <: Abstr return nothing end type_checks(convex_set::AbstractConvexSet) = nothing -# type_checks(convex_set::Union{PsdCone{BigFloat}, PsdConeTriangle{BigFloat}}) = throw(ArgumentError("COSMO currently does not support the combination of PSD constraints and BigFloat.")) - function check_A_dim(A::Union{AbstractVector{<:Real},AbstractMatrix{<:Real}}, n::Int) From 75dc3a0911f73cc520794863b8546464efd9ebc3 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Tue, 29 Jun 2021 17:30:13 +0200 Subject: [PATCH 3/9] Better dispatch --- src/convexset.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/convexset.jl b/src/convexset.jl index 026e5ab5..478b6ff7 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -183,7 +183,7 @@ for (syevr, elty) in end #@eval end #for -function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T} +function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T <: Union{Float32,Float64}} #computes the upper triangular part of the projection of X onto the PSD cone @@ -203,7 +203,7 @@ function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T} rank_k_update!(X, ws) end -function _project!(X::AbstractMatrix{BigFloat}, ws::PsdBlasWorkspace{BigFloat}) +function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T} w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) D = LinearAlgebra.Diagonal(max.(w,0)) X .= Z'D'Z From 2421625e13b8d5ebaee454c1001050a7c85c42c4 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Tue, 29 Jun 2021 17:44:38 +0200 Subject: [PATCH 4/9] Remove D --- src/convexset.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/convexset.jl b/src/convexset.jl index 478b6ff7..eda70f5a 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -205,8 +205,7 @@ end function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T} w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) - D = LinearAlgebra.Diagonal(max.(w,0)) - X .= Z'D'Z + X .= Z'LinearAlgebra.Diagonal(max.(w,0))*Z end function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where {T} From 337a36547e95cf74d8300ce31e685c81e4e8ca89 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Wed, 30 Jun 2021 14:21:41 +0200 Subject: [PATCH 5/9] formula was false --- src/convexset.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/convexset.jl b/src/convexset.jl index eda70f5a..4b2bca1d 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -205,7 +205,7 @@ end function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T} w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) - X .= Z'LinearAlgebra.Diagonal(max.(w,0))*Z + X .= Z*LinearAlgebra.Diagonal(max.(w,0))*Z' end function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where {T} From 3e25ee54655b57748e87d6ffb26dcffa6e687245 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Wed, 30 Jun 2021 14:35:17 +0200 Subject: [PATCH 6/9] Tring MutableArithmetics + small change. --- src/convexset.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/convexset.jl b/src/convexset.jl index 4b2bca1d..ded47547 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -205,6 +205,7 @@ end function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T} w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) + #X .= MutableArithmetics.mul!(Z,LinearAlgebra.Diagonal(max.(w,0)),Z') " does not work. X .= Z*LinearAlgebra.Diagonal(max.(w,0))*Z' end @@ -279,9 +280,11 @@ function project!(x::AbstractVector{T}, cone::Union{PsdCone{T}, DensePsdCone{T}} symmetrize_upper!(X) _project!(X, cone.work) - #fill in the lower triangular part - for j=1:n, i=1:(j-1) - X[j,i] = X[i,j] + #fill in the lower triangular part, only needed when calling BLAS. + if T <: Union{Float32,Float64} + for j=1:n, i=1:(j-1) + X[j,i] = X[i,j] + end end end return nothing From c5b1f70bf059685449a1a3a8934300694c807ff6 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Wed, 30 Jun 2021 16:53:11 +0200 Subject: [PATCH 7/9] Move to mutable arithmetics for aribitrary type PSD projections --- Project.toml | 1 + src/COSMO.jl | 3 ++- src/convexset.jl | 15 +++++++++++++-- 3 files changed, 16 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 82c0d905..c5544afb 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QDLDL = "bfc457fd-c171-5ab7-bd9e-d5dbfc242d63" diff --git a/src/COSMO.jl b/src/COSMO.jl index 88a2559a..dbeb534e 100644 --- a/src/COSMO.jl +++ b/src/COSMO.jl @@ -1,7 +1,7 @@ __precompile__() module COSMO -using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg, DataStructures, Requires, Printf, IterTools, GenericLinearAlgebra +using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg, DataStructures, Requires, Printf, IterTools, GenericLinearAlgebra, MutableArithmetics using Reexport using COSMOAccelerators @@ -12,6 +12,7 @@ export assemble!, warm_start!, empty_model!, update! const DefaultFloat = Float64 const DefaultInt = LinearAlgebra.BlasInt +const MA = MutableArithmetics include("./kktsolver.jl") diff --git a/src/convexset.jl b/src/convexset.jl index ded47547..1bf7381d 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -205,8 +205,19 @@ end function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T} w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) - #X .= MutableArithmetics.mul!(Z,LinearAlgebra.Diagonal(max.(w,0)),Z') " does not work. - X .= Z*LinearAlgebra.Diagonal(max.(w,0))*Z' + + # The follwoing lines performe the operation : X = Z*LinearAlgebra.Diagonal(max.(w,0))*Z' + # using MutableArithmetics. + X .= T(0) + buffer = zero(X) + for i in eachindex(w) + w[i] <= 0 && continue + z = view(Z, :, i) + MA.mutable_operate_to!(buffer, *, z, z') + for I in eachindex(X) + X[I] = MA.add_mul!(X[I], w[i], buffer[I]) + end + end end function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where {T} From edfcbe9d719fa442a8d24e3fe62d2d69ebe35827 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Thu, 1 Jul 2021 10:49:08 +0200 Subject: [PATCH 8/9] Change comment --- src/convexset.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/convexset.jl b/src/convexset.jl index 1bf7381d..03c8c7c4 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -205,9 +205,7 @@ end function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T} w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) - - # The follwoing lines performe the operation : X = Z*LinearAlgebra.Diagonal(max.(w,0))*Z' - # using MutableArithmetics. + # The follwoing lines uses MutableArithmetics to perform the operation : X .= Z*LinearAlgebra.Diagonal(max.(w,0))*Z' X .= T(0) buffer = zero(X) for i in eachindex(w) From ae81709b41ef8230c86ae24131f8f3e4ae41c7c3 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Thu, 1 Jul 2021 11:14:27 +0200 Subject: [PATCH 9/9] Temporary fix --- src/convexset.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/convexset.jl b/src/convexset.jl index 03c8c7c4..dfadf11b 100644 --- a/src/convexset.jl +++ b/src/convexset.jl @@ -206,14 +206,16 @@ end function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T} w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X)) # The follwoing lines uses MutableArithmetics to perform the operation : X .= Z*LinearAlgebra.Diagonal(max.(w,0))*Z' - X .= T(0) + X = zero(X) buffer = zero(X) for i in eachindex(w) - w[i] <= 0 && continue + w[i] <= T(0) && continue z = view(Z, :, i) MA.mutable_operate_to!(buffer, *, z, z') - for I in eachindex(X) - X[I] = MA.add_mul!(X[I], w[i], buffer[I]) + for j in eachindex(X) + #The following line should be : X[j] = MA.add_mul!(X[j], w[i], buffer[j]) + # However it does not work for BigFLoats yet, so i keep the allocating working line : + X[j] += w[i]*buffer[j] end end end