diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index e1250f9..5159fdf 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -17,13 +17,13 @@ jobs:
           - version: '1'
             os: ubuntu-latest
             arch: x64
-          - version: '1.0'
+          - version: '1.6'
             os: ubuntu-latest
             arch: x64
           - version: '1'
             os: windows-latest
             arch: x64
-          - version: '1.0'
+          - version: '1.6'
             os: windows-latest
             arch: x64
     steps:
diff --git a/Project.toml b/Project.toml
index b49b699..5774b2f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ProxSDP"
 uuid = "65e78d25-6039-50a4-9445-38022e3d2eb3"
 repo = "https://github.com/mariohsouto/ProxSDP.jl.git"
-version = "1.6.2"
+version = "1.7.0"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
@@ -18,7 +18,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 julia = "1"
 Arpack = "0.3.2, 0.5.1"
 KrylovKit = "0.5.2"
-MathOptInterface = "0.9.14"
+MathOptInterface = "0.10.6"
 TimerOutputs = "0.5.0"
 
 [extras]
diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
index ad419c1..892e86e 100644
--- a/src/MOI_wrapper.jl
+++ b/src/MOI_wrapper.jl
@@ -1,92 +1,89 @@
-
 const MOI = MathOptInterface
-const CI = MOI.ConstraintIndex
-const VI = MOI.VariableIndex
-
-const MOIU = MOI.Utilities
-
-Base.@kwdef mutable struct MOISolution
-    ret_val::Int = 0
-    raw_status::String = "Problem not solved"
-    primal::Vector{Float64} = Float64[] # primal of variables
-    dual_eq::Vector{Float64} = Float64[] # dual of constraints
-    dual_in::Vector{Float64} = Float64[] # dual of constraints
-    dual_cone::Vector{Float64} = Float64[] # dual of constraints
-    slack_eq::Vector{Float64} = Float64[]
-    slack_in::Vector{Float64} = Float64[]
-    # slack_sd::Vector{Vector{Float64}}
-    # slack_so::Vector{Vector{Float64}}
-    primal_residual::Float64 = NaN
-    dual_residual::Float64 = NaN
-    objval::Float64 = NaN
-    dual_objval::Float64 = NaN
-    # obj_cte::Float64 = NaN # used in rays - see SCS
-    gap::Float64 = NaN
-    time::Float64 = NaN
-    final_rank::Int = -1
-    primal_feasible_user_tol::Bool = false
-    dual_feasible_user_tol::Bool = false
-    certificate_found::Bool = false
-end
-
-# `ModelData` struct is discarded
-mutable struct ModelData
-    m::Int # Number of rows/constraints
-    n::Int # Number of cols/variables
-
-    # equality matrix
-    I::Vector{Int} # List of rows
-    J::Vector{Int} # List of cols
-    V::Vector{Float64} # List of coefficients
-    # equality rhs
-    b::Vector{Float64} # constants
-
-    # inequality matrix
-    Ii::Vector{Int} # List of rows
-    Ji::Vector{Int} # List of cols
-    Vi::Vector{Float64} # List of coefficients
-    # inequality rhs
-    h::Vector{Float64} # constants
-
-    objective_constant::Float64 # The objective is min c'x + objective_constant
-    c::Vector{Float64}
-end
-
-Base.@kwdef mutable struct ConeData
-
-    cone_cols::Int = 0
 
-    # matrix lines
-    eq_tot::Int = 0 # number of linear equality constraints
-    eq_start::Vector{Int} = Int[]
-    eq_len::Vector{Int} = Int[]
-    in_tot::Int = 0 # length of LP cone
-    in_start::Vector{Int} = Int[]
-    in_len::Vector{Int} = Int[]
+MOI.Utilities.@product_of_sets(Zeros, MOI.Zeros)
+
+MOI.Utilities.@product_of_sets(Nonpositives, MOI.Nonpositives)
+
+MOI.Utilities.@struct_of_constraints_by_set_types(
+    StructCache,
+    MOI.Zeros,
+    MOI.Nonpositives,
+    MOI.SecondOrderCone,
+    MOI.PositiveSemidefiniteConeTriangle,
+)
+
+const OptimizerCache{T} = MOI.Utilities.GenericModel{
+    T,
+    MOI.Utilities.ObjectiveContainer{T},
+    MOI.Utilities.VariablesContainer{T},
+    StructCache{T}{
+        MOI.Utilities.MatrixOfConstraints{
+            T,
+            MOI.Utilities.MutableSparseMatrixCSC{
+                T,
+                Int64,
+                MOI.Utilities.OneBasedIndexing,
+            },
+            Vector{T},
+            Zeros{T},
+        },
+        MOI.Utilities.MatrixOfConstraints{
+            T,
+            MOI.Utilities.MutableSparseMatrixCSC{
+                T,
+                Int64,
+                MOI.Utilities.OneBasedIndexing,
+            },
+            Vector{T},
+            Nonpositives{T},
+        },
+        MOI.Utilities.VectorOfConstraints{
+            MOI.VectorOfVariables,
+            MOI.SecondOrderCone,
+        },
+        MOI.Utilities.VectorOfConstraints{
+            MOI.VectorOfVariables,
+            MOI.PositiveSemidefiniteConeTriangle,
+        },
+    },
+}
 
-    # cones
-    sdc::Vector{Vector{Int}} = Vector{Int}[] # semidefinite
+Base.@kwdef mutable struct ConeData
+    psc::Vector{Vector{Int}} = Vector{Int}[] # semidefinite
     soc::Vector{Vector{Int}} = Vector{Int}[] # second order
 end
 
 mutable struct Optimizer <: MOI.AbstractOptimizer
-    cone::ConeData
-    maxsense::Bool
-    data::Union{Nothing, ModelData} # only non-Void between MOI.copy_to and MOI.optimize!
-    sol::MOISolution
+    cones::Union{Nothing, ConeData}
+    zeros::Union{Nothing, Zeros{Float64}}
+    nonps::Union{Nothing, Nonpositives{Float64}}
+    sol::Result
     options::Options
 end
+
 function Optimizer(;kwargs...)
-    optimizer = Optimizer(ConeData(), false, nothing, MOISolution(), Options())
+    optimizer = Optimizer(
+        nothing,
+        nothing,
+        nothing,
+        Result(),
+        Options()
+        )
     for (key, value) in kwargs
-        MOI.set(optimizer, MOI.RawParameter(key), value)
+        MOI.set(optimizer, MOI.RawOptimizerAttribute(string(key)), value)
     end
     return optimizer
 end
 
+#=
+    Basic Attributes
+=#
+
 MOI.get(::Optimizer, ::MOI.SolverName) = "ProxSDP"
 
-function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
+MOI.get(::Optimizer, ::MOI.SolverVersion) = "1.7.0"
+
+function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
     fields = fieldnames(Options)
     name = Symbol(param.name)
     if name in fields
@@ -96,7 +93,8 @@ function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
     end
     return value
 end
-function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
+
+function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
     fields = fieldnames(Options)
     name = Symbol(param.name)
     if name in fields
@@ -107,12 +105,14 @@ function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
 end
 
 MOI.supports(::Optimizer, ::MOI.Silent) = true
+
 function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
     if value == true
         optimizer.options.timer_verbose = false
     end
     optimizer.options.log_verbose = !value
 end
+
 function MOI.get(optimizer::Optimizer, ::MOI.Silent)
     if optimizer.options.log_verbose
         return false
@@ -124,446 +124,180 @@ function MOI.get(optimizer::Optimizer, ::MOI.Silent)
 end
 
 MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
+
 function MOI.set(optimizer::Optimizer, ::MOI.TimeLimitSec, value)
     optimizer.options.time_limit = value
 end
+
 function MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec)
     return optimizer.options.time_limit
 end
 
-
 MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = false
 
 function MOI.is_empty(optimizer::Optimizer)
-    !optimizer.maxsense && optimizer.data === nothing
+    return optimizer.cones === nothing &&
+    optimizer.zeros === nothing &&
+    optimizer.nonps === nothing &&
+    optimizer.sol.status == 0
 end
+
 function MOI.empty!(optimizer::Optimizer)
-    optimizer.maxsense = false
-    optimizer.data = nothing # It should already be nothing except if an error is thrown inside copy_to
-    optimizer.cone = ConeData()
-    optimizer.sol.ret_val = 0
+    optimizer.cones = nothing
+    optimizer.zeros = nothing
+    optimizer.nonps = nothing
+    optimizer.sol = Result()
+    return
 end
 
-MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names
+function _rows(
+    optimizer::Optimizer,
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{T},MOI.Zeros},
+) where T
+    return MOI.Utilities.rows(optimizer.zeros, ci)
+end
 
-function MOI.supports(::Optimizer,
-                      ::Union{MOI.ObjectiveSense,
-                              MOI.ObjectiveFunction{MOI.SingleVariable},
-                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
-    return true
+function _rows(
+    optimizer::Optimizer,
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{T},MOI.Nonpositives},
+) where T
+    return MOI.Utilities.rows(optimizer.nonps, ci)
 end
 
-function MOI.supports_constraint(
+# MOI.supports
+
+function MOI.supports(
     ::Optimizer,
-    ::Type{<:MOI.VectorAffineFunction{Float64}},
-    ::Type{<:Union{MOI.Zeros, MOI.Nonpositives}})
+    ::Union{
+        MOI.ObjectiveSense,
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
+    },
+) where T
     return true
 end
 
 function MOI.supports_constraint(
     ::Optimizer,
-    ::Type{<:MOI.VectorOfVariables},
-    ::Type{<:Union{MOI.SecondOrderCone,
-                   MOI.PositiveSemidefiniteConeTriangle}})
+    ::Type{MOI.VectorAffineFunction{T}},
+    ::Type{<:Union{MOI.Zeros, MOI.Nonpositives}}
+) where T
     return true
 end
 
-function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
-    return MOIU.automatic_copy_to(dest, src; kws...)
-end
-
-const SupportedSets = Union{MOI.PositiveSemidefiniteConeTriangle,
-                            MOI.SecondOrderCone}
-function get_new_variable_indexes_unique(optimizer::Optimizer, set::SupportedSets)
-    len = MOI.dimension(set)
-    first = optimizer.cone.cone_cols + 1
-    last = optimizer.cone.cone_cols + len
-    optimizer.cone.cone_cols += len
-    return collect(first:last)
-end
-function get_cone_list(cone::ConeData, set::MOI.PositiveSemidefiniteConeTriangle)
-    return cone.sdc
-end
-function get_cone_list(cone::ConeData, set::MOI.SecondOrderCone)
-    return cone.soc
-end
-function get_cone_list(optimizer::Optimizer, set::SupportedSets)
-    get_cone_list(optimizer.cone, set)
-end
-
-function MOIU.allocate_constrained_variables(optimizer::Optimizer, set::SupportedSets)
-    inds = get_new_variable_indexes_unique(optimizer,
-    # func,
-    set)
-    vars = MOI.VariableIndex.(inds)
-    list = get_cone_list(optimizer, set)
-    push!(list, inds)
-    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(length(list))
-    return vars, ci
-end
-
-function MOIU.load_constrained_variables(
-    optimizer::Optimizer, vis::Vector{MOI.VariableIndex},
-    ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
-    set::SupportedSets)
-end
-
-function _allocate_constraint(cone::ConeData, f, set::SupportedSets)
-    list = get_cone_list(cone, set)
-    push!(list, Int[])
-    return length(list)
-end
-
-function _allocate_constraint(cone::ConeData, f, s::MOI.Zeros)
-    len = MOI.dimension(s)
-    push!(cone.eq_start, cone.eq_tot + 1)
-    push!(cone.eq_len, len)
-    cone.eq_tot += len
-    return length(cone.eq_len)
-end
-
-function _allocate_constraint(cone::ConeData, f, s::MOI.Nonpositives)
-    len = MOI.dimension(s)
-    push!(cone.in_start, cone.in_tot + 1)
-    push!(cone.in_len, len)
-    cone.in_tot += len
-    return length(cone.in_len)
-end
-
-function MOIU.load_constraint(optimizer::Optimizer, ci::CI,
-    f::MOI.VectorOfVariables,
-    set::SupportedSets)
-
-    vec = get_cone_list(optimizer.cone, set)[ci.value]
-
-    append!(vec, [i.value for i in f.variables])
-    nothing
+function MOI.supports_add_constrained_variables(
+    ::Optimizer,
+    ::Type{<:Union{
+        MOI.SecondOrderCone,
+        MOI.PositiveSemidefiniteConeTriangle,
+        }
+    }
+)
+    return true
 end
 
-function MOIU.allocate_constraint(optimizer::Optimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
-    return CI{F, S}(_allocate_constraint(optimizer.cone, f, s))
+function cone_data(src::OptimizerCache, ::Type{T}) where T
+    cache = MOI.Utilities.constraints(
+        src.constraints,
+        MOI.VectorOfVariables,
+        T,
+    )
+    indices = MOI.get(
+        cache,
+        MOI.ListOfConstraintIndices{MOI.VectorOfVariables, T}(),
+    )
+    funcs = MOI.get.(cache, MOI.ConstraintFunction(), indices)
+    return Vector{Int64}[Int64[v.value for v in f.variables] for f in funcs]
 end
 
 # Vectorized length for matrix dimension n
-sympackedlen(n) = div(n*(n+1), 2)
+sympackedlen(n) = div(n * (n + 1), 2)
 # Matrix dimension for vectorized length n
-sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
-
-output_index(t::MOI.VectorAffineTerm) = t.output_index
-variable_index_value(v::MOI.VariableIndex) = v.value
-variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
-variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
-coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
-coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)
-
-function MOIU.load_constraint(optimizer::Optimizer, ci::CI,
-                              f::MOI.VectorAffineFunction,
-                              set::MOI.AbstractVectorSet)
-
-    n = length(f.terms)
-    if n != 1
-        A = sparse(output_index.(f.terms), variable_index_value.(f.terms),
-        coefficient.(f.terms))
-        # sparse combines duplicates with + but does
-        # not remove zeros created so we call dropzeros!
-        dropzeros!(A)
-        I, J, V = findnz(A)
-    else
-        I, J, V = output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms)
-    end
-#=
-    n = length(f.terms)
-    I = Int[]
-    J = Int[]
-    V = Float64[]
-    sizehint!(I, n)
-    sizehint!(J, n)
-    sizehint!(V, n)
-    for i in 1:n
-        if !iszero(coefficient(f.terms[i]))
-            push!(I, output_index(f.terms[i]))
-            push!(J, variable_index_value(f.terms[i]))
-            push!(V, coefficient(f.terms[i]))
-        end
-    end
-=#
-    #=
-    I,J,V = output_index.(f.terms), variable_index_value.(f.terms), coefficient.(f.terms)
-    =#
-
-    offset = ctr_offset(optimizer, set, ci)
-    rows = 1:len(optimizer, set, ci)
-    i = (offset-1) .+ rows
-
-    rhs = matrix_rhs_vec(optimizer, set)
-    rhs[i] .= -f.constants
-
-    # MOI:  Ax + b {==, <=} 0
-    # psdp: Ax {==, <=} -b
-    I_, J_, V_ = matrix_triplets(optimizer, set)
-    append!(I_, offset-1 .+ I)
-    append!(J_, J)
-    append!(V_, V)
-
-    nothing
-end
-function ctr_offset(optimizer::Optimizer, set::MOI.Zeros, ci)
-    return optimizer.cone.eq_start[ci.value]
-end
-function len(optimizer::Optimizer, set::MOI.Zeros, ci)
-    return optimizer.cone.eq_len[ci.value]
-end
-function ctr_offset(optimizer::Optimizer, set::MOI.Nonpositives, ci)
-    return optimizer.cone.in_start[ci.value]
-end
-function len(optimizer::Optimizer, set::MOI.Nonpositives, ci)
-    return optimizer.cone.in_len[ci.value]
-end
-function matrix_triplets(optimizer::Optimizer, set::MOI.Zeros)
-    return optimizer.data.I, optimizer.data.J, optimizer.data.V
-end
-function matrix_rhs_vec(optimizer::Optimizer, set::MOI.Zeros)
-    return optimizer.data.b
-end
-function matrix_triplets(optimizer::Optimizer, set::MOI.Nonpositives)
-    return optimizer.data.Ii, optimizer.data.Ji, optimizer.data.Vi
-end
-function matrix_rhs_vec(optimizer::Optimizer, set::MOI.Nonpositives)
-    return optimizer.data.h
-end
-
-# called after ctr var - only aditonal
-function MOIU.allocate_variables(optimizer::Optimizer, nvars::Integer)
-    allocated = optimizer.cone.cone_cols
-    VI.(allocated .+ (1:nvars))
-end
-
-function MOIU.load_variables(optimizer::Optimizer, nvars::Integer)
-    # @show "load var"
-    # @show nvars
-    cone = optimizer.cone
-    m = cone.eq_tot + cone.in_tot # + cone.q + cone.s + 3cone.ep + cone.ed
-    I = Int[]
-    J = Int[]
-    V = Float64[]
-    sizehint!(I, cone.eq_tot)
-    sizehint!(J, cone.eq_tot)
-    sizehint!(V, cone.eq_tot)
-    b = zeros(cone.eq_tot)
-    Ii = Int[]
-    Ji = Int[]
-    Vi = Float64[]
-    sizehint!(Ii, cone.in_tot)
-    sizehint!(Ji, cone.in_tot)
-    sizehint!(Vi, cone.in_tot)
-    h = zeros(cone.in_tot)
-
-    tot_vars = nvars# + cone.cone_cols
-    c = zeros(tot_vars)
-    optimizer.data = ModelData(m, tot_vars, I, J, V, b, Ii, Ji, Vi, h, 0., c)
-    # `optimizer.sol` contains the result of the previous optimization.
-    # It is used as a warm start if its length is the same, e.g.
-    # probably because no variable and/or constraint has been added.
-    if length(optimizer.sol.primal) != tot_vars
-        optimizer.sol.primal = zeros(tot_vars)
-    end
-    # TODO @joaquim
-    # @assert length(optimizer.sol.dual) == length(optimizer.sol.slack)
-    # if length(optimizer.sol.dual) != m
-    #     optimizer.sol.dual = zeros(m)
-    #     optimizer.sol.slack = zeros(m)
-    # end
-    return nothing
-end
-
-function MOIU.allocate(::Optimizer, ::MOI.VariablePrimalStart,
-                       ::MOI.VariableIndex, ::Union{Nothing, Float64})
-end
-function MOIU.allocate(::Optimizer, ::MOI.ConstraintPrimalStart,
-                       ::MOI.ConstraintIndex, ::Float64)
-end
-function MOIU.allocate(::Optimizer, ::MOI.ConstraintDualStart,
-                       ::MOI.ConstraintIndex, ::Float64)
-end
-function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense,
-                       sense::MOI.OptimizationSense)
-    optimizer.maxsense = sense == MOI.MAX_SENSE
-end
-function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction,
-    ::MOI.Union{MOI.SingleVariable,
-                MOI.ScalarAffineFunction{Float64}})
-end
+sympackeddim(n) = div(isqrt(1 + 8n) - 1, 2)
 
