Merge 21a1387 into a46b5f6

osqp · Mar 22, 2019 · c5fce9b · c5fce9b
2 parents a46b5f6 + 21a1387
commit c5fce9b
Show file tree

Hide file tree

Showing 35 changed files with 369 additions and 358 deletions.
diff --git a/docs/conf.py b/docs/conf.py
@@ -51,7 +51,7 @@
 
 # General information about the project.
 project = 'OSQP'
-copyright = '2018, Bartolomeo Stellato, Goran Banjac'
+copyright = '2019, Bartolomeo Stellato, Goran Banjac'
 author = 'Bartolomeo Stellato, Goran Banjac'
 
 # The version info for the project you're documenting, acts as replacement for

diff --git a/docs/contributing/index.rst b/docs/contributing/index.rst
@@ -75,7 +75,7 @@ The linear system solver object is defined in :code:`mysolver.h` as follows
         };
 
         // Initialize mysolver solver
-        mysolver_solver *init_linsys_solver_mysolver(const csc * P, const csc * A, c_float sigma, c_float * rho_vec, c_int polish);
+        mysolver_solver *init_linsys_solver_mysolver(const csc * P, const csc * A, c_float sigma, c_float * rho_vec, c_int polish, c_int * exitflag);
 
         // Solve linear system and store result in b
         c_int solve_linsys_mysolver(mysolver_solver * s, c_float * b, const OSQPSettings * settings);

diff --git a/docs/examples/demo.rst b/docs/examples/demo.rst
@@ -101,17 +101,17 @@ C
 
     int main(int argc, char **argv) {
         // Load problem data
-        c_float P_x[4] = {4.00, 1.00, 1.00, 2.00, };
-        c_int P_nnz = 4;
-        c_int P_i[4] = {0, 1, 0, 1, };
-        c_int P_p[3] = {0, 2, 4, };
-        c_float q[2] = {1.00, 1.00, };
-        c_float A_x[4] = {1.00, 1.00, 1.00, 1.00, };
+        c_float P_x[3] = {4.0, 1.0, 2.0, };
+        c_int P_nnz = 3;
+        c_int P_i[3] = {0, 0, 1, };
+        c_int P_p[3] = {0, 1, 3, };
+        c_float q[2] = {1.0, 1.0, };
+        c_float A_x[4] = {1.0, 1.0, 1.0, 1.0, };
         c_int A_nnz = 4;
         c_int A_i[4] = {0, 1, 0, 2, };
         c_int A_p[3] = {0, 2, 4, };
-        c_float l[3] = {1.00, 0.00, 0.00, };
-        c_float u[3] = {1.00, 0.70, 0.70, };
+        c_float l[3] = {1.0, 0.0, 0.0, };
+        c_float u[3] = {1.0, 0.7, 0.7, };
         c_int n = 2;
         c_int m = 3;
 
@@ -132,13 +132,12 @@ C
         data->l = l;
         data->u = u;
 
-
         // Define Solver settings as default
         osqp_set_default_settings(settings);
         settings->alpha = 1.0; // Change alpha parameter
 
         // Setup workspace
-        work = osqp_setup(data, settings);
+        osqp_setup(&work, data, settings);
 
         // Solve Problem
         osqp_solve(work);

diff --git a/docs/interfaces/CC++.rst b/docs/interfaces/CC++.rst
@@ -80,7 +80,7 @@ Data
 .. doxygenstruct:: OSQPData
    :members:
 
-The matrices are defined in `Compressed Sparse Column (CSC) format <https://people.sc.fsu.edu/~jburkardt/data/cc/cc.html>`_.
+The matrices are defined in `Compressed Sparse Column (CSC) format <https://people.sc.fsu.edu/~jburkardt/data/cc/cc.html>`_ using zero-based indexing.
 
 .. doxygenstruct:: csc
    :members:

diff --git a/examples/osqp_demo.c b/examples/osqp_demo.c
@@ -4,27 +4,20 @@
 
 int main(int argc, char **argv) {
   // Load problem data
-  c_float P_x[4] =
-  { 4.00000000000000000000, 1.00000000000000000000, 1.00000000000000000000,
-    2.00000000000000000000, };
-  c_int   P_nnz  = 4;
-  c_int   P_i[4] = { 0, 1, 0, 1, };
-  c_int   P_p[3] = { 0, 2, 4, };
-  c_float q[2]   = { 1.00000000000000000000, 1.00000000000000000000, };
-  c_float A_x[4] =
-  { 1.00000000000000000000, 1.00000000000000000000, 1.00000000000000000000,
-    1.00000000000000000000, };
-  c_int   A_nnz  = 4;
-  c_int   A_i[4] = { 0, 1, 0, 2, };
-  c_int   A_p[3] = { 0, 2, 4, };
-  c_float l[3]   =
-  { 1.00000000000000000000, 0.00000000000000000000, 0.00000000000000000000, };
-  c_float u[3] =
-  { 1.00000000000000000000, 0.69999999999999995559, 0.69999999999999995559, };
+  c_float P_x[3] = { 4.0, 1.0, 2.0, };
+  c_int P_nnz  = 3;
+  c_int P_i[3] = { 0, 0, 1, };
+  c_int P_p[3] = { 0, 1, 3, };
+  c_float q[2]   = { 1.0, 1.0, };
+  c_float A_x[4] = { 1.0, 1.0, 1.0, 1.0, };
+  c_int A_nnz  = 4;
+  c_int A_i[4] = { 0, 1, 0, 2, };
+  c_int A_p[3] = { 0, 2, 4, };
+  c_float l[3]   = { 1.0, 0.0, 0.0, };
+  c_float u[3] = { 1.0, 0.7, 0.7, };
   c_int n = 2;
   c_int m = 3;
 
-
   // Problem settings
   OSQPSettings *settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings));
 
