diff --git a/Project.toml b/Project.toml
index 1735d92e..c5544afb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,9 +7,11 @@ 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"
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QDLDL = "bfc457fd-c171-5ab7-bd9e-d5dbfc242d63"
@@ -29,8 +31,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..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
+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 29a31498..dfadf11b 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,6 +203,23 @@ function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T}
 		rank_k_update!(X, ws)
 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 = zero(X)
+  buffer = zero(X)
+  for i in eachindex(w)
+    w[i] <= T(0) && continue
+    z = view(Z, :, i)
+    MA.mutable_operate_to!(buffer, *, z, z')
+    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
+
 function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where {T}
   n = size(X, 1)
   @. X = zero(T)
@@ -274,9 +291,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
diff --git a/src/interface.jl b/src/interface.jl
index 302dbe44..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)