Highs mods

cupdlp/README.md

@@ -0,0 +1,66 @@
+# cuPDLP-C observations
+
+This directory contains files from [cuPDLP-C v0.3.0](https://github.com/COPT-Public/cuPDLP-C/tree/v0.3.0). Below are some issues experienced when integrating them into HiGHS.
+
+## Preprocessing issue
+
+The following line is not recognised by g++, 
+
+> #if !(CUPDLP_CPU)
+
+so I've had to replace all ocurrences by
+
+> #ifndef CUPDLP_CPU
+
+This yields a compiler warning about "extra tokens at end of #ifndef
+directive" in the case of the following, but it's not a problem for
+now, as CUPDLP_CPU is set
+
+> #ifndef CUPDLP_CPU & USE_KERNELS
+
+## cmake issues
+
+CUPDLP_CPU and CUPDLP_DEBUG should both set when building. However, they are not recognised so are forced by the following lines in cupdlp_defs.h
+
+#define CUPDLP_CPU
+
+#define CUPDLP_DEBUG (1)
+
+## Use of macro definitions within C
+
+Although the macro definitions in [glbopts.h](https://github.com/ERGO-Code/HiGHS/blob/add-pdlp/src/pdlp/cupdlp/glbopts.h) are fine when used in C under Linux, they cause the following compiler errors on Windows.
+
+> error C2146: syntax error: missing ';' before identifier 'calloc' (or 'malloc')
+
+In HiGHS, all the macros using `typeof` have been replaced by multiple type-specific macros
+
+## Problem with sys/time.h
+
+The HiGHS branch add-pdlp compiles and runs fine on @jajhall's Linux machine, but CI tests on GitHub fail utterly due to `sys/time.h` not being found. Until this is fixed, or HiGHS passes its own timer for use within `cuPDLP-c`, timing within `cuPDLP-c` can be disabled using the compiler directive `CUPDLP_TIMER`. By default this is defined, so the `cuPDLP-c` is retained.
+
+## Termination of cuPDLP-C
+
+cuPDLP-C terminates when either the current or averaged iterates satisfy primal/dual feasibility (and a duality gap criterion), using a 2-norm measure relative to the size of the RHS/costs. HiGHS assesses primal/dual feasibility using a infinity-norm absolute measure. Thus the cuPDLP-C result frequently fails to satisfy HiGHS primal/dual feasibility. To get around this, `iInfNormAbsLocalTermination` has been introduced into cuPDLP-C. 
+
+By default, `iInfNormAbsLocalTermination` is false, so that the original cuPDLP-C termination criteria are used.
+
+When `iInfNormAbsLocalTermination` is true, cuPDLP-C terminates only when primal/dual feasibility is satisfied for the infinity-norm absolute measure of the current iterate, so that HiGHS primal/dual feasibility is satisfied. 
+
+## Controlling the `cuPDLP-c` logging
+
+As a research code, `cuPDLP-c` naturally produces a lot of logging output. HiGHS must be able to run with less logging output, or completely silently. This is achieved using the `nLogLevel` parameter in `cuPDLP-c`. 
+
+By default, `nLogLevel` is 2, so all the original `cuPDLP-c` logging is produced.
+
+* If `nLogLevel` is 1, then the `cuPDLP-c` logging is less verbose 
+* If `nLogLevel` is 0, then there is no `cuPDLP-c` logging
+
+A related issue is the use of `fp` and `fp_sol`. HiGHS won't be using these, so sets them to null pointers. `cuPDLP-c` already doesn't print the solution if `fp_sol` is a null pointer, so the call to `writeJson(fp, pdhg);` is now conditional on `if (fp)`. 
+
+## Returning the iteration count
+
+The `cuPDLP-c` iteration count is held in `pdhg->timers->nIter`, but `pdhg` is destroyed in `LP_SolvePDHG`, so `cupdlp_int* num_iter` has been added to the parameter list of this method.
+
+
+
+

cupdlp/cupdlp.h

@@ -5,7 +5,7 @@
 #include "cupdlp_cs.h"
 #include "cupdlp_defs.h"
 #include "cupdlp_linalg.h"
-#include "cupdlp_mps.h"
+//#include "cupdlp_mps.h" Is this still needed?
 #include "cupdlp_proj.h"
 #include "cupdlp_restart.h"
 #include "cupdlp_scaling_cuda.h"

cupdlp/cupdlp_defs.h

@@ -1,6 +1,8 @@
 #ifndef CUPDLP_H_GUARD
 #define CUPDLP_H_GUARD
 
+#define CUPDLP_TIMER (1)
+
 #if !(CUPDLP_CPU)
 #include "cuda/cupdlp_cuda_kernels.cuh"
 #include "cuda/cupdlp_cudalinalg.cuh"
@@ -96,19 +98,23 @@ typedef enum {
   IF_RUIZ_SCALING,
   IF_L2_SCALING,
   IF_PC_SCALING,
+  N_LOG_LEVEL,
   N_LOG_INTERVAL,
   IF_PRESOLVE,
+  I_INF_NORM_ABS_LOCAL_TERMINATION,
+  N_INT_USER_PARAM
 } CUPDLP_INT_USER_PARAM_INDEX;
-#define N_INT_USER_PARAM 10
+  //#define N_INT_USER_PARAM 12
 typedef enum {
   D_SCALING_LIMIT = 0,
   D_PRIMAL_TOL,
   D_DUAL_TOL,
   D_GAP_TOL,
   D_FEAS_TOL,
   D_TIME_LIM,
+  N_FLOAT_USER_PARAM
 } CUPDLP_FLOAT_USER_PARAM_INDEX;
-#define N_FLOAT_USER_PARAM 6
+  //#define N_FLOAT_USER_PARAM 6
 
 // used in sparse matrix-dense vector multiplication
 struct CUPDLP_CUDA_DENSE_VEC {
@@ -174,10 +180,14 @@ struct CUPDLP_SETTINGS {
   cupdlp_float dPrimalTol;
   cupdlp_float dDualTol;
   cupdlp_float dGapTol;
+  cupdlp_int iInfNormAbsLocalTermination;
 
   // max iter and time
   cupdlp_int nIterLim;
   cupdlp_float dTimeLim;
+
+  // Logging
+  cupdlp_int nLogLevel;
   cupdlp_int nLogInterval;
 
   // restart
@@ -408,4 +418,4 @@ struct CUPDLP_WORK {
 #ifdef __cplusplus
 }
 #endif
-#endif
+#endif

cupdlp/cupdlp_linalg.c

@@ -141,6 +141,27 @@ double nrminf(cupdlp_int n, const double *x, cupdlp_int incx) {
 #endif
 }
 
+cupdlp_int nrminfindex(cupdlp_int n, const double *x, cupdlp_int incx) {
+#ifdef USE_MY_BLAS
+  assert(incx == 1);
+
+  double nrm = 0.0;
+  cupdlp_int index = 0;
+
+  for (int i = 0; i < n; ++i) {
+    double tmp = fabs(x[i]);
+    if (tmp > nrm) {
+      nrm = tmp;
+      index = i;
+    }
+  }
+
+  return index;
+#else
+  return dnrminfindex(n, x, incx);
+#endif
+}
+
 double twoNorm(double *x, cupdlp_int n) { return nrm2(n, x, 1); }
 
 double twoNormSquared(double *x, cupdlp_int n) { return pow(twoNorm(x, n), 2); }
@@ -290,7 +311,7 @@ cupdlp_float diffDotDiff(cupdlp_float *x1, cupdlp_float *x2, cupdlp_float *y1,
 
 /*------------------------ new added --------------------*/
 
-double dot(cupdlp_int n, cupdlp_float *x, cupdlp_int incx, cupdlp_float *y,
+double dot(cupdlp_int n, const cupdlp_float *x, cupdlp_int incx, const cupdlp_float *y,
            cupdlp_int incy) {
 #ifdef USE_MY_BLAS
   assert(incx == 1 && incy == 1);
@@ -307,11 +328,11 @@ double dot(cupdlp_int n, cupdlp_float *x, cupdlp_int incx, cupdlp_float *y,
 #endif
 }
 
-double Dotprod(cupdlp_float *x, cupdlp_float *y, cupdlp_int n) {
+double Dotprod(const cupdlp_float *x, const cupdlp_float *y, cupdlp_int n) {
   return dot(n, x, 1, y, 1);
 }
 
-double Dotprod_Neumaier(cupdlp_float *x, cupdlp_float *y, cupdlp_int n) {
+double Dotprod_Neumaier(const cupdlp_float *x, const cupdlp_float *y, cupdlp_int n) {
   return dot(n, x, 1, y, 1);
 }
 
@@ -571,6 +592,20 @@ cupdlp_int cupdlp_twoNorm(CUPDLPwork *w, const cupdlp_int n,
   return 0;
 }
 
+cupdlp_int cupdlp_infNormIndex(CUPDLPwork *w, const cupdlp_int n,
+                               const cupdlp_float *x, cupdlp_int *res) {
+#ifndef CUPDLP_CPU
+#ifndef SFLOAT
+  CHECK_CUBLAS(cublasIdamax(w->cublashandle, n, x, 1, res));
+#else
+  CHECK_CUBLAS(cublasIsamax(w->cublashandle, n, x, 1, res));
+#endif
+#else
+  *res = nrminfindex(n, x, 1);
+#endif
+  return 0;
+}
+
 cupdlp_int cupdlp_scaleVector(CUPDLPwork *w, const cupdlp_float weight,
                               cupdlp_float *x, const cupdlp_int n) {
 #if !(CUPDLP_CPU)

cupdlp/cupdlp_linalg.h

@@ -27,6 +27,8 @@ double twoNormSquared(double *x, cupdlp_int n);
 
 double infNorm(double *x, cupdlp_int n);
 
+cupdlp_int infNormIndex(double *x, cupdlp_int n);
+
 /*------------------------ new added --------------------*/
 
 double GenNorm(double *x, cupdlp_int n, cupdlp_float p);
@@ -67,13 +69,13 @@ void cupdlp_projNegative(cupdlp_float *x, const cupdlp_int len);
 
 /*------------------------ new added --------------------*/
 
-extern double dot(cupdlp_int n, cupdlp_float *x, cupdlp_int incx,
-                  cupdlp_float *y, cupdlp_int incy);
+extern double dot(cupdlp_int n, const cupdlp_float *x, cupdlp_int incx,
+                  const cupdlp_float *y, cupdlp_int incy);
 
-extern double Dotprod(cupdlp_float *x, cupdlp_float *y, cupdlp_int n);
+extern double Dotprod(const cupdlp_float *x, const cupdlp_float *y, cupdlp_int n);
 
 // todo, add this
-extern double Dotprod_Neumaier(cupdlp_float *x, cupdlp_float *y, cupdlp_int n);
+extern double Dotprod_Neumaier(const cupdlp_float *x, const cupdlp_float *y, cupdlp_int n);
 
 /* x = x + weight * y */
 void AddToVector(cupdlp_float *x, const cupdlp_float weight,
@@ -111,6 +113,12 @@ cupdlp_int cupdlp_dot(CUPDLPwork *w, const cupdlp_int n, const cupdlp_float *x,
 cupdlp_int cupdlp_twoNorm(CUPDLPwork *w, const cupdlp_int n,
                           const cupdlp_float *x, cupdlp_float *res);
 
+cupdlp_int cupdlp_infNorm(CUPDLPwork *w, const cupdlp_int n,
+                          const cupdlp_float *x, cupdlp_float *res);
+
+cupdlp_int cupdlp_infNormIndex(CUPDLPwork *w, const cupdlp_int n,
+                               const cupdlp_float *x, cupdlp_int *res);
+
 cupdlp_int cupdlp_scaleVector(CUPDLPwork *w, const cupdlp_float weight,
                               cupdlp_float *x, const cupdlp_int n);
 
@@ -172,4 +180,4 @@ void cupdlp_initvec(cupdlp_float *x, const cupdlp_float val,
 void cupdlp_compute_interaction_and_movement(CUPDLPwork *w,
                                              cupdlp_float *dMovement,
                                              cupdlp_float *dIteraction);
-#endif  // CUPDLP_CUPDLP_LINALG_H
+#endif  // CUPDLP_CUPDLP_LINALG_H

cupdlp/cupdlp_proj.c

@@ -143,4 +143,4 @@ void PDHG_Restart_Iterate_GPU(CUPDLPwork *pdhg) {
   PDHG_Compute_Residuals(pdhg);
   // cupdlp_printf("Recomputed stepsize ratio: %e,  sqrt(ratio)=%e",
   // stepsize->dBeta, sqrt(stepsize->dBeta));
-}
+}

cupdlp/cupdlp_restart.c

@@ -84,10 +84,12 @@ PDHG_restart_choice PDHG_Check_Restart_GPU(CUPDLPwork *work) {
 
   if (restart_choice != PDHG_NO_RESTART) {
     if (muCurrent < muAverage) {
-      cupdlp_printf("Last restart was iter %d: %s", iterates->iLastRestartIter,
+        if (work->settings->nLogLevel > 1)
+	  cupdlp_printf("Last restart was iter %d: %s", iterates->iLastRestartIter,
                     "current\n");
     } else {
-      cupdlp_printf("Last restart was iter %d: %s", iterates->iLastRestartIter,
+        if (work->settings->nLogLevel > 1)
+	  cupdlp_printf("Last restart was iter %d: %s", iterates->iLastRestartIter,
                     "average\n");
     }
   }
@@ -117,4 +119,4 @@ cupdlp_float PDHG_Restart_Score_GPU(cupdlp_float weightSquared,
            dDualFeas * dDualFeas / weightSquared + dDualityGap * dDualityGap);
 
   return dScoreGPU;
-}
+}