aãdd lin_sys_solver field to info, update docs

docs/src/algorithm/scale.rst

@@ -17,8 +17,8 @@ algorithm then becomes:
   w^{k+1} &= w^k + u^{k+1} - \tilde u^{k+1} \\
   \end{align}
 
-which yields :math:`w^k \rightarrow u^\star + R^{-1} \mathcal{Q}(u^\star)` where :math:`0
-\in \mathcal{Q}(u^\star) + N_{\mathcal{C}_+}(u^\star)`.
+which yields :math:`w^k \rightarrow u^\star + R^{-1} \mathcal{Q}(u^\star)` where
+:math:`0 \in \mathcal{Q}(u^\star) + N_{\mathcal{C}_+}(u^\star)`.
 
 This changes the first two steps of the procedure. The :ref:`linear projection
 <linear_solver>` and the cone projection (explained next).
@@ -77,6 +77,40 @@ The quantity :math:`\rho_x` is determined by the :ref:`setting <settings>` value
 described in :ref:`Dynamic scale updating <updating_scale>`. Finally,  :math:`d`
 is determined by :code:`TAU_FACTOR` in the code defined in :code:`glbopts.h`.
 
+Root-plus function
+------------------
+
+Finally, the :code:`root_plus` function is modified to be the solution
+of the following quadratic equation:
+
+.. math::
+  \tau^2 (d + r^\top R_{-1} r) + \tau (r^\top R_{-1} \mu^k - 2 r^\top R_{-1} p^k - d \eta^k) + p^k R_{-1} (p^k - \mu^k) = 0,
+
+where :math:`R_{-1}` corresponds to the first :math:`n+m` entries of :math:`R`.
+Other than when computing :math:`\kappa` (which does not affect the algorithm)
+this is the *only* place where :math:`d` appears, so we have a lot of
+flexibility in how to choose it and it can even change from iteration to
+iteration. It is an open question on how best to select this parameter.
+
+Moreau decomposition
+--------------------
+The projection onto a cone under a diagonal scaling also satisfies a
+Moreau-style decomposition identity, as follows:
+
+.. math::
+   x + R^{-1} \Pi_\mathcal{C^*}^{R^{-1}} ( - R x ) = \Pi_\mathcal{C}^R ( x )
+
+where :math:`\Pi_\mathcal{C}^R` denotes the projection onto convex cone
+:math:`\mathcal{C}` under the :math:`R`-norm, which is defined as
+
+.. math::
+  \|x\|_R = \sqrt{x^\top R x}.
+
+This identity is useful when deriving cone projection routines, though most
+cones are either invariant to this or we enforce that :math:`R` is constant
+across them. Note that the two components of the decomposition are
+:math:`R`-orthogonal.
+
 Dual vector
 -----------
 
@@ -91,39 +125,14 @@ and we have
 .. math::
   \begin{align}
   v^{k+1} &= R( u^{k+1} + w^k - 2 \tilde u^{k+1} ) \\
-          &= R( \Pi_{\mathcal{C}_+} (2 \tilde u^{k+1} - w^k) + w^k - 2 \tilde u^{k+1}) \\
-          &= R( \Pi_{\mathcal{C}^*_+} (-2 \tilde u^{k+1} + w^k)) \\
-          &= R \tilde v^{k+1} \\
+          &= R( \Pi^R_{\mathcal{C}_+} (2 \tilde u^{k+1} - w^k) + w^k - 2 \tilde u^{k+1}) \\
+          &= R( R^{-1} \Pi^{R^{-1}}_{\mathcal{C}^*_+} (R(w^k -2 \tilde u^{k+1}))) \\
+          &= \Pi^{R^{-1}}_{\mathcal{C}^*_+} (R(w^k -2 \tilde u^{k+1})) \\
           &\in \mathcal{C}^*_+
   \end{align}
 
-by Moreau and the fact that :math:`R` is chosen to be constant
-within each sub-cone :math:`\mathcal{K}`. Finally note that :math:`v^k \perp
-u^k`, since
-
-.. math::
-  \begin{align}
-  (u^k)^\top v^k &= (u^k)^\top R \tilde v^{k+1}  \\
-   &= \sum_{\mathcal{K} \in \mathcal{C}_+} r_\mathcal{K} (u^k_\mathcal{K})^\top \tilde v^k_\mathcal{K} \\
-   &= 0
-  \end{align}
-
-since :math:`\tilde v^k \perp u^k` and :math:`R` is chosen to be constant
-within each sub-cone :math:`\mathcal{K}`.
-
-Root-plus function
-------------------
-
-Finally, the :code:`root_plus` function is modified to be the solution
-of the following quadratic equation:
-
-.. math::
-  \tau^2 (d + r^\top R r) + \tau (r^\top R \mu^k - 2 r^\top R p^k - d \eta^k) + p^k R (p^k - \mu^k) = 0.
-
-Other than when computing :math:`\kappa` (which does not affect the algorithm)
-this is the *only* place where :math:`d` appears, so we have a lot of
-flexibility in how to choose it and it can even change from iteration to
-iteration. It is an open question on how best to select this parameter.
+by Moreau, and finally note that :math:`v^k \perp
+u^k` from the fact that the Moreau decomposition is :math:`R`-orthogonal.
 
 .. _updating_scale:
 

