diff --git a/Project.toml b/Project.toml
index d9c68ce..d099594 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,11 +17,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 DataFrames = "1.6"
 LinearAlgebra = "1.9, 1.10"
 LinearOperators = "2.8"
+Logging = "1.9"
 Printf = "1.9, 1.10"
-QRupdate = "1.0"
 Random = "1.9, 1.10"
 SparseArrays = "1.9, 1.10"
-Logging = "1.9"
 julia = "1.9, 1.10"
 
 [extras]
diff --git a/src/BPDual.jl b/src/BPDual.jl
index e8690f8..ee39ffe 100644
--- a/src/BPDual.jl
+++ b/src/BPDual.jl
@@ -100,6 +100,13 @@ function bpdual(
     # ------------------------------------------------------------------
     m, n = size(A)
 
+
+    work = rand(n)
+    work2 = rand(n)
+    work3 = rand(n)
+    work4 = rand(n)
+    work5 = rand(m)
+
     tracer = ASPTracer(
         Int[],                  # iteration
         Float64[],              # lambda
@@ -110,8 +117,8 @@ function bpdual(
         active = Vector{Int}([])
         state = Vector{Int}(zeros(Int, n))
         y = Vector{Float64}(zeros(Float64, m))
-        S = Matrix{Float64}(zeros(Float64, m, 0))
-        R = Matrix{Float64}(zeros(Float64, 0, 0))
+        S = Matrix{Float64}(zeros(Float64, m, n))
+        R = Matrix{Float64}(zeros(Float64, n, n))
     end
 
     if homotopy
@@ -171,6 +178,8 @@ function bpdual(
     numtrim = 0
     nprodA = 0
     nprodAt = 0
+    cur_r_size = 0
+
 
     # ------------------------------------------------------------------
     # Cold/warm-start initialization.
@@ -211,25 +220,25 @@ function bpdual(
         sL, sU = infeasibilities(bl, bu, z)
         g = b - λ*y  # Steepest-descent direction
 
-        if isempty(R)
+        if isempty(R[1:cur_r_size,1:cur_r_size])
             condS = 1
         else
-            rmin = minimum(diag(R))
-            rmax = maximum(diag(R))
+            rmin = minimum(diag(R[1:cur_r_size, 1:cur_r_size]))
+            rmax = maximum(diag(R[1:cur_r_size, 1:cur_r_size]))
             condS = rmax / rmin
         end
 
         if condS > 1e+10
             eFlag = :EXIT_SINGULAR_LS
             # Pad x with enough zeros to make it compatible with S.
-            npad = size(S, 2) - size(x, 1)
+            npad = size(S[:, 1:cur_r_size], 2) - size(x, 1)
             x = [x; zeros(npad)]
         else
-            dx, dy = newtonstep(S, R, g, x, λ)
+            dx, dy = newtonstep(S[:,1:cur_r_size], R[1:cur_r_size, 1:cur_r_size], g, x, λ)
             x .+= dx
         end
 
-        r = b - S*x
+        r = b - S[:,1:cur_r_size]*x
 
         # Print to log.
         yNorm = norm(y, 2)
@@ -273,7 +282,7 @@ function bpdual(
                 @info "\nOptimal solution found. Trimming multipliers..."
             end
             g = b - λin*y
-            trimx(x, S, R, active, state, g, b, λ, feaTol, optTol, loglevel)
+            trimx(x, S[:, 1:cur_r_size], R[1:cur_r_size, 1:cur_r_size], active, state, g, b, λ, feaTol, optTol, loglevel)
             numtrim = nact - length(active)
             nact = length(active)
         end
@@ -288,7 +297,7 @@ function bpdual(
         p = q = 0
 
         if homotopy
-            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R, S, x, y, sL, sU, λ, lamFinal)
+            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R[1:cur_r_size, 1:cur_r_size], S[:,1:cur_r_size], x, y, sL, sU, λ, lamFinal)
             nprodAt += 1
         else
             if norm(dy, Inf) < eps()        
@@ -339,8 +348,9 @@ function bpdual(
             a = A * zerovec
             nprodA += 1
             zerovec[p] = 0
-            R = qraddcol(S, R, a)
-            S = [S a]
+            qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5)  # Update R
+            cur_r_size +=1             
+            # S = [S a]
             push!(active, p)
             push!(x, 0)
         else
@@ -358,10 +368,12 @@ function bpdual(
                 _, qa = findmax(abs.(x .* dropa))
                 q = active[qa]
                 state[q] = 0
-                S = S[:, 1:nact .!= qa]
+                # S = S[:, 1:nact .!= qa]
                 deleteat!(active, qa)
                 deleteat!(x, qa)
-                R = qrdelcol(R, qa)
+                # R = qrdelcol(R, qa)
+                qrdelcol!(S, R, qa)
+                cur_r_size -=1           
             else
                 eFlag = :EXIT_OPTIMAL
             end
@@ -390,4 +402,4 @@ function bpdual(
         @info "Solution time (sec): $tottime"
     end
     return tracer
-end
+end
\ No newline at end of file
diff --git a/src/helpers.jl b/src/helpers.jl
index f13df18..a71364d 100644
--- a/src/helpers.jl
+++ b/src/helpers.jl
@@ -387,4 +387,4 @@ function htpynewlam(active, state, A, R, S, x, y, s1, s2, λ, lamFinal)
     end
     
     return x, dy, dz, step, λ, p
-end
+end
\ No newline at end of file
diff --git a/src/omp.jl b/src/omp.jl
index 14a7f32..d315c5d 100644
--- a/src/omp.jl
+++ b/src/omp.jl
@@ -46,28 +46,30 @@ function asp_omp(
     optTol::Real = 1e-05,
     gapTol::Real = 1e-06,
     pivTol::Real = 1e-12,
-    actMax::Real = Inf)
+    actMax::Union{Real, Nothing} = nothing) 
     
-    # Start the clock and size up the problem.
     time0 = time()
 
     z = A' * b
     
     m = length(b)
     n = length(z)
+
+    work = rand(size(A,2))
+    work2 = rand(size(A,2))
+    work3 = rand(size(A,2))
+    work4 = rand(size(A,2))
+    work5 = rand(size(A,1))
+
     nprodA = 0
     nprodAt = 1
 
-    # Initialize the tracer
-
     tracer = OMPTracer(
-        Int[],                  # iteration
-        Float64[],              # lambda
-        Vector{SparseVector{Float64}}() # now stores full sparse solutions
+        Int[],                 
+        Float64[],              
+        Vector{SparseVector{Float64}}() 
     )
 
-    # Print log header.
-
     if loglevel > 0
         @info "-"^124
         @info @sprintf("%-30s : %-10d    %-30s : %-10.4e", "No. rows", m, "λ", λin)
@@ -76,6 +78,7 @@ function asp_omp(
         @info "-"^124
     end
 
+
     # Initialize local variables.
     EXIT_INFO = Dict(
         :EXIT_OPTIMAL => "Optimal solution found -- full Newton step",
@@ -84,6 +87,7 @@ function asp_omp(
         :EXIT_LAMBDA => "Reached minimum value of lambda",
         :EXIT_RHS_ZERO => "b = 0. The solution is x = 0",
         :EXIT_UNCONSTRAINED => "Unconstrained solution r = b is optimal",
+        :EXIT_ACTMAX => "Max no. of active constraints reached",
         :EXIT_UNKNOWN => "unknown exit"
     )
 
@@ -92,8 +96,8 @@ function asp_omp(
     x = zeros(Float64, 0)
     zerovec = zeros(Float64, n)
     p = 0
+    cur_r_size = 0
 
-    # Quick exit if the RHS is zero.
     if norm(b, Inf) == 0
         r = zeros(m)
         eFlag = :EXIT_RHS_ZERO
@@ -113,11 +117,17 @@ function asp_omp(
         state = zeros(Int, n)
     end
     if R === nothing
-        R = Matrix{Float64}(undef, 0, 0)
+        R = Matrix{Float64}(undef, size(A,2), size(A,2))
+        S = Matrix{Float64}(undef, size(A,1), size(A,2))
     end
 
-    @info @sprintf("%4s  %8s %12s %12s %12s", "Itn", "Var", "λ", "rNorm", "xNorm")
+    if actMax === nothing
+        actMax = size(A, 2)
+    end
 
+    if loglevel>0
+        @info @sprintf("%4s  %8s %12s %12s %12s", "Itn", "Var", "λ", "rNorm", "xNorm")
+    end
 
     # Main loop.
     while true
@@ -129,20 +139,22 @@ function asp_omp(
             nprodAt += 1
             zmax = norm(z, Inf)
         else
-            x,y = csne(R, S, vec(b))
+            x,y = csne(R[1:cur_r_size, 1:cur_r_size], S[:,1:cur_r_size], vec(b))
             if norm(x, Inf) > 1e12
                 eFlag = :EXIT_SINGULAR_LS
                 break
             end
 
-            Sx = S * x
+            Sx = S[:,1:cur_r_size] * x
             r = b - Sx
         end
 
         rNorm = norm(r, 2)
         xNorm = norm(x, 1)
 
-        @info @sprintf("%4i  %8i %12.5e %12.5e %12.5e", itn, p, zmax, rNorm, xNorm)
+        if loglevel>0
+            @info @sprintf("%4i  %8i %12.5e %12.5e %12.5e", itn, p, zmax, rNorm, xNorm)
+        end
 
         # Check exit conditions.
         if eFlag != :EXIT_UNKNOWN
@@ -153,8 +165,10 @@ function asp_omp(
             eFlag = :EXIT_OPTIMAL
         elseif itn >= itnMax
             eFlag = :EXIT_TOO_MANY_ITNS
+        elseif itn == actMax
+            eFlag = :EXIT_ACTMAX
         end
-        
+
         if eFlag != :EXIT_UNKNOWN
             break
         end
@@ -180,14 +194,13 @@ function asp_omp(
         nprodA += 1
         zerovec[p] = 0
 
-        R = qraddcol(S, R, a)  # Update R
-        S = hcat(S, a)  # Expand S, active
-
+        qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5)  # Update R
+        # S = hcat(S, a)  # Expand S, active
+        cur_r_size +=1 
         push!(tracer.iteration, itn)
         push!(tracer.lambda, zmax)
         sparse_x_full = SparseVector(n, copy(active), copy(x))
         push!(tracer.solution, copy(sparse_x_full))
-        
         push!(active, p)
 
     end #while true
diff --git a/test/test_bpdn.jl b/test/test_bpdn.jl
index 106b597..257ba41 100644
--- a/test/test_bpdn.jl
+++ b/test/test_bpdn.jl
@@ -1,11 +1,11 @@
 
 # ------------------------------------------------------------------
-#  Test Basis pursuit
+#  Test BPDN & OMP
 # ------------------------------------------------------------------ 
 
 using Test, LinearAlgebra, Random, SparseArrays, ActiveSetPursuit, LinearOperators
 
-function test_bpdn()
+function test_pursuits()
     m = 600
     n = 2560
     k = 20
@@ -21,12 +21,16 @@ function test_bpdn()
     b = A * x
     
     # Solve the basis pursuit problem
-    tracer = asp_bpdn(A, b, 0.0, loglevel =0);
+    tracer_bpdn = asp_bpdn(A, b, 0.0, loglevel =0);
+    tracer_omp = asp_omp(A, b, 0.0, loglevel =0);
+    xx_bpdn, _ = tracer_bpdn[end]
+    xx_omp, _ = tracer_omp[end]
 
-    xx, λ = tracer[end]
-    pFeas = norm(A * xx - b, Inf) / max(1, norm(xx, Inf))
-    @test pFeas <= 1e-6
+    pFeas_bpdn = norm(A * xx_bpdn - b, Inf) / max(1, norm(xx_bpdn, Inf))
+    pFeas_omp = norm(A * xx_omp - b, Inf) / max(1, norm(xx_omp, Inf))
 
+    @test pFeas_bpdn <= 1e-6
+    @test pFeas_omp <= 1e-6
     # ------------------------
     # Linear operator
     # ------------------------
@@ -35,15 +39,20 @@ function test_bpdn()
     b_op = OP * x
         
     # Solve the basis pursuit problem
-    tracer_op = asp_bpdn(OP, b_op, 0.0, loglevel =0);
+    tracer_op_bpdn = asp_bpdn(OP, b_op, 0.0, loglevel =0);
+    tracer_op_omp = asp_omp(OP, b_op, 0.0, loglevel =0);
 
-    xx_op, λ_op = tracer_op[end]
-    pFeas_op = norm(OP * xx_op - b_op, Inf) / max(1, norm(xx_op, Inf))
-    @test pFeas_op <= 1e-6
+    xx_op_bpdn, _ = tracer_op_bpdn[end]
+    xx_op_omp, _ = tracer_op_omp[end]
+
+    pFeas_bpdn_op = norm(OP * xx_op_bpdn - b_op, Inf) / max(1, norm(xx_op_bpdn, Inf))
+    pFeas_omp_op = norm(OP * xx_op_omp - b_op, Inf) / max(1, norm(xx_op_omp, Inf))
+    @test pFeas_bpdn_op <= 1e-6
+    @test pFeas_omp_op <= 1e-6
 end
 
 Random.seed!(1234)
 
 for ntest = 1:10
-    test_bpdn()
+    test_pursuits()
 end
\ No newline at end of file
diff --git a/test/test_recover_decaying.jl b/test/test_recover_decaying.jl
index 687a632..ad43f98 100644
--- a/test/test_recover_decaying.jl
+++ b/test/test_recover_decaying.jl
@@ -22,7 +22,7 @@ function test_recover_decaying()
     bu = +ones(n)
 
     # Solve the basis pursuit problem
-    tracer = asp_homotopy(A, b, min_lambda = 0.0, itnMax = 400, loglevel =0) 
+    tracer = asp_homotopy(A, b, min_lambda = 0.0, actMax = 400, loglevel =0) 
     xs, λ = tracer[end]
 
     cumulative_norm = cumsum(abs.(x[sortperm(abs.(x), rev=true)]))