-function MOIU.load(::Optimizer, ::MOI.VariablePrimalStart,
-    ::MOI.VariableIndex, ::Nothing)
-end
-function MOIU.load(optimizer::Optimizer, ::MOI.VariablePrimalStart,
-                   vi::MOI.VariableIndex, value::Float64)
-    optimizer.sol.primal[vi.value] = value
-end
-function MOIU.load(::Optimizer, ::MOI.ConstraintPrimalStart,
-    ::MOI.ConstraintIndex, ::Nothing)
-end
-function MOIU.load(::Optimizer, ::MOI.ConstraintDualStart,
-    ::MOI.ConstraintIndex, ::Nothing)
-end
-function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense)
-end
-function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
-                   f::MOI.SingleVariable)
-    MOIU.load(optimizer,
-          MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
-          MOI.ScalarAffineFunction{Float64}(f))
-    return nothing
-end
-function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
-               f::MOI.ScalarAffineFunction)
-    c0 = Vector(sparsevec(variable_index_value.(f.terms), coefficient.(f.terms),
-                        optimizer.data.n))
-    optimizer.data.objective_constant = f.constant
-    optimizer.data.c = optimizer.maxsense ? -c0 : c0
-    return nothing
-end
-
-matindices(n::Integer) = (LinearIndices(tril(trues(n,n))))[findall(tril(trues(n,n)))]
-
-function MOI.optimize!(optimizer::Optimizer)
-    if optimizer.data === nothing
-        # optimize! has already been called and no new model has been copied
-        return
-    end
-
-    # @show "in opt"
+matindices(n::Integer) = LinearIndices(trues(n,n))[findall(LinearAlgebra.tril(trues(n,n)))]
 
+function _optimize!(dest::Optimizer, src::OptimizerCache)
+    MOI.empty!(dest)
     TimerOutputs.reset_timer!()
 
     @timeit "pre-processing" begin
 
-    # parse options
-    options = optimizer.options
-
-    cone = optimizer.cone
-
-    #= 
-        Build linear sets
+    #=
+        Affine
     =#
-
-    m = optimizer.data.m #rows
-    n = optimizer.data.n #cols
-
-    objective_constant = optimizer.data.objective_constant
-    # c = optimizer.data.c
-
-    # Build Prox SDP Affine Sets
-    # b = optimizer.data.b
-    # A = sparse(optimizer.data.I, optimizer.data.J, optimizer.data.V, length(b), n)
-    # h = optimizer.data.h
-    # G = sparse(optimizer.data.Ii, optimizer.data.Ji, optimizer.data.Vi, length(h), n)
-
-    # Dimensions (of affine sets)
-
-    aff = AffineSets(0, 0, 0, 0,
-    sparse(optimizer.data.I, optimizer.data.J, optimizer.data.V, length(optimizer.data.b), n),
-    sparse(optimizer.data.Ii, optimizer.data.Ji, optimizer.data.Vi, length(optimizer.data.h), n),
-    optimizer.data.b,
-    optimizer.data.h,
-    optimizer.data.c)
-    # A, G, b, h, c)
-
-    n_variables = size(aff.A)[2] # primal
-    n_eqs = size(aff.A)[1]
-    n_ineqs = size(aff.G)[1]
-    aff.n = n_variables
-    aff.p = n_eqs
-    aff.m = n_ineqs
-
-    #= 
-        Build conic sets
+    Ab = MOI.Utilities.constraints(
+        src.constraints,
+        MOI.VectorAffineFunction{Float64},
+        MOI.Zeros,
+    )
+    A = Ab.coefficients
+    b = -Ab.constants
+    Gh = MOI.Utilities.constraints(
+        src.constraints,
+        MOI.VectorAffineFunction{Float64},
+        MOI.Nonpositives,
+    )
+    G = Gh.coefficients
+    h = -Gh.constants
+    @assert A.n == G.n
+    #=
+        Objective
+    =#
+    max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
+    obj =
+        MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
+    objective_constant = MOI.constant(obj)
+    c = zeros(A.n)
+    obj_sign = ifelse(max_sense, -1.0, 1.0)
+    for term in obj.terms
+        c[term.variable.value] += obj_sign * term.coefficient
+    end
+    dest.zeros = deepcopy(Ab.sets) # TODO copy(Ab.sets)
+    dest.nonps = deepcopy(Gh.sets) # TODO copy(Gh.sets)
+    # TODO: simply this after
+    # https://github.com/jump-dev/MathOptInterface.jl/issues/1711
+    _A = if A.m == 0
+        SparseArrays.sparse(Int64[], Int64[], Float64[], 0, A.n)
+    else
+        convert(SparseArrays.SparseMatrixCSC{Float64,Int64}, A)
+    end
+    _G = if G.m == 0
+        SparseArrays.sparse(Int64[], Int64[], Float64[], 0, G.n)
+    else
+        convert(SparseArrays.SparseMatrixCSC{Float64,Int64}, G)
+    end
+    aff = AffineSets(A.n, A.m, G.m, 0,
+        _A, _G, b, h, c)
+    #=
+        Cones
     =#
-
-    # Build SDP Sets
     con = ConicSets(
         SDPSet[],
         SOCSet[]
         )
-
-    # create extra variables
-    n_tot_variables = n_variables
-    In, Jn, Vn = Int[], Int[], Float64[]
-
-    # this way there is a single elements per column
-    # because we assume VOV in SET and NOT AFF in SET
-    for i in eachindex(cone.soc)
-        n_vars = length(cone.soc[i])
-        vec_inds = cone.soc[i]
-        n_tot_variables += fix_duplicates!(vec_inds, n_tot_variables, In, Jn, Vn)
-        push!(con.socone, SOCSet(vec_inds, n_vars))
-        # first_ind_local += n_vars
-    end
-
-    for i in eachindex(cone.sdc)
-        n_vars = length(cone.sdc[i])
-        vec_inds = cone.sdc[i]
-        n_tot_variables += fix_duplicates!(vec_inds, n_tot_variables, In, Jn, Vn)
-        sq_side = sympackeddim(n_vars)
-        sq_len = sq_side*sq_side
-        tri_len = n_vars
+    soc_s = cone_data(src, MOI.SecondOrderCone)
+    for soc in soc_s
+        push!(con.socone, SOCSet(soc, length(soc)))
+    end
+    psc_s = cone_data(src, MOI.PositiveSemidefiniteConeTriangle)
+    for psc in psc_s
+        tri_len = length(psc)
+        sq_side = sympackeddim(tri_len)
         mat_inds = matindices(sq_side)
-        push!(con.sdpcone, SDPSet(vec_inds, mat_inds, tri_len, sq_len, sq_side))
-        # first_ind_local += n_vars
-    end
-
-    optimizer.data = nothing # Allows GC to free optimizer.data before A is loaded
-
-    #= 
-        Pre-process to remove duplicates from differente cones
-    =#
-
-    for i in 1:length(con.socone)
-        for j in i+1:length(con.socone)
-            vec1 = con.socone[i].idx
-            vec2 = con.socone[j].idx
-            n_tot_variables += fix_duplicates!(vec1, vec2, n_tot_variables, In, Jn, Vn)
-        end
+        push!(con.sdpcone, SDPSet(psc, mat_inds, tri_len, sq_side))
     end
-
-    for i in 1:length(con.sdpcone)
-        for j in i+1:length(con.sdpcone)
-            vec1 = con.sdpcone[i].vec_i
-            vec2 = con.sdpcone[j].vec_i
-            n_tot_variables += fix_duplicates!(vec1, vec2, n_tot_variables, In, Jn, Vn)
-        end
-    end
-
-    for i in 1:length(con.sdpcone)
-        for j in 1:length(con.socone)
-            vec1 = con.sdpcone[i].vec_i
-            vec2 = con.socone[j].idx
-            n_tot_variables += fix_duplicates!(vec1, vec2, n_tot_variables, In, Jn, Vn)
-        end
-    end
-
-    n_new_variables = n_tot_variables - n_variables
-
-    A2 = sparse(In .- n_variables, Jn, Vn, n_new_variables, n_tot_variables)
-    b2 = zeros(n_new_variables)
-
-    append!(aff.b, b2)
-    aff.A = hvcat((2,1), aff.A, spzeros(size(aff.A)[1], n_new_variables), A2)
-    aff.G = hcat(aff.G, spzeros(size(aff.G)[1], n_new_variables))
-    append!(aff.c, zeros(n_new_variables))
-    aff.n += n_new_variables
-    aff.p += n_new_variables
-    aff.extra = n_new_variables
-
-    In, Jn, Vn = Int[], Int[], Float64[]
-    A2 = sparse(In, Jn, Vn, 0, 0)
+    dest.cones = ConeData(
+        psc_s,
+        soc_s,
+    )
+    #
+    end # timeit
 
     #= 
         Solve modified problem
     =#
 
-    end
+    options = dest.options
 
     # warm = WarmStart()
 
@@ -580,17 +314,6 @@ function MOI.optimize!(optimizer::Optimizer)
         Logging.global_logger(global_log)
     end
 
-    #= 
-        Unload solution
-    =#
-
-    ret_val = sol.status
-    primal = sol.primal[1:aff.n-aff.extra]
-    dual_cone = sol.dual_cone[1:aff.n-aff.extra]
-    dual_eq = sol.dual_eq[1:aff.p-aff.extra]
-    slack_eq = sol.slack_eq[1:aff.p-aff.extra]
-    objval = sol.objval
-
     if options.timer_verbose
         TimerOutputs.print_timer(TimerOutputs.DEFAULT_TIMER)
         print("\n")
@@ -606,137 +329,55 @@ function MOI.optimize!(optimizer::Optimizer)
         close(f)
     end
 
-    optimizer.sol = MOISolution(
-        ret_val,
-        sol.status_string,
-        primal,
-        dual_eq,
-        sol.dual_in,
-        dual_cone,
-        slack_eq,
-        sol.slack_in,
-        sol.primal_residual,
-        sol.dual_residual,
-        (optimizer.maxsense ? -1 : 1) * objval+objective_constant,
-        (optimizer.maxsense ? -1 : 1) * sol.dual_objval+objective_constant,
-        sol.gap,
-        sol.time,
-        sol.final_rank,
-        sol.primal_feasible_user_tol,
-        sol.dual_feasible_user_tol,
-        sol.certificate_found,
-        )
-end
+    #= 
+        Fix solution
+    =#
 
-function get_indices_cone(A, rows, n_vars, first_ind_local)
-    vec_inds = zeros(Int, n_vars)
-    for (idx, col) in enumerate(first_ind_local:first_ind_local+n_vars-1)
-        # columns -> pos in cone (because is tranposed)
-        for j in nzrange(A, col)
-            row = rows[j]
-            position = idx
-            var_idx = row # global
-            vec_inds[position] = var_idx
-        end
-    end
-    return vec_inds
-end
+    sol.objval = obj_sign * sol.objval + objective_constant
+    sol.dual_objval = obj_sign * sol.dual_objval + objective_constant
 
-function fix_duplicates!(vec1::Vector{Int}, vec2::Vector{Int}, n::Int,
-                         In::Vector{Int}, Jn::Vector{Int}, Vn::Vector{Float64})
-    duplicates = intersect(vec2, vec1)
-    n_dups = length(duplicates)
-    if n_dups == 0
-        return 0
-    end
+    dest.sol = sol
 
-    append!(Jn, duplicates)
-    append!(Jn, collect(1:n_dups) .+ n)
-    append!(In, collect(1:n_dups) .+ n)
-    append!(In, collect(1:n_dups) .+ n)
-    append!(Vn,  ones(n_dups))
-    append!(Vn, -ones(n_dups))
-    cont = 1
-    for i in eachindex(vec2)
-        if vec2[i] == duplicates[cont]
-            vec2[i] = n + cont
-            cont += 1
-        end
-    end
-    @assert cont-1 == n_dups
-    return n_dups
-end
-
-function fix_duplicates!(vec::Vector{Int}, n::Int,
-                         In::Vector{Int}, Jn::Vector{Int}, Vn::Vector{Float64})
-    seen = Set{eltype(vec)}()
-    dups = Vector{eltype(vec)}()
-    new_vars = 0
-    for i in eachindex(vec)
-        x = vec[i]
-        if in(x, seen)
-            push!(dups, x)
-            new_vars += 1
-            vec[i] = new_vars + n
-        else
-            push!(seen, x)
-        end
-    end
-    n_dups = length(dups)
-    if n_dups == 0
-        return 0
-    end
-    append!(Jn, dups)
-    append!(Jn, collect(1:n_dups) .+ n)
-    append!(In, collect(1:n_dups) .+ n)
-    append!(In, collect(1:n_dups) .+ n)
-    append!(Vn,  ones(n_dups))
-    append!(Vn, -ones(n_dups))
-    return n_dups
-end
-
-function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
-    n = sympackeddim(length(v))
-    n1, n2 = size(out)
-    @assert n == n1 == n2
-    c = 0
-    for j in 1:n, i in 1:j
-        c += 1
-        out[i,j] = v[c]
-    end
-    return out
+    return
 end
-function ivech(v::AbstractVector{T}) where T
-    n = sympackeddim(length(v))
-    out = zeros(n, n)
-    ivech!(out, v)
-    
-    return out
+
+function MOI.optimize!(dest::Optimizer, src::OptimizerCache)
+    _optimize!(dest, src)
+    return MOI.Utilities.identity_index_map(src), false
 end
 
-ivec(X) = Matrix(Symmetric(ivech(X),:U))
+function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
+    cache = OptimizerCache{Float64}()
+    index_map = MOI.copy_to(cache, src)
+    _optimize!(dest, cache)
+    return index_map, false
+end
 
 #=
-    Status
+    Attributes
 =#
 
-function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
-    return optimizer.sol.time
+MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) = optimizer.sol.time
+
+"""
+    PDHGIterations()
+
+The number of PDHG iterations completed during the solve.
+"""
+struct PDHGIterations <: MOI.AbstractModelAttribute end
+
+MOI.is_set_by_optimize(::PDHGIterations) = true
+
+function MOI.get(optimizer::Optimizer, ::PDHGIterations)
+    return Int64(optimizer.sol.iter)
 end
 
 function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
-    return optimizer.sol.raw_status
+    return optimizer.sol.status_string
 end
 
-# Implements getter for result value and statuses
-# ProxSDP returns one of the following integers:
-# 0 - MOI.OPTIMIZE_NOT_CALLED
-# 1 - MOI.OPTIMAL
-# 2 - MOI.TIME_LIMIT
-# 3 - MOI.ITERATION_LIMIT
-# 4 - MOI.INFEASIBLE_OR_UNBOUNDED
 function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
-    s = optimizer.sol.ret_val
+    s = optimizer.sol.status
     @assert 0 <= s <= 6
     if s == 0
         return MOI.OPTIMIZE_NOT_CALLED