docs/src/api/c.rst

@@ -9,7 +9,8 @@ C / C++
 Main solver API
 ---------------
 
-The main C API is imported from the header :code:`scs.h`.
+The C API is imported from the header :code:`scs.h`, available `here
+<https://github.com/cvxgrp/scs/blob/master/include/scs.h>`_.
 
 .. doxygenfunction:: scs
 

docs/src/api/info.rst

@@ -19,6 +19,9 @@ the following fields.
    * - :code:`status`
      - :code:`char *`
      - Status string (e.g., 'solved')
+   * - :code:`lin_sys_solver`
+     - :code:`char *`
+     - Linear system solver used
    * - :code:`status_val`
      - :code:`scs_int`
      - Status integer :ref:`exit flag <exit_flags>`

docs/src/examples/c.rst

@@ -11,15 +11,16 @@ C code to solve this is below.
    :language: c
 
 After following the CMake :ref:`install instructions <c_install>`, we can
-compile the code using:
+compile the code using (assuming the library was installed in
+:code:`/usr/local/`):
 
 .. code::
 
 	  gcc -I/usr/local/include/scs -L/usr/local/lib/ -lscsdir qp.c -o qp.out
 
 .. ./qp.out > qp.c.out
 
-Then running the code yields output
+Then running the binary yields output
 
 .. literalinclude:: qp.c.out
    :language: none

docs/src/examples/qp.c