@@ -33,7 +26,7 @@ int main(int argc, char **argv) {
   OSQPData *data;      // OSQPData
 
   // Populate data
-  data    = (OSQPData *)c_malloc(sizeof(OSQPData));
+  data = (OSQPData *)c_malloc(sizeof(OSQPData));
   data->n = n;
   data->m = m;
   data->P = csc_matrix(data->n, data->n, P_nnz, P_x, P_i, P_p);
@@ -42,12 +35,11 @@ int main(int argc, char **argv) {
   data->l = l;
   data->u = u;
 
-
-  // Define Solver settings as default
+  // Define solver settings as default
   osqp_set_default_settings(settings);
 
   // Setup workspace
-  work = osqp_setup(data, settings);
+  osqp_setup(&work, data, settings);
 
   // Solve Problem
   osqp_solve(work);
@@ -59,6 +51,5 @@ int main(int argc, char **argv) {
   c_free(data);
   c_free(settings);
 
-
   return 0;
 }
diff --git a/include/constants.h b/include/constants.h
@@ -12,6 +12,16 @@ extern "C" {
 # define OSQP_VERSION ("0.5.0") /* string literals automatically null-terminated
                                    */
 
+/******************
+* Setup Status  *
+******************/
+# define OSQP_DATA_VALIDATION_ERROR (1)
+# define OSQP_SETTINGS_VALIDATION_ERROR (2)
+# define OSQP_MEMORY_ALLOCATION_ERROR (3)
+# define OSQP_LOAD_LINSYS_SOLVER_ERROR (4)
+# define OSQP_INIT_LINSYS_SOLVER_ERROR (5)
+# define OSQP_INIT_LINSYS_SOLVER_NONCVX (6)
+
 
 /******************
 * Solver Status  *

diff --git a/include/lin_sys.h b/include/lin_sys.h
@@ -32,22 +32,24 @@ c_int unload_linsys_solver(enum linsys_solver_type linsys_solver);
 
 /**
  * Initialize linear system solver structure
- * @param P             Cost function matrix
+ * @param   P             Cost function matrix
  * @param	A             Constraint matrix
  * @param	sigma         Algorithm parameter
  * @param	rho_vec       Algorithm parameter
+ * @param	linsys_solver Linear system solver
  * @param	polish        0/1 depending whether we are allocating for
  *polishing or not
- * @param	linsys_solver Linear system solver
+ * @param	exitflag      Exitflag for error (0 if no errors)
  * @return Pointer to linear system solver structure on success, OSQP_NULL on
  *failure.
  */
 LinSysSolver* init_linsys_solver(const csc              *P,
                                  const csc              *A,
                                  c_float                 sigma,
-                                 c_float                *rho_vec,
+                                 const c_float          *rho_vec,
                                  enum linsys_solver_type linsys_solver,
-                                 c_int                   polish);
+                                 c_int                   polish,
+                                 c_int                  *exitflag);
 
 
 # endif // EMBEDDED

diff --git a/include/osqp.h b/include/osqp.h
@@ -50,12 +50,12 @@ void osqp_set_default_settings(OSQPSettings *settings);
  * NB: This is the only function that allocates dynamic memory and is not used
  *during code generation
  *
+ * @param  work         Solver workspace
  * @param  data         Problem data
  * @param  settings     Solver settings
- * @return              Solver environment
+ * @return              Exitflag for errors (0 if no errors)
  */
-OSQPWorkspace* osqp_setup(const OSQPData *data,
-                          OSQPSettings   *settings);
+c_int osqp_setup(OSQPWorkspace** work, const OSQPData* data, const OSQPSettings* settings);
 
 # endif // #ifndef EMBEDDED
 

diff --git a/include/types.h b/include/types.h
@@ -126,9 +126,8 @@ typedef struct {
 typedef struct {
   c_int    n; ///< number of variables n
   c_int    m; ///< number of constraints m
-  csc     *P; ///< quadratic part of the cost P in csc format (size n x n). It
-              ///  can be either the full P or only the upper triangular part. The
-              ///  workspace stores only the upper triangular part
+  csc     *P; ///< the upper triangular part of the quadratic cost matrix
+              ///  P in csc format (size n x n).
   csc     *A; ///< linear constraints matrix A in csc format (size m x n)
   c_float *q; ///< dense array for linear part of cost function (size n)
   c_float *l; ///< dense array for lower bound (size m)
@@ -314,25 +313,25 @@ typedef struct {
  *      on the choice
  */
 struct linsys_solver {
-  enum linsys_solver_type type; ///< Linear system solver type (see type.h)
+  enum linsys_solver_type type; ///< Linear system solver type
   // Functions
   c_int (*solve)(LinSysSolver       *self,
                  c_float            *b,
                  const OSQPSettings *settings); ///< Solve linear system
 
-    # ifndef EMBEDDED
+# ifndef EMBEDDED
   void (*free)(LinSysSolver *self);             ///< Free linear system solver
                                                 // (only in desktop version)
-    # endif // ifndef EMBEDDED
+# endif // ifndef EMBEDDED
 
-    # if EMBEDDED != 1
+# if EMBEDDED != 1
   c_int (*update_matrices)(LinSysSolver *self, const csc *P, const csc *A,
                            const OSQPSettings *settings); ///< Update matrices P
                                                           // and A in the solver
   c_int (*update_rho_vec)(LinSysSolver  *s,
                           const c_float *rho_vec,
                           const c_int    m);              ///< Update rho
-    # endif // if EMBEDDED != 1
+# endif // if EMBEDDED != 1
 
 # ifndef EMBEDDED
   c_int nthreads; ///< Number of threads active

diff --git a/include/util.h b/include/util.h
@@ -30,7 +30,7 @@ const char* osqp_version(void);
  * @param  settings Settings to be copied
  * @return          New settings structure
  */
-OSQPSettings* copy_settings(OSQPSettings *settings);
+OSQPSettings* copy_settings(const OSQPSettings *settings);
 
 # endif // #ifndef EMBEDDED
 
@@ -94,9 +94,9 @@ void print_footer(OSQPInfo *info,
   // Some R packages clash with elements
   // of the windows.h header, so use a
   // slimmer version for conflict avoidance
-  # ifdef R_LANG
-    #define NOGDI
-  # endif
+# ifdef R_LANG
+#define NOGDI
+# endif
 
 #   include <windows.h>
 

diff --git a/lin_sys/direct/pardiso/pardiso_interface.c b/lin_sys/direct/pardiso/pardiso_interface.c
@@ -48,7 +48,7 @@ void free_linsys_solver_pardiso(pardiso_solver *s) {
 
 
 // Initialize factorization structure
-pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float sigma, c_float * rho_vec, c_int polish){
+pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float sigma, const c_float * rho_vec, c_int polish, c_int * exitflag){
     c_int i;                     // loop counter
     c_int nnzKKT;                // Number of nonzeros in KKT
     // Define Variables
@@ -93,9 +93,10 @@ pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float
 
     // Check if matrix has been created
     if (!(s->KKT)) {
-    #ifdef PRINTING
+#ifdef PRINTING
 	    c_eprint("Error in forming KKT matrix");
-    #endif
+#endif
+        *exitflag = OSQP_INIT_LINSYS_SOLVER_ERROR;
 	    return OSQP_NULL;
     } else {
 	    // Adjust indexing for Pardiso
@@ -113,11 +114,11 @@ pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float
     }
 
     // Set MKL interface layer (Long integers if activated)
-    #ifdef DLONG
+#ifdef DLONG
     mkl_set_interface_layer(MKL_INTERFACE_ILP64);
-    #else
+#else
     mkl_set_interface_layer(MKL_INTERFACE_LP64);
-    #endif
+#endif
 
     // Set Pardiso variables
     s->mtype = -2;        // Real symmetric indefinite matrix
@@ -148,10 +149,11 @@ pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float
              &(s->n), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
              s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
     if ( s->error != 0 ){
-        #ifdef PRINTING
-            c_eprint("Error during symbolic factorization: %d", (int)s->error);
-        #endif
+#ifdef PRINTING
+        c_eprint("Error during symbolic factorization: %d", (int)s->error);
+#endif
         free_linsys_solver_pardiso(s);
+        *exitflag = OSQP_INIT_LINSYS_SOLVER_ERROR;
         return OSQP_NULL;
     }
 
@@ -161,10 +163,11 @@ pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float
              &(s->n), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
              s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
     if ( s->error != 0 ){
-        #ifdef PRINTING
-            c_eprint("Error during numerical factorization: %d", (int)s->error);
-        #endif
+#ifdef PRINTING
+        c_eprint("Error during numerical factorization: %d", (int)s->error);
+#endif
         free_linsys_solver_pardiso(s);
+        *exitflag = OSQP_INIT_LINSYS_SOLVER_ERROR;
         return OSQP_NULL;
     }
 
@@ -177,6 +180,8 @@ pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float
     // Assign type
     s->type = MKL_PARDISO_SOLVER;
 
+    // No error
+    *exitflag = 0;
     return s;
 }
 
@@ -188,9 +193,9 @@ c_int solve_linsys_pardiso(pardiso_solver * s, c_float * b, const OSQPSettings *
              &(s->n), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
              s->iparm, &(s->msglvl), b, s->bp, &(s->error));
     if ( s->error != 0 ){
-        #ifdef PRINTING
+#ifdef PRINTING
         c_eprint("Error during solution: %d", (int)s->error);
-        #endif
+#endif
         return 1;
     }
 

diff --git a/lin_sys/direct/pardiso/pardiso_interface.h b/lin_sys/direct/pardiso/pardiso_interface.h
@@ -39,7 +39,7 @@ struct pardiso {
      */
     // Attributes
     csc *KKT;         ///< KKT matrix (in CSR format!)
-    c_int *KKT_i;     ///< KKT column ondeces in 1-indexing for Pardiso
+    c_int *KKT_i;     ///< KKT column indices in 1-indexing for Pardiso
     c_int *KKT_p;     ///< KKT row pointers in 1-indexing for Pardiso
     c_float *bp;      ///< workspace memory for solves (rhs)
 
@@ -69,14 +69,15 @@ struct pardiso {
 /**
  * Initialize Pardiso Solver
  *
- * @param  P      Cost function matrix (upper triangular form)
- * @param  A      Constraints matrix
- * @param  sigma   Algorithm parameter. If polish, then sigma = delta.
- * @param  rho_vec Algorithm parameter. If polish, then rho_vec = OSQP_NULL.
- * @param  polish Flag whether we are initializing for polish or not
- * @return        Initialized private structure
+ * @param  P         Cost function matrix (upper triangular form)
+ * @param  A         Constraints matrix
+ * @param  sigma     Algorithm parameter. If polish, then sigma = delta.
+ * @param  rho_vec   Algorithm parameter. If polish, then rho_vec = OSQP_NULL.
+ * @param  polish    Flag whether we are initializing for polish or not
+ * @param  exitflag  Exitflag for error (0 if no errors)
+ * @return           Initialized private structure
  */
-pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float sigma, c_float * rho_vec, c_int polish);
+pardiso_solver *init_linsys_solver_pardiso(const csc * P, const csc * A, c_float sigma, const c_float * rho_vec, c_int polish, c_int * exitflag);
 
 
 /**