diff --git a/src/MOIWrapper.jl b/src/MOIWrapper.jl
index 1733ee99..7533f260 100644
--- a/src/MOIWrapper.jl
+++ b/src/MOIWrapper.jl
@@ -359,7 +359,7 @@ function processconstraints!(optimizer::Optimizer, src::MOI.ModelLike, idxmap, r
     set_constant = zeros(Float64, m)
     # loop over constraints and modify A, l, u and constants
     for (F, S) in MOI.get(src, MOI.ListOfConstraints())
-        processconstraints!((I, J, V), b, COSMOconvexSets, constant, set_constant, src, idxmap, rowranges, F, S)
+        processconstraints!(optimizer, (I, J, V), b, COSMOconvexSets, constant, set_constant, src, idxmap, rowranges, F, S)
     end
     optimizer.set_constant = set_constant
     # subtract constant from right hand side
@@ -379,7 +379,7 @@ function processconstraints!(optimizer::Optimizer, src::MOI.ModelLike, idxmap, r
     return nothing
 end
 
-function processconstraints!(triplets::SparseTriplets, b::Vector, COSMOconvexSets::Array{COSMO.AbstractConvexSet{Float64}}, constant::Vector{Float64}, set_constant::Vector{Float64},
+function processconstraints!(optimizer, triplets::SparseTriplets, b::Vector, COSMOconvexSets::Array{COSMO.AbstractConvexSet{Float64}}, constant::Vector{Float64}, set_constant::Vector{Float64},
         src::MOI.ModelLike, idxmap, rowranges::Dict{Int, UnitRange{Int}},
         F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
     cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
@@ -394,7 +394,7 @@ function processconstraints!(triplets::SparseTriplets, b::Vector, COSMOconvexSet
             set_constant[rows] =  MOI.constant(s)
         end
         processConstraint!(triplets, f, rows, idxmap, s)
-        processSet!(b, rows, COSMOconvexSets, s)
+        processSet!(b, rows, COSMOconvexSets, s, optimizer.inner.settings)
     end
     nothing
 end
@@ -425,7 +425,7 @@ end
 
 
 # process function like f(x) = coeff*x_1
-function processConstraint!(triplets::SparseTriplets, f::MOI.ScalarAffineFunction, row::Int, idxmap::MOIU.IndexMap, s::MOI.AbstractSet)
+function processConstraint!(triplets::SparseTriplets, f::MOI.ScalarAffineFunction, row::Int, idxmap, s::MOI.AbstractSet)
     (I, J, V) = triplets
     for term in f.terms
         push!(I, row)
@@ -434,7 +434,7 @@ function processConstraint!(triplets::SparseTriplets, f::MOI.ScalarAffineFunctio
     end
 end
 
-function processConstraint!(triplets::SparseTriplets, f::MOI.VectorOfVariables, rows::UnitRange{Int}, idxmap::MOIU.IndexMap, s::MOI.AbstractSet)
+function processConstraint!(triplets::SparseTriplets, f::MOI.VectorOfVariables, rows::UnitRange{Int}, idxmap, s::MOI.AbstractSet)
     (I, J, V) = triplets
     cols = zeros(length(rows))
     for (i, var) in enumerate(f.variables)
@@ -470,75 +470,75 @@ end
 ##############################
 
 # process the following sets Union{LessThan, GreaterThan, EqualTo}
-function processSet!(b::Vector, row::Int, cs, s::LessThan)
+function processSet!(b::Vector, row::Int, cs, s::LessThan, settings::Settings)
     b[row] += s.upper
     push!(cs, COSMO.Nonnegatives{Float64}(1))
     nothing
 end
-function processSet!(b::Vector, row::Int, cs, s::GreaterThan)
+function processSet!(b::Vector, row::Int, cs, s::GreaterThan, settings::Settings)
     b[row] -= s.lower
     push!(cs, COSMO.Nonnegatives{Float64}(1))
     nothing
 end
-function processSet!(b::Vector, row::Int, cs, s::EqualTo)
+function processSet!(b::Vector, row::Int, cs, s::EqualTo, settings::Settings)
     b[row] += s.value
     push!(cs, COSMO.ZeroSet{Float64}(1))
     nothing
 end
 
-function processSet!(b::Vector, row::Int, cs, s::MOI.Interval)
+function processSet!(b::Vector, row::Int, cs, s::MOI.Interval, settings::Settings)
     push!(cs, COSMO.Box{Float64}([s.lower], [s.upper]))
     nothing
 end
 
 # process the following sets Zeros
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::Zeros)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::Zeros, settings::Settings)
     push!(cs, COSMO.ZeroSet{Float64}(length(rows)))
     nothing
 end
 
 # process the following sets Union{Zeros, Nonnegatives, Nonpositives}
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonnegatives)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonnegatives, settings::Settings)
     push!(cs, COSMO.Nonnegatives{Float64}(length(rows)))
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonpositives)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonpositives, settings::Settings)
     push!(cs, COSMO.Nonnegatives{Float64}(length(rows)))
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::SOC)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::SOC, settings::Settings)
     push!(cs, COSMO.SecondOrderCone{Float64}(length(rows)))
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeSquare)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeSquare, settings::Settings)
     push!(cs, COSMO.PsdCone{Float64}(length(rows)))
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeTriangle)
-    push!(cs, COSMO.PsdConeTriangle{Float64}(length(rows)))
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeTriangle, settings::Settings)
+    push!(cs, settings.psd_projector(length(rows)))
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.ExponentialCone)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.ExponentialCone, settings::Settings)
     push!(cs, COSMO.ExponentialCone{Float64}())
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualExponentialCone)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualExponentialCone, settings::Settings)
     push!(cs, COSMO.DualExponentialCone{Float64}())
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PowerCone)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PowerCone, settings::Settings)
     push!(cs, COSMO.PowerCone{Float64}(s.exponent))
     nothing
 end
 
-function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualPowerCone)
+function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualPowerCone, settings::Settings)
     push!(cs, COSMO.DualPowerCone{Float64}(s.exponent))
     nothing
 end
diff --git a/src/algebra.jl b/src/algebra.jl
index 0dfd21f4..8853fadb 100644
--- a/src/algebra.jl
+++ b/src/algebra.jl
@@ -1,4 +1,6 @@
 using  LinearAlgebra
+import LinearAlgebra.lmul!
+import LinearAlgebra.rmul!
 const  IdentityMatrix = UniformScaling{Bool}
 
 
@@ -111,11 +113,11 @@ end
 
 lmul!(L::IdentityMatrix, M::AbstractMatrix) = L.λ ? M : M .= zero(eltype(M))
 
-function lmul!(L::Diagonal, x::AbstractVector)
-	(length(L.diag) == length(x)) || throw(DimensionMismatch())
-	@. x = x * L.diag
-	return nothing
-end
+function lmul!(L::Diagonal, x::AbstractVector)	
+	(length(L.diag) == length(x)) || throw(DimensionMismatch())	
+	@. x = x * L.diag	
+	return nothing	
+end	
 
 lmul!(L::IdentityMatrix, x::AbstractVector) = L.λ ? x : x .= zero(eltype(x))
 
@@ -217,4 +219,48 @@ function is_neg_def!(X::AbstractMatrix{T}, tol::T=zero(T)) where T
 end
 
 is_pos_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(copy(X), tol)
