diff --git a/codegen/codegen.c b/codegen/codegen.c
index 2cfa4d1a..4288e01d 100644
--- a/codegen/codegen.c
+++ b/codegen/codegen.c
@@ -203,7 +203,7 @@ void write_daqp_workspace_src(FILE* f, DAQPWorkspace* work, const char* prefix){
             prefix,prefix,prefix,prefix, 0); // reuse_ind
     fprintf(f, "%sWS, %d,\n", prefix, 0); //n_active
     fprintf(f, "%d,%d,\n",0,-1); //iterations + sing_id
-    fprintf(f, "%sprox_mask, %d,\n", prefix, n); // prox_mask
+    fprintf(f, "%sprox_mask, %d,\n", prefix, work->n_prox); // prox_mask
     fprintf(f, "%f,\n",0.0); // Soft slack
     fprintf(f, "&%ssettings, \n", prefix);
     // BnB
diff --git a/include/api.h b/include/api.h
index b130ef05..87e93cd5 100644
--- a/include/api.h
+++ b/include/api.h
@@ -30,9 +30,10 @@ void daqp_quadprog(DAQPResult* res, DAQPProblem* qp,DAQPSettings* settings);
 void daqp_avi(DAQPResult *res, DAQPProblem* problem, DAQPSettings *settings);
 
 int setup_daqp(DAQPProblem *qp, DAQPWorkspace* work, c_float* setup_time);
-int setup_daqp_ldp(DAQPWorkspace *work, DAQPProblem* qp);
+int setup_daqp_main(DAQPProblem *qp, DAQPWorkspace* work, c_float* setup_time, int check_unc);
+int setup_daqp_ldp(DAQPWorkspace *work, DAQPProblem* qp, const int check_unc);
 void setup_daqp_hiqp(DAQPWorkspace *work, int* break_points, int nh);
-int setup_daqp_bnb(DAQPWorkspace* work, int nb, int ns);
+int setup_daqp_bnb(DAQPWorkspace* work, int* sense, int nb, int ns);
 int setup_daqp_avi(DAQPAVI* avi, DAQPProblem* p, DAQPWorkspace* work, c_float* setup_time);
 
 void allocate_daqp_settings(DAQPWorkspace *work);