@@ -756,19 +397,20 @@ function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
 end
 
 function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue)
-    result_count_err(optimizer, attr)
+    MOI.check_result_index_bounds(optimizer, attr)
     value = optimizer.sol.objval
     return value
 end
+
 function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue)
-    result_count_err(optimizer, attr)
+    MOI.check_result_index_bounds(optimizer, attr)
     value = optimizer.sol.dual_objval
     return value
 end
 
 function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
-    s = optimizer.sol.ret_val
-    if attr.N > 1 || s == 0
+    s = optimizer.sol.status
+    if attr.result_index > 1 || s == 0
         return MOI.NO_SOLUTION
     end
     if s == 5 && optimizer.sol.certificate_found
@@ -780,9 +422,10 @@ function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus)
         return MOI.INFEASIBLE_POINT
     end
 end
+
 function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
-    s = optimizer.sol.ret_val
-    if attr.N > 1 || s ==0
+    s = optimizer.sol.status
+    if attr.result_index > 1 || s ==0
         return MOI.NO_SOLUTION
     end
     if s == 6 && optimizer.sol.certificate_found
@@ -795,95 +438,91 @@ function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus)
     end
 end
 
-function result_count_err(optimizer::Optimizer, attr)
-    if optimizer.sol.ret_val == 0
-        throw(MOI.ResultIndexBoundsError{typeof(attr)}(attr, 0))
-    elseif result_index(attr) > 1
-        throw(MOI.ResultIndexBoundsError{typeof(attr)}(attr, 1))
-    end
-    return nothing
+function MOI.get(optimizer::Optimizer, ::MOI.ResultCount)
+    return optimizer.sol.result_count
 end
 
-result_index(a::MOI.ObjectiveValue) = a.result_index
-result_index(a::MOI.DualObjectiveValue) = a.result_index
-result_index(a::MOI.ConstraintDual) = a.N
-result_index(a::MOI.ConstraintPrimal) = a.N
-result_index(a::MOI.VariablePrimal) = a.N
-
 #=
-    Solution
+    Results
 =#
 
-function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::VI)
-    result_count_err(optimizer, attr)
-    optimizer.sol.primal[vi.value]
-end
-function MOI.get(optimizer::Optimizer, a::MOI.VariablePrimal, vi::Vector{VI})
-    result_count_err(optimizer, a)
-    return MOI.get.(optimizer, a, vi)
-end
-
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
-                 ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
-    result_count_err(optimizer, attr)
-    cone = optimizer.cone
-    inds = (cone.eq_start[ci.value]-1) .+ (1:cone.eq_len[ci.value])
-    return optimizer.sol.slack_eq[inds]
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.VariablePrimal,
+    vi::MOI.VariableIndex,
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.sol.primal[vi.value]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintPrimal,
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Zeros},
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.sol.slack_eq[_rows(optimizer, ci)]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintPrimal,
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Nonpositives},
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.sol.slack_in[_rows(optimizer, ci)]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintPrimal,
+    ci::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.sol.primal[optimizer.cones.psc[ci.value]]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintPrimal,
+    ci::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.SecondOrderCone}
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.sol.primal[optimizer.cones.soc[ci.value]]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},MOI.Zeros},
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return -optimizer.sol.dual_eq[_rows(optimizer, ci)]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},MOI.Nonpositives},
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return -optimizer.sol.dual_in[_rows(optimizer, ci)]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.sol.dual_cone[optimizer.cones.psc[ci.value]]
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.SecondOrderCone}
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.sol.dual_cone[optimizer.cones.soc[ci.value]]
 end
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
-    ci::CI{<:MOI.AbstractFunction, MOI.Nonpositives})
-    result_count_err(optimizer, attr)
-    cone = optimizer.cone
-    inds = (cone.in_start[ci.value]-1) .+ (1:cone.in_len[ci.value])
-    return optimizer.sol.slack_in[inds]
-end
-
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
-        ci::CI{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle})
-    result_count_err(optimizer, attr)
-    optimizer.sol.primal[optimizer.cone.sdc[ci.value]]
-end
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
-    ci::CI{MOI.VectorOfVariables, MOI.SecondOrderCone})
-    result_count_err(optimizer, attr)
-    optimizer.sol.primal[optimizer.cone.soc[ci.value]]
-end
-
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
-    ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
-    result_count_err(optimizer, attr)
-    cone = optimizer.cone
-    inds = (cone.eq_start[ci.value]-1) .+ (1:cone.eq_len[ci.value])
-    return -optimizer.sol.dual_eq[inds]
-end
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
-ci::CI{<:MOI.AbstractFunction, MOI.Nonpositives})
-    result_count_err(optimizer, attr)
-    cone = optimizer.cone
-    inds = (cone.in_start[ci.value]-1) .+ (1:cone.in_len[ci.value])
-    return -optimizer.sol.dual_in[inds]
-end
-
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
-    ci::CI{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle})
-    result_count_err(optimizer, attr)
-    optimizer.sol.dual_cone[optimizer.cone.sdc[ci.value]]
-end
-function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
-    ci::CI{MOI.VectorOfVariables, MOI.SecondOrderCone})
-    result_count_err(optimizer, attr)
-    optimizer.sol.dual_cone[optimizer.cone.soc[ci.value]]
-end
-
-function MOI.get(optimizer::Optimizer, ::MOI.ResultCount)
-    if MOI.get(optimizer, MOI.TerminationStatus()) in [
-                # MOI.INFEASIBLE_OR_UNBOUNDED,
-                # MOI.INFEASIBLE,
-                # MOI.DUAL_INFEASIBLE,
-                MOI.OPTIMIZE_NOT_CALLED
-        ]
-        return 0
-    else
-        return 1
-    end
-end
\ No newline at end of file
diff --git a/src/ProxSDP.jl b/src/ProxSDP.jl
index 629a83d..4fcd5e2 100644
--- a/src/ProxSDP.jl
+++ b/src/ProxSDP.jl
@@ -1,32 +1,30 @@
 module ProxSDP
 
-    using Arpack
-    using KrylovKit
-    using MathOptInterface
-    using TimerOutputs
-
-    using Printf
-    using SparseArrays
-    using LinearAlgebra
-    import Random
-    import Logging
-
-    import LinearAlgebra: BlasInt
-
-    include("structs.jl")
-    include("util.jl")
-    include("printing.jl")
-    include("scaling.jl")
-    include("equilibration.jl")
-    include("pdhg.jl")
-    include("residuals.jl")
-    include("eigsolver.jl")
-    include("prox_operators.jl")
-
-    include("MOI_wrapper.jl")
-
-    MOIU.@model _ProxSDPModelData () (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan) (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.PositiveSemidefiniteConeTriangle) () (MOI.SingleVariable,) (MOI.ScalarAffineFunction,) (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)
-
-    Solver(;args...) = MOIU.CachingOptimizer(_ProxSDPModelData{Float64}(), ProxSDP.Optimizer(args))
+import Arpack
+import KrylovKit
+import MathOptInterface
+import TimerOutputs
+import TimerOutputs: @timeit
+
+import Printf
+import SparseArrays
+import LinearAlgebra
+import LinearAlgebra: BlasInt
+import Random
+import Logging
+
+
+include("structs.jl")
+include("options.jl")
+include("util.jl")
+include("printing.jl")
+include("scaling.jl")
+include("equilibration.jl")
+include("pdhg.jl")
+include("residuals.jl")
+include("eigsolver.jl")
+include("prox_operators.jl")
+
+include("MOI_wrapper.jl")
 
 end
\ No newline at end of file
diff --git a/src/eigsolver.jl b/src/eigsolver.jl
index bfbd313..059417d 100644
--- a/src/eigsolver.jl
+++ b/src/eigsolver.jl
@@ -354,7 +354,7 @@ mutable struct EigSolverAlloc{T}
     bmat::String
     which::String
     mode::Int
-    Amat::Symmetric{T,Matrix{T}}
+    Amat::LinearAlgebra.Symmetric{T,Matrix{T}}
     lworkl::Int
     TOL::Base.RefValue{T}
     v::Matrix{T}
@@ -384,7 +384,8 @@ hasconverged(arc::EigSolverAlloc) = arc.converged
 
 function EigSolverAlloc(T::DataType, n::Int64, opt::Options)::EigSolverAlloc
     arc = EigSolverAlloc{T}()
-    @timeit "init_arc" _init_eigsolver!(arc, Symmetric(Matrix{T}(I, n, n), :U), 1, n, opt)
+    @timeit "init_arc" _init_eigsolver!(
+        arc, LinearAlgebra.Symmetric(Matrix{T}(LinearAlgebra.I, n, n), :U), 1, n, opt)
     return arc
 end
 
@@ -409,7 +410,7 @@ function eigsolver_update_resid!(arc, opt, init)
     return nothing
 end
 
-function _init_eigsolver!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T
+function _init_eigsolver!(arc::EigSolverAlloc{T}, A::LinearAlgebra.Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T
     if opt.eigsolver == 1
         arpack_init!(arc, A, nev, n, opt)
     else
@@ -426,7 +427,7 @@ function arpack_getvectors(arc::EigSolverAlloc)
     return arc.v
 end
 
-function arpack_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T
+function arpack_init!(arc::EigSolverAlloc{T}, A::LinearAlgebra.Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T
 
     # IDO - Integer.  (INPUT/OUTPUT)  (in julia its is a integer vector)
     # Reverse communication flag.
@@ -568,7 +569,7 @@ function arpack_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::In
     return nothing
 end
 
-function arpack_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, opt::Options, up_ncv::Bool) where T
+function arpack_update!(arc::EigSolverAlloc{T}, A::LinearAlgebra.Symmetric{T,Matrix{T}}, nev::Int64, opt::Options, up_ncv::Bool) where T
 
     need_resize = up_ncv
     if !need_resize
@@ -674,14 +675,14 @@ function _saupd!(arc::EigSolverAlloc{T})::Nothing where T
             @simd for i in 1:arc.n
                 @inbounds arc.x[i] = arc.workd[i-1+arc.ipntr[1]]
             end
-            mul!(arc.y, arc.Amat, arc.x)
+            LinearAlgebra.mul!(arc.y, arc.Amat, arc.x)
             @simd for i in 1:arc.n
                  @inbounds arc.workd[i-1+arc.ipntr[2]] = arc.y[i]
             end
             # using views
             # x = view(arc.workd, arc.ipntr[1] .+ arc.zernm1)
             # y = view(arc.workd, arc.ipntr[2] .+ arc.zernm1)
-            # mul!(y, arc.Amat, x)
+            # LinearAlgebra.mul!(y, arc.Amat, x)
         elseif arc.ido[] == 99
             # In this case, don't call _EUPD! (according to https://help.scilab.org/docs/5.3.3/en_US/dseupd.html)
             break
@@ -744,7 +745,7 @@ function _seupd!(arc::EigSolverAlloc{T})::Nothing where T
     return nothing
 end
 
-function arpack_eig!(arc::EigSolverAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T
+function arpack_eig!(arc::EigSolverAlloc, A::LinearAlgebra.Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T
 
     up_ncv = false
     for i in 1:1
@@ -776,7 +777,7 @@ function krylovkit_getvector(arc::EigSolverAlloc, i)
     return arc.vecs[i]
 end
 
-function krylovkit_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T
+function krylovkit_init!(arc::EigSolverAlloc{T}, A::LinearAlgebra.Symmetric{T,Matrix{T}}, nev::Int64, n::Int64, opt::Options) where T
     arc.n = n
     arc.nev = nev
     eigsolver_init_resid!(arc, opt, n)
@@ -785,7 +786,7 @@ function krylovkit_init!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev:
     return nothing
 end
 
-function krylovkit_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, nev::Int64, opt::Options) where T
+function krylovkit_update!(arc::EigSolverAlloc{T}, A::LinearAlgebra.Symmetric{T,Matrix{T}}, nev::Int64, opt::Options) where T
     arc.nev = nev
     if opt.krylovkit_reset_resid
         eigsolver_update_resid!(arc, opt, opt.krylovkit_resid_init)
@@ -794,7 +795,7 @@ function krylovkit_update!(arc::EigSolverAlloc{T}, A::Symmetric{T,Matrix{T}}, ne
     return nothing
 end
 
-function krylovkit_eig!(arc::EigSolverAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T
+function krylovkit_eig!(arc::EigSolverAlloc, A::LinearAlgebra.Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T
     arc.converged = true
     @timeit "update_arc" krylovkit_update!(arc, A, nev, opt)::Nothing
     @timeit "krylovkit" begin
diff --git a/src/equilibration.jl b/src/equilibration.jl
index dd06ec7..be390b2 100644
--- a/src/equilibration.jl
+++ b/src/equilibration.jl
@@ -1,6 +1,3 @@
-
-using LinearAlgebra
-
 function equilibrate!(M, aff, opt)
     max_iters=opt.equilibration_iters
     lb = opt.equilibration_lb
@@ -16,11 +13,11 @@ function equilibrate!(M, aff, opt)
         u_, v_         = zeros(aff.m + aff.p), zeros(aff.n)
         u_grad, v_grad = zeros(aff.m + aff.p), zeros(aff.n)
         row_norms, col_norms = zeros(aff.m + aff.p), zeros(aff.n)
-        E = Diagonal(u)
-        D = Diagonal(v)
+        E = LinearAlgebra.Diagonal(u)
+        D = LinearAlgebra.Diagonal(v)
         M_ = copy(M)
 
-        rows_M_ = rowvals(M_)
+        rows_M_ = SparseArrays.rowvals(M_)
     end
 
     for iter in 1:max_iters
@@ -29,8 +26,8 @@ function equilibrate!(M, aff, opt)
             D.diag .= exp.(v)
         end
         @timeit "M_" begin 
-            mul!(M_, M, D)
-            mul!(M_, E, M_)
+            LinearAlgebra.mul!(M_, M, D)
+            LinearAlgebra.mul!(M_, E, M_)
         end
 
         step_size = 2. / (γ * (iter + 1.))
@@ -40,7 +37,7 @@ function equilibrate!(M, aff, opt)
             fill!(row_norms, 0.0)
             fill!(col_norms, 0.0)
             for col in 1:aff.n
-                for line in nzrange(M_, col)
+                for line in SparseArrays.nzrange(M_, col)
                     row_norms[rows_M_[line]] += abs2(M_[rows_M_[line], col])
                     col_norms[col] += abs2(M_[rows_M_[line], col])
                 end
diff --git a/src/options.jl b/src/options.jl
new file mode 100644
index 0000000..4c4da4a
--- /dev/null
+++ b/src/options.jl
@@ -0,0 +1,132 @@
+Base.@kwdef mutable struct Options
+
+    # Printing options
+    log_verbose::Bool = false
+    log_freq::Int = 1000
+    timer_verbose::Bool = false
+    timer_file::Bool = false
+    disable_julia_logger::Bool = true
+
+    # time options
+    time_limit::Float64 = 3600_00. #100 hours
+
+    warn_on_limit::Bool = false
+    extended_log::Bool = false
+    extended_log2::Bool = false
+    log_repeat_header::Bool = false
+
+    # Default tolerances
+    tol_gap::Float64 = 1e-4
+    tol_feasibility::Float64 = 1e-4
+    tol_feasibility_dual::Float64 = 1e-4
+    tol_primal::Float64 = 1e-4
+    tol_dual::Float64 = 1e-4
+    tol_psd::Float64 = 1e-7
+    tol_soc::Float64 = 1e-7
+
+    check_dual_feas::Bool = false
+    check_dual_feas_freq::Int = 1000
+
+    max_obj::Float64 = 1e20
+    min_iter_max_obj::Int = 10
+
+    # infeasibility check
+    min_iter_time_infeas::Int = 1000
+    infeas_gap_tol::Float64 = 1e-4
+    infeas_limit_gap_tol::Float64 = 1e-1
+    infeas_stable_gap_tol::Float64 = 1e-4
+    infeas_feasibility_tol::Float64 = 1e-4
+    infeas_stable_feasibility_tol::Float64 = 1e-8
+
+    certificate_search::Bool = true
+    certificate_obj_tol::Float64 = 1e-1
+    certificate_fail_tol::Float64 = 1e-8
+
+    # Bounds on beta (dual_step / primal_step) [larger bounds may lead to numerical inaccuracy]
+    min_beta::Float64 = 1e-5
+    max_beta::Float64 = 1e+5
+    initial_beta::Float64 = 1.
+
+    # Adaptive primal-dual steps parameters [adapt_decay above .7 may lead to slower convergence]
+    initial_adapt_level::Float64 = .9
+    adapt_decay::Float64 = .8
+    adapt_window::Int64 = 50
+
+    # PDHG parameters
+    convergence_window::Int = 200
+    convergence_check::Int = 50
+    max_iter::Int = 0
+    min_iter::Int = 40
+    divergence_min_update::Int = 50
+    max_iter_lp::Int = 10_000_000
+    max_iter_conic::Int = 1_000_000
+    max_iter_local::Int = 0 #ignores user setting
+
+    advanced_initialization::Bool = true
+
+    # Linesearch parameters
+    line_search_flag::Bool = true
+    max_linsearch_steps::Int = 5000
+    delta::Float64 = .9999
+    initial_theta::Float64 = 1.
+    linsearch_decay::Float64 = .75
+
+    # Spectral decomposition parameters
+    full_eig_decomp::Bool = false
+    max_target_rank_krylov_eigs::Int = 16
+    min_size_krylov_eigs::Int = 100
+    warm_start_eig::Bool = true
+    rank_increment::Int = 1 # 0=multiply, 1 = add
+    rank_increment_factor::Int = 1 # 0 multiply, 1 = add
+
+    # eigsolver selection
+    #=
+        1: Arpack [dsaupd] (tipically non-deterministic)
+        2: KrylovKit [eigsolve/Lanczos] (DEFAULT)
+    =#
+    eigsolver::Int = 2
+    eigsolver_min_lanczos::Int = 25
+    eigsolver_resid_seed::Int = 1234
+
+    # Arpack
+    # note that Arpack is Non-deterministic
+    # (https://github.com/mariohsouto/ProxSDP.jl/issues/69)
+    arpack_tol::Float64 = 1e-10
+    #=
+        0: arpack random [usually faster - NON-DETERMINISTIC - slightly]
+        1: all ones [???]
+        2: julia random uniform (eigsolver_resid_seed) [medium for DETERMINISTIC]
+        3: julia normalized random normal (eigsolver_resid_seed) [best for DETERMINISTIC]
+    =#
+    arpack_resid_init::Int = 3
+    arpack_reset_resid::Bool = true # true for determinism
+    # larger is more stable to converge and more deterministic
+    arpack_max_iter::Int = 10_000
+    # see remark for of dsaupd
+
+    # KrylovKit
+    krylovkit_reset_resid::Bool = false
+    krylovkit_resid_init::Int = 3
+    krylovkit_tol::Float64 = 1e-12
+    krylovkit_max_iter::Int = 100
+    krylovkit_eager::Bool = false
+    krylovkit_verbose::Int = 0
+
+    # Reduce rank [warning: heuristics]
+    reduce_rank::Bool = false
+    rank_slack::Int = 3
+
+    full_eig_freq::Int = 10_000_000
+    full_eig_len::Int = 0
+
+    # equilibration parameters
+    equilibration::Bool = false
+    equilibration_iters::Int = 1000
+    equilibration_lb::Float64 = -10.0
+    equilibration_ub::Float64 = +10.0
+    equilibration_limit::Float64 = 0.9
+    equilibration_force::Bool = false
+
+    # spectral norm [using exact norm via svds may result in nondeterministic behavior]
+    approx_norm::Bool = true
+end
diff --git a/src/pdhg.jl b/src/pdhg.jl
index a94f171..6e8a446 100644
--- a/src/pdhg.jl
+++ b/src/pdhg.jl
@@ -1,4 +1,8 @@
-function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CPResult
+function chambolle_pock(
+    affine_sets::AffineSets,
+    conic_sets::ConicSets,
+    opt
+)::Result
 
     # Initialize parameters
     p = Params()
@@ -7,9 +11,9 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
     p.window = opt.convergence_window
     p.beta = opt.initial_beta
     p.time0 = time()
-    p.norm_b = norm(affine_sets.b, 2)
-    p.norm_h = norm(affine_sets.h, 2)
-    p.norm_c = norm(affine_sets.c, 2)
+    p.norm_b = LinearAlgebra.norm(affine_sets.b, 2)
+    p.norm_h = LinearAlgebra.norm(affine_sets.h, 2)
+    p.norm_c = LinearAlgebra.norm(affine_sets.c, 2)
     p.rank_update, p.stop_reason, p.update_cont = 0, 0, 0
     p.stop_reason_string = "Not optimized"
     p.target_rank = 2 * ones(length(conic_sets.sdpcone))
@@ -20,8 +24,10 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
     p.certificate_search = false
     p.certificate_search_min_iter = 0
     p.certificate_found = false
-    sol = Array{CPResult}(undef, 0)
-    arc_list = [EigSolverAlloc(Float64, sdp.sq_side, opt) for (idx, sdp) in enumerate(conic_sets.sdpcone)]
+    sol = Array{Result}(undef, 0)
+    arc_list = [
+        EigSolverAlloc(Float64, sdp.sq_side, opt)
+            for (idx, sdp) in enumerate(conic_sets.sdpcone)]
     ada_count = 0
 
     if opt.max_iter <= 0
@@ -69,8 +75,8 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
                 opt.equilibration = true
             end
             if opt.equilibration
-                @timeit "equilibrate inner" E, D = equilibrate!(M, affine_sets, opt)
-                @timeit "equilibrate scaling" begin
+                @timeit "equilib inner" E, D = equilibrate!(M, affine_sets, opt)
+                @timeit "equilib scaling" begin
                     M = E * M * D
                     affine_sets.A = M[1:affine_sets.p, :]
                     affine_sets.G = M[affine_sets.p + 1:end, :]
@@ -80,8 +86,8 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
                     affine_sets.c = D * affine_sets.c
                 end
             else
-                E = Diagonal(zeros(affine_sets.m + affine_sets.p))
-                D = Diagonal(zeros(affine_sets.n))
+                E = LinearAlgebra.Diagonal(zeros(affine_sets.m + affine_sets.p))
+                D = LinearAlgebra.Diagonal(zeros(affine_sets.n))
             end
         end
         
@@ -105,14 +111,17 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
                     spectral_norm = Arpack.svds(M, nsv=1)[1].S[1]
                 catch
                     println("    WARNING: Failed to compute spectral norm of M, shifting to Frobenius norm")
-                    spectral_norm = norm(M)
+                    spectral_norm = LinearAlgebra.norm(M)
                 end
             else
                 F = LinearAlgebra.svd!(Matrix(M))
                 spectral_norm = maximum(F.S)
             end
         else
-            spectral_norm = norm(M)
+            spectral_norm = LinearAlgebra.norm(M)
+        end
+        if spectral_norm < 1e-10
+            spectral_norm = 1.0
         end
 
         # Build struct for storing matrices
@@ -128,8 +137,8 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
     # Initialization
     if opt.advanced_initialization
         pair.x .= p.primal_step .* mat.c
-        mul!(a.Mx, mat.M, pair.x)
-        mul!(a.Mx_old, mat.M, pair.x_old)
+        LinearAlgebra.mul!(a.Mx, mat.M, pair.x)
+        LinearAlgebra.mul!(a.Mx_old, mat.M, pair.x_old)
     end
 
     # Fixed-point loop
@@ -419,7 +428,7 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
             residuals.feasibility[p.iter] > opt.infeas_feasibility_tol &&
             max_abs_diff(residuals.feasibility) < opt.infeas_stable_feasibility_tol
             )
-            
+
             p.stop_reason = 6 # Infeasible
             p.stop_reason_string = "Infeasible: feasibility stalled at $(residuals.feasibility[p.iter])"
 
@@ -499,7 +508,7 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
     if opt.certificate_search && p.certificate_search
 
         @assert length(sol) == 1
-        
+
         if p.certificate_found
 
             if p.stop_reason == 6
@@ -508,7 +517,7 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
 
             pop!(sol)
             push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
-                c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))    
+                c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))
         end
 
     else
