Modularize unconstrained shortcut, and disable it when setting up workspace outside of single solve

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

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);

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

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);

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)
 

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
+

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)

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;

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:

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.

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;
         }