-is_neg_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(-X, tol)
\ No newline at end of file
+is_neg_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(-X, tol)
+
+function populate_upper_triangle!(A::AbstractMatrix{T}, x::AbstractVector{T}, scaling_factor::T=T(1/sqrt(2))) where T
+	k = 0
+	 for j in 1:size(A, 2)
+		for i in 1:j
+		   k += 1
+		   if i != j
+			   A[i, j] = scaling_factor * x[k]
+		   else
+			   A[i, j] = x[k]
+		   end
+		 end
+	 end
+	 nothing
+end
+
+function extract_upper_triangle!(A::AbstractMatrix{T}, x::AbstractVector{T}, scaling_factor::T=T(sqrt(2))) where T
+   k = 0
+	 for j in 1:size(A, 2)
+		for i in 1:j
+		   k += 1
+		   if i != j
+			   x[k] = scaling_factor * A[i, j]
+		   else
+			   x[k] = A[i, j]
+		   end
+		 end
+	 end
+   nothing
+end
+
+function populate_upper_triangle(x::AbstractVector{T}, scaling_factor::T=T(1/sqrt(2))) where T
+	n = Int(1/2*(sqrt(8*length(x) + 1) - 1)) # Solution of (n^2 + n)/2 = length(x) obtained by WolframAlpha
+	A = zeros(n, n)
+	populate_upper_triangle!(A, x, scaling_factor)
+	return Symmetric(A)
+end
+
+function extract_upper_triangle(A::AbstractMatrix{T}, scaling_factor::T=T(sqrt(2))) where T
+	n = size(A, 2)
+	x = zeros(Int(n*(n + 1)/2))
+	extract_upper_triangle!(A, x, scaling_factor)
+	return x
+end
\ No newline at end of file
diff --git a/src/convexset.jl b/src/convexset.jl
index ea8ca536..4125cce6 100644
--- a/src/convexset.jl
+++ b/src/convexset.jl
@@ -317,11 +317,11 @@ mutable struct PsdConeTriangle{T} <: AbstractConvexCone{T}
     clique_ind::Int64
 
     function PsdConeTriangle{T}(dim::Int, tree_ind::Int64, clique_ind::Int64) where{T}
-        dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
-        side_dimension = Int(sqrt(0.25 + 2 * dim) - 0.5);
-        new(dim, side_dimension, zeros(side_dimension, side_dimension), PsdBlasWorkspace{T}(side_dimension), tree_ind, clique_ind)
+      dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
+      side_dimension = Int(sqrt(0.25 + 2 * dim) - 0.5);
+      new(dim, side_dimension, zeros(side_dimension, side_dimension), PsdBlasWorkspace{T}(side_dimension), tree_ind, clique_ind)
 
-    end
+  end
 end
 PsdConeTriangle(dim) = PsdConeTriangle{DefaultFloat}(dim)
 PsdConeTriangle{T}(dim::Int64) where{T} = PsdConeTriangle{T}(dim, 0, 0)