@@ -520,7 +529,15 @@ function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CP
     return sol[1]
 end
 
-function linesearch!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets, mat::Matrices, opt::Options, p::Params)::Nothing
+function linesearch!(
+    pair::PrimalDual,
+    a::AuxiliaryData,
+    affine_sets::AffineSets,
+    mat::Matrices,
+    opt::Options,
+    p::Params
+)::Nothing
+
     p.primal_step = p.primal_step * sqrt(1. + p.theta)
 
     for i in 1:opt.max_linsearch_steps
@@ -536,14 +553,14 @@ function linesearch!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets
             a.y_temp .-= (p.beta * p.primal_step) .* a.y_half
         end
 
-        @timeit "linesearch 3" mul!(a.Mty, mat.Mt, a.y_temp)
+        @timeit "linesearch 3" LinearAlgebra.mul!(a.Mty, mat.Mt, a.y_temp)
         
         # In-place norm
         @timeit "linesearch 4" begin
             a.Mty .-= a.Mty_old
             a.y_temp .-= pair.y_old
-            y_norm = norm(a.y_temp)
-            Mty_norm = norm(a.Mty)
+            y_norm = LinearAlgebra.norm(a.y_temp)
+            Mty_norm = LinearAlgebra.norm(a.Mty)
         end
 
         if sqrt(p.beta) * p.primal_step * Mty_norm <= opt.delta * y_norm
@@ -564,7 +581,14 @@ function linesearch!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets
     return nothing
 end
 
-function dual_step!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets, mat::Matrices, opt::Options, p::Params)::Nothing
+function dual_step!(
+    pair::PrimalDual,
+    a::AuxiliaryData,
+    affine_sets::AffineSets,
+    mat::Matrices,
+    opt::Options,
+    p::Params
+)::Nothing
 
     @timeit "dual step 1" begin
         a.y_half .= pair.y .+ p.dual_step * (2. * a.Mx .- a.Mx_old)
@@ -576,7 +600,7 @@ function dual_step!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets,
         a.y_temp .-= p.dual_step * a.y_half
     end
 
-    @timeit "linesearch 3" mul!(a.Mty, mat.Mt, a.y_temp)
+    @timeit "linesearch 3" LinearAlgebra.mul!(a.Mty, mat.Mt, a.y_temp)
 
     copyto!(pair.y, a.y_temp)
     p.primal_step_old = p.primal_step
@@ -584,7 +608,16 @@ function dual_step!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets,
     return nothing
 end
 
-function primal_step!(pair::PrimalDual, a::AuxiliaryData, cones::ConicSets, mat::Matrices, opt::Options, p::Params, arc_list, iter::Int64)::Nothing
+function primal_step!(
+    pair::PrimalDual,
+    a::AuxiliaryData,
+    cones::ConicSets,
+    mat::Matrices,
+    opt::Options,
+    p::Params,
+    arc_list,
+    iter::Int64,
+)::Nothing
 
     pair.x .-= p.primal_step .* (a.Mty .+ mat.c)
 