diff --git a/include/constants.h b/include/constants.h
index 80cdade6..c59384c2 100644
--- a/include/constants.h
+++ b/include/constants.h
@@ -51,6 +51,7 @@ extern "C" {
 #define DAQP_UPDATE_d 8
 #define DAQP_UPDATE_sense 16 
 #define DAQP_UPDATE_hierarchy 32
+#define DAQP_UPDATE_unconstrained 64
 
 // CONSTRAINT MASKS 
 #define DAQP_ACTIVE 1
diff --git a/include/utils.h b/include/utils.h
index 0e922bc2..ecd415a0 100644
--- a/include/utils.h
+++ b/include/utils.h
@@ -15,6 +15,7 @@ int daqp_update_d(DAQPWorkspace *work, c_float *bupper, c_float *blower);
 int daqp_check_bounds(DAQPWorkspace* work, c_float* bupper, c_float* blower);
 void daqp_normalize_Rinv(DAQPWorkspace *work);
 int daqp_normalize_M(DAQPWorkspace *work);
+int daqp_check_unconstrained(DAQPWorkspace* work, const int mask);
 
 int daqp_update_avi(DAQPAVI *avi, DAQPProblem *problem);
 int daqp_lu(c_float* A, int* P, int n);
diff --git a/interfaces/daqp-julia/src/api.jl b/interfaces/daqp-julia/src/api.jl
index a65610a1..6fe975a6 100644
--- a/interfaces/daqp-julia/src/api.jl
+++ b/interfaces/daqp-julia/src/api.jl
@@ -55,7 +55,7 @@ function quadprog(H::Union{Matrix{Float64}, Cholesky},f::Vector{Float64},
     d = DAQPBase.Model() 
     !isnothing(settings) && DAQPBase.settings(d,settings)
     exitflag,setup_time = DAQPBase.setup(d,QPj(H,f,A,bupper,blower,sense;A_rowmaj);
-                                         primal_start,dual_start)
+                                         primal_start,dual_start, check_unconstrained=true)
     return DAQPBase.solve(d;setup_time);
 end
 
@@ -177,7 +177,6 @@ is interpreted as
 function avi(H::Matrix{Float64},f::Vector{Float64}, 
         A::Matrix{Float64},bupper::Vector{Float64},blower::Vector{Float64}=Float64[],sense::Vector{Cint}=Cint[];A_rowmaj=false,settings=nothing, primal_start::Vector{Cdouble}=Cdouble[],
         dual_start::Vector{Cdouble}=Cdouble[])
-    #return quadprog(QPj(H,f,A,bupper,blower,sense;A_rowmaj,is_avi=true);settings)
     d = DAQPBase.Model() 
     exitflag,setup_time = DAQPBase.setup(d,QPj(H,f,A,bupper,blower,sense;A_rowmaj,is_avi=true);
                                          primal_start,dual_start)
@@ -228,8 +227,8 @@ function delete!(daqp::DAQPBase.Model)
     end
 end
 
-function setup(daqp::DAQPBase.Model, qp::DAQPBase.QPj;
-        primal_start::Vector{Cdouble}=Cdouble[],dual_start::Vector{Cdouble}=Cdouble[])
+function setup(daqp::DAQPBase.Model, qp::DAQPBase.QPj;primal_start::Vector{Cdouble}=Cdouble[],
+        dual_start::Vector{Cdouble}=Cdouble[],check_unconstrained=false)
     daqp.qpj = qp
     daqp.qpc = DAQPBase.QPc(daqp.qpj)
     old_settings = settings(daqp); # in case setup fails
@@ -252,7 +251,7 @@ function setup(daqp::DAQPBase.Model, qp::DAQPBase.QPj;
     end
     
     # Handle AVI with other setup 
-    exitflag = ccall((:setup_daqp,DAQPBase.libdaqp),Cint,(Ptr{DAQPBase.QPc}, Ptr{DAQPBase.Workspace}, Ptr{Cdouble}), daqp.qpc_ptr, daqp.work, setup_time)
+    exitflag = ccall((:setup_daqp_main,DAQPBase.libdaqp),Cint,(Ptr{DAQPBase.QPc}, Ptr{DAQPBase.Workspace}, Ptr{Cdouble},Cint), daqp.qpc_ptr, daqp.work, setup_time, Cint(check_unconstrained))
     if(exitflag < 0)
         # XXX: if setup fails DAQP currently clears settings
         ccall((:allocate_daqp_settings,DAQPBase.libdaqp),Nothing,(Ptr{DAQPBase.Workspace},),daqp.work)
@@ -329,7 +328,8 @@ function settings(p::Ptr{DAQPBase.Workspace},changes::Dict{Symbol,<:Any})
     return new_settings;
 end
 
-function update(daqp::DAQPBase.Model, H,f,A,bupper,blower,sense=nothing,break_points=nothing) 
+function update(daqp::DAQPBase.Model, H,f,A,bupper,blower,sense=nothing,break_points=nothing,
+        check_unconstrained = true)
     update_mask = Cint(0);
     work = unsafe_load(daqp.work);
     if(!isnothing(H) && work.n == size(H,1) && work.n == size(H,2))
@@ -362,6 +362,7 @@ function update(daqp::DAQPBase.Model, H,f,A,bupper,blower,sense=nothing,break_po
         daqp.qpj.break_points .= break_points
         update_mask+=32
     end
+    #check_unconstrained && (update_mask += 64) # Enable shortcut for unconstrained optimum
     daqp.qpc = QPc(daqp.qpj);
     unsafe_store!(daqp.qpc_ptr, daqp.qpc)
 
diff --git a/interfaces/daqp-julia/test/core_tests.jl b/interfaces/daqp-julia/test/core_tests.jl
index 5c71efaf..c8819900 100644
--- a/interfaces/daqp-julia/test/core_tests.jl
+++ b/interfaces/daqp-julia/test/core_tests.jl
@@ -494,3 +494,19 @@ end
     @test norm(x3 .- (-5.0)) < tol  # all at lower bound
 end
 
+@testset "Unconstrained shortcut" begin
+    H = diagm(ones(10))
+    f = zeros(10)
+    A = randn(10000,10)
+    b = rand(10000)
+    tquad = @elapsed x,fval,exitflag,info = quadprog(H,f,A,b)
+    @test norm(x) < tol
+    tsetup = @elapsed begin
+        d = DAQPBase.Model()
+        setup(d,H,f,A,b)
+        x,fval,exitflag,info2 = solve(d)
+        @test norm(x) < tol
+    end
+    @test tquad < 10*tsetup #Should be orders of magnitude faster
+end
+
diff --git a/interfaces/daqp-matlab/daqp.m b/interfaces/daqp-matlab/daqp.m
index 6c1116a3..066d803a 100644
--- a/interfaces/daqp-matlab/daqp.m
+++ b/interfaces/daqp-matlab/daqp.m
@@ -15,7 +15,7 @@
             if nargin < 7, primal_start = []; end
             if nargin < 8, dual_start = []; end
             d = daqp();
-            [exitflag,setup_time] = d.setup(H,f,A,bupper,blower,sense,[],0,primal_start,dual_start);
+            [exitflag,setup_time] = d.setup(H,f,A,bupper,blower,sense,[],0,primal_start,dual_start,1);
             if(exitflag <0)
                 x = [];fval=[];info=[];
                 return; 
@@ -104,7 +104,7 @@ function delete(this)
             [x,fval,exitflag,info] = daqpmex('solve', this.work_ptr,...
                 this.H,this.f,this.A,this.bupper,this.blower,this.sense,this.break_points);
         end
-        function [exitflag,setup_time] = setup(this,H,f,A,bupper,blower,sense,break_points,problem_type,primal_start,dual_start)
+        function [exitflag,setup_time] = setup(this,H,f,A,bupper,blower,sense,break_points,problem_type,primal_start,dual_start,check_unconstrained)
             if(nargin < 8)
                 break_points = [];
             end
@@ -117,6 +117,9 @@ function delete(this)
             if(nargin < 11)
                 dual_start = [];
             end
+            if(nargin < 11)
+                check_unconstrained = 0;
+            end
 
             % TODO Check validity
             % TODO match double/single with c_float... 
@@ -143,7 +146,7 @@ function delete(this)
             this.break_points= int32(break_points);
             [exitflag,setup_time] = daqpmex('setup', this.work_ptr,this.H,this.f,...
                 this.A,this.bupper,this.blower,this.sense,this.break_points,int32(problem_type),...
-                primal_start,dual_start);
+                primal_start,dual_start,check_unconstrained);
         end
 
         function settings = settings(this,varargin)
diff --git a/interfaces/daqp-matlab/daqpmex.c b/interfaces/daqp-matlab/daqpmex.c
index 91522316..9857d646 100644
--- a/interfaces/daqp-matlab/daqpmex.c
+++ b/interfaces/daqp-matlab/daqpmex.c
@@ -101,9 +101,10 @@ void mexFunction( int nlhs, mxArray *plhs[],
 	  } else if (nrhs > 10 && !mxIsEmpty(prhs[10])) {
 		daqp_primal_init_active(qp, (c_float *)mxGetPr(prhs[10]));
 	  }
+      int check_unconstrained = mxGetM(prhs[12]);
 
 	  c_float solve_time;
-	  error_flag = setup_daqp(qp,work,&solve_time);
+	  error_flag = setup_daqp_main(qp,work,&solve_time,check_unconstrained);
 	  if(error_flag < 0){
 		free(work->qp);
 		work->qp = NULL;
diff --git a/interfaces/daqp-python/daqp.pxd b/interfaces/daqp-python/daqp.pxd
index 08e05391..4e9c505f 100644
--- a/interfaces/daqp-python/daqp.pxd
+++ b/interfaces/daqp-python/daqp.pxd
@@ -73,7 +73,7 @@ cdef extern from "api.h":
     cdef extern nogil:
         void daqp_set_primal_start(DAQPWorkspace *work, double *x)
     cdef extern nogil:
-        int setup_daqp(DAQPProblem *qp, DAQPWorkspace *work, double *setup_time)
+        int setup_daqp_main(DAQPProblem *qp, DAQPWorkspace *work, double *setup_time, int check_unc)
     cdef extern nogil:
         void daqp_solve(DAQPResult *res, DAQPWorkspace *work)
     cdef extern nogil:
diff --git a/interfaces/daqp-python/daqp.pyx b/interfaces/daqp-python/daqp.pyx
index f252791a..90e2e7d8 100644
--- a/interfaces/daqp-python/daqp.pyx
+++ b/interfaces/daqp-python/daqp.pyx
@@ -43,7 +43,7 @@ cdef void _solve_warm_start(DAQPProblem* problem, DAQPSettings* settings,
         work.settings = settings
         setup_time_c = 0.0
         with nogil:
-            setup_flag = setup_daqp(problem, work, &setup_time_c)
+            setup_flag = setup_daqp_main(problem, work, &setup_time_c,1)
         res.setup_time = setup_time_c
         res.exitflag  = setup_flag
         if setup_flag >= 0:
@@ -416,7 +416,7 @@ cdef class Model:
 
         # ---- call setup_daqp ----
         with nogil:
-            setup_flag = setup_daqp(&self._qp, self._work, &setup_time_c)
+            setup_flag = setup_daqp_main(&self._qp, self._work, &setup_time_c,0)
 
         if setup_flag < 0:
             # setup_daqp freed settings on failure; re-allocate and restore.
diff --git a/src/api.c b/src/api.c
index 5a6fe7ff..2bf75187 100644
--- a/src/api.c
+++ b/src/api.c
@@ -58,7 +58,7 @@ void daqp_quadprog(DAQPResult *res, DAQPProblem* qp, DAQPSettings *settings){
 
     DAQPWorkspace work;
     work.settings = settings;
-    setup_flag = setup_daqp(qp,&work,&(res->setup_time));
+    setup_flag = setup_daqp_main(qp,&work,&(res->setup_time),1);
     res->exitflag = setup_flag;
 
     if(setup_flag >= 0){
@@ -70,14 +70,19 @@ void daqp_quadprog(DAQPResult *res, DAQPProblem* qp, DAQPSettings *settings){
     }
 }
 
-// XXX should be very similar to quadprog now
 void daqp_avi(DAQPResult *res, DAQPProblem* problem, DAQPSettings *settings){
     // Set the flag correctly
+    problem->problem_type = 1;
     daqp_quadprog(res,problem,settings);
 }
 
 // Setup workspace and transform QP to LDP
 int setup_daqp(DAQPProblem* qp, DAQPWorkspace *work, c_float* setup_time){
+    return setup_daqp_main(qp,work,setup_time,0);//  
+}
+
+// Internal function for setting workspace and transform QP to LDP
+int setup_daqp_main(DAQPProblem* qp, DAQPWorkspace *work, c_float* setup_time, int check_unc){
     int errorflag;
     int own_settings=1;
     (void)setup_time; // avoids warning when compiling without profiling
@@ -122,14 +127,15 @@ int setup_daqp(DAQPProblem* qp, DAQPWorkspace *work, c_float* setup_time){
         work->avi= malloc(sizeof(DAQPAVI));
         allocate_daqp_avi(work->avi,qp->n);
     }
-    errorflag = setup_daqp_ldp(work,qp);
+
+    errorflag = setup_daqp_bnb(work, qp->sense, nb, ns);
     if(errorflag < 0){
         if(own_settings==0) work->settings = NULL;
         free_daqp_workspace(work);
         return errorflag;
     }
 
-    errorflag = setup_daqp_bnb(work, nb, ns);  
+    errorflag = setup_daqp_ldp(work,qp,check_unc);
     if(errorflag < 0){
         if(own_settings==0) work->settings = NULL;
         free_daqp_workspace(work);
@@ -146,7 +152,7 @@ int setup_daqp(DAQPProblem* qp, DAQPWorkspace *work, c_float* setup_time){
 }
 
 //  Setup LDP from QP  
-int setup_daqp_ldp(DAQPWorkspace *work, DAQPProblem *qp){
+int setup_daqp_ldp(DAQPWorkspace *work, DAQPProblem *qp, const int check_unc){
     // Always update M, d and sense
     int update_mask = DAQP_UPDATE_M+DAQP_UPDATE_d+DAQP_UPDATE_sense; 
     int error_flag;
@@ -164,22 +170,26 @@ int setup_daqp_ldp(DAQPWorkspace *work, DAQPProblem *qp){
         update_mask+=DAQP_UPDATE_v;
     }
 
+    // For LPs (no Hessian), mark all directions as needing proximal regularisation.
+    // This lets daqp_solve dispatch to daqp_prox based on n_prox rather than eps_prox.
+    if(qp->H == NULL && qp->f!=NULL)
+        work->n_prox = work->n;
+
     // Allocate memory for LDP
     allocate_daqp_ldp(work, qp->n, qp->m, qp->ms, alloc_R, alloc_v);
 
     // Update hierarchy if hqp
     if(qp->nh > 1) update_mask += DAQP_UPDATE_hierarchy;
 
+    if(check_unc) update_mask += DAQP_UPDATE_unconstrained; 
+
     // Form LDP
     error_flag = daqp_update_ldp(update_mask, work, qp);
     if(error_flag<0){
         free_daqp_ldp(work);
         return error_flag;
     }
-    // For LPs (no Hessian), mark all directions as needing proximal regularisation.
-    // This lets daqp_solve dispatch to daqp_prox based on n_prox rather than eps_prox.
-    if(qp->H == NULL && qp->f!=NULL)
-        work->n_prox = work->n;
+
     return 1;
 }
 
@@ -191,7 +201,7 @@ void setup_daqp_hiqp(DAQPWorkspace* work, int* break_points, int nh){
     }
 }
 
-int setup_daqp_bnb(DAQPWorkspace* work, int nb, int ns){
+int setup_daqp_bnb(DAQPWorkspace* work, int* sense, int nb, int ns){
     int i, nadded;
     if(nb > work->n) return DAQP_EXIT_OVERDETERMINED_INITIAL;
     if((work->bnb == NULL) && (nb >0)){
@@ -201,7 +211,7 @@ int setup_daqp_bnb(DAQPWorkspace* work, int nb, int ns){
         // Detect which constraints are binary
         work->bnb->bin_ids = malloc(nb*sizeof(int));
         for(i = 0, nadded = 0; nadded < nb; i++){
-            if(work->qp->sense[i] & DAQP_BINARY)
+            if(sense[i] & DAQP_BINARY)
                 work->bnb->bin_ids[nadded++] = i;
         }
 
diff --git a/src/utils.c b/src/utils.c
index 4b6119bd..1cfa3908 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -30,7 +30,7 @@ int daqp_update_ldp(const int mask, DAQPWorkspace *work, DAQPProblem* qp){
     }
 
     // Check bounds early
-    if(mask&DAQP_UPDATE_Rinv||mask&DAQP_UPDATE_M||mask&DAQP_UPDATE_v||mask&DAQP_UPDATE_d){
+    if(mask&DAQP_UPDATE_M||mask&DAQP_UPDATE_v||mask&DAQP_UPDATE_d){
         error_flag = daqp_check_bounds(work,qp->bupper,qp->blower);
         if(error_flag<0) return error_flag;
         if(error_flag==1) do_activate = 1;
@@ -53,78 +53,8 @@ int daqp_update_ldp(const int mask, DAQPWorkspace *work, DAQPProblem* qp){
         daqp_update_v(qp->f,work,mask);
     }
 
-    // The unconstrained shortcut is invalid when:
-    //   - The problem is an LP (H == NULL, f != NULL): Rinv and RinvD are both NULL
-    //   - Equality constraints are present
-    const int has_hessian = (work->Rinv != NULL || work->RinvD != NULL || work->v == NULL);
-    const int check_unconstrained_base =
-        (mask&DAQP_UPDATE_Rinv||mask&DAQP_UPDATE_M||mask&DAQP_UPDATE_v||mask&DAQP_UPDATE_d) &&
-        work->bnb == NULL && work->nh <= 1 && work->avi == NULL;
-    int has_equalities = 0;
-    if(check_unconstrained_base && has_hessian){
-        for(i = 0; i < work->m; i++)
-            if(work->sense[i]&(DAQP_ACTIVE + DAQP_IMMUTABLE)){ has_equalities = 1; break; }
-    }
-
-    const int check_unconstrained = check_unconstrained_base && has_hessian && !has_equalities;
-    /** Check if unconstrained optimum is primal feasible.
-     *  Compute x_unc = Rinv * (-v) (using the raw, un-normalized Rinv) and
-     *  verify dupper_unnorm = bupper - A*x_unc >= 0 and
-     *         dlower_unnorm = blower - A*x_unc <= 0 for every constraint.
-     *  If so, x_unc is optimal and we can skip the expensive M = A*Rinv step.
-     *  Only applicable for standard QPs (no BnB, no hierarchy, no AVI, no prox,
-     *  no equality constraints, and with a positive-definite Hessian). **/
-    if(check_unconstrained){
-        int j, disp;
-        c_float sum;
-        c_float* swp_ptr;
-        int feasible = 1;
-        const int n = work->n;
-        const c_float primal_tol = work->settings->primal_tol;
-
-        // Compute x_unc = Rinv_unnorm * (-v) stored temporarily in work->x.
-        swp_ptr = work->x; work->u = work->xold; work->x = work->xold; work->xold = swp_ptr;
-        if(work->v != NULL){
-            if(work->Rinv != NULL){
-                // Upper-triangular back-substitution: u[i] = sum_{j>=i} Rinv[i,j]*(-v[j])
-                for(i = 0, disp = 0; i < n; i++){
-                    for(j = i, sum = 0; j < n; j++)
-                        sum += work->Rinv[disp++] * work->v[j];
-                    work->x[i] = -sum;
-                }
-            } else if(work->RinvD != NULL){
-                for(i = 0; i < n; i++) work->x[i] = -work->RinvD[i] * work->v[i];
-            } else {
-                for(i = 0; i < n; i++) work->x[i] = -work->v[i];
-            }
-        }
-        // If v == NULL (f == 0), x_unc = 0 and work->x is already 0.
-
-        // Check simple bounds: blower[i] <= x_unc[i] <= bupper[i]
-        for(i = 0; i < work->ms; i++){
-            work->dupper[i] = qp->bupper[i] - work->x[i];
-            work->dlower[i] = qp->blower[i] - work->x[i];
-            if(work->dupper[i] < -primal_tol || work->dlower[i] > primal_tol)
-                feasible = 0;
-        }
-
-        // Check general constraints: blower[i] <= A[i,:]*x_unc <= bupper[i]
-        for(i = work->ms, disp = 0; i < work->m; i++){
-            for(j = 0, sum = 0.0; j < n; j++)
-                sum += qp->A[disp++] * work->x[j];
-            work->dupper[i] = qp->bupper[i] - sum;
-            work->dlower[i] = qp->blower[i] - sum;
-            if(work->dupper[i] < -primal_tol || work->dlower[i] > primal_tol)
-                feasible = 0;
-        }
-        if(feasible){
-            reset_daqp_workspace(work);
-            work->sing_ind = DAQP_UNCONSTRAINED_OPTIMAL;
-            return 0;
-        }
-        // Switch back such that any warm starts a preserved 
-        swp_ptr = work->x; work->u = work->xold; work->x = work->xold; work->xold = swp_ptr;
-    }
+    int unconstrained_flag = daqp_check_unconstrained(work,mask);
+    if(unconstrained_flag == DAQP_UNCONSTRAINED_OPTIMAL) return 0;
 
     /** Update M **/
     if(mask&DAQP_UPDATE_Rinv||mask&DAQP_UPDATE_M){
@@ -140,7 +70,7 @@ int daqp_update_ldp(const int mask, DAQPWorkspace *work, DAQPProblem* qp){
 
     /** Update d **/
     if(mask&DAQP_UPDATE_Rinv||mask&DAQP_UPDATE_M||mask&DAQP_UPDATE_v||mask&DAQP_UPDATE_d){
-        if(check_unconstrained){ // Already computed d, just need to normalize
+        if(unconstrained_flag == 1){ // Already computed d, just need to normalize
             if(work->scaling != NULL){
                 for(i = 0; i < work->m; i++){
                     work->dupper[i]*=work->scaling[i];
@@ -481,6 +411,71 @@ int daqp_normalize_M(DAQPWorkspace* work){
     return 0;
 }
 
+// Returns 0 if it did not even compute the unconstrained solution 
+// Returns 1 if it computed the unconstrained solution, but it was not optimal
+// Returns DAQP_UNCONSTRAINED_OPTIMAL if the unconstrained solution is optimal
+int daqp_check_unconstrained(DAQPWorkspace* work, const int mask){
+    int i;
+    if ((mask&DAQP_UPDATE_unconstrained)==0) return 0;
+    if ((mask&(DAQP_UPDATE_Rinv+DAQP_UPDATE_M+DAQP_UPDATE_v+DAQP_UPDATE_d)) == 0) return 0; // Nothing to update
+    if (work->bnb != NULL || work->nh > 1 || work->avi != NULL || work->n_prox >0) return 0; // Not a QP 
+    for(i = 0; i < work->m; i++) if(work->sense[i]&(DAQP_ACTIVE + DAQP_IMMUTABLE)) return 0; // No equalities
+
+    // Check if unconstrained optimum is primal feasible.
+    // Compute x_unc = Rinv * (-v) (using the raw, un-normalized Rinv) and
+    // verify dupper_unnorm = bupper - A*x_unc >= 0 and
+    //        dlower_unnorm = blower - A*x_unc <= 0 for every constraint.
+    int j, disp;
+    c_float sum;
+    c_float* swp_ptr;
+    int feasible = 1;
+    const int n = work->n;
+    const c_float primal_tol = work->settings->primal_tol;
+
+    // Compute x_unc = Rinv_unnorm * (-v) stored temporarily in work->x.
+    swp_ptr = work->x; work->u = work->xold; work->x = work->xold; work->xold = swp_ptr;
+    if(work->v != NULL){
+        if(work->Rinv != NULL){
+            // Upper-triangular back-substitution: u[i] = sum_{j>=i} Rinv[i,j]*(-v[j])
+            for(i = 0, disp = 0; i < n; i++){
+                for(j = i, sum = 0; j < n; j++)
+                    sum += work->Rinv[disp++] * work->v[j];
+                work->x[i] = -sum;
+            }
+        } else if(work->RinvD != NULL){
+            for(i = 0; i < n; i++) work->x[i] = -work->RinvD[i] * work->v[i];
+        } else {
+            for(i = 0; i < n; i++) work->x[i] = -work->v[i];
+        }
+    }
+    //
+    // Check simple bounds: blower[i] <= x_unc[i] <= bupper[i]
+    for(i = 0; i < work->ms; i++){
+        work->dupper[i] = work->qp->bupper[i] - work->x[i];
+        work->dlower[i] = work->qp->blower[i] - work->x[i];
+        if(work->dupper[i] < -primal_tol || work->dlower[i] > primal_tol)
+            feasible = 0;
+    }
+
+    // Check general constraints: blower[i] <= A[i,:]*x_unc <= bupper[i]
+    for(i = work->ms, disp = 0; i < work->m; i++){
+        for(j = 0, sum = 0.0; j < n; j++)
+            sum += work->qp->A[disp++] * work->x[j];
+        work->dupper[i] = work->qp->bupper[i] - sum;
+        work->dlower[i] = work->qp->blower[i] - sum;
+        if(work->dupper[i] < -primal_tol || work->dlower[i] > primal_tol)
+            feasible = 0;
+    }
+    if(feasible){
+        reset_daqp_workspace(work);
+        work->sing_ind = DAQP_UNCONSTRAINED_OPTIMAL;
+        return DAQP_UNCONSTRAINED_OPTIMAL;
+    }
+    // Switch back such that any warm starts a preserved
+    swp_ptr = work->x; work->u = work->xold; work->x = work->xold; work->xold = swp_ptr;
+    return 1;
+}
+
 int daqp_update_avi(DAQPAVI* avi, DAQPProblem* p){
     const int n = p->n;
     // Setup matrices Hsym, Hs_rho, and H_rho, LU_H