@@ -14,9 +14,9 @@ int main(int argc, char **argv) {
   int A_p[3] = {0, 2, 4};
   double b[3] = {-1., 0.3, -0.5};
   double c[2] = {-1., -1.};
-  /* data shape */
-  int n = 2;
-  int m = 3;
+  /* data shapes */
+  int n = 2; /* number of variables */
+  int m = 3; /* number of constraints */
 
   /* Allocate SCS structs */
   ScsCone *k = (ScsCone *)calloc(1, sizeof(ScsCone));
@@ -25,7 +25,7 @@ int main(int argc, char **argv) {
   ScsSolution *sol = (ScsSolution *)calloc(1, sizeof(ScsSolution));
   ScsInfo *info = (ScsInfo *)calloc(1, sizeof(ScsInfo));
 
-  /* Fill in structs */
+  /* Fill in data struct */
   d->m = m;
   d->n = n;
   d->b = b;
@@ -36,7 +36,7 @@ int main(int argc, char **argv) {
   /* Cone */
   k->l = m;
 
-  /* Utility to set some default settings */
+  /* Utility to set default settings */
   scs_set_default_settings(stgs);
 
   /* Modify tolerances */
@@ -49,8 +49,9 @@ int main(int argc, char **argv) {
   /* Verify that SCS solved the problem */
   printf("SCS solved successfully: %i\n", exitflag == SCS_SOLVED);
 
-  /* Print iterations taken */
-  printf("SCS took %i iters\n", info->iter);
+  /* Print some info about the solve */
+  printf("SCS took %i iters, using the %s linear solver.\n", info->iter,
+         info->lin_sys_solver);
 
   /* Print solution x */
   printf("Optimal solution vector x*:\n");

docs/src/examples/qp.c.out

@@ -13,21 +13,21 @@ lin-sys:  sparse-direct
 ------------------------------------------------------------------
  iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
 ------------------------------------------------------------------
-     0| 1.78e+00  1.00e+00  1.42e+00 -1.04e+00  1.00e-01  2.14e-05 
-   175| 3.34e-11  4.17e-10  4.07e-10  1.24e+00  2.08e+01  1.25e-04 
+     0| 1.78e+00  1.00e+00  1.42e+00 -1.04e+00  1.00e-01  2.09e-05 
+   175| 4.30e-11  5.94e-10  5.70e-10  1.24e+00  2.08e+01  1.27e-04 
 ------------------------------------------------------------------
 status:  solved
-timings: total: 2.01e-04s = setup: 7.52e-05s + solve: 1.26e-04s
-	 lin-sys: 2.56e-05s, cones: 1.43e-05s, accel: 2.20e-05s
+timings: total: 2.11e-04s = setup: 8.29e-05s + solve: 1.28e-04s
+	 lin-sys: 2.44e-05s, cones: 1.23e-05s, accel: 2.31e-05s
 ------------------------------------------------------------------
 objective = 1.235000
 ------------------------------------------------------------------
 SCS solved successfully: 1
-SCS took 175 iters
+SCS took 175 iters, using the sparse-direct linear solver.
 Optimal solution vector x*:
 x[0] = 0.300000
 x[1] = -0.700000
 Optimal dual vector y*:
 y[0] = 2.700000
 y[1] = 2.100000
-y[2] = 0.000000
+y[2] = -0.000000

docs/src/install/c.rst

@@ -54,9 +54,9 @@ project
   make
 
 The CMake build-system exports two CMake targets called :code:`scs::scsdir` and
-:code:`scs::scsindir` as well as a header file `scs.h` that defines the API. The
-libraries can be imported using the find_package CMake command and used by
-calling target_link_libraries as in the following example:
+:code:`scs::scsindir` as well as a header file :code:`scs.h` that defines the
+API. The libraries can be imported using the find_package CMake command and used
+by calling target_link_libraries as in the following example:
 
 .. code:: bash
 
@@ -73,7 +73,7 @@ calling target_link_libraries as in the following example:
 
 Makefile
 ^^^^^^^^
-Alternatively you can use the Makefile and manage the libraries yourself 
+Alternatively you can use the Makefile and manage the libraries yourself
 
 .. code:: bash
 
@@ -92,7 +92,7 @@ To compile and run the tests execute
 If make completes successfully, it will produce two static library files,
 :code:`libscsdir.a`, :code:`libscsindir.a`, and two dynamic library files
 :code:`libscsdir.ext`, :code:`libscsindir.ext` (where :code:`.ext` extension is
-platform dependent) in the :code:`out` folder.  
+platform dependent) in the :code:`out` folder.
 
 If you have a GPU and have CUDA installed, you can also execute make gpu to
 compile SCS to run on the GPU which will create additional libraries and demo
@@ -108,5 +108,5 @@ binaries in the out folder corresponding to the GPU version.  Note that the GPU
 To use the libraries in your own source code, compile your code with the linker
 option :code:`-L(PATH_TO_SCS_LIBS)` and :code:`-lscsdir` or :code:`-lscsindir`
 (as needed). The API and required data structures are defined in the file
-:code:`include/scs.h`. The four main API functions are:
+:code:`include/scs.h`.
 

include/cones.h

@@ -22,8 +22,8 @@ struct SCS_CONE_WORK {
   scs_int *cone_boundaries;
   scs_int cone_boundaries_len;
   scs_int scaled_cones; /* boolean, whether the cones have been scaled */
-  scs_float *s; /* used for Moreau decomposition in projection */
-  scs_int m;    /* total length of cone */
+  scs_float *s;         /* used for Moreau decomposition in projection */
+  scs_int m;            /* total length of cone */
   /* box cone quantities */
   scs_float *bl, *bu, box_t_warm_start;
 #ifdef USE_LAPACK

include/linsys.h

@@ -17,10 +17,9 @@ extern "C" {
 /**
  * Initialize `ScsLinSysWork` structure and perform any necessary preprocessing.
  *
- *  @param  A          A data matrix.
- *  @param  P          P data matrix.
- *  @param  rho_y_vec  `rho_y > 0` diagonal entries.
- *  @param  rho_x      `rho_x > 0` float.
+ *  @param  A          `A` data matrix, `m x n`.
+ *  @param  P          `P` data matrix, `n x n`.
+ *  @param  diag_r     `R > 0` diagonal entries of length `m + n`.
  *  @return            Linear system solver workspace.
  *
  */
@@ -35,15 +34,15 @@ ScsLinSysWork *SCS(init_lin_sys_work)(const ScsMatrix *A, const ScsMatrix *P,
 void SCS(free_lin_sys_work)(ScsLinSysWork *w);
 
 /**
- * Solves the linear system required by SCS at each iteration:
+ * Solves the linear system as required by SCS at each iteration:
  * \f[
  *    \begin{bmatrix}
- *    (\rho_x I + P) & A^\top \\
- *     A   &  -\mathrm{diag}(\rho_y) \\
+ *    (R_x + P) & A^\top \\
+ *     A   &  -R_y \\
  *    \end{bmatrix} x = b
  *  \f]
  *
- *  for `x`. Overwrites `b` with result.
+ *  for `x`, where `diag(R_x, R_y) = R`. Overwrites `b` with result.
  *
  *  @param  w    Linear system private workspace.
  *  @param  b    Right hand side, contains solution at the end.
@@ -59,11 +58,11 @@ scs_int SCS(solve_lin_sys)(ScsLinSysWork *w, scs_float *b, const scs_float *s,
  *  direct method for solving the linear system might need to update the
  *  factorization of the matrix.
  *
- *  @param  w          Linear system private workspace.
- *  @param  diag_r     `diag_r` diagonal entries of R (assumed to be diagonal).
+ *  @param  w             Linear system private workspace.
+ *  @param  new_diag_r    Updated `diag_r`, diagonal entries of R.
  *
  */
-void SCS(update_lin_sys_diag_r)(ScsLinSysWork *w, scs_float *diag_r);
+void SCS(update_lin_sys_diag_r)(ScsLinSysWork *w, scs_float *new_diag_r);
 
 /**
  * Name of the linear solver.

include/scs.h

@@ -153,6 +153,8 @@ typedef struct {
   scs_int iter;
   /** Status string, e.g. 'solved'. */
   char status[128];
+  /** Linear system solver used. */
+  char lin_sys_solver[128];
   /** Status as scs_int, defined in glbopts.h. */
   scs_int status_val;
   /** Number of updates to scale. */

linsys/gpu/indirect/private.c

@@ -374,8 +374,7 @@ ScsLinSysWork *SCS(init_lin_sys_work)(const ScsMatrix *A, const ScsMatrix *P,
 }
 
 /* solves (R_x + P + A' R_y^{-1} A)x = b, s warm start, solution stored in
- * b */
-/* on GPU */
+ * b, on GPU */
 static scs_int pcg(ScsLinSysWork *pr, const scs_float *s, scs_float *bg,
                    scs_int max_its, scs_float tol) {
   scs_int i, n = pr->n;

src/cones.c

@@ -596,11 +596,8 @@ static void normalize_box_cone(ScsConeWork *c, scs_float *D, scs_int bsize) {
   }
 }
 
-/* project onto { (t, s) | t * l <= s <= t * u, t >= 0 }, Newton's method on t
-   tx = [t; s], total length = bsize
-
-  under Euclidean metric r_box
-
+/* Project onto { (t, s) | t * l <= s <= t * u, t >= 0 }, Newton's method on t
+   tx = [t; s], total length = bsize, under Euclidean metric 1/r_box.
 */
 static scs_float proj_box_cone(scs_float *tx, const scs_float *bl,
                                const scs_float *bu, scs_int bsize,
@@ -735,6 +732,9 @@ static void proj_power_cone(scs_float *v, scs_float a) {
 }
 
 /* project onto the primal K cone in the paper */
+/* the r_y vector determines the INVERSE metric, ie, project under the
+ * diag(r_y)^-1 norm.
+ */
 static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c,
                          scs_int normalize, scs_float *r_y) {
   scs_int i, status;
@@ -894,13 +894,16 @@ void scale_box_cone(const ScsCone *k, ScsConeWork *c, ScsScaling *scal) {
   }
 }
 
-/* outward facing cone projection routine
-   performs projection in-place
-   if normalize > 0 then will use normalized (equilibrated) cones if applicable.
+/* Outward facing cone projection routine, performs projection in-place.
+   If normalize > 0 then will use normalized (equilibrated) cones if applicable.
+
+   Moreau decomposition for R-norm projections:
+
+    `x + R^{-1} \Pi_{C^*}^{R^{-1}} ( - R x ) = \Pi_C^R ( x )`
 
-   x + R^{-1} \Pi_{C^*}^{R^{-1}} ( - R x ) = \Pi_C^R ( x )
+   where \Pi^R_C is the projection onto C under the R-norm:
 
-   where \Pi^R_C is the projection onto C under the R-Euclidean norm.
+    `||x||_R = \sqrt{x ' R x}`.
 
 */
 scs_int SCS(proj_dual_cone)(scs_float *x, ScsConeWork *c, ScsScaling *scal,

src/scs.c

@@ -975,6 +975,7 @@ scs_int scs_solve(ScsWork *w, ScsSolution *sol, ScsInfo *info) {
   /* initialize ctrl-c support */
   scs_start_interrupt_listener();
   SCS(tic)(&solve_timer);
+  strcpy(info->lin_sys_solver, SCS(get_lin_sys_method)());
   info->status_val = SCS_UNFINISHED; /* not yet converged */
   update_work(d, w, sol);