@@ -598,7 +631,7 @@ function primal_step!(pair::PrimalDual, a::AuxiliaryData, cones::ConicSets, mat:
         @timeit "soc proj" soc_projection!(pair.x, a, cones, opt, p)
     end
 
-    @timeit "linesearch -1" mul!(a.Mx, mat.M, pair.x)
+    @timeit "linesearch -1" LinearAlgebra.mul!(a.Mx, mat.M, pair.x)
 
     return nothing
 end
@@ -606,9 +639,9 @@ end
 function certificate_dual_infeasibility(affine_sets, p, opt)
 
     if opt.log_verbose
-        println("---------------------------------------------------------------------------------------")
+        println("-"^87)
         println("    Begin search for dual infeasibility certificate")
-        println("---------------------------------------------------------------------------------------")
+        println("-"^87)
     end
 
     fill!(affine_sets.b, 0.0)
@@ -622,9 +655,9 @@ end
 function certificate_infeasibility(affine_sets, p, opt)
 
     if opt.log_verbose
-        println("---------------------------------------------------------------------------------------")
+        println("-"^87)
         println("    Begin search for infeasibility certificate")
-        println("---------------------------------------------------------------------------------------")
+        println("-"^87)
     end
 
     fill!(affine_sets.c, 0.0)
@@ -649,7 +682,7 @@ function cone_feas(v, cones, a, num = sqrt(2))
         if sdp.sq_side == 1
             sdp_viol = max(sdp_viol, -min(0.0, a.m[idx][1]))
         else
-            fact = eigen!(a.m[idx])
+            fact = LinearAlgebra.eigen!(a.m[idx])
             sdp_viol = max(sdp_viol, -min(0.0, minimum(fact.values)))
         end
     end
@@ -659,7 +692,7 @@ function cone_feas(v, cones, a, num = sqrt(2))
         len = soc.len
         push!(a.soc_s, view(v, cont + 1))
         s = v[cont+1]
-        sdp_viol = max(sdp_viol, -min(0.0, s - norm(view(v, cont + 2:cont + len))))
+        sdp_viol = max(sdp_viol, -min(0.0, s - LinearAlgebra.norm(view(v, cont + 2:cont + len))))
         cont += len
     end
     return max(sdp_viol, soc_viol), cont
@@ -729,7 +762,7 @@ function cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
 
     dual_feasibility = dual_feas(dual_in, dual_cone, conic_sets, a)
 
-    return CPResult(
+    return Result(
         p.stop_reason,
         p.stop_reason_string,
         pair.x[var_ordering],
@@ -744,9 +777,11 @@ function cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
         residuals.dual_obj[p.iter],
         residuals.dual_gap[p.iter],
         time() - p.time0,
+        p.iter,
         sum(p.current_rank),
         residuals.feasibility[p.iter] <= opt.tol_feasibility,
         dual_feasibility <= opt.tol_feasibility_dual,
         p.certificate_found,
+        1,
     )
 end
\ No newline at end of file
diff --git a/src/printing.jl b/src/printing.jl
index 079069a..446ea68 100644
--- a/src/printing.jl
+++ b/src/printing.jl
@@ -96,21 +96,21 @@ end
 function print_progress(residuals::Residuals, p::Params, opt, val = -1.0)
     primal_res = residuals.primal_residual[p.iter]
     dual_res = residuals.dual_residual[p.iter]
-    s_k = @sprintf("%d", p.iter)
+    s_k = Printf.@sprintf("%d", p.iter)
     s_k *= " |"
-    s_s = @sprintf("%.5f", residuals.dual_gap[p.iter])
+    s_s = Printf.@sprintf("%.5f", residuals.dual_gap[p.iter])
     s_s *= " |"
-    s_o = @sprintf("%.3f", residuals.prim_obj[p.iter])
+    s_o = Printf.@sprintf("%.3f", residuals.prim_obj[p.iter])
     s_o *= " |"
-    s_f = @sprintf("%.5f", residuals.feasibility[p.iter])
+    s_f = Printf.@sprintf("%.5f", residuals.feasibility[p.iter])
     s_f *= " |"
-    s_p = @sprintf("%.5f", primal_res)
+    s_p = Printf.@sprintf("%.5f", primal_res)
     s_p *= " |"
-    s_d = @sprintf("%.5f", dual_res)
+    s_d = Printf.@sprintf("%.5f", dual_res)
     s_d *= " |"
-    s_target_rank = @sprintf("%.0f", sum(p.target_rank))
+    s_target_rank = Printf.@sprintf("%.0f", sum(p.target_rank))
     s_target_rank *= " |"
-    s_time = @sprintf("%.4f", time() - p.time0)
+    s_time = Printf.@sprintf("%.4f", time() - p.time0)
     s_time *= " |"
     a = "|"
     a *= " "^max(0, 9 - length(s_k))
@@ -131,12 +131,12 @@ function print_progress(residuals::Residuals, p::Params, opt, val = -1.0)
     a *= s_time
 
     if opt.extended_log || opt.extended_log2
-        str = @sprintf("%.3f", residuals.dual_obj[p.iter]) * " |"
+        str = Printf.@sprintf("%.3f", residuals.dual_obj[p.iter]) * " |"
         str = " "^max(0, 11 - length(str)) * str
         a *= str
     end
     if opt.extended_log2
-        str = @sprintf("%.5f", val) * " |"
+        str = Printf.@sprintf("%.5f", val) * " |"
         str = " "^max(0, 11 - length(str)) * str
         a *= str
     end
diff --git a/src/prox_operators.jl b/src/prox_operators.jl
index 6d6637a..79992fb 100644
--- a/src/prox_operators.jl
+++ b/src/prox_operators.jl
@@ -110,7 +110,7 @@ end
 
 function full_eig!(a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing
     p.current_rank[idx] = 0
-    fact = eigen!(a.m[idx])
+    fact = LinearAlgebra.eigen!(a.m[idx])
     p.min_eig[idx] = 0.0 #minimum(fact.values)
     fill!(a.m[idx].data, 0.)
     for i in 1:length(fact.values)
@@ -128,7 +128,7 @@ end
 function min_eig(a::AuxiliaryData, idx::Int, p::Params)
     # if p.min_eig[idx] == -Inf
     #     @timeit "bad min eig" begin
-    #         fact = eigen!(a.m[idx])
+    #         fact = LinearAlgebra.eigen!(a.m[idx])
     #         p.min_eig[idx] = minimum(fact.values)
     #     end
     # end
@@ -143,7 +143,7 @@ function soc_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets,
 end
 
 function soc_projection!(v::ViewVector, s::ViewScalar)::Nothing
-    nv = norm(v, 2)
+    nv = LinearAlgebra.norm(v, 2)
     if nv <= -s[]
         s[] = 0.
         v .= 0.
diff --git a/src/residuals.jl b/src/residuals.jl
index bdbd971..89ed782 100644
--- a/src/residuals.jl
+++ b/src/residuals.jl
@@ -19,13 +19,13 @@ function compute_gap!(residuals::Residuals, pair::PrimalDual, a::AuxiliaryData,
     residuals.feasibility[p.iter] = max(residuals.equa_feasibility, residuals.ineq_feasibility)
 
     # Primal-dual gap
-    residuals.prim_obj[p.iter] = dot(aff.c, pair.x)
+    residuals.prim_obj[p.iter] = LinearAlgebra.dot(aff.c, pair.x)
     residuals.dual_obj[p.iter] = 0.
     if aff.p > 0
-        residuals.dual_obj[p.iter] -= dot(aff.b, @view pair.y[1:aff.p])
+        residuals.dual_obj[p.iter] -= LinearAlgebra.dot(aff.b, @view pair.y[1:aff.p])
     end
     if aff.m > 0
-        residuals.dual_obj[p.iter] -= dot(aff.h, @view pair.y[aff.p+1:end])
+        residuals.dual_obj[p.iter] -= LinearAlgebra.dot(aff.h, @view pair.y[aff.p+1:end])
     end
     residuals.dual_gap[p.iter] =
         abs(residuals.prim_obj[p.iter] - residuals.dual_obj[p.iter]) /
@@ -43,7 +43,9 @@ function compute_residual!(residuals::Residuals, pair::PrimalDual, a::AuxiliaryD
     pair.x_old .= pair.x .- p.primal_step .* a.Mty
     # Px - Px_old
     pair.x_old .-= a.Mty_old 
-    residuals.primal_residual[p.iter] = sqrt(aff.n) * norm(pair.x_old, Inf) / max(norm(a.Mty_old, Inf), p.norm_b, p.norm_h, 1.)
+    residuals.primal_residual[p.iter] =
+        sqrt(aff.n) * LinearAlgebra.norm(pair.x_old, Inf) /
+            max(LinearAlgebra.norm(a.Mty_old, Inf), p.norm_b, p.norm_h, 1.)
 
     # Dual residual
     # Py_old
@@ -52,7 +54,9 @@ function compute_residual!(residuals::Residuals, pair::PrimalDual, a::AuxiliaryD
     pair.y_old .= pair.y .- p.dual_step .* a.Mx
     # Py - Py_old
     pair.y_old .-= a.Mx_old
-    residuals.dual_residual[p.iter] = sqrt(aff.m + aff.p) * norm(pair.y_old, Inf) / max(norm(a.Mx_old, Inf), p.norm_c, 1.)
+    residuals.dual_residual[p.iter] =
+        sqrt(aff.m + aff.p) * LinearAlgebra.norm(pair.y_old, Inf) /
+            max(LinearAlgebra.norm(a.Mx_old, Inf), p.norm_c, 1.)
 
     # Compute combined residual
     residuals.comb_residual[p.iter] = max(residuals.primal_residual[p.iter], residuals.dual_residual[p.iter])
@@ -78,7 +82,7 @@ end
 
 function soc_gap(v::ViewVector, s::ViewScalar)
 
-    return norm(v, 2) - s[]
+    return LinearAlgebra.norm(v, 2) - s[]
 end
 
 function convergedrank(a, p::Params, cones::ConicSets, opt::Options)::Bool
diff --git a/src/scaling.jl b/src/scaling.jl
index 6a7996b..a18a6bc 100644
--- a/src/scaling.jl
+++ b/src/scaling.jl
@@ -10,7 +10,7 @@ function preprocess!(aff::AffineSets, conic_sets::ConicSets)
             for i in eachindex(iv)
                 M[im[i]] = iv[i]
             end
-            X = Symmetric(M, :L)
+            X = LinearAlgebra.Symmetric(M, :L)
 
             n = size(X)[1] # columns or line
             cont = 1
@@ -41,21 +41,21 @@ end
 
 function norm_scaling(affine_sets::AffineSets, cones::ConicSets)
     cte = (sqrt(2.) / 2.)
-    rows = rowvals(affine_sets.A)
+    rows = SparseArrays.rowvals(affine_sets.A)
     cont = 1
     for sdp in cones.sdpcone, j in 1:sdp.sq_side, i in 1:j
         if i != j
-            for line in nzrange(affine_sets.A, cont)
+            for line in SparseArrays.nzrange(affine_sets.A, cont)
                 affine_sets.A[rows[line], cont] *= cte
             end
         end
         cont += 1
     end
-    rows = rowvals(affine_sets.G)
+    rows = SparseArrays.rowvals(affine_sets.G)
     cont = 1
     for sdp in cones.sdpcone, j in 1:sdp.sq_side, i in 1:j
         if i != j
-            for line in nzrange(affine_sets.G, cont)
+            for line in SparseArrays.nzrange(affine_sets.G, cont)
                 affine_sets.G[rows[line], cont] *= cte
             end
         end
diff --git a/src/structs.jl b/src/structs.jl
index d75b0a1..b52ccf7 100644
--- a/src/structs.jl
+++ b/src/structs.jl
@@ -29,146 +29,13 @@ function Base.length(V::CircularVector{T}) where T
     return V.l
 end
 
-Base.@kwdef mutable struct Options
-
-    # Printing options
-    log_verbose::Bool = false
-    log_freq::Int = 1000
-    timer_verbose::Bool = false
-    timer_file::Bool = false
-    disable_julia_logger = true
-
-    # time options
-    time_limit::Float64 = 3600_00. #100 hours
-
-    warn_on_limit::Bool = false
-    extended_log::Bool = false
-    extended_log2::Bool = false
-    log_repeat_header::Bool = false
-
-    # Default tolerances
-    tol_gap::Float64 = 1e-4
-    tol_feasibility::Float64 = 1e-4
-    tol_feasibility_dual::Float64 = 1e-4
-    tol_primal::Float64 = 1e-4
-    tol_dual::Float64 = 1e-4
-    tol_psd::Float64 = 1e-7
-    tol_soc::Float64 = 1e-7
-
-    check_dual_feas::Bool = false
-    check_dual_feas_freq::Int = 1000
-
-    max_obj::Float64 = 1e20
-    min_iter_max_obj::Int = 10
-
-    # infeasibility check
-    min_iter_time_infeas::Int = 1000
-    infeas_gap_tol::Float64 = 1e-4
-    infeas_limit_gap_tol::Float64 = 1e-1
-    infeas_stable_gap_tol::Float64 = 1e-4
-    infeas_feasibility_tol::Float64 = 1e-4
-    infeas_stable_feasibility_tol::Float64 = 1e-8
-
-    certificate_search::Bool = true
-    certificate_obj_tol::Float64 = 1e-1
-    certificate_fail_tol::Float64 = 1e-8
-
-    # Bounds on beta (dual_step / primal_step) [larger bounds may lead to numerical inaccuracy]
-    min_beta::Float64 = 1e-5
-    max_beta::Float64 = 1e+5
-    initial_beta::Float64 = 1.
-
-    # Adaptive primal-dual steps parameters [adapt_decay above .7 may lead to slower convergence]
-    initial_adapt_level::Float64 = .9
-    adapt_decay::Float64 = .8
-    adapt_window::Int64 = 50
-
-    # PDHG parameters
-    convergence_window::Int = 200
-    convergence_check::Int = 50
-    max_iter::Int = 0
-    min_iter::Int = 40
-    divergence_min_update::Int = 50
-    max_iter_lp::Int = 10_000_000
-    max_iter_conic::Int = 1_000_000
-    max_iter_local::Int = 0 #ignores user setting
-
-    advanced_initialization::Bool = true
-
-    # Linesearch parameters
-    line_search_flag::Bool = true
-    max_linsearch_steps::Int = 5000
-    delta::Float64 = .9999
-    initial_theta::Float64 = 1.
-    linsearch_decay::Float64 = .75
-
-    # Spectral decomposition parameters
-    full_eig_decomp::Bool = false
-    max_target_rank_krylov_eigs::Int = 16
-    min_size_krylov_eigs::Int = 100
-    warm_start_eig::Bool = true
-    rank_increment::Int = 1 # 0=multiply, 1 = add
-    rank_increment_factor::Int = 1 # 0 multiply, 1 = add
-
-    # eigsolver selection
-    #=
-        1: Arpack [dsaupd] (tipically non-deterministic)
-        2: KrylovKit [eigsolve/Lanczos] (DEFAULT)
-    =#
-    eigsolver::Int = 2
-    eigsolver_min_lanczos::Int = 25
-    eigsolver_resid_seed::Int = 1234
-
-    # Arpack
-    # note that Arpack is Non-deterministic
-    # (https://github.com/mariohsouto/ProxSDP.jl/issues/69)
-    arpack_tol::Float64 = 1e-10
-    #=
-        0: arpack random [usually faster - NON-DETERMINISTIC - slightly]
-        1: all ones [???]
-        2: julia random uniform (eigsolver_resid_seed) [medium for DETERMINISTIC]
-        3: julia normalized random normal (eigsolver_resid_seed) [best for DETERMINISTIC]
-    =#
-    arpack_resid_init::Int = 3
-    arpack_reset_resid::Bool = true # true for determinism
-    # larger is more stable to converge and more deterministic
-    arpack_max_iter::Int = 10_000
-    # see remark for of dsaupd
-    
-    # KrylovKit
-    krylovkit_reset_resid::Bool = false
-    krylovkit_resid_init::Int = 3
-    krylovkit_tol::Float64 = 1e-12
-    krylovkit_max_iter::Int = 100
-    krylovkit_eager::Bool = false
-    krylovkit_verbose::Int = 0
-
-    # Reduce rank [warning: heuristics]
-    reduce_rank::Bool = false
-    rank_slack::Int = 3
-
-    full_eig_freq::Int = 10000000
-    full_eig_len::Int = 0
-
-    # equilibration parameters
-    equilibration::Bool = false
-    equilibration_iters::Int = 1000
-    equilibration_lb::Float64 = -10.0
-    equilibration_ub::Float64 = +10.0
-    equilibration_limit::Float64 = 0.9
-    equilibration_force::Bool = false
-
-    # spectral norm [using exact norm via svds may result in nondeterministic behavior]
-    approx_norm::Bool = true
-end
-
 mutable struct AffineSets
     n::Int  # Size of primal variables
     p::Int  # Number of linear equalities
     m::Int  # Number of linear inequalities
     extra::Int  # Number of adition linear equalities (for disjoint cones)
-    A::SparseMatrixCSC{Float64,Int64} 
-    G::SparseMatrixCSC{Float64,Int64}
+    A::SparseArrays.SparseMatrixCSC{Float64,Int64}
+    G::SparseArrays.SparseMatrixCSC{Float64,Int64}
     b::Vector{Float64}
     h::Vector{Float64}
     c::Vector{Float64}
@@ -178,7 +45,6 @@ mutable struct SDPSet
     vec_i::Vector{Int}
     mat_i::Vector{Int}
     tri_len::Int
-    sq_len::Int
     sq_side::Int
 end
 
@@ -192,25 +58,27 @@ mutable struct ConicSets
     socone::Vector{SOCSet}
 end
 
-mutable struct CPResult
-    status::Int
-    status_string::String
-    primal::Vector{Float64}
-    dual_cone::Vector{Float64}
-    dual_eq::Vector{Float64}
-    dual_in::Vector{Float64}
-    slack_eq::Vector{Float64}
-    slack_in::Vector{Float64}
-    primal_residual::Float64
-    dual_residual::Float64
-    objval::Float64
-    dual_objval::Float64
-    gap::Float64
-    time::Float64
-    final_rank::Int
-    primal_feasible_user_tol::Bool
-    dual_feasible_user_tol::Bool
-    certificate_found::Bool
+Base.@kwdef mutable struct Result
+    status::Int = 0
+    status_string::String = "Problem not solved"
+    primal::Vector{Float64} = Float64[]
+    dual_cone::Vector{Float64} = Float64[]
+    dual_eq::Vector{Float64} = Float64[]
+    dual_in::Vector{Float64} = Float64[]
+    slack_eq::Vector{Float64} = Float64[]
+    slack_in::Vector{Float64} = Float64[]
+    primal_residual::Float64 = NaN
+    dual_residual::Float64 = NaN
+    objval::Float64 = NaN
+    dual_objval::Float64 = NaN
+    gap::Float64 = NaN
+    time::Float64 = NaN
+    iter::Int = -1
+    final_rank::Int = -1
+    primal_feasible_user_tol::Bool = false
+    dual_feasible_user_tol::Bool = false
+    certificate_found::Bool = false
+    result_count::Int = 0
 end
 
 mutable struct PrimalDual
@@ -258,7 +126,7 @@ const ViewVector = SubArray#{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int}},
 const ViewScalar = SubArray#{Float64, 1, Vector{Float64}, Tuple{Int}, true}
 
 mutable struct AuxiliaryData
-    m::Vector{Symmetric{Float64,Matrix{Float64}}}
+    m::Vector{LinearAlgebra.Symmetric{Float64,Matrix{Float64}}}
 
     Mty::Vector{Float64}
     Mty_old::Vector{Float64}
@@ -274,7 +142,7 @@ mutable struct AuxiliaryData
 
     function AuxiliaryData(aff::AffineSets, cones::ConicSets) 
         new(
-            [Symmetric(zeros(sdp.sq_side, sdp.sq_side), :U) for sdp in cones.sdpcone], 
+            [LinearAlgebra.Symmetric(zeros(sdp.sq_side, sdp.sq_side), :U) for sdp in cones.sdpcone], 
             zeros(aff.n), zeros(aff.n),
             zeros(aff.p+aff.m), zeros(aff.p+aff.m),
             zeros(aff.p+aff.m), zeros(aff.p+aff.m),
@@ -284,11 +152,9 @@ mutable struct AuxiliaryData
 end
 
 mutable struct Matrices
-    M::SparseMatrixCSC{Float64,Int64}
-    Mt::SparseMatrixCSC{Float64,Int64}
+    M::SparseArrays.SparseMatrixCSC{Float64,Int64}
+    Mt::SparseArrays.SparseMatrixCSC{Float64,Int64}
     c::Vector{Float64}
-
-    Matrices(M, Mt, c) = new(M, Mt, c)
 end
 
 mutable struct Params
diff --git a/src/util.jl b/src/util.jl
index ec4b734..e96621e 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -13,4 +13,26 @@ function map_socs!(v::Vector{Float64}, cones::ConicSets, a::AuxiliaryData)
         cont += len
     end
     return nothing
-end
\ No newline at end of file
+end
+
+function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
+    n = sympackeddim(length(v))
+    n1, n2 = size(out)
+    @assert n == n1 == n2
+    c = 0
+    for j in 1:n, i in 1:j
+        c += 1
+        out[i,j] = v[c]
+    end
+    return out
+end
+
+function ivech(v::AbstractVector{T}) where T
+    n = sympackeddim(length(v))
+    out = zeros(n, n)
+    ivech!(out, v)
+    
+    return out
+end
+
+ivec(X) = Matrix(LinearAlgebra.Symmetric(ivech(X),:U))
diff --git a/test/base_mimo.jl b/test/base_mimo.jl
index 7e54afb..36064d0 100644
--- a/test/base_mimo.jl
+++ b/test/base_mimo.jl
@@ -1,4 +1,5 @@
-using Random
+import Random
+
 function mimo_data(seed, m, n)
     rng = Random.MersenneTwister(seed)
 
@@ -17,11 +18,11 @@ end
 
 function mimo_eval(s, H, y, L, XX)
     x_hat = sign.(XX[1:end-1, end])
-    rank = length([eig for eig in eigen(XX).values if eig > 1e-7])
+    rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig > 1e-7])
     @show decode_error = sum(abs.(x_hat - s))
     @show rank
-    @show norm(y - H * x_hat)
-    @show norm(y - H * s)
+    @show LinearAlgebra.norm(y - H * x_hat)
+    @show LinearAlgebra.norm(y - H * s)
     @show tr(L * XX)
     return nothing
 end
\ No newline at end of file
diff --git a/test/base_randsdp.jl b/test/base_randsdp.jl
index 7e53f47..4740224 100644
--- a/test/base_randsdp.jl
+++ b/test/base_randsdp.jl
@@ -1,5 +1,5 @@
-using Random
-using LinearAlgebra
+import Random
+import LinearAlgebra
 
 function randsdp_data(seed, m, n)
     rng = Random.MersenneTwister(seed)
@@ -22,9 +22,9 @@ function randsdp_data(seed, m, n)
     return A, b, C
 end
 function randsdp_eval(A,b,C,n,m,XX)
-    @show minus_rank = length([eig for eig in eigen(XX).values if eig < -1e-10])
+    @show minus_rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig < -1e-10])
 
-    @show rank = length([eig for eig in eigen(XX).values if eig > 1e-10])
+    @show rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig > 1e-10])
 
     @show tr(C * XX)
     for i in 1:m
diff --git a/test/base_sdplib.jl b/test/base_sdplib.jl
index bf08e31..33e10a9 100644
--- a/test/base_sdplib.jl
+++ b/test/base_sdplib.jl
@@ -1,6 +1,6 @@
 function sdplib_data(path)
     # Read data from file
-    data = readdlm(path, use_mmap=true)
+    data = DelimitedFiles.readdlm(path, use_mmap=true)
 
     # Parse SDPLIB data
     m = data[1, 1]
@@ -23,7 +23,7 @@ function sdplib_data(path)
     end
     n = abs(cum_blks[end])
     n = length(c)
-    F = Dict(i => spzeros(n, n) for i = 0:m)
+    F = Dict(i => SparseArrays.spzeros(n, n) for i = 0:m)
     for k=5:size(data)[1]
         idx = cum_blks[data[k, 2]]
         # if data[k, 2] == 1
@@ -44,7 +44,7 @@ function sdplib_data(path)
     return n, m, F, c
 end
 function sdplib_eval(F,c,n,m,XX)
-    rank = length([eig for eig in eigen(XX).values if eig > 1e-10])
+    rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig > 1e-10])
     @show rank
     @show tr(F[0] * XX)
     nothing