@@ -357,54 +357,28 @@ function project!(x::AbstractArray, cone::Union{PsdConeTriangle{T}, DensePsdCone
 end
 
 # Notice that we are using a (faster) in-place version that modifies the input
-function in_dual!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where{T}
+function in_dual!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where{T}
     n = cone.sqrt_dim
     populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
     return COSMO.is_pos_def!(cone.X, tol)
 end
-in_dual(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where {T} = in_dual!(x, cone, tol)
+in_dual(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where {T} = in_dual!(x, cone, tol)
 
-function in_pol_recc!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where{T}
+function in_pol_recc!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where{T}
     n = cone.sqrt_dim
     populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
     Xs = Symmetric(cone.X)
     return COSMO.is_neg_def!(cone.X, tol)
 end
-in_pol_recc(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where {T} = in_pol_recc!(x, cone, tol)
+in_pol_recc(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where {T} = in_pol_recc!(x, cone, tol)
 
 
 function allocate_memory!(cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}) where {T}
   cone.X = zeros(cone.sqrt_dim, cone.sqrt_dim)
 end
 
-function populate_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
- 	k = 0
-  	for j in 1:size(A, 2)
-     	for i in 1:j
-        	k += 1
-        	if i != j
-        		A[i, j] = scaling_factor * x[k]
-        	else
-        		A[i, j] = x[k]
-        	end
-      	end
-  	end
-  	nothing
-end
-
-function extract_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
-	k = 0
-  	for j in 1:size(A, 2)
-     	for i in 1:j
-        	k += 1
-        	if i != j
-        		x[k] = scaling_factor * A[i, j]
-        	else
-        		x[k] = A[i, j]
-        	end
-      	end
-  	end
-	nothing
+function reset_iteration_counters!(cone::AbstractConvexSet)
+  nothing
 end
 
 """
@@ -848,7 +822,7 @@ function scale!(cone::AbstractConvexCone{T}, ::AbstractVector{T}) where{T}
   return nothing
 end
 
-function rectify_scaling!(E, work, set::Union{SecondOrderCone{T}, PsdCone{T}, DensePsdCone{T}, PsdConeTriangle{T}, DensePsdConeTriangle{T}, PowerCone{T}, DualPowerCone{T}, ExponentialCone{T}, DualExponentialCone{T}}) where{T}
+function rectify_scaling!(E, work, set::Union{SecondOrderCone{T}, PsdCone{T}, DensePsdCone{T}, PsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}, DensePsdConeTriangle{T}, PowerCone{T}, DualPowerCone{T}, ExponentialCone{T}, DualExponentialCone{T}}) where{T}
   return rectify_scalar_scaling!(E, work)
 end
 rectify_scaling!(E, work, set::Union{ZeroSet{<:Real}, Nonnegatives{<:Real}, Box{<:Real}}) = false
diff --git a/src/lobpcg.jl b/src/lobpcg.jl
new file mode 100644
index 00000000..d455adfe
--- /dev/null
+++ b/src/lobpcg.jl
@@ -0,0 +1,353 @@
+using LinearAlgebra, Statistics
+using Printf, Random
+
+mutable struct LOBPCGIterable{T}
+    n::Int
+    m::Int
+    m_initial::Int
+    iteration::Int
+    A::AbstractMatrix{T}
+    λ_barrier::T
+    X::Matrix{T}
+    P::Matrix{T}
+
+    AX::Matrix{T}
+    AP::Matrix{T}
+    residual_norms::Vector{T}
+
+    λ::Vector{T}
+    condition_estimates::Vector{T}
+    active_indices::BitVector
+    tol::T
+    verbosity::Int
+    buffer_size::Int
+    delta::Int
+    which::Symbol
+    min_full_iters::Int
+    max_dim::Int
+    print_interval::Int
+
+    start_time::Float64
+
+    function LOBPCGIterable(A::AbstractMatrix{T};
+        tol::T=T(1e-4), verbosity=1,
+        buffer_size::Int=max(Int(floor(size(A, 1) / 50)), 3),
+        which=:largest,
+        λ_barrier=T(NaN),
+        max_dim=floor(size(A, 1)/4),
+        min_full_iters=0,
+        print_interval=1
+    ) where T
+        n = size(A, 1)
+
+        @assert issymmetric(A) "Matrix must be symmetric"
+        @assert which == :largest || which == :smallest "The argument 'which' can take the values :largest or :smallest" 
+        new{T}(n, 0, 0, 0, A,
+            which == :largest ? λ_barrier : -λ_barrier,
+            Matrix{T}(undef, n, 0), Matrix{T}(undef, n, 0), # X, P
+            Matrix{T}(undef, n, 0), Matrix{T}(undef, n, 0), # ΑΧ, AP
+            zeros(T, 0), # residual_norms
+            zeros(T, 0), # λ
+            [-log10(eps(T))/2], # condition_estimates
+            BitArray(undef, 0), # active_indices
+            tol, verbosity,
+            buffer_size, # buffer
+            0,
+            which,
+            min_full_iters, max_dim,
+            print_interval,
+            0.0
+        )
+    end
+end
+
+function initialize!(data::LOBPCGIterable{T}, A::AbstractMatrix{T}, X::Matrix{T}; kwargs...) where T
+    @assert data.n == size(A, 1) == size(A, 2)
+    @assert issymmetric(A) "Matrix must be symmetric"
+    data.A = A
+    initialize!(data, X; kwargs...)
+end
+
+function initialize!(data::LOBPCGIterable{T}, X::Matrix{T}; is_orthonormal=false, which=:unchanged) where T
+    data.m = size(X, 2)
+    data.m_initial = size(X, 2)
+    @assert size(X, 1) == data.n
+    if which != :unchanged
+        @assert which == :largest || which == :smallest "The argument 'which' can take the values :largest or :smallest" 
+        data.which = which
+    end
+    data.m > data.max_dim && throw(ArgumentError("Too many columns in initial matrix.")) # m_max?
+
+    if data.m < data.buffer_size
+        data.buffer_size = data.m
+    end
+    data.X = X
+    data.AX = mul_A(data.X, data)
+    restart_subspace!(data)
+    data.λ, _ = rayleigh_ritz!(data.X, data.AX; is_orthonormal=is_orthonormal)
+    data.active_indices = BitArray(undef, data.m)
+    data.active_indices .= true
+
+    append!(data.residual_norms, T(NaN))
+    data.iteration = 0
+    data.condition_estimates = [-log10(eps(T))/2]
+end
+
+function print_statistics(data::LOBPCGIterable)
+    if data.iteration == 1 && data.verbosity > 0
+        println("-----------------------")
+        println("Computing ", string(data.which), " eigenpairs of a matrix of size: ", size(data.A))
+        println("Number of initial Ritz pairs: ", size(data.X, 2), ". Tolerance: ", data.tol)
+        println("Maximum allowed subspace size: ", data.max_dim, ". Minimum full iterations: ", data.min_full_iters)
+        if !isnan(data.λ_barrier)
+            range_string = data.which == :largest ? string("(", data.λ_barrier,", Inf)") : string("(-Inf, ", -data.λ_barrier,")")
+            println("All the eigenvalues in ", range_string, " were requested. Buffer size: ", data.buffer_size)
+        end
+        println("-----------------------")
+    end
+    if data.verbosity > 1
+        if mod(data.iteration - 1, data.print_interval*10) == 0
+            @printf("Iter\t Non-conv pairs\t Residuals(min|max)\t Ritz Values(min|λmax) \t Time(s) \t Restarted\n")
+        end
+        if data.delta == 0
+            size_string = string(sum(data.active_indices), "      ")
+        else
+            size_string = string(sum(data.active_indices), " (", data.delta, ")  ")
+        end
+        min_λ = data.which == :largest ? minimum(data.λ) : -maximum(data.λ)
+        max_λ = data.which == :largest ? maximum(data.λ) : -minimum(data.λ)
+        if mod(data.iteration - 1, data.print_interval) == 0
+            @printf("%d\t %s\t %.2e | %.2e\t %+.2e | %+.2e \t %.3e \t %s\n",
+                data.iteration, size_string, minimum(data.residual_norms), maximum(data.residual_norms),
+                min_λ, max_λ, time() - data.start_time,
+                is_restarted(data) ? "YES" : "-")
+        end
+    end
+end
+
+function mul_A!(Y, X, data::LOBPCGIterable)
+    mul!(Y, data.A, X)
+    if data.which == :smallest
+        @. Y = -Y
+    end
+    return Y
+end
+mul_A(X, data::LOBPCGIterable) = return mul_A!(similar(X), X, data)
+
+function lobpcg!(data::LOBPCGIterable{Float64}, max_iter=10)
+    data.start_time = time()
+    status = :not_converged
+    while data.iteration < max_iter && status == :not_converged
+        data.iteration += 1
+        status = iterate!(data)
+    end
+
+    if data.verbosity > 0
+        if status == :excessive_column_num
+            println("Too many desired eigenvalues. Use direct eigendecomposition instead.")
+        elseif status == :excessive_column_num
+        end
+        @printf("%d eigenvalues converged in %d iterations, %.3e seconds\n", sum(.!data.active_indices), data.iteration - 1, time() - data.start_time)
+    end
+    if data.which == :smallest
+        @. data.λ = -data.λ
+    end
+    return data, status
+end
+
+function lobpcg(A::AbstractMatrix{Float64}, X::Matrix{Float64}, max_iter=10;  kwargs...)
+    data = LOBPCGIterable(A; kwargs...)
+    initialize!(data, A, X)
+    lobpcg!(data, max_iter)
+end
+
+function iterate!(data::LOBPCGIterable{T}) where T
+    R = data.AX - data.X*Diagonal(data.λ)
+    R, data.P, data.AP, converged = drop_converged_indices(data, R)
+    converged && return :converged
+    R -= data.X*(data.X'*R) # Orthogonalize w.r.t. X
+    _, success = orthonormalize!(R)
+    if !success 
+        @warn "LOBPCG: Failed to orthonormalize residual"
+        return :ill_conditioned
+    end
+
+    ### Main computational time in the following time
+    AR = mul_A(R, data)
+    ###
+    F, success = orthonormalize!(data.P)
+    if success
+        rdiv!(data.AP, F.U)
+    else
+        data.verbosity > 0 && @warn string("Orthonormalization of P on iteration ", data.iteration, " failed. Restarting subspace.")
+        restart_subspace!(data)
+    end
+    
+    # The explicit_computation flag is not consistent with Knyzev's MATLAB's implementation.
+    # As a matter of fact, Knyazev has two different MATLAB versions one of which has explicit_computation always true. See
+    # https://github.com/lobpcg/blopex/blob/3e5ad076eb41aad0d7d4204bc7c70c151dedaaf9/blopex_tools/matlab/lobpcg/lobpcg_Rev_4_13.m
+    # https://github.com/lobpcg/blopex/blob/33b5e9f7b9de3a9c28edbc0e22381f3f2edce315/blopex_matlab/driver/lobpcg.m
+    explicit_computation = any(data.residual_norms .<= eps(T)^T(0.5))
+    G, Q = compute_gramm_matrix(data, R, AR, explicit_computation)
+    success, condition_estimate = check_conditioning(data, Q, R)
+    if !success && !is_restarted(data)
+        data.verbosity > 1 && @info string("Restarting due to condition estimate: ", condition_estimate)
+        G, Q = restart_subspace!(data, G, Q)
+        success, condition_estimate = check_conditioning(data, Q, R)
+    end
+    !success && data.verbosity > 0 && @warn "Gramm matrix ill-conditioned: results unpredictable"
+    push!(data.condition_estimates, condition_estimate)
+
+    ### Main computational time in the following block
+    if !explicit_computation && is_restarted(data)
+        λ, W = eigen!(G) # Q is simply identity in this case
+    else
+        λ, W = eigen!(G, Q)
+    end
+    ###
+    print_statistics(data)
+    return update_subspaces!(data, R, AR, λ, W)
+end
+
+function check_conditioning(data, Q, R)
+    condition_estimate = log10(cond(Q)) + 1;
+    condition_estimate_mean = mean(data.condition_estimates[
+        max(1, data.iteration - 10 - Int(round(log(size(R, 2))))):data.iteration])
+    failure = (condition_estimate / condition_estimate_mean > 2 && condition_estimate > 2) || condition_estimate > 8
+    return !failure, condition_estimate
+end
+
+is_restarted(data) = length(data.P) == 0
+
+function restart_subspace!(data::LOBPCGIterable{T}) where T
+    data.P = zeros(T, data.n, 0)
+    data.AP = zeros(T, data.n, 0)
+end
+
+function restart_subspace!(data::LOBPCGIterable{T}, G, Q) where T
+    new_dim = size(G, 1) - size(data.P, 2)
+    restart_subspace!(data)
+    return Symmetric(view(G.data, 1:new_dim, 1:new_dim)), Symmetric(view(Q.data, 1:new_dim, 1:new_dim))
+end
+
+function update_subspaces!(data::LOBPCGIterable{T}, R, AR, λ, W) where T
+    subspace_diff = count(λi -> λi > data.λ_barrier, λ) - data.m
+
+    # Contract subspace if necessary
+    data.delta = 0
+    if subspace_diff < -data.buffer_size && !isnan(data.λ_barrier)
+        m_contracted = max(data.m + subspace_diff, data.buffer_size)
+        data.delta = m_contracted - data.m
+        data.m = m_contracted
+        data.active_indices = data.active_indices[-data.delta+1:end]
+    end
+    data.λ = λ[end-data.m+1:end]
+    W1 = view(W, :, size(W, 2) - data.m + 1:size(W,2))
+
+    W11 = view(W1, 1:size(data.X, 2), :); W21 = view(W1, size(data.X, 2)+1:size(W, 1), :)
+    data.P = [R data.P] * W21
+    data.X = data.X * W11 + data.P
+    data.AP = [AR data.AP] * W21
+    data.AX = data.AX * W11 + data.AP
+
+    # Expand subspace if necessary
+    if subspace_diff > -data.buffer_size/2 && !isnan(data.λ_barrier) 
+        data.delta = subspace_diff + data.buffer_size
+        data.m += data.delta
+        data.m > data.max_dim && return :excessive_column_num
+        X_ = randn(data.n, data.m - size(data.X, 2))
+        data.X = [X_ data.X]
+        # @show svdvals(data.X)
+        data.AX = [mul_A(X_, data) data.AX]
+        data.λ, _ = rayleigh_ritz!(data.X, data.AX)
+        # @show norm(data.X'*data.X - I)
+        prepend!(data.active_indices, [true for i = 1:data.delta])
+    end
+
+    return :not_converged
+end
+
+function compute_gramm_matrix(data::LOBPCGIterable{T}, R, AR, explicit_computation) where T
+    #= Computes the following two matrices efficiently
+    G = Symmetric([data.X R data.P]' * mul_A([data.X R data.P], data))
+    Q = Symmetric([data.X R data.P]' * [data.X R data.P])
+    =#
+    X, P, AX, AP = data.X, data.P, data.AX, data.AP
+    if !explicit_computation
+        XAX = diagm(0 => data.λ)
+        XX, RR, PP = I, I, I
+        XR = zeros(T, size(X, 2), size(R, 2))
+    else
+        XAX = X'*AX
+        XX, RR, PP = X'*X, R'*R, P'*P
+        XR = X'*R
+    end
+    XAR = X'*AR;    XAP = X'*AP
+    RAR = R'*AR;    RAP = R'*AP
+                    PAP = P'*AP
+    XP = X'*P; RP = R'*P
+
+    G = [XAX    XAR     XAP;
+         XAR'   RAR     RAP;
+         XAP'   RAP'    PAP]
+
+    Q = [XX     XR      XP;
+         XR'    RR      RP;
+         XP'    RP'     PP]
+    #@assert norm(Symmetric(Q) - Symmetric([data.X R data.P]' * [data.X R data.P])) <= 1e-6 norm(Symmetric(Q) - Symmetric([data.X R data.P]' * [data.X R data.P]))
+    return Symmetric(G), Symmetric(Q)
+end
+
+function drop_converged_indices(data::LOBPCGIterable{T}, R) where T
+    data.residual_norms = column_norms(R)
+    # Technically, the &= in the next command should be =
+    # We have it like this to match the original implementation in MATLAB (BLOPEX).
+    @. data.active_indices &= (data.residual_norms > data.tol)
+    converged = false
+    if data.iteration > data.min_full_iters
+        R = R[:, data.active_indices]
+        if !is_restarted(data)
+            data.P = data.P[:, data.active_indices[1:size(data.P, 2)]]
+            data.AP = data.AP[:, data.active_indices[1:size(data.AP, 2)]]
+        end
+        num_converged = length(data.active_indices) - count(data.active_indices)
+        if num_converged >= data.m_initial
+            # Terminate if all eigeivnvalues in the desired region have converged
+            converged = !any(data.active_indices[data.λ .+ data.residual_norms .> data.λ_barrier])
+        end
+        converged |= num_converged == data.m
+    end
+    
+    return R, data.P, data.AP, converged
+end
+
+function rayleigh_ritz!(X::AbstractMatrix{T}, AX::AbstractMatrix{T}; is_orthonormal=false) where T
+    if is_orthonormal
+        λ, V = eigen!(Symmetric(X' * AX))
+    else
+        λ, V = eigen!(Symmetric(X' * AX), Symmetric(X'*X))
+    end
+    copyto!(X, X * V)
+    copyto!(AX, AX * V)
+    return λ, V
+end
+
+function orthonormalize!(A::AbstractMatrix{T}) where {T}
+    # @show svdvals(A)
+    m = size(A, 2)
+    F = cholesky!(Symmetric(A'*A), Val(false); check=false)
+    if issuccess(F)
+        rdiv!(A, F.U)
+        return F, true
+    else
+        return F, false
+    end
+end
+
+function column_norms(A::AbstractMatrix{T}) where T
+    v = zeros(T, size(A, 2))
+	@inbounds @simd for i = 1:size(A, 2)
+		v[i] = norm(view(A, :, i))
+	end
+	return v
+end
\ No newline at end of file
diff --git a/src/lobpcg_projection.jl b/src/lobpcg_projection.jl
new file mode 100644
index 00000000..2cff7fbf
--- /dev/null
+++ b/src/lobpcg_projection.jl
@@ -0,0 +1,115 @@
+mutable struct PsdConeTriangleLOBPCG{T} <: AbstractConvexCone{T}
+    dim::Int
+    sqrt_dim::Int
+    iteration::Int
+    is_subspace_positive::Bool
+    X::Matrix{T}
+    U::Matrix{T}
+    lobpcg::LOBPCGIterable{T}
+    exact_projections::Int
+    lobpcg_iterations::Int
+    max_iter::Int
+    subspace_dim_history::Vector{Int}
+
+    function PsdConeTriangleLOBPCG{T}(dim::Int; buffer_size::Int=-1, max_iter::Int=10, max_subspace_dim::Int=-1, verbosity::Int=0) where {T}
+        dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
+        sqrt_dim = Int(1 / 2 * (sqrt(8 * dim + 1) - 1)) # Solution of (sqrt_dim^2 + sqrt_dim)/2 = length(x) obtained by WolframAlpha
+        sqrt_dim * (sqrt_dim + 1) / 2 == dim || throw(DomainError(dim, "dimension must be a square"))
+
+        if buffer_size < 0
+            buffer_size = Int(floor(min(max(sqrt_dim / 50, 3), 20)))
+        end
+        if max_subspace_dim < 0
+            max_subspace_dim = Int(floor(sqrt_dim/4))
+        end
+
+        lobpcg = LOBPCGIterable(zeros(T, sqrt_dim, sqrt_dim), verbosity=verbosity,
+            buffer_size=buffer_size, λ_barrier=T(0), max_dim=max_subspace_dim, min_full_iters=1)
+
+        new(dim, sqrt_dim, 0, true,
+            zeros(T, sqrt_dim, sqrt_dim), # X
+            zeros(T, sqrt_dim, 0), # U
+            lobpcg,
+            0, 0, max_iter, # exact_projections, lobpcg_iterations
+            zeros(Int, 0) # subspace_dim_history
+        )
+    end
+end
+PsdConeTriangleLOBPCG(args...; kwargs...) = PsdConeTriangleLOBPCG{DefaultFloat}(args...; kwargs...)
+
+function get_tolerance(cone::PsdConeTriangleLOBPCG{T}) where T
+    return T(max(sqrt(cone.sqrt_dim) / cone.iteration^(1.01) * 10, 1e-7))
+end
+
+function project_exact!(x::AbstractArray{T}, cone::PsdConeTriangleLOBPCG{T}) where {T}
+    cone.exact_projections += 1
+    populate_upper_triangle!(cone.X, x)
+    λ, U  = eigen!(Symmetric(cone.X))
+    zero_tol = 1e-9
+    cone.is_subspace_positive = sum(λ .> zero_tol) < cone.sqrt_dim/2
+    if cone.is_subspace_positive
+        range = max(sum(λ .<= zero_tol) - cone.lobpcg.buffer_size + 1, 1):cone.sqrt_dim
+    else
+        range = 1:min(sum(λ .< -zero_tol) + cone.lobpcg.buffer_size, cone.sqrt_dim)
+        λ .= -λ
+        populate_upper_triangle!(cone.X, x) # Restore cone.X altered by eigen!
+    end
+    cone.U = view(U, :, range)
+    construct_projection!(x, cone, view(λ, range))
+end
+
+function construct_projection!(x, cone::PsdConeTriangleLOBPCG{T}, λ) where {T}
+    append!(cone.subspace_dim_history, size(cone.U, 2))
+    scaled_U = rmul!(copy(cone.U), Diagonal(sqrt.(max.(λ, 0))))
+    if cone.is_subspace_positive
+        BLAS.syrk!('U', 'N', one(T), scaled_U, zero(T), cone.X)
+    else
+        BLAS.syrk!('U', 'N', one(T), scaled_U, one(T), cone.X)
+    end
+    extract_upper_triangle!(cone.X, x)
+end
+
+function project!(x::AbstractArray, cone::PsdConeTriangleLOBPCG{T}) where {T}
+    sqrt_dim = cone.sqrt_dim
+    cone.iteration += 1
+
+    if size(cone.U, 2) <= cone.lobpcg.max_dim && cone.iteration > 1
+        populate_upper_triangle!(cone.X, x)
+        initialize!(cone.lobpcg, Symmetric(cone.X), cone.U, is_orthonormal=true,
+            which = cone.is_subspace_positive ? :largest : :smallest)
+        cone.lobpcg.tol = get_tolerance(cone)
+        #try
+            cone.lobpcg, status = lobpcg!(cone.lobpcg, cone.max_iter)
+            cone.lobpcg_iterations += cone.lobpcg.iteration - 1
+            #=
+        catch e
+            if isa(e, PosDefException)
+                status = :error
+                @warn "LOBPCG failed (PosDefException); switching to exact eigendecompositon."
+            else
+                throw(e)
+            end
+        end
+        =#
+    else
+        status = :not_called
+    end
+
+    if status != :converged
+        return project_exact!(x, cone)
+    else
+        cone.U = cone.lobpcg.X
+        if !cone.is_subspace_positive
+            cone.lobpcg.λ .= -cone.lobpcg.λ
+        end
+        construct_projection!(x, cone, cone.lobpcg.λ)
+    end
+end
+
+function reset_iteration_counters!(cone::PsdConeTriangleLOBPCG)
+    cone.exact_projections = 0
+    cone.lobpcg_iterations = 0
+    cone.iteration = 1
+end
+
+allocate_memory!(cone::PsdConeTriangleLOBPCG) = nothing
\ No newline at end of file
diff --git a/src/printing.jl b/src/printing.jl
index 14e61625..cb10bc59 100644
--- a/src/printing.jl
+++ b/src/printing.jl
@@ -12,7 +12,12 @@ function print_header(ws::COSMO.Workspace)
 	println("Problem:  x ∈ R^{$(n)},\n          constraints: A ∈ R^{$(m)x$(n)} ($(count(!iszero, ws.p.A)) nnz),\n          matrix size to factor: $(n + m)x$(n + m) ($(nnz_in_M) nnz)")
 	for (iii, set) in enumerate(sort(ws.p.C.sets, by = x -> -x.dim))
 		set_name = split(string(typeof(set).name), ".")[end]
-		iii == 1 ? println("Sets:"*" "^5*"$(set_name) of dim: $(set.dim)") : println(" "^10*"$(set_name) of dim: $(set.dim)")
+		if isa(set, PsdCone) || isa(set, PsdConeTriangle) || isa(set, PsdConeTriangleLOBPCG)
+			set_dim = string(set.sqrt_dim, "x", set.sqrt_dim)
+		else
+			set_dim = string(set.dim)
+		end
+		iii == 1 ? println("Sets:"*" "^5*"$(set_name) of dim: $(set_dim)") : println(" "^10*"$(set_name) of dim: $(set_dim)")
 		if iii > 5
 			println(" "^10*"... and $(length(ws.p.C.sets)-5) more")
 			break
@@ -22,21 +27,33 @@ function print_header(ws::COSMO.Workspace)
 		println("Decomp:   Num of original PSD cones: $(ws.ci.num_psd_cones)\n" * " "^10*"Num decomposable PSD cones: $(ws.ci.num_decomposable)\n" * " "^10*"Num PSD cones after decomposition: $(ws.ci.num_decom_psd_cones)\n" * " "^10*"Merge Strategy: $(stringify(ws.settings.merge_strategy))")
 	end
 	println("Settings: ϵ_abs = $(@sprintf("%.1e", settings.eps_abs)), ϵ_rel = $(@sprintf("%.1e", settings.eps_rel)),\n" * " "^10 * "ϵ_prim_inf = $(@sprintf("%.1e", settings.eps_prim_inf)), ϵ_dual_inf = $(@sprintf("%.1e", settings.eps_dual_inf)),\n" * " "^10 * "ρ = $(settings.rho), σ = $(settings.sigma), α = $(settings.alpha),\n" * " "^10 * "max_iter = $(settings.max_iter),\n" * " "^10 * "scaling iter = $(settings.scaling) ($(scaling_status)),\n" * " "^10 * "check termination every $(settings.check_termination) iter,\n" * " "^10 * "check infeasibility every $(settings.check_infeasibility) iter,\n" *	 " "^10 * "KKT system solver: $(print_lin_sys(settings.kkt_solver.ObjectType))")
-
+	
 	println("Setup Time: $(round.(ws.times.setup_time*1000; digits=2))ms\n")
-	println("Iter:\tObjective:\tPrimal Res:\tDual Res:\tRho:")
+	if settings.psd_projector == PsdConeTriangleLOBPCG || (isa(settings.psd_projector, OptionsFactory) && settings.psd_projector.ObjectType == PsdConeTriangleLOBPCG)
+		println("Iter:\tObjective:\tPrimal Res:\tDual Res:\tRho:\t\tTime(s):\tSubspaces:\tLOBPCG avg it:\tExact avg it:")
+	else
+		println("Iter:\tObjective:\tPrimal Res:\tDual Res:\tRho:\t\tTime(s):")
+	end
+	flush(stdout)
 	nothing
 end
 
-function print_iteration(ws::COSMO.Workspace, iter::Int64, cost::Float64, r_prim::Float64, r_dual::Float64)
+function print_iteration(ws::COSMO.Workspace, iter::Int64, cost::Float64, r_prim::Float64, r_dual::Float64, running_time::Float64)
 	settings = ws.settings
 	if mod(iter, 1) == 0 || iter == 1 || iter == 2 || iter == settings.max_iter
 		if mod(iter, settings.check_termination) == 0
-			@printf("%d\t%.4e\t%.4e\t%.4e\t%.4e\n", iter, cost, r_prim, r_dual, ws.ρ)
+			lobpcg_dims_str, lobpcg_iterations, exact_projections = extract_lobpcg_statistics(ws)
+			if isnan(lobpcg_iterations)
+				@printf("%d\t%.4e\t%.4e\t%.4e\t%.4e\t%.2e\n", iter, cost, r_prim, r_dual, ws.ρ, running_time)
+			else
+				@printf("%d\t%.4e\t%.4e\t%.4e\t%.4e\t%.2e\t%s\t%.2f\t\t%.2f\n", iter, cost, r_prim, r_dual, ws.ρ, running_time,
+					lobpcg_dims_str, lobpcg_iterations/settings.check_termination, exact_projections/settings.check_termination)
+			end
 		else
-			@printf("%d\t%.4e\t ---\t\t\t---\n", iter, cost)
+			@printf("%d\t%.4e\t ---\t\t\t---\t%.2e\n", iter, cost, running_time)
 		end
 	end
+	flush(stdout)
 	nothing
 end
 
@@ -50,6 +67,37 @@ function stringify(merge_strategy::Union{Type{<: AbstractMergeStrategy}, Options
 	end
 end
 
+function extract_lobpcg_statistics(ws; max_printed_dimensions=3)
+	lobpcg_iterations, exact_projections = 0, 0
+	lobpcg_dims = Int[]
+	sets_counter = 0
+	for set in ws.p.C.sets
+		if isa(set, COSMO.PsdConeTriangleLOBPCG)
+			if length(lobpcg_dims) < max_printed_dimensions
+				lobpcg_dim = set.subspace_dim_history[end]
+				if set.is_subspace_positive
+					append!(lobpcg_dims, lobpcg_dim)
+				else
+					append!(lobpcg_dims, -lobpcg_dim)
+				end
+			end
+
+			lobpcg_iterations += set.lobpcg_iterations
+			exact_projections += set.exact_projections
+			set.lobpcg_iterations = 0
+			set.exact_projections = 0
+			sets_counter += 1
+		end
+	end
+	lobpcg_dims_str = string(lobpcg_dims)
+	if length(lobpcg_dims) < sets_counter
+		lobpcg_dims_str = string(lobpcg_dims_str[1:end-1], ", ...]")
+	end
+	lobpcg_dims_str = rpad(lobpcg_dims_str, 12, " ")
+	return lobpcg_dims_str, lobpcg_iterations/sets_counter, exact_projections/sets_counter
+end
+
+
 function print_result(status::Symbol, iter::Int64, cost::Float64, rt::Float64)
 	println("\n" * "-"^66 * "\n>>> Results\nStatus: $(status)\nIterations: $(iter)\nOptimal objective: $(round.(cost; digits = 4))\nRuntime: $(round.(rt; digits = 3))s ($(round.(rt * 1000; digits = 2))ms)\n")
 	nothing
diff --git a/src/projections.jl b/src/projections.jl
index aa5288d3..74c2e29a 100644
--- a/src/projections.jl
+++ b/src/projections.jl
@@ -31,4 +31,6 @@ struct CompositeConvexSet{T} <: AbstractConvexSet{T}
 end
 
 include("./splitvector.jl")
+include("./lobpcg.jl")
+include("./lobpcg_projection.jl")
 include("./convexset.jl")
diff --git a/src/settings.jl b/src/settings.jl
index a4ab7df7..d0808a19 100644
--- a/src/settings.jl
+++ b/src/settings.jl
@@ -55,6 +55,7 @@ mutable struct Settings
 	max_iter::Int64
 	verbose::Bool
 	kkt_solver::Union{Type{<: AbstractKKTSolver}, OptionsFactory{<: AbstractKKTSolver}}
+	psd_projector::Union{Type{<: AbstractConvexSet}, OptionsFactory{<: AbstractConvexSet}}
 	check_termination::Int64
 	check_infeasibility::Int64
 	scaling::Int64
@@ -86,6 +87,7 @@ mutable struct Settings
 		max_iter=2500,
 		verbose=false,
 		kkt_solver=QdldlKKTSolver,
+		psd_projector=PsdConeTriangle,
 		check_termination=40,
 		check_infeasibility=40,
 		scaling=10,
@@ -99,7 +101,7 @@ mutable struct Settings
 		RHO_MAX = 1e6,
 		RHO_TOL = 1e-4,
 		decompose = true,
-    complete_dual = false,
+    	complete_dual = false,
 		time_limit = 0.0,
 		obj_true = NaN,
 		obj_true_tol = 1e-3,
@@ -113,6 +115,6 @@ mutable struct Settings
 	if !isa(merge_strategy, OptionsFactory)
 		merge_strategy = with_options(merge_strategy)
 	end
-	new(rho, sigma, alpha, eps_abs, eps_rel, eps_prim_inf, eps_dual_inf, max_iter, verbose, kkt_solver, check_termination, check_infeasibility, scaling, MIN_SCALING, MAX_SCALING, adaptive_rho, adaptive_rho_interval, adaptive_rho_tolerance, verbose_timing, RHO_MIN, RHO_MAX, RHO_TOL, decompose, complete_dual, time_limit, obj_true, obj_true_tol, merge_strategy, compact_transformation)
+	new(rho, sigma, alpha, eps_abs, eps_rel, eps_prim_inf, eps_dual_inf, max_iter, verbose, kkt_solver, psd_projector, check_termination, check_infeasibility, scaling, MIN_SCALING, MAX_SCALING, adaptive_rho, adaptive_rho_interval, adaptive_rho_tolerance, verbose_timing, RHO_MIN, RHO_MAX, RHO_TOL, decompose, complete_dual, time_limit, obj_true, obj_true_tol, merge_strategy, compact_transformation)
 end
 end
diff --git a/src/solver.jl b/src/solver.jl
index 83535c1e..3418cbf4 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -125,7 +125,7 @@ function optimize!(ws::COSMO.Workspace)
 			end
 
 			# print iteration steps
-			settings.verbose && print_iteration(ws, iter, cost, r_prim, r_dual)
+			settings.verbose && print_iteration(ws, iter, cost, r_prim, r_dual, time() - solver_time_start)
 
 			if has_converged(ws, r_prim, r_dual)
 				status = :Solved
diff --git a/test/UnitTests/LOBPCG/lobpcg_projection_tests.jl b/test/UnitTests/LOBPCG/lobpcg_projection_tests.jl
new file mode 100644
index 00000000..85baf256
--- /dev/null
+++ b/test/UnitTests/LOBPCG/lobpcg_projection_tests.jl
@@ -0,0 +1,105 @@
+using Random, Test, LinearAlgebra
+
+verbosity = 0 # Change to 2 to see LOBPCG's output
+Random.seed!(1)
+λ = [1:6; -50.00:-0.00]
+n = length(λ)
+Q = Matrix(qr(randn(n, n)).Q)
+X = Q*Diagonal(λ)*Q'
+
+n_ = Int(n*(n + 1)/2)
+
+function exact_projection(A::Matrix)
+    l, V = eigen(Symmetric(A))
+    return Symmetric(V*Diagonal(max.(l, 0))*V')
+end
+
+function lobpcg_unit_test!(cone, X::Matrix)
+    cone.exact_projections = 0
+    cone.lobpcg_iterations = 0
+    x = COSMO.extract_upper_triangle(X)
+    X_proj_exact = exact_projection(X)
+    COSMO.project!(x, cone)
+    X_proj = COSMO.populate_upper_triangle(x)
+
+    @test norm(X_proj - X_proj_exact) < 1e-5 # This is tolerance is hand-wavy
+end
+
+function set_tolerance!(cone, tol)
+    # Currently, the tolerance of LOBPCG is determined as a
+    # monotonically decreasing function of LOBPCG's iteration number.
+    # This function (i.e. set_tolerance!) uses binary search to find
+    # the iteration number that gives a close tolerance to the specified tol.
+    iter_min = 0
+    iter_max = 1e12
+    for i = 1:100
+        cone.iteration = Int(round((iter_min + iter_max)/2))
+        if COSMO.get_tolerance(cone) < tol
+            iter_max = cone.iteration
+        else
+            iter_min = cone.iteration
+        end
+    end
+    # @show COSMO.get_tolerance(cone), tol
+end
+
+function lobpcg_unit_tests(X::Matrix)
+    cone = COSMO.PsdConeTriangleLOBPCG(n_, verbosity=verbosity, max_iter=200)
+    @testset "Initial call (exact eigendecomposition)" begin
+        lobpcg_unit_test!(cone, X)
+        @test cone.exact_projections == 1
+        @test cone.lobpcg_iterations == 0
+    end
+
+    @testset "Warm starting - already solved to tolerance" begin
+        set_tolerance!(cone, 1e-1)
+        lobpcg_unit_test!(cone, X + 1e-6*Symmetric(randn(n, n)))
+        @test cone.exact_projections == 0
+        @test cone.lobpcg_iterations == 1
+    end
+    @testset "Warm starting - a few iterations required" begin
+        set_tolerance!(cone, 1e-7)
+        # Perturb sligtly the matrix under projection
+        lobpcg_unit_test!(cone, X + 1e-6*Symmetric(randn(n, n)))
+        @test cone.exact_projections == 0
+        @test cone.lobpcg_iterations >= 2
+        @test cone.lobpcg_iterations < 10
+    end
+
+    @testset "Warm starting - expand subspace" begin
+        set_tolerance!(cone, 1e-7)
+        # Shift the spectrum
+        if cone.is_subspace_positive
+            lobpcg_unit_test!(cone, X + 3.1*I)
+        else
+            lobpcg_unit_test!(cone, X - 3.1*I)
+        end
+        @test size(cone.U, 2) > 6
+        @test cone.exact_projections == 0
+        @test cone.lobpcg_iterations < 35
+    end
+
+    @testset "Warm starting - contract subspace" begin
+        set_tolerance!(cone, 1e-7)
+        # Shift the spectrum
+        if cone.is_subspace_positive
+            lobpcg_unit_test!(cone, X - 1.1*I)
+        else
+            lobpcg_unit_test!(cone, X + 1.1*I)
+        end
+        @test size(cone.U, 2) < 6
+        @test cone.exact_projections == 0
+        @test cone.lobpcg_iterations < 35
+    end
+end
+
+@testset "LOBPCG PSD Projection" begin
+    @testset "Track Positive Eigenspace" begin
+        lobpcg_unit_tests(X)
+    end
+    @testset "Track Negative Eigenspace" begin
+        lobpcg_unit_tests(-X)
+    end
+end
+
+nothing
diff --git a/test/UnitTests/LOBPCG/lobpcg_sanity_checks.jl b/test/UnitTests/LOBPCG/lobpcg_sanity_checks.jl
new file mode 100644
index 00000000..d23982f6
--- /dev/null
+++ b/test/UnitTests/LOBPCG/lobpcg_sanity_checks.jl
@@ -0,0 +1,68 @@
+using LinearAlgebra
+using Random
+using Test
+
+function generate_test_matrix(n, n_pos, n_small, gap)
+    # Generates the n x n symmetric matrix A with
+    # n_pos eigenvalues logarithmically on [0, 1]
+    # n_small eigenvalues uniformly on [-gap, gap]
+    # and the rest uniformly on [-1, 0].
+    n_neg = n - n_pos - n_small;
+    @assert n_neg >=0 "Size of matrix is too small"
+
+    eigenvalues_positive = 10.0.^range(-10, 0, length=n_pos)
+    eigenvalues_negative = range(-1, 0, length=n_neg);
+    eigenvalues_small = gap*range(-1, 1, length=n_small);
+
+    Q = Matrix(qr(randn(n, n)).Q)
+    l = [eigenvalues_positive; eigenvalues_negative; eigenvalues_small]
+    return Symmetric(Q'*Diagonal(l)*Q)
+end
+
+Random.seed!(1)
+n = 2000; n_pos = 3; n_small = 2;
+A = generate_test_matrix(n, n_pos, n_small, 1e-10);
+X0 = randn(n, n_pos)
+max_iter = 1000
+tol=1e-7
+verbosity=1
+function check_residuals(A, data, barrier = 0.0)
+    data.X = data.X[:, data.λ .>= barrier]
+    data.λ = data.λ[data.λ .>= barrier]
+    # The code does not really guarrantie the requested tolerance. See how active_indices are defined in the Code.
+    @test all(COSMO.column_norms(A*data.X - data.X*Diagonal(data.λ)) .< 4*tol)
+    @test size(data.X, 2) >= sum(eigvals(A) .> tol)
+end
+
+@testset "LOBPCG - Basic Tests" begin
+    @testset "Matrix with predefined eigenvalues" begin
+        data, flag = COSMO.lobpcg(A, copy(X0), max_iter, tol=tol,
+            λ_barrier=0.0, verbosity=verbosity)
+        check_residuals(A, data, 0.0)
+        # Test which=:smallest
+        data, flag = COSMO.lobpcg(-A, copy(X0), max_iter, tol=tol,
+            λ_barrier=0.0, verbosity=verbosity, which=:smallest)
+        data.λ = -data.λ
+        check_residuals(A, data, 0.0)
+    end
+    @testset "Excessive desired eigenvalues" begin
+        data, flag = COSMO.lobpcg(-A, copy(X0), max_iter, tol=tol,
+            λ_barrier=0.0, verbosity=verbosity)
+        @test flag == :excessive_column_num
+    end
+
+    @testset "Matrix with no positive eigenvalues" begin
+        A = randn(n, 50); A = Symmetric(-A*A');
+        data, flag = COSMO.lobpcg(A, copy(X0), max_iter, tol=tol,
+            λ_barrier=0.0,
+            verbosity=verbosity)
+        check_residuals(A, data, 0.0)
+    end
+
+    @testset "Matrix of all zeros" begin
+        A = zeros(n, n)
+        data, flag = COSMO.lobpcg(A, copy(X0), max_iter, tol=tol,
+            λ_barrier=0.0, verbosity=verbosity)
+        check_residuals(A, data, 0.0)
+    end
+end
\ No newline at end of file
diff --git a/test/UnitTests/LOBPCG/lobpcg_stripped.m b/test/UnitTests/LOBPCG/lobpcg_stripped.m
new file mode 100755
index 00000000..8649e39a
--- /dev/null
+++ b/test/UnitTests/LOBPCG/lobpcg_stripped.m
@@ -0,0 +1,236 @@
+function [X,lambda,varargout] = lobpcg_stripped(X,A, tol, varargin)
+    [n,blockSize]=size(X);
+    verbosityLevel = 1;
+    failureFlag = 1;
+    
+    if nargin < 4
+        maxIterations = min(n, 20);
+    else
+        maxIterations = varargin{1};
+    end
+
+    %Make X orthonormal
+    % [X, flag] = cholQ(X);
+    % assert(flag == 0, 'The initial approximation is not full rank')
+    
+    % Preallocation
+    condestGhistory=zeros(1,maxIterations+1);
+    P=zeros(n,blockSize, 'like', X);
+    AP=zeros(n,blockSize, 'like', X);
+    
+    [X, lambda, AX, R] = RR(X, A);
+    assert(norm(X'*X - eye(size(X, 2)), 'fro') <= 1e-8)
+    condestGhistory(1)=-log10(eps)/2;  %if too small cause unnecessary restarts
+    active_indices = true(blockSize,1);
+    
+    for iterationNumber=1:maxIterations
+        residualNorms = vecnorm(R)';
+        
+        active_indices = (residualNorms > tol) & active_indices;
+        % sum(active_indices)
+        R = R(:, active_indices);
+        P = P(:, active_indices);
+        AP = AP(:, active_indices);
+        
+        if sum(active_indices) == 0
+            failureFlag=0; %all eigenpairs converged
+            break
+        end
+        % Making active residuals orthogonal to X
+        % assert(norm(X*(X'*R)) < 1e-8)
+        R = R - X*(X'*R);
+        
+        %Making active residuals orthonormal
+        [R, flag] = cholQ(R);
+        if flag ~= 0
+            warning('BLOPEX:lobpcg:ResidualNotFullRank',...
+                'The residual is not full rank.');
+            break
+        end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        [P, AP, X, AX, lambda, condestGhistory] = RayeleighRitz(A, ...
+            X, AX, lambda, R, P, AP, residualNorms, condestGhistory, ...
+            iterationNumber, blockSize);
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%end RR
+        
+        if false % verbosityLevel
+            fprintf('Iteration %i current block size %i \n',...
+                iterationNumber, sum(active_indices));
+            fprintf('Eigenvalues lambda %17.16e \n',lambda);
+            fprintf('Residual Norms %e \n',residualNorms');
+        end
+
+        if iterationNumber == maxIterations
+            fprintf('Warning: maximum iterations reached without converging\n')
+        else
+            R = AX - X.*lambda';
+        end
+    end
+    %The main step of the method was the CG cycle: end
+
+    %Postprocessing
+    %Making sure X's are "exactly" othonormalized by final "exact" RR
+    [X, lambda, ~, R] = RR(X, A);
+    residualNorms=vecnorm(R)';
+
+    varargout(1)={failureFlag};
+    varargout(2)={residualNorms};
+end
+
+
+function [P, AP, X, AX, lambda, condestGhistory] = RayeleighRitz(A,...
+    X, AX, lambda, R, P, AP, residualNorms, condestGhistory,...
+    iterationNumber, blockSize)
+
+    AR = A*R;
+    
+    % The Raileight-Ritz method for [X R P]
+    if residualNorms <= eps^0.6  %suggested by Garrett Moran, private
+        explicitGramFlag = 1;
+    else
+        explicitGramFlag = 0;
+    end
+
+    if iterationNumber == 1
+        restart = 1;
+    else
+        %Making active conjugate directions orthonormal
+        [P, AP, restart] = cholQ2(P, AP);
+        if restart
+            warning('BLOPEX:lobpcg:DirectionNotFullRank',...
+                'The direction matrix is not full rank.');
+        end
+    end
+
+    for cond_try=1:2           %cond_try == 2 when restart
+        [oA, oB] = getGram(X, R, P, AX, AR, AP, lambda,...
+            blockSize, explicitGramFlag, restart);
+        
+        condestG = log10(cond(oB))+1;
+        condestGmean = mean(condestGhistory(max(1,iterationNumber-10-...
+            round(log(size(R,2)))):iterationNumber));
+        if (condestG/condestGmean > 2 && condestG > 2 )|| condestG > 8
+            %black magic - need to guess the restart
+            if true %verbosityLevel
+                fprintf('Restart on step %i as condestG %5.4e \n',...
+                    iterationNumber,condestG);
+            end
+            if cond_try == 1 && ~restart
+                restart=1; %steepest descent restart for stability
+            else
+                warning('BLOPEX:lobpcg:IllConditioning',...
+                    'Gramm matrix ill-conditioned: results unpredictable');
+            end
+        else
+            break
+        end
+    end
+    [W, L]=eig(oA,oB);
+    % Shouldn't we sort the eigenvalues??
+    % assert(norm(sort(diag(L)) - diag(L)) == 0)
+    W = W(:, 1:blockSize);
+    lambda=diag(L(1:blockSize,1:blockSize));
+    
+    if ~restart
+        % X = [X R P]*W;
+        P = [R P]*W(blockSize+1:end, :);
+        AP = [AR AP]*W(blockSize+1:end, :);
+    else %use block steepest descent
+        % X = [X R]*W;
+        P = R*W(blockSize+1:end, :);
+        AP = AR*W(blockSize+1:end, :);
+    end
+   
+    X = X*W(1:blockSize,:) + P;
+    AX = AX*W(1:blockSize,:) + AP;
+    clear W;
+    %%end RR
+    condestGhistory(iterationNumber+1) = condestG;
+end
+
+function [A, flag] = cholQ(A)
+    [C, flag] = chol(A'*A);
+    if flag == 0
+        A = A/C;
+    end
+end
+
+function [A, B, flag] = cholQ2(A, B)
+    [C, flag] = chol(A'*A);
+    if flag == 0
+        A = A/C;
+        B = B/C;
+    end
+end
+
+function [X, l, AX, R] = RR(X, A)
+    XX=X'*X;
+    XX=(XX' + XX)*0.5;
+    AX = A*X;
+    XAX = X'*AX;
+    XAX = (XAX + XAX')*0.5;
+    [V, L] = eig(XAX,XX);
+    l=diag(L);
+    X = X*V;
+    AX = AX*V;
+    R = AX - X.*l';
+end
+
+function [oA, oB] = getGram(X, R, P, AX, AR, AP, lambda,...
+    blockSize, explicitGramFlag, restart)
+
+    oXAR=AX'*R;
+    oRAR=AR'*R;
+    oRAR=(oRAR'+oRAR)*0.5;
+    
+    if ~restart
+        oXAP=AX'*P;
+        oRAP=AR'*P;
+        oPAP=AP'*P;
+        oPAP=(oPAP'+oPAP)*0.5;
+        oXP=X'*P;
+        oRP=R'*P;
+    end
+
+    if explicitGramFlag
+        oXAX=AX'*X;
+        oXAX=(oXAX'+oXAX)*0.5;
+        oXX=X'*X;
+        oXX=(oXX'+oXX)*0.5;
+        oRR=R'*R;
+        oXR=X'*R;
+        oRR=(oRR'+oRR)*0.5;
+
+        if ~restart
+            oA = [oXAX      oXAR     oXAP
+                  oXAR'     oRAR     oRAP
+                  oXAP'     oRAP'    oPAP];
+        
+            oPP=P'*P;
+            oPP=(oPP'+oPP)*0.5;
+            oB = [oXX    oXR    oXP
+                  oXR'   oRR    oRP
+                  oXP'   oRP'   oPP];
+        else
+            oA = [oXAX     oXAR
+                  oXAR'    oRAR];
+              
+            oB = [oXX      oXR
+                  oXR'     eye(size(R, 2))];
+        end
+    else
+        if ~restart
+            oA = [diag(lambda)   oXAR     oXAP
+                  oXAR'          oRAR     oRAP
+                  oXAP'          oRAP'    oPAP];
+            oB = [eye(blockSize)                  zeros(blockSize,size(R, 2))   oXP
+                  zeros(blockSize,size(R, 2))'    eye(size(R, 2))               oRP
+                  oXP'                            oRP'                          eye(size(P, 2)) ];
+        else
+            oA = [diag(lambda)   oXAR
+                  oXAR'          oRAR];
+            oB = eye(blockSize+size(R, 2));
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/UnitTests/LOBPCG/test_against_matlab.jl b/test/UnitTests/LOBPCG/test_against_matlab.jl
new file mode 100644
index 00000000..e481214c
--- /dev/null
+++ b/test/UnitTests/LOBPCG/test_against_matlab.jl
@@ -0,0 +1,32 @@
+using MATLAB
+using LinearAlgebra, Random
+using Test
+include("../../../src/lobpcg.jl")
+
+Random.seed!(1)
+n = 2500 # Matrix size
+k = 100  # Number of deisred eigenvalues
+tol = 1e-4
+A = Matrix(Symmetric(randn(n, n)))
+X0 = randn(n, k)
+max_iter = 10
+verbosity = 2
+
+function compute_eigenvector_error(X1, X2)
+    error_sum = 0
+    for i = 1:k
+        error_sum += min(norm(X1[:, i] - X2[:, i]), norm(X1[:, i] + X2[:, i]))
+    end
+    return error_sum
+end
+
+@testset "Random Matrix: Compare with MATLAB's implementation" begin
+    data, converged = lobpcg(A, copy(X0), tol=tol, max_iter, which=:smallest, verbosity=verbosity)
+    blockVectorX, lambda, failureFlag = mxcall(:lobpcg_stripped, 3, X0, A, tol, max_iter)
+    # show(stdout, "text/plain", data.X); println()
+    # show(stdout, "text/plain", blockVectorX); println()
+    @test norm(lambda - data.λ[end:-1:1]) < 1e-8
+    @test compute_eigenvector_error(data.X[:, end:-1:1], blockVectorX) < 1e-8
+end
+
+nothing
\ No newline at end of file
diff --git a/test/UnitTests/print.jl b/test/UnitTests/print.jl
index 9942908f..4b736946 100644
--- a/test/UnitTests/print.jl
+++ b/test/UnitTests/print.jl
@@ -7,11 +7,12 @@ r_prim = 1.5e-3
 r_dual = 1.2e-2
 status = :Solved
 rt = 0.7
+time = 0.1
 
 
 @testset "Printing" begin
-   @test COSMO.print_iteration(ws, iter, cost, r_prim, r_dual) == nothing
-   @test COSMO.print_iteration(ws, ws.settings.check_termination, cost, r_prim, r_dual) == nothing
+   @test COSMO.print_iteration(ws, iter, cost, r_prim, r_dual, time) == nothing
+   @test COSMO.print_iteration(ws, ws.settings.check_termination, cost, r_prim, r_dual, time) == nothing
    @test COSMO.print_result(status, iter, cost, rt) == nothing
 end
 nothing
diff --git a/test/run_cosmo_tests.jl b/test/run_cosmo_tests.jl
index b5f0b60a..da407d10 100644
--- a/test/run_cosmo_tests.jl
+++ b/test/run_cosmo_tests.jl
@@ -24,6 +24,8 @@ include("./UnitTests/COSMOTestUtils.jl")
   include("./UnitTests/psd_completion.jl")
   include("./UnitTests/psd_completion_and_merging.jl")
   include("./UnitTests/clique_merging_example.jl")
+  include("./UnitTests/LOBPCG/lobpcg_sanity_checks.jl")
+  include("./UnitTests/LOBPCG/lobpcg_projection_tests.jl")
 
     # optional unittests
   if in("Pardiso",keys(Pkg.installed()))