diff --git a/test/base_sensorloc.jl b/test/base_sensorloc.jl
index 07ad1c2..60f0c34 100644
--- a/test/base_sensorloc.jl
+++ b/test/base_sensorloc.jl
@@ -10,18 +10,18 @@ function sensorloc_data(seed, n)
     # Sensor true position (2 dimensional)
     x_true = Random.rand(rng, Float64, (2, n))
     # Distances from sensors to sensors
-    d = Dict((i, j) => norm(x_true[:, i] - x_true[:, j]) for i in 1:n for j in 1:i)
+    d = Dict((i, j) => LinearAlgebra.norm(x_true[:, i] - x_true[:, j]) for i in 1:n for j in 1:i)
 
     # Anchor positions
     a = Dict(i => Random.rand(rng, Float64, (2, 1)) for i in 1:m)
     # Distances from anchor to sensors
-    d_bar = Dict((k, j) => norm(x_true[:, j] - a[k]) for k in 1:m for j in 1:n)
+    d_bar = Dict((k, j) => LinearAlgebra.norm(x_true[:, j] - a[k]) for k in 1:m for j in 1:n)
 
     return m, x_true, a, d, d_bar
 end
 
 function sensorloc_eval(n, m, x_true, XX)
-    @show norm(x_true - XX[1:2, 3:n + 2])
-    @show rank = length([eig for eig in eigen(XX).values if eig > 1e-7])
+    @show LinearAlgebra.norm(x_true - XX[1:2, 3:n + 2])
+    @show rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig > 1e-7])
     return nothing
 end
\ No newline at end of file
diff --git a/test/jump_mimo.jl b/test/jump_mimo.jl
index 97be138..0115d02 100644
--- a/test/jump_mimo.jl
+++ b/test/jump_mimo.jl
@@ -28,7 +28,7 @@ function jump_mimo(solver, seed, n; verbose = false, test = false)
     verbose && mimo_eval(s,H,y,L,XX)
 
     objval = objective_value(model)
-    stime = MOI.get(model, MOI.SolveTime())
+    stime = MOI.get(model, MOI.SolveTimeSec())
 
     rank = -1
     try
diff --git a/test/jump_randsdp.jl b/test/jump_randsdp.jl
index 2646f75..a405db3 100644
--- a/test/jump_randsdp.jl
+++ b/test/jump_randsdp.jl
@@ -15,7 +15,7 @@ function jump_randsdp(solver, seed, n, m, verbose = false)
     verbose && randsdp_eval(A,b,C,n,m,XX)
     
     objval = objective_value(model)
-    stime = MOI.get(model, MOI.SolveTime())
+    stime = MOI.get(model, MOI.SolveTimeSec())
 
     # @show tp = typeof(model.moi_backend.optimizer.model.optimizer)
     # @show fieldnames(tp)
diff --git a/test/jump_sdplib.jl b/test/jump_sdplib.jl
index 4b48def..015ca7a 100644
--- a/test/jump_sdplib.jl
+++ b/test/jump_sdplib.jl
@@ -11,12 +11,12 @@ function jump_sdplib(solver, path; verbose = false, test = false)
 
     # Objective function
     @objective(model, Min, sum(F[0][idx...] * X[idx...] 
-        for idx in zip(findnz(F[0])[1:end-1]...)))
+        for idx in zip(SparseArrays.findnz(F[0])[1:end-1]...)))
 
     # Linear equality constraints
     for k = 1:m
         @constraint(model, sum(F[k][idx...] * X[idx...]
-            for idx in zip(findnz(F[k])[1:end-1]...)) == c[k])
+            for idx in zip(SparseArrays.findnz(F[k])[1:end-1]...)) == c[k])
     end
     
     teste = @time optimize!(model)
@@ -26,7 +26,7 @@ function jump_sdplib(solver, path; verbose = false, test = false)
     verbose && sdplib_eval(F,c,n,m,XX)
 
     objval = objective_value(model)
-    stime = MOI.get(model, MOI.SolveTime())
+    stime = MOI.get(model, MOI.SolveTimeSec())
 
 
     # @show tp = typeof(model.moi_backend.optimizer.model.optimizer)
@@ -48,7 +48,7 @@ function jump_sdplib(solver, path; verbose = false, test = false)
     max_lin_viol = 0.0
     for k = 1:m
         val = abs(sum(F[k][idx...] * XX[idx...]
-            for idx in zip(findnz(F[k])[1:end-1]...)) - c[k])
+            for idx in zip(SparseArrays.findnz(F[k])[1:end-1]...)) - c[k])
         if val > 0.0
             if val > max_lin_viol
                 max_lin_viol = val
diff --git a/test/jump_sensorloc.jl b/test/jump_sensorloc.jl
index 0235131..f4ac504 100644
--- a/test/jump_sensorloc.jl
+++ b/test/jump_sensorloc.jl
@@ -1,4 +1,4 @@
-using Random
+import Random
 function jump_sensorloc(solver, seed, n; verbose = false, test = false)
 
     m, x_true, a, d, d_bar = sensorloc_data(seed, n)
@@ -59,7 +59,7 @@ function jump_sensorloc(solver, seed, n; verbose = false, test = false)
     verbose && sensorloc_eval(n, m, x_true, XX)
 
     objval = objective_value(model)
-    stime = MOI.get(model, MOI.SolveTime())
+    stime = MOI.get(model, MOI.SolveTimeSec())
 
     # @show tp = typeof(model.moi_backend.optimizer.model.optimizer)
     # @show fieldnames(tp)
diff --git a/test/moi_init.jl b/test/moi_init.jl
deleted file mode 100644
index 86d0572..0000000
--- a/test/moi_init.jl
+++ /dev/null
@@ -1,107 +0,0 @@
-const MOI = MathOptInterface
-const MOIT = MOI.Test
-const MOIB = MOI.Bridges
-const MOIU = MOI.Utilities
-
-sympackedlen(n) = div(n*(n+1), 2)
-sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
-function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
-    n = sympackeddim(length(v))
-    n1, n2 = size(out)
-    @assert n == n1 == n2
-    c = 0
-    for j in 1:n, i in 1:j
-        c += 1
-        out[i,j] = v[c]
-    end
-    return out
-end
-function ivech(v::AbstractVector{T}) where T
-    n = sympackeddim(length(v))
-    out = zeros(n, n)
-    ivech!(out, v)
-    
-    return out
-end
-
-solvers = Tuple{String, Any}[]
-
-#MOIU.@model ProxSDPModelData () (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan) (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.PositiveSemidefiniteConeTriangle) () (MOI.SingleVariable,) (MOI.ScalarAffineFunction,) (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)
-
-#=
-    CSDP & SDPA
-=#
-
-# ]activate soi
-# using CSDP
-# using SDPA
-
-# MOIU.@model(SOI_ModelData,
-#             (),
-#             (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
-#             (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
-#              MOI.PositiveSemidefiniteConeTriangle),
-#             (),
-#             (MOI.SingleVariable,),
-#             (MOI.ScalarAffineFunction,),
-#             (MOI.VectorOfVariables,),
-#             (MOI.VectorAffineFunction,))
-
-# const CSDP_optimizer = 
-# MOIB.full_bridge_optimizer(
-#     MOIU.CachingOptimizer(
-#         SOI_ModelData{Float64}(), CSDP.Optimizer()),#printlevel=0))
-#     Float64
-#     )
-
-# push!(solvers, ("CSDP", CSDP_optimizer)) # eps = ???
-
-# const SDPA_optimizer = 
-# MOIB.full_bridge_optimizer(
-#     MOIU.CachingOptimizer(
-#         SOI_ModelData{Float64}(), SDPA.Optimizer()),#printlevel=0))
-#     Float64
-#     )
-
-# push!(solvers, ("SDPA", SDPA_optimizer)) # eps = ???
-
-#=
-    COSMO
-=#
-
-# using COSMO
-
-# MOIU.@model(COSMOModelData,
-#         (),
-#         (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
-#         (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
-#          MOI.PositiveSemidefiniteConeSquare, MOI.PositiveSemidefiniteConeTriangle, MOI.ExponentialCone, MOI.DualExponentialCone),
-#         (MOI.PowerCone, MOI.DualPowerCone),
-#         (MOI.SingleVariable,),
-#         (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
-#         (MOI.VectorOfVariables,),
-#         (MOI.VectorAffineFunction,),);
-
-# const COSMO_optimizer = MOIU.CachingOptimizer(MOIU.UniversalFallback(COSMOModelData{Float64}()),
-#     COSMO.Optimizer(eps_abs = 1e-4, eps_rel = 1e-4 ));
-
-# push!(solvers, ("COSMO", COSMO_optimizer))
-
-#=
-    SeDuMi
-=#
-
-import SeDuMi
-const optimizer = SeDuMi.Optimizer(fid=0)
-
-MOIU.@model(SeDuMi_ModelData, (), (),
-            (MOI.Zeros, MOI.Nonnegatives, MOI.SecondOrderCone,
-             MOI.RotatedSecondOrderCone, MOI.PositiveSemidefiniteConeTriangle),
-            (), (), (), (MOI.VectorOfVariables,), (MOI.VectorAffineFunction,))
-
-const SeDuMi_optimizer = MOIB.full_bridge_optimizer(MOIU.CachingOptimizer(
-    MOIU.UniversalFallback(SeDuMi_ModelData{Float64}()),
-    SeDuMi.Optimizer(fid=0)),
-    Float64)
-    
-push!(solvers, ("SeDuMi", SeDuMi_optimizer))
diff --git a/test/moi_mimo.jl b/test/moi_mimo.jl
index a91408c..a05d525 100644
--- a/test/moi_mimo.jl
+++ b/test/moi_mimo.jl
@@ -32,8 +32,8 @@ function moi_mimo(optimizer, seed, n; verbose = false, test = false, scalar = fa
     end
 
     Xsq = Matrix{MOI.VariableIndex}(undef, n+1,n+1)
-    ivech!(Xsq, X)
-    Xsq = Matrix(Symmetric(Xsq,:U))
+    ProxSDP.ivech!(Xsq, X)
+    Xsq = Matrix(LinearAlgebra.Symmetric(Xsq,:U))
 
     vov = MOI.VectorOfVariables(X)
     cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(n+1))
@@ -61,7 +61,7 @@ function moi_mimo(optimizer, seed, n; verbose = false, test = false, scalar = fa
 
     stime = -1.0
     try
-        stime = MOI.get(optimizer, MOI.SolveTime())
+        stime = MOI.get(optimizer, MOI.SolveTimeSec())
     catch
         println("could not query time")
     end
diff --git a/test/moi_proxsdp_unit.jl b/test/moi_proxsdp_unit.jl
new file mode 100644
index 0000000..876a09f
--- /dev/null
+++ b/test/moi_proxsdp_unit.jl
@@ -0,0 +1,362 @@
+function simple_lp(optimizer)
+
+    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
+    MOI.empty!(bridged)
+    @test MOI.is_empty(bridged)
+
+    # add 10 variables - only diagonal is relevant
+    X = MOI.add_variables(bridged, 2)
+
+    # add sdp constraints - only ensuring positivenesse of the diagonal
+    vov = MOI.VectorOfVariables(X)
+
+    c1 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(2.0, X[1]),
+            MOI.ScalarAffineTerm(1.0, X[2])
+        ], 0.0), MOI.EqualTo(4.0))
+
+    c2 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(1.0, X[1]),
+            MOI.ScalarAffineTerm(2.0, X[2])
+        ], 0.0), MOI.EqualTo(4.0))
+
+    b1 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(1.0, X[1])
+        ], 0.0), MOI.GreaterThan(0.0))
+
+    b2 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(1.0, X[2])
+        ], 0.0), MOI.GreaterThan(0.0))
+
+    MOI.set(bridged, 
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
+        )
+    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.optimize!(bridged)
+
+    obj = MOI.get(bridged, MOI.ObjectiveValue())
+
+    @test obj ≈ -9.33333 atol = 1e-2
+
+    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)
+
+    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2
+
+end
+
+function simple_lp_2_1d_sdp(optimizer)
+
+    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
+    MOI.empty!(bridged)
+    @test MOI.is_empty(bridged)
+
+    # add 10 variables - only diagonal is relevant
+    X = MOI.add_variables(bridged, 2)
+
+    # add sdp constraints - only ensuring positivenesse of the diagonal
+    vov = MOI.VectorOfVariables(X)
+
+    c1 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(2.0, X[1]),
+            MOI.ScalarAffineTerm(1.0, X[2])
+        ], 0.0), MOI.EqualTo(4.0))
+
+    c2 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(1.0, X[1]),
+            MOI.ScalarAffineTerm(2.0, X[2])
+        ], 0.0), MOI.EqualTo(4.0))
+
+    b1 = MOI.add_constraint(bridged, 
+        MOI.VectorOfVariables([X[1]]), MOI.PositiveSemidefiniteConeTriangle(1))
+
+    b2 = MOI.add_constraint(bridged, 
+        MOI.VectorOfVariables([X[2]]), MOI.PositiveSemidefiniteConeTriangle(1))
+
+    MOI.set(bridged, 
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
+        )
+    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.optimize!(bridged)
+
+    obj = MOI.get(bridged, MOI.ObjectiveValue())
+
+    @test obj ≈ -9.33333 atol = 1e-2
+
+    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)
+
+    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2
+
+end
+
+function lp_in_SDP_equality_form(optimizer)
+
+    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
+    MOI.empty!(bridged)
+    @test MOI.is_empty(bridged)
+
+    # add 10 variables - only diagonal is relevant
+    X = MOI.add_variables(bridged, 10)
+
+    # add sdp constraints - only ensuring positivenesse of the diagonal
+    vov = MOI.VectorOfVariables(X)
+    cX = MOI.add_constraint(bridged, vov, MOI.PositiveSemidefiniteConeTriangle(4))
+
+    c1 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(2.0, X[1]),
+            MOI.ScalarAffineTerm(1.0, X[3]),
+            MOI.ScalarAffineTerm(1.0, X[6])
+        ], 0.0), MOI.EqualTo(4.0))
+
+    c2 = MOI.add_constraint(bridged, 
+        MOI.ScalarAffineFunction([
+            MOI.ScalarAffineTerm(1.0, X[1]),
+            MOI.ScalarAffineTerm(2.0, X[3]),
+            MOI.ScalarAffineTerm(1.0, X[10])
+        ], 0.0), MOI.EqualTo(4.0))
+
+    MOI.set(bridged, 
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[3]]), 0.0)
+        )
+    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.optimize!(bridged)
+
+    obj = MOI.get(bridged, MOI.ObjectiveValue())
+
+    @test obj ≈ -9.33333 atol = 1e-2
+
+    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)
+
+    @test Xr ≈ [1.3333, .0, 1.3333, .0, .0, .0, .0, .0, .0, .0] atol = 1e-2
+
+end
+
+function lp_in_SDP_inequality_form(optimizer)
+
+    MOI.empty!(optimizer)
+    @test MOI.is_empty(optimizer)
+
+    # add 10 variables - only diagonal is relevant
+    X = MOI.add_variables(optimizer, 3)
+
+    # add sdp constraints - only ensuring positivenesse of the diagonal
+    vov = MOI.VectorOfVariables(X)
+    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
+
+    c1 = MOI.add_constraint(optimizer, 
+        MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,1],[
+            MOI.ScalarAffineTerm(2.0, X[1]),
+            MOI.ScalarAffineTerm(1.0, X[3]),
+        ]), [-4.0]), MOI.Nonpositives(1))
+
+    c2 = MOI.add_constraint(optimizer, 
+        MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,1],[
+            MOI.ScalarAffineTerm(1.0, X[1]),
+            MOI.ScalarAffineTerm(2.0, X[3]),
+        ]), [-4.0]), MOI.Nonpositives(1))
+
+    MOI.set(optimizer, 
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 3.0], [X[1], X[3]]), 0.0)
+        )
+    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    MOI.optimize!(optimizer)
+
+    obj = MOI.get(optimizer, MOI.ObjectiveValue())
+
+    @test obj ≈ 9.33333 atol = 1e-2
+
+    Xr = MOI.get(optimizer, MOI.VariablePrimal(), X)
+
+    @test Xr ≈ [1.3333, .0, 1.3333] atol = 1e-2
+
+    c1_d = MOI.get(optimizer, MOI.ConstraintDual(), c1)
+    c2_d = MOI.get(optimizer, MOI.ConstraintDual(), c2)
+
+end
+
+function sdp_from_moi(optimizer)
+    # min X[1,1] + X[2,2]    max y
+    #     X[2,1] = 1         [0   y/2     [ 1  0
+    #                         y/2 0    <=   0  1]
+    #     X >= 0              y free
+    # Optimal solution:
+    #
+    #     ⎛ 1   1 ⎞
+    # X = ⎜       ⎟           y = 2
+    #     ⎝ 1   1 ⎠
+    MOI.empty!(optimizer)
+    @test MOI.is_empty(optimizer)
+
+    X = MOI.add_variables(optimizer, 3)
+
+    vov = MOI.VectorOfVariables(X)
+    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
+
+    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))
+
+    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[3]]), 0.0))
+
+    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.optimize!(optimizer)
+
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
+    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT
+
+    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2 atol=1e-2
+
+    Xv = ones(3)
+    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
+    # @test MOI.get(optimizer, MOI.ConstraintPrimal(), cX) ≈ Xv atol=1e-2
+
+    # @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ 2 atol=1e-2
+    # @show MOI.get(optimizer, MOI.ConstraintDual(), c)
+
+end
+
+function double_sdp_from_moi(optimizer)
+    # solve simultaneously two of these:
+    # min X[1,1] + X[2,2]    max y
+    #     X[2,1] = 1         [0   y/2     [ 1  0
+    #                         y/2 0    <=   0  1]
+    #     X >= 0              y free
+    # Optimal solution:
+    #
+    #     ⎛ 1   1 ⎞
+    # X = ⎜       ⎟           y = 2
+    #     ⎝ 1   1 ⎠
+    MOI.empty!(optimizer)
+    @test MOI.is_empty(optimizer)
+
+    X = MOI.add_variables(optimizer, 3)
+    Y = MOI.add_variables(optimizer, 3)
+
+    vov = MOI.VectorOfVariables(X)
+    vov2 = MOI.VectorOfVariables(Y)
+    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
+    cY = MOI.add_constraint(optimizer, vov2, MOI.PositiveSemidefiniteConeTriangle(2))
+
+    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))
+    c2 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, Y[2]))], [-1.0]), MOI.Zeros(1))
+
+    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[end], Y[1], Y[end]]), 0.0))
+
+    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.optimize!(optimizer)
+
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
+    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT
+
+    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2*2 atol=1e-2
+
+    Xv = ones(3)
+    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
+    Yv = ones(3)
+    @test MOI.get(optimizer, MOI.VariablePrimal(), Y) ≈ Yv atol=1e-2
+    # @test MOI.get(optimizer, MOI.ConstraintPrimal(), cX) ≈ Xv atol=1e-2
+
+    # @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ 2 atol=1e-2
+    # @show MOI.get(optimizer, MOI.ConstraintDual(), c)
+
+end
+
+function double_sdp_with_duplicates(optimizer)
+
+    cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
+    optimizer0 = ProxSDP.Optimizer(
+        log_verbose = true,
+        log_freq = 10,
+        check_dual_feas = true,
+        check_dual_feas_freq = 10,
+    )
+    MOI.empty!(cache);
+    optimizer1 = MOIU.CachingOptimizer(cache, optimizer0);
+    optimizer = MOIB.full_bridge_optimizer(optimizer1, Float64);
+
+    MOI.empty!(optimizer)
+
+    x = MOI.add_variable(optimizer)
+    X = [x, x, x]
+
+    vov = MOI.VectorOfVariables(X)
+    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
+
+    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))
+
+    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[3]]), 0.0))
+
+    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.optimize!(optimizer)
+
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+
+    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
+    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT
+
+    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2 atol=1e-2
+
+    Xv = ones(3)
+    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
+
+end
+
+function sdp_wiki(optimizer)
+    # https://en.wikipedia.org/wiki/Semidefinite_programming
+    MOI.empty!(optimizer)
+    @test MOI.is_empty(optimizer)
+
+    X = MOI.add_variables(optimizer, 6)
+
+    vov = MOI.VectorOfVariables(X)
+    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(3))
+
+    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[1]))], [-1.0]), MOI.Zeros(1))
+    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[3]))], [-1.0]), MOI.Zeros(1))
+    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[6]))], [-1.0]), MOI.Zeros(1))
+
+    c12_ub = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [0.1]), MOI.Nonpositives(1))   # x <= -0.1 -> x + 0.1 <= 0
+    c12_lb = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(-1.0, X[2]))], [-0.2]), MOI.Nonpositives(1)) # x >= -0.2 -> -x + -0.2 <= 0
+
+    c23_ub = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[5]))], [-0.5]), MOI.Nonpositives(1)) # x <= 0.5 ->  x - 0.5 <= 0
+    c23_lb = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(-1.0, X[5]))], [0.4]), MOI.Nonpositives(1)) # x >= 0.4 -> -x + 0.4 <= 0
+
+    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[4]]), 0.0))
+
+    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+
+    MOI.optimize!(optimizer)
+
+    obj = MOI.get(optimizer, MOI.ObjectiveValue())
+    @test obj ≈ -0.978 atol=1e-2
+
+    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+
+    MOI.optimize!(optimizer)
+
+    obj = MOI.get(optimizer, MOI.ObjectiveValue())
+    @test obj ≈ 0.872 atol=1e-2
+
+end
+
+
+simple_lp(optimizer)
+simple_lp_2_1d_sdp(optimizer)
+lp_in_SDP_equality_form(optimizer)
+lp_in_SDP_inequality_form(optimizer)
+sdp_from_moi(optimizer)
+double_sdp_from_moi(optimizer)
+double_sdp_with_duplicates(optimizer)
+sdp_wiki(optimizer)
\ No newline at end of file
diff --git a/test/moi_randsdp.jl b/test/moi_randsdp.jl
index ba8ce9d..6329358 100644
--- a/test/moi_randsdp.jl
+++ b/test/moi_randsdp.jl
@@ -12,8 +12,8 @@ function moi_randsdp(optimizer, seed, n, m; verbose = false, test = false, atol
     X = MOI.add_variables(optimizer, nvars)
 
     Xsq = Matrix{MOI.VariableIndex}(undef,n,n)
-    ivech!(Xsq, X)
-    Xsq = Matrix(Symmetric(Xsq,:U))
+    ProxSDP.ivech!(Xsq, X)
+    Xsq = Matrix(LinearAlgebra.Symmetric(Xsq,:U))
 
     for k in 1:m
         ctr_k = vec([MOI.ScalarAffineTerm(A[k][i,j], Xsq[i,j]) for j in 1:n, i in 1:n])
@@ -58,7 +58,7 @@ function moi_randsdp(optimizer, seed, n, m; verbose = false, test = false, atol
     objval = Inf
     stime = -1.0
     try
-        stime = MOI.get(optimizer, MOI.SolveTime())
+        stime = MOI.get(optimizer, MOI.SolveTimeSec())
     catch
         println("could not query time")
     end
@@ -67,11 +67,11 @@ function moi_randsdp(optimizer, seed, n, m; verbose = false, test = false, atol
 
         Xsq_s = MOI.get.(optimizer, MOI.VariablePrimal(), Xsq)
 
-        minus_rank = length([eig for eig in eigen(Xsq_s).values if eig < -1e-4])
+        minus_rank = length([eig for eig in LinearAlgebra.eigen(Xsq_s).values if eig < -1e-4])
         if test
             @test minus_rank == 0
         end
-        # rank = length([eig for eig in eigen(XX).values if eig > 1e-10])
+        # rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig > 1e-10])
         # @show rank
         if test
             @test tr(C * Xsq_s) - objval < atol
diff --git a/test/moi_sdplib.jl b/test/moi_sdplib.jl
index 59d8b95..fa6e571 100644
--- a/test/moi_sdplib.jl
+++ b/test/moi_sdplib.jl
@@ -17,19 +17,19 @@ function moi_sdplib(optimizer, path; verbose = false, test = false, scalar = fal
     cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(n))
 
     Xsq = Matrix{MOI.VariableIndex}(undef, n,n)
-    ivech!(Xsq, X)
-    Xsq = Matrix(Symmetric(Xsq,:U))
+    ProxSDP.ivech!(Xsq, X)
+    Xsq = Matrix(LinearAlgebra.Symmetric(Xsq,:U))
 
     # Objective function
     objf_t = [MOI.ScalarAffineTerm(F[0][idx...], Xsq[idx...])
-        for idx in zip(findnz(F[0])[1:end-1]...)]
+        for idx in zip(SparseArrays.findnz(F[0])[1:end-1]...)]
     MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(objf_t, 0.0))
     MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
 
     # Linear equality constraints
     for k in 1:m
         ctr_k = [MOI.ScalarAffineTerm(F[k][idx...], Xsq[idx...]) 
-            for idx in zip(findnz(F[k])[1:end-1]...)]
+            for idx in zip(SparseArrays.findnz(F[k])[1:end-1]...)]
         if scalar
             MOI.add_constraint(optimizer, MOI.ScalarAffineFunction(ctr_k, 0.0), MOI.EqualTo(c[k]))
         else
@@ -44,13 +44,13 @@ function moi_sdplib(optimizer, path; verbose = false, test = false, scalar = fal
 
     stime = -1.0
     try
-        stime = MOI.get(optimizer, MOI.SolveTime())
+        stime = MOI.get(optimizer, MOI.SolveTimeSec())
     catch
         println("could not query time")
     end
 
     Xsq_s = MOI.get.(optimizer, MOI.VariablePrimal(), Xsq)
-    minus_rank = length([eig for eig in eigen(Xsq_s).values if eig < -1e-4])
+    minus_rank = length([eig for eig in LinearAlgebra.eigen(Xsq_s).values if eig < -1e-4])
     if test
         @test minus_rank == 0
     end
diff --git a/test/moi_sensorloc.jl b/test/moi_sensorloc.jl
index 522c373..6324bbf 100644
--- a/test/moi_sensorloc.jl
+++ b/test/moi_sensorloc.jl
@@ -13,8 +13,8 @@ function moi_sensorloc(optimizer, seed, n; verbose = false, test = false, scalar
     nvars = ProxSDP.sympackedlen(n + 2) 
     X = MOI.add_variables(optimizer, nvars)
     Xsq = Matrix{MOI.VariableIndex}(undef, n + 2, n + 2)
-    ivech!(Xsq, X)
-    Xsq = Matrix(Symmetric(Xsq, :U))
+    ProxSDP.ivech!(Xsq, X)
+    Xsq = Matrix(LinearAlgebra.Symmetric(Xsq, :U))
     vov = MOI.VectorOfVariables(X)
     cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(n + 2))
 
@@ -103,7 +103,7 @@ function moi_sensorloc(optimizer, seed, n; verbose = false, test = false, scalar
 
     stime = -1.0
     try
-        stime = MOI.get(optimizer, MOI.SolveTime())
+        stime = MOI.get(optimizer, MOI.SolveTimeSec())
     catch
         println("could not query time")
     end
diff --git a/test/moitest.jl b/test/moitest.jl
index 61e02a5..a5f6aba 100644
--- a/test/moitest.jl
+++ b/test/moitest.jl
@@ -1,14 +1,15 @@
-push!(Base.LOAD_PATH,joinpath(dirname(@__FILE__),"..",".."))
+push!(Base.LOAD_PATH, joinpath(dirname(@__FILE__), "..", ".."))
 
-using ProxSDP, MathOptInterface, Test, LinearAlgebra, Random, SparseArrays, DelimitedFiles
-
-using LinearAlgebra
-LinearAlgebra.symmetric_type(::Type{MathOptInterface.VariableIndex}) = MathOptInterface.VariableIndex
-LinearAlgebra.symmetric(v::MathOptInterface.VariableIndex, ::Symbol) = v
-LinearAlgebra.transpose(v::MathOptInterface.VariableIndex) = v
+using Test
+import ProxSDP
+import MathOptInterface
+import LinearAlgebra
+import Random
+import SparseArrays
+import DelimitedFiles
 
 const MOI = MathOptInterface
-const MOIT = MOI.Test
+const MOIT = MOI.DeprecatedTest
 const MOIB = MOI.Bridges
 const MOIU = MOI.Utilities
 
@@ -50,20 +51,94 @@ const optimizer_lowacc_arpack = MOIU.CachingOptimizer(cache,
     ProxSDP.Optimizer(eigsolver = 1, tol_gap = 1e-3, tol_feasibility = 1e-3, log_verbose = false))
 const optimizer_lowacc_krylovkit = MOIU.CachingOptimizer(cache,
     ProxSDP.Optimizer(eigsolver = 2, tol_gap = 1e-3, tol_feasibility = 1e-3, log_verbose = false))
-const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, infeas_certificates = true)
-const config_conic = MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals = true, infeas_certificates = true)
-# const config_conic_nodual = MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals = false, infeas_certificates = true)
+const config = MOIT.Config{Float64}(atol=1e-3, rtol=1e-3, infeas_certificates = true)
+const config_conic = MOIT.Config{Float64}(atol=1e-3, rtol=1e-3, duals = true, infeas_certificates = true)
 
 @testset "SolverName" begin
     @test MOI.get(optimizer, MOI.SolverName()) == "ProxSDP"
 end
 
-@testset "supports_default_copy_to" begin
-    @test MOIU.supports_allocate_load(ProxSDP.Optimizer(), false)
-    @test !MOIU.supports_allocate_load(ProxSDP.Optimizer(), true)
+@testset "SolverVersion" begin
+    ver = readlines(joinpath(@__DIR__, "..", "Project.toml"))[4][12:16]
+    @test MOI.get(optimizer, MOI.SolverVersion()) == ver
+end
+
+function test_runtests()
+
+    config = MOI.Test.Config(
+        atol = 1e-4,
+        rtol = 1e-3,
+        exclude = Any[
+            MOI.ConstraintBasisStatus,
+            MOI.VariableBasisStatus,
+            MOI.ConstraintName,
+            MOI.VariableName,
+            MOI.ObjectiveBound,
+            MOI.ScalarFunctionConstantNotZero, # ignored by UniversalFallback
+        ],
+    )
+
+    opt = ProxSDP.Optimizer(
+        tol_gap = 1e-6, tol_feasibility= 1e-6,
+        # max_iter = 100_000,
+        time_limit = 1.0, #seconds FAST
+        warn_on_limit = true,
+        # log_verbose = true, log_freq = 100000
+        )
+    MOI.set(opt, MOI.Silent(), true)
+    model = MOI.Bridges.full_bridge_optimizer(
+        MOI.Utilities.CachingOptimizer(
+            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+            opt,
+        ),
+        Float64,
+    )
+    MOI.Test.runtests(
+        model,
+        config,
+        exclude = String[
+            # unexpected failure. But this is probably in the bridge
+            # layer, not ProxSDP.
+            # see: https://github.com/jump-dev/MathOptInterface.jl/issues/1665
+            "test_model_UpperBoundAlreadySet",
+            "test_model_LowerBoundAlreadySet",
+            # # TODO(joaquimg): good catch, but very pathological
+            # "test_objective_ObjectiveFunction_blank",
+            # poorly scaled problem (solved bellow with higher accuracy)
+            "test_linear_add_constraints",
+        ],
+    )
+
+    opt = ProxSDP.Optimizer(
+        tol_primal = 1e-7, tol_dual = 1e-7,
+        tol_gap = 1e-7, tol_feasibility = 1e-7,
+        time_limit = 5.0,
+        warn_on_limit = true,
+        )
+    MOI.set(opt, MOI.Silent(), true)
+    model = MOI.Bridges.full_bridge_optimizer(
+        MOI.Utilities.CachingOptimizer(
+            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+            opt,
+        ),
+        Float64,
+    )
+    MOI.empty!(model)
+    MOI.Test.test_linear_add_constraints(
+        model,
+        config,
+    )
+
+    @test MOI.get(model, ProxSDP.PDHGIterations()) >= 0
+
+    return
+end
+
+@testset "MOI Unit" begin
+    test_runtests()
 end
 
-@testset "Unit" begin
+@testset "Old MOI Unit" begin
     bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
     MOIT.unittest(bridged, config,[
         # not supported attributes
@@ -87,13 +162,13 @@ end
         ]
     )
     # TODO:
-    bridged_slow = MOIB.full_bridge_optimizer(optimizer_slow, Float64)
+    # bridged_slow = MOIB.full_bridge_optimizer(optimizer_slow, Float64)
     # MOIT.solve_farkas_interval_upper(bridged_slow, config)
     # MOIT.solve_farkas_interval_lower(bridged, config)
     # MOIT.solve_farkas_equalto_upper(bridged_slow, config)
     # MOIT.solve_farkas_equalto_lower(bridged, config)
     # MOIT.solve_farkas_variable_lessthan_max(bridged_slow, config)
-    MOIT.solve_farkas_variable_lessthan(bridged_slow, config)
+    # MOIT.solve_farkas_variable_lessthan(bridged_slow, config)
     # MOIT.solve_farkas_lessthan(bridged_slow, config)
     # MOIT.solve_farkas_greaterthan(bridged, config)
 end
@@ -109,17 +184,20 @@ end
         "linear8a",
         #"linear8b", "linear8c", "linear12",
         # poorly conditioned
-        "linear10",
         "linear5",
         "linear9",
+        "linear10",
         # primalstart not accepted
         "partial_start",
+        # certificate
+        "linear12",
         ]
     )
     MOIT.linear8atest(MOIB.full_bridge_optimizer(optimizer_high_acc, Float64), config)
-    MOIT.linear9test(MOIB.full_bridge_optimizer(optimizer_high_acc, Float64), config)
     MOIT.linear5test(MOIB.full_bridge_optimizer(optimizer_high_acc, Float64), config)
+    MOIT.linear9test(MOIB.full_bridge_optimizer(optimizer_high_acc, Float64), config)
     MOIT.linear10test(MOIB.full_bridge_optimizer(optimizer_high_acc, Float64), config)
+    MOIT.linear12test(MOIB.full_bridge_optimizer(optimizer_high_acc, Float64), config)
 end
 
 @testset "MOI Continuous Conic" begin
@@ -142,362 +220,18 @@ end
         # "normone2", "norminf2", "rotatedsoc2"#
         # slow to find certificate
         "normone2",
+        "norminf2", # new
         ]
     )
     # # these fail due to infeasibility certificate not being disabled
     # MOIT.norminf2test(MOIB.full_bridge_optimizer(optimizer, Float64), config_conic_nodual)
-    MOIT.normone2test(MOIB.full_bridge_optimizer(optimizer_slow, Float64), config_conic)
+    # MOIT.normone2test(MOIB.full_bridge_optimizer(optimizer_slow, Float64), config_conic)
     # # requires certificates always
     # MOIT.rotatedsoc2test(MOIB.full_bridge_optimizer(optimizer, Float64), config_conic_nodual)
 end
 
-@testset "Simple LP" begin
-
-    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
-    MOI.empty!(bridged)
-    @test MOI.is_empty(bridged)
-
-    # add 10 variables - only diagonal is relevant
-    X = MOI.add_variables(bridged, 2)
-
-    # add sdp constraints - only ensuring positivenesse of the diagonal
-    vov = MOI.VectorOfVariables(X)
-
-    c1 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(2.0, X[1]),
-            MOI.ScalarAffineTerm(1.0, X[2])
-        ], 0.0), MOI.EqualTo(4.0))
-
-    c2 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(1.0, X[1]),
-            MOI.ScalarAffineTerm(2.0, X[2])
-        ], 0.0), MOI.EqualTo(4.0))
-
-    b1 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(1.0, X[1])
-        ], 0.0), MOI.GreaterThan(0.0))
-
-    b2 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(1.0, X[2])
-        ], 0.0), MOI.GreaterThan(0.0))
-
-    MOI.set(bridged, 
-        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
-        )
-    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.optimize!(bridged)
-
-    obj = MOI.get(bridged, MOI.ObjectiveValue())
-
-    @test obj ≈ -9.33333 atol = 1e-2
-
-    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)
-
-    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2
-
-end
-
-@testset "Simple LP with 2 1D SDP" begin
-
-    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
-    MOI.empty!(bridged)
-    @test MOI.is_empty(bridged)
-
-    # add 10 variables - only diagonal is relevant
-    X = MOI.add_variables(bridged, 2)
-
-    # add sdp constraints - only ensuring positivenesse of the diagonal
-    vov = MOI.VectorOfVariables(X)
-
-    c1 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(2.0, X[1]),
-            MOI.ScalarAffineTerm(1.0, X[2])
-        ], 0.0), MOI.EqualTo(4.0))
-
-    c2 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(1.0, X[1]),
-            MOI.ScalarAffineTerm(2.0, X[2])
-        ], 0.0), MOI.EqualTo(4.0))
-
-    b1 = MOI.add_constraint(bridged, 
-        MOI.VectorOfVariables([X[1]]), MOI.PositiveSemidefiniteConeTriangle(1))
-
-    b2 = MOI.add_constraint(bridged, 
-        MOI.VectorOfVariables([X[2]]), MOI.PositiveSemidefiniteConeTriangle(1))
-
-    MOI.set(bridged, 
-        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
-        )
-    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.optimize!(bridged)
-
-    obj = MOI.get(bridged, MOI.ObjectiveValue())
-
-    @test obj ≈ -9.33333 atol = 1e-2
-
-    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)
-
-    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2
-
-end
-
-@testset "LP in SDP EQ form" begin
-
-    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
-    MOI.empty!(bridged)
-    @test MOI.is_empty(bridged)
-
-    # add 10 variables - only diagonal is relevant
-    X = MOI.add_variables(bridged, 10)
-
-    # add sdp constraints - only ensuring positivenesse of the diagonal
-    vov = MOI.VectorOfVariables(X)
-    cX = MOI.add_constraint(bridged, vov, MOI.PositiveSemidefiniteConeTriangle(4))
-
-    c1 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(2.0, X[1]),
-            MOI.ScalarAffineTerm(1.0, X[3]),
-            MOI.ScalarAffineTerm(1.0, X[6])
-        ], 0.0), MOI.EqualTo(4.0))
-
-    c2 = MOI.add_constraint(bridged, 
-        MOI.ScalarAffineFunction([
-            MOI.ScalarAffineTerm(1.0, X[1]),
-            MOI.ScalarAffineTerm(2.0, X[3]),
-            MOI.ScalarAffineTerm(1.0, X[10])
-        ], 0.0), MOI.EqualTo(4.0))
-
-    MOI.set(bridged, 
-        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[3]]), 0.0)
-        )
-    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.optimize!(bridged)
-
-    obj = MOI.get(bridged, MOI.ObjectiveValue())
-
-    @test obj ≈ -9.33333 atol = 1e-2
-
-    Xr = MOI.get(bridged, MOI.VariablePrimal(), X)
-
-    @test Xr ≈ [1.3333, .0, 1.3333, .0, .0, .0, .0, .0, .0, .0] atol = 1e-2
-
-end
-
-@testset "LP in SDP INEQ form" begin
-
-    MOI.empty!(optimizer)
-    @test MOI.is_empty(optimizer)
-
-    # add 10 variables - only diagonal is relevant
-    X = MOI.add_variables(optimizer, 3)
-
-    # add sdp constraints - only ensuring positivenesse of the diagonal
-    vov = MOI.VectorOfVariables(X)
-    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
-
-    c1 = MOI.add_constraint(optimizer, 
-        MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,1],[
-            MOI.ScalarAffineTerm(2.0, X[1]),
-            MOI.ScalarAffineTerm(1.0, X[3]),
-        ]), [-4.0]), MOI.Nonpositives(1))
-
-    c2 = MOI.add_constraint(optimizer, 
-        MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,1],[
-            MOI.ScalarAffineTerm(1.0, X[1]),
-            MOI.ScalarAffineTerm(2.0, X[3]),
-        ]), [-4.0]), MOI.Nonpositives(1))
-
-    MOI.set(optimizer, 
-        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 3.0], [X[1], X[3]]), 0.0)
-        )
-    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.optimize!(optimizer)
-
-    obj = MOI.get(optimizer, MOI.ObjectiveValue())
-
-    @test obj ≈ 9.33333 atol = 1e-2
-
-    Xr = MOI.get(optimizer, MOI.VariablePrimal(), X)
-
-    @test Xr ≈ [1.3333, .0, 1.3333] atol = 1e-2
-
-    c1_d = MOI.get(optimizer, MOI.ConstraintDual(), c1)
-    c2_d = MOI.get(optimizer, MOI.ConstraintDual(), c2)
-
-end
-
-@testset "SDP from MOI" begin
-    # min X[1,1] + X[2,2]    max y
-    #     X[2,1] = 1         [0   y/2     [ 1  0
-    #                         y/2 0    <=   0  1]
-    #     X >= 0              y free
-    # Optimal solution:
-    #
-    #     ⎛ 1   1 ⎞
-    # X = ⎜       ⎟           y = 2
-    #     ⎝ 1   1 ⎠
-    MOI.empty!(optimizer)
-    @test MOI.is_empty(optimizer)
-
-    X = MOI.add_variables(optimizer, 3)
-
-    vov = MOI.VectorOfVariables(X)
-    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
-
-    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))
-
-    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[3]]), 0.0))
-
-    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.optimize!(optimizer)
-
-    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
-
-    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
-    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT
-
-    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2 atol=1e-2
-
-    Xv = ones(3)
-    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
-    # @test MOI.get(optimizer, MOI.ConstraintPrimal(), cX) ≈ Xv atol=1e-2
-
-    # @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ 2 atol=1e-2
-    # @show MOI.get(optimizer, MOI.ConstraintDual(), c)
-
-end
-
-@testset "Double SDP from MOI" begin
-    # solve simultaneously two of these:
-    # min X[1,1] + X[2,2]    max y
-    #     X[2,1] = 1         [0   y/2     [ 1  0
-    #                         y/2 0    <=   0  1]
-    #     X >= 0              y free
-    # Optimal solution:
-    #
-    #     ⎛ 1   1 ⎞
-    # X = ⎜       ⎟           y = 2
-    #     ⎝ 1   1 ⎠
-    MOI.empty!(optimizer)
-    @test MOI.is_empty(optimizer)
-
-    X = MOI.add_variables(optimizer, 3)
-    Y = MOI.add_variables(optimizer, 3)
-
-    vov = MOI.VectorOfVariables(X)
-    vov2 = MOI.VectorOfVariables(Y)
-    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
-    cY = MOI.add_constraint(optimizer, vov2, MOI.PositiveSemidefiniteConeTriangle(2))
-
-    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))
-    c2 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, Y[2]))], [-1.0]), MOI.Zeros(1))
-
-    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[end], Y[1], Y[end]]), 0.0))
-
-    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.optimize!(optimizer)
-
-    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
-
-    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
-    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT
-
-    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2*2 atol=1e-2
-
-    Xv = ones(3)
-    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
-    Yv = ones(3)
-    @test MOI.get(optimizer, MOI.VariablePrimal(), Y) ≈ Yv atol=1e-2
-    # @test MOI.get(optimizer, MOI.ConstraintPrimal(), cX) ≈ Xv atol=1e-2
-
-    # @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ 2 atol=1e-2
-    # @show MOI.get(optimizer, MOI.ConstraintDual(), c)
-
-end
-
-@testset "SDP with duplicates from MOI" begin
-
-    cache = MOIU.UniversalFallback(MOIU.Model{Float64}());
-    #optimizer0 = SCS.Optimizer(linear_solver=SCS.Direct, eps=1e-8);
-    optimizer0 = ProxSDP.Optimizer()#linear_solver=SCS.Direct, eps=1e-8);
-    MOI.empty!(cache);
-    optimizer1 = MOIU.CachingOptimizer(cache, optimizer0);
-    optimizer = MOIB.full_bridge_optimizer(optimizer1, Float64);
-
-    MOI.empty!(optimizer)
-
-    x = MOI.add_variable(optimizer)
-    X = [x, x, x]
-
-    vov = MOI.VectorOfVariables(X)
-    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(2))
-
-    c = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [-1.0]), MOI.Zeros(1))
-
-    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[1], X[3]]), 0.0))
-
-    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.optimize!(optimizer)
-
-    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
-
-    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
-    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.FEASIBLE_POINT
-
-    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2 atol=1e-2
-
-    Xv = ones(3)
-    @test MOI.get(optimizer, MOI.VariablePrimal(), X) ≈ Xv atol=1e-2
-
-end
-
-@testset "SDP from Wikipedia" begin
-    # https://en.wikipedia.org/wiki/Semidefinite_programming
-    MOI.empty!(optimizer)
-    @test MOI.is_empty(optimizer)
-
-    X = MOI.add_variables(optimizer, 6)
-
-    vov = MOI.VectorOfVariables(X)
-    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(3))
-
-    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[1]))], [-1.0]), MOI.Zeros(1))
-    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[3]))], [-1.0]), MOI.Zeros(1))
-    cd1 = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[6]))], [-1.0]), MOI.Zeros(1))
-
-    c12_ub = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[2]))], [0.1]), MOI.Nonpositives(1))   # x <= -0.1 -> x + 0.1 <= 0
-    c12_lb = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(-1.0, X[2]))], [-0.2]), MOI.Nonpositives(1)) # x >= -0.2 -> -x + -0.2 <= 0
-
-    c23_ub = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(1.0, X[5]))], [-0.5]), MOI.Nonpositives(1)) # x <= 0.5 ->  x - 0.5 <= 0
-    c23_lb = MOI.add_constraint(optimizer, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1,MOI.ScalarAffineTerm(-1.0, X[5]))], [0.4]), MOI.Nonpositives(1)) # x >= 0.4 -> -x + 0.4 <= 0
-
-    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [X[4]]), 0.0))
-
-    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-
-    MOI.optimize!(optimizer)
-
-    obj = MOI.get(optimizer, MOI.ObjectiveValue())
-    @test obj ≈ -0.978 atol=1e-2
-
-    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-
-    MOI.optimize!(optimizer)
-
-    obj = MOI.get(optimizer, MOI.ObjectiveValue())
-    @test obj ≈ 0.872 atol=1e-2
-
+@testset "ProxSDP MOI Units tests" begin
+    include("moi_proxsdp_unit.jl")
 end
 
 @testset "MIMO Sizes" begin
@@ -557,17 +291,17 @@ end
 @testset "Full eig" begin
     MOIT.psdt0vtest(
         MOIB.full_bridge_optimizer(optimizer_full, Float64),
-        MOIT.TestConfig(atol=1e-3, rtol=1e-3, duals = false)
+        MOIT.Config{Float64}(atol=1e-3, rtol=1e-3, duals = false)
         )
 end
 
 @testset "Print" begin
-    MOIT.linear15test(optimizer_print, MOIT.TestConfig(atol=1e-3, rtol=1e-3))
+    MOIT.linear15test(optimizer_print, MOIT.Config{Float64}(atol=1e-3, rtol=1e-3))
 end
 
 @testset "Unsupported argument" begin
     MOI.empty!(cache)
-    @test_throws ErrorException  optimizer_unsupportedarg = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(unsupportedarg = 10))
+    @test_throws ErrorException optimizer_unsupportedarg = MOIU.CachingOptimizer(cache, ProxSDP.Optimizer(unsupportedarg = 10))
     # @test_throws ErrorException MOI.optimize!(optimizer_unsupportedarg)
 end
 
diff --git a/test/quad_knapsack.jl b/test/quad_knapsack.jl
index 930a2e3..9d63db9 100644
--- a/test/quad_knapsack.jl
+++ b/test/quad_knapsack.jl
@@ -1,4 +1,4 @@
-using StatsBase
+import StatsBase
 
 function quad_knapsack(solver, seed)
     rng = Random.MersenneTwister(seed)
@@ -28,7 +28,7 @@ function quad_knapsack(solver, seed)
     C = zeros((n, n))
     for i in 1:n
         for j in 1:n
-            if sample(rng, [1, 0],Weights([delta, 1.0 - delta])) == 1
+            if StatsBase.sample(rng, [1, 0], StatsBase.Weights([delta, 1.0 - delta])) == 1
                 c_ = - Random.rand(rng, 1:100)
                 C[i, j] = c_
                 C[j, i] = c_
@@ -64,7 +64,7 @@ function quad_knapsack(solver, seed)
     else
         XX = getvalue.(X)
     end
-    rank = length([eig for eig in eigen(XX).values if eig > 1e-10])
+    rank = length([eig for eig in LinearAlgebra.eigen(XX).values if eig > 1e-10])
     @show rank
     @show diag(XX)
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index e3a97ca..f98ee96 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,27 +2,6 @@ path = joinpath(dirname(@__FILE__), "..", "..")
 push!(Base.LOAD_PATH, path)
 datapath = joinpath(dirname(@__FILE__), "data")
 
-using ProxSDP
-
-sympackedlen(n) = div(n*(n+1), 2)
-sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
-function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
-    n = sympackeddim(length(v))
-    n1, n2 = size(out)
-    @assert n == n1 == n2
-    c = 0
-    for j in 1:n, i in 1:j
-        c += 1
-        out[i,j] = v[c]
-    end
-    return out
-end
-function ivech(v::AbstractVector{T}) where T
-    n = sympackeddim(length(v))
-    out = zeros(n, n)
-    ivech!(out, v)
-    
-    return out
-end
+import ProxSDP
 
 include("moitest.jl")
\ No newline at end of file