Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpSolve.cpp

    r2333 r2385  
    1616#include <math.h>
    1717#ifdef _MSC_VER
    18 #include <windows.h>   // for Sleep()
     18#include <windows.h> // for Sleep()
    1919#ifdef small
    2020#undef small
    2121#endif
    2222#else
    23 #include <unistd.h>    // for usleep()
     23#include <unistd.h> // for usleep()
    2424#endif
    2525
     
    5353#endif
    5454#include "CoinStructuredModel.hpp"
    55 double zz_slack_value=0.0;
     55double zz_slack_value = 0.0;
    5656#ifdef CLP_USEFUL_PRINTOUT
    5757double debugDouble[10];
     
    8383class lpHook : public VOL_user_hooks {
    8484private:
    85      lpHook(const lpHook&);
    86      lpHook& operator=(const lpHook&);
     85  lpHook(const lpHook &);
     86  lpHook &operator=(const lpHook &);
     87
    8788private:
    88      /// Pointer to dense vector of structural variable upper bounds
    89      double *colupper_;
    90      /// Pointer to dense vector of structural variable lower bounds
    91      double *collower_;
    92      /// Pointer to dense vector of objective coefficients
    93      double *objcoeffs_;
    94      /// Pointer to dense vector of right hand sides
    95      double *rhs_;
    96      /// Pointer to dense vector of senses
    97      char    *sense_;
     89  /// Pointer to dense vector of structural variable upper bounds
     90  double *colupper_;
     91  /// Pointer to dense vector of structural variable lower bounds
     92  double *collower_;
     93  /// Pointer to dense vector of objective coefficients
     94  double *objcoeffs_;
     95  /// Pointer to dense vector of right hand sides
     96  double *rhs_;
     97  /// Pointer to dense vector of senses
     98  char *sense_;
    9899
    99      /// The problem matrix in a row ordered form
    100      CoinPackedMatrix rowMatrix_;
    101      /// The problem matrix in a column ordered form
    102      CoinPackedMatrix colMatrix_;
     100  /// The problem matrix in a row ordered form
     101  CoinPackedMatrix rowMatrix_;
     102  /// The problem matrix in a column ordered form
     103  CoinPackedMatrix colMatrix_;
    103104
    104105public:
    105      lpHook(double* clb, double* cub, double* obj,
    106             double* rhs, char* sense, const CoinPackedMatrix& mat);
    107      virtual ~lpHook();
     106  lpHook(double *clb, double *cub, double *obj,
     107    double *rhs, char *sense, const CoinPackedMatrix &mat);
     108  virtual ~lpHook();
    108109
    109110public:
    110      // for all hooks: return value of -1 means that volume should quit
    111      /** compute reduced costs
     111  // for all hooks: return value of -1 means that volume should quit
     112  /** compute reduced costs
    112113         @param u (IN) the dual variables
    113114         @param rc (OUT) the reduced cost with respect to the dual values
    114115     */
    115      virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc);
     116  virtual int compute_rc(const VOL_dvector &u, VOL_dvector &rc);
    116117
    117      /** Solve the subproblem for the subgradient step.
     118  /** Solve the subproblem for the subgradient step.
    118119         @param dual (IN) the dual variables
    119120         @param rc (IN) the reduced cost with respect to the dual values
     
    123124         @param pcost (OUT) the primal objective value of <code>x</code>
    124125     */
    125      virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
    126                                   double& lcost, VOL_dvector& x, VOL_dvector& v,
    127                                   double& pcost);
    128      /** Starting from the primal vector x, run a heuristic to produce
     126  virtual int solve_subproblem(const VOL_dvector &dual, const VOL_dvector &rc,
     127    double &lcost, VOL_dvector &x, VOL_dvector &v,
     128    double &pcost);
     129  /** Starting from the primal vector x, run a heuristic to produce
    129130         an integer solution
    130131         @param x (IN) the primal vector
     
    132133         <code>DBL_MAX</code> here if no feas sol was found
    133134     */
    134      virtual int heuristics(const VOL_problem& p,
    135                             const VOL_dvector& x, double& heur_val) {
    136           return 0;
    137      }
     135  virtual int heuristics(const VOL_problem &p,
     136    const VOL_dvector &x, double &heur_val)
     137  {
     138    return 0;
     139  }
    138140};
    139141
    140142//#############################################################################
    141143
    142 lpHook::lpHook(double* clb, double* cub, double* obj,
    143                double* rhs, char* sense,
    144                const CoinPackedMatrix& mat)
     144lpHook::lpHook(double *clb, double *cub, double *obj,
     145  double *rhs, char *sense,
     146  const CoinPackedMatrix &mat)
    145147{
    146      colupper_ = cub;
    147      collower_ = clb;
    148      objcoeffs_ = obj;
    149      rhs_ = rhs;
    150      sense_ = sense;
    151      assert (mat.isColOrdered());
    152      colMatrix_.copyOf(mat);
    153      rowMatrix_.reverseOrderedCopyOf(mat);
     148  colupper_ = cub;
     149  collower_ = clb;
     150  objcoeffs_ = obj;
     151  rhs_ = rhs;
     152  sense_ = sense;
     153  assert(mat.isColOrdered());
     154  colMatrix_.copyOf(mat);
     155  rowMatrix_.reverseOrderedCopyOf(mat);
    154156}
    155157
     
    162164//#############################################################################
    163165
    164 int
    165 lpHook::compute_rc(const VOL_dvector& u, VOL_dvector& rc)
     166int lpHook::compute_rc(const VOL_dvector &u, VOL_dvector &rc)
    166167{
    167      rowMatrix_.transposeTimes(u.v, rc.v);
    168      const int psize = rowMatrix_.getNumCols();
     168  rowMatrix_.transposeTimes(u.v, rc.v);
     169  const int psize = rowMatrix_.getNumCols();
    169170
    170      for (int i = 0; i < psize; ++i)
    171           rc[i] = objcoeffs_[i] - rc[i];
    172      return 0;
     171  for (int i = 0; i < psize; ++i)
     172    rc[i] = objcoeffs_[i] - rc[i];
     173  return 0;
    173174}
    174175
    175176//-----------------------------------------------------------------------------
    176177
    177 int
    178 lpHook::solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
    179                          double& lcost, VOL_dvector& x, VOL_dvector& v,
    180                          double& pcost)
     178int lpHook::solve_subproblem(const VOL_dvector &dual, const VOL_dvector &rc,
     179  double &lcost, VOL_dvector &x, VOL_dvector &v,
     180  double &pcost)
    181181{
    182      int i;
    183      const int psize = x.size();
    184      const int dsize = v.size();
     182  int i;
     183  const int psize = x.size();
     184  const int dsize = v.size();
    185185
    186      // compute the lagrangean solution corresponding to the reduced costs
    187      for (i = 0; i < psize; ++i)
    188           x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
     186  // compute the lagrangean solution corresponding to the reduced costs
     187  for (i = 0; i < psize; ++i)
     188    x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
    189189
    190      // compute the lagrangean value (rhs*dual + primal*rc)
    191      lcost = 0;
    192      for (i = 0; i < dsize; ++i)
    193           lcost += rhs_[i] * dual[i];
    194      for (i = 0; i < psize; ++i)
    195           lcost += x[i] * rc[i];
     190  // compute the lagrangean value (rhs*dual + primal*rc)
     191  lcost = 0;
     192  for (i = 0; i < dsize; ++i)
     193    lcost += rhs_[i] * dual[i];
     194  for (i = 0; i < psize; ++i)
     195    lcost += x[i] * rc[i];
    196196
    197      // compute the rhs - lhs
    198      colMatrix_.times(x.v, v.v);
    199      for (i = 0; i < dsize; ++i)
    200           v[i] = rhs_[i] - v[i];
     197  // compute the rhs - lhs
     198  colMatrix_.times(x.v, v.v);
     199  for (i = 0; i < dsize; ++i)
     200    v[i] = rhs_[i] - v[i];
    201201
    202      // compute the lagrangean primal objective
    203      pcost = 0;
    204      for (i = 0; i < psize; ++i)
    205           pcost += x[i] * objcoeffs_[i];
     202  // compute the lagrangean primal objective
     203  pcost = 0;
     204  for (i = 0; i < psize; ++i)
     205    pcost += x[i] * objcoeffs_[i];
    206206
    207      return 0;
     207  return 0;
    208208}
    209209
     
    213213inline void
    214214convertBoundToSense(const double lower, const double upper,
    215                     char& sense, double& right,
    216                     double& range)
     215  char &sense, double &right,
     216  double &range)
    217217{
    218      range = 0.0;
    219      if (lower > -1.0e20) {
    220           if (upper < 1.0e20) {
    221                right = upper;
    222                if (upper == lower) {
    223                     sense = 'E';
    224                } else {
    225                     sense = 'R';
    226                     range = upper - lower;
    227                }
    228           } else {
    229                sense = 'G';
    230                right = lower;
    231           }
    232      } else {
    233           if (upper < 1.0e20) {
    234                sense = 'L';
    235                right = upper;
    236           } else {
    237                sense = 'N';
    238                right = 0.0;
    239           }
    240      }
     218  range = 0.0;
     219  if (lower > -1.0e20) {
     220    if (upper < 1.0e20) {
     221      right = upper;
     222      if (upper == lower) {
     223        sense = 'E';
     224      } else {
     225        sense = 'R';
     226        range = upper - lower;
     227      }
     228    } else {
     229      sense = 'G';
     230      right = lower;
     231    }
     232  } else {
     233    if (upper < 1.0e20) {
     234      sense = 'L';
     235      right = upper;
     236    } else {
     237      sense = 'N';
     238      right = 0.0;
     239    }
     240  }
    241241}
    242242
    243243static int
    244 solveWithVolume(ClpSimplex * model, int numberPasses, int doIdiot)
     244solveWithVolume(ClpSimplex *model, int numberPasses, int doIdiot)
    245245{
    246      VOL_problem volprob;
    247      volprob.parm.gap_rel_precision = 0.00001;
    248      volprob.parm.maxsgriters = 3000;
    249      if(numberPasses > 3000) {
    250           volprob.parm.maxsgriters = numberPasses;
    251           volprob.parm.primal_abs_precision = 0.0;
    252           volprob.parm.minimum_rel_ascent = 0.00001;
    253      } else if (doIdiot > 0) {
    254           volprob.parm.maxsgriters = doIdiot;
    255      }
    256      if (model->logLevel() < 2)
    257           volprob.parm.printflag = 0;
    258      else
    259           volprob.parm.printflag = 3;
    260      const CoinPackedMatrix* mat = model->matrix();
    261      int psize = model->numberColumns();
    262      int dsize = model->numberRows();
    263      char * sense = new char[dsize];
    264      double * rhs = new double [dsize];
     246  VOL_problem volprob;
     247  volprob.parm.gap_rel_precision = 0.00001;
     248  volprob.parm.maxsgriters = 3000;
     249  if (numberPasses > 3000) {
     250    volprob.parm.maxsgriters = numberPasses;
     251    volprob.parm.primal_abs_precision = 0.0;
     252    volprob.parm.minimum_rel_ascent = 0.00001;
     253  } else if (doIdiot > 0) {
     254    volprob.parm.maxsgriters = doIdiot;
     255  }
     256  if (model->logLevel() < 2)
     257    volprob.parm.printflag = 0;
     258  else
     259    volprob.parm.printflag = 3;
     260  const CoinPackedMatrix *mat = model->matrix();
     261  int psize = model->numberColumns();
     262  int dsize = model->numberRows();
     263  char *sense = new char[dsize];
     264  double *rhs = new double[dsize];
    265265
    266      // Set the lb/ub on the duals
    267      volprob.dsize = dsize;
    268      volprob.psize = psize;
    269      volprob.dual_lb.allocate(dsize);
    270      volprob.dual_ub.allocate(dsize);
    271      int i;
    272      const double * rowLower = model->rowLower();
    273      const double * rowUpper = model->rowUpper();
    274      for (i = 0; i < dsize; ++i) {
    275           double range;
    276           convertBoundToSense(rowLower[i], rowUpper[i],
    277                               sense[i], rhs[i], range);
    278           switch (sense[i]) {
    279           case 'E':
    280                volprob.dual_lb[i] = -1.0e31;
    281                volprob.dual_ub[i] = 1.0e31;
    282                break;
    283           case 'L':
    284                volprob.dual_lb[i] = -1.0e31;
    285                volprob.dual_ub[i] = 0.0;
    286                break;
    287           case 'G':
    288                volprob.dual_lb[i] = 0.0;
    289                volprob.dual_ub[i] = 1.0e31;
    290                break;
    291           default:
    292                printf("Volume Algorithm can't work if there is a non ELG row\n");
    293                return 1;
    294           }
    295      }
    296      // Check bounds
    297      double * saveLower = model->columnLower();
    298      double * saveUpper = model->columnUpper();
    299      bool good = true;
    300      for (i = 0; i < psize; i++) {
    301           if (saveLower[i] < -1.0e20 || saveUpper[i] > 1.0e20) {
    302                good = false;
    303                break;
    304           }
    305      }
    306      if (!good) {
    307           saveLower = CoinCopyOfArray(model->columnLower(), psize);
    308           saveUpper = CoinCopyOfArray(model->columnUpper(), psize);
    309           for (i = 0; i < psize; i++) {
    310                if (saveLower[i] < -1.0e20)
    311                     saveLower[i] = -1.0e20;
    312                if(saveUpper[i] > 1.0e20)
    313                     saveUpper[i] = 1.0e20;
    314           }
    315      }
    316      lpHook myHook(saveLower, saveUpper,
    317                    model->objective(),
    318                    rhs, sense, *mat);
     266  // Set the lb/ub on the duals
     267  volprob.dsize = dsize;
     268  volprob.psize = psize;
     269  volprob.dual_lb.allocate(dsize);
     270  volprob.dual_ub.allocate(dsize);
     271  int i;
     272  const double *rowLower = model->rowLower();
     273  const double *rowUpper = model->rowUpper();
     274  for (i = 0; i < dsize; ++i) {
     275    double range;
     276    convertBoundToSense(rowLower[i], rowUpper[i],
     277      sense[i], rhs[i], range);
     278    switch (sense[i]) {
     279    case 'E':
     280      volprob.dual_lb[i] = -1.0e31;
     281      volprob.dual_ub[i] = 1.0e31;
     282      break;
     283    case 'L':
     284      volprob.dual_lb[i] = -1.0e31;
     285      volprob.dual_ub[i] = 0.0;
     286      break;
     287    case 'G':
     288      volprob.dual_lb[i] = 0.0;
     289      volprob.dual_ub[i] = 1.0e31;
     290      break;
     291    default:
     292      printf("Volume Algorithm can't work if there is a non ELG row\n");
     293      return 1;
     294    }
     295  }
     296  // Check bounds
     297  double *saveLower = model->columnLower();
     298  double *saveUpper = model->columnUpper();
     299  bool good = true;
     300  for (i = 0; i < psize; i++) {
     301    if (saveLower[i] < -1.0e20 || saveUpper[i] > 1.0e20) {
     302      good = false;
     303      break;
     304    }
     305  }
     306  if (!good) {
     307    saveLower = CoinCopyOfArray(model->columnLower(), psize);
     308    saveUpper = CoinCopyOfArray(model->columnUpper(), psize);
     309    for (i = 0; i < psize; i++) {
     310      if (saveLower[i] < -1.0e20)
     311        saveLower[i] = -1.0e20;
     312      if (saveUpper[i] > 1.0e20)
     313        saveUpper[i] = 1.0e20;
     314    }
     315  }
     316  lpHook myHook(saveLower, saveUpper,
     317    model->objective(),
     318    rhs, sense, *mat);
    319319
    320      volprob.solve(myHook, false /* no warmstart */);
     320  volprob.solve(myHook, false /* no warmstart */);
    321321
    322      if (saveLower != model->columnLower()) {
    323           delete [] saveLower;
    324           delete [] saveUpper;
    325      }
    326      //------------- extract the solution ---------------------------
     322  if (saveLower != model->columnLower()) {
     323    delete[] saveLower;
     324    delete[] saveUpper;
     325  }
     326  //------------- extract the solution ---------------------------
    327327
    328      //printf("Best lagrangean value: %f\n", volprob.value);
     328  //printf("Best lagrangean value: %f\n", volprob.value);
    329329
    330      double avg = 0;
    331      for (i = 0; i < dsize; ++i) {
    332           switch (sense[i]) {
    333           case 'E':
    334                avg += CoinAbs(volprob.viol[i]);
    335                break;
    336           case 'L':
    337                if (volprob.viol[i] < 0)
    338                     avg += (-volprob.viol[i]);
    339                break;
    340           case 'G':
    341                if (volprob.viol[i] > 0)
    342                     avg += volprob.viol[i];
    343                break;
    344           }
    345      }
     330  double avg = 0;
     331  for (i = 0; i < dsize; ++i) {
     332    switch (sense[i]) {
     333    case 'E':
     334      avg += CoinAbs(volprob.viol[i]);
     335      break;
     336    case 'L':
     337      if (volprob.viol[i] < 0)
     338        avg += (-volprob.viol[i]);
     339      break;
     340    case 'G':
     341      if (volprob.viol[i] > 0)
     342        avg += volprob.viol[i];
     343      break;
     344    }
     345  }
    346346
    347      //printf("Average primal constraint violation: %f\n", avg/dsize);
     347  //printf("Average primal constraint violation: %f\n", avg/dsize);
    348348
    349      // volprob.dsol contains the dual solution (dual feasible)
    350      // volprob.psol contains the primal solution
    351      //              (NOT necessarily primal feasible)
    352      CoinMemcpyN(volprob.dsol.v, dsize, model->dualRowSolution());
    353      CoinMemcpyN(volprob.psol.v, psize, model->primalColumnSolution());
    354      return 0;
     349  // volprob.dsol contains the dual solution (dual feasible)
     350  // volprob.psol contains the primal solution
     351  //              (NOT necessarily primal feasible)
     352  CoinMemcpyN(volprob.dsol.v, dsize, model->dualRowSolution());
     353  CoinMemcpyN(volprob.psol.v, psize, model->primalColumnSolution());
     354  return 0;
    355355}
    356356#endif
    357 static ClpInterior * currentModel2 = NULL;
     357static ClpInterior *currentModel2 = NULL;
    358358#endif
    359359//#############################################################################
     
    362362
    363363#include "CoinSignal.hpp"
    364 static ClpSimplex * currentModel = NULL;
     364static ClpSimplex *currentModel = NULL;
    365365#ifdef ABC_INHERIT
    366 static AbcSimplex * currentAbcModel = NULL;
     366static AbcSimplex *currentAbcModel = NULL;
    367367#endif
    368368
    369369extern "C" {
    370      static void
     370static void
    371371#if defined(_MSC_VER)
    372      __cdecl
     372  __cdecl
    373373#endif // _MSC_VER
    374      signal_handler(int /*whichSignal*/)
    375      {
    376           if (currentModel != NULL)
    377                currentModel->setMaximumIterations(0); // stop at next iterations
     374  signal_handler(int /*whichSignal*/)
     375{
     376  if (currentModel != NULL)
     377    currentModel->setMaximumIterations(0); // stop at next iterations
    378378#ifdef ABC_INHERIT
    379           if (currentAbcModel != NULL)
    380                currentAbcModel->setMaximumIterations(0); // stop at next iterations
     379  if (currentAbcModel != NULL)
     380    currentAbcModel->setMaximumIterations(0); // stop at next iterations
    381381#endif
    382382#ifndef SLIM_CLP
    383           if (currentModel2 != NULL)
    384                currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
    385 #endif
    386           return;
    387      }
     383  if (currentModel2 != NULL)
     384    currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
     385#endif
     386  return;
    388387}
    389 #if ABC_INSTRUMENT>1
     388}
     389#if ABC_INSTRUMENT > 1
    390390int abcPricing[20];
    391391int abcPricingDense[20];
     
    395395#define MAX_COUNT 20
    396396#define MAX_FRACTION 101
    397 static char * types[MAX_TYPES];
     397static char *types[MAX_TYPES];
    398398static double counts[MAX_TYPES][MAX_COUNT];
    399399static double countsFraction[MAX_TYPES][MAX_FRACTION];
    400 static double * currentCounts;
    401 static double * currentCountsFraction;
     400static double *currentCounts;
     401static double *currentCountsFraction;
    402402static int currentType;
    403403static double workMultiplier[MAX_TYPES];
     
    410410void instrument_initialize(int numberRows)
    411411{
    412   trueNumberRows=numberRows;
    413   numberTypes=0;
    414   memset(counts,0,sizeof(counts));
    415   currentCounts=NULL;
    416   memset(countsFraction,0,sizeof(countsFraction));
    417   currentCountsFraction=NULL;
    418   memset(workMultiplier,0,sizeof(workMultiplier));
    419   memset(work,0,sizeof(work));
    420   memset(otherWork,0,sizeof(otherWork));
    421   memset(timesCalled,0,sizeof(timesCalled));
    422   memset(timesStarted,0,sizeof(timesStarted));
    423   currentType=-1;
    424   fractionDivider=(numberRows+MAX_FRACTION-2)/(MAX_FRACTION-1);
     412  trueNumberRows = numberRows;
     413  numberTypes = 0;
     414  memset(counts, 0, sizeof(counts));
     415  currentCounts = NULL;
     416  memset(countsFraction, 0, sizeof(countsFraction));
     417  currentCountsFraction = NULL;
     418  memset(workMultiplier, 0, sizeof(workMultiplier));
     419  memset(work, 0, sizeof(work));
     420  memset(otherWork, 0, sizeof(otherWork));
     421  memset(timesCalled, 0, sizeof(timesCalled));
     422  memset(timesStarted, 0, sizeof(timesStarted));
     423  currentType = -1;
     424  fractionDivider = (numberRows + MAX_FRACTION - 2) / (MAX_FRACTION - 1);
    425425}
    426 void instrument_start(const char * type,int numberRowsEtc)
     426void instrument_start(const char *type, int numberRowsEtc)
    427427{
    428   if (currentType>=0)
     428  if (currentType >= 0)
    429429    instrument_end();
    430   currentType=-1;
    431   currentWork=0.0;
    432   for (int i=0;i<numberTypes;i++) {
    433     if (!strcmp(types[i],type)) {
    434       currentType=i;
     430  currentType = -1;
     431  currentWork = 0.0;
     432  for (int i = 0; i < numberTypes; i++) {
     433    if (!strcmp(types[i], type)) {
     434      currentType = i;
    435435      break;
    436436    }
    437437  }
    438   if (currentType==-1) {
    439     assert (numberTypes<MAX_TYPES);
    440     currentType=numberTypes;
    441     types[numberTypes++]=strdup(type);
     438  if (currentType == -1) {
     439    assert(numberTypes < MAX_TYPES);
     440    currentType = numberTypes;
     441    types[numberTypes++] = strdup(type);
    442442  }
    443443  currentCounts = &counts[currentType][0];
    444444  currentCountsFraction = &countsFraction[currentType][0];
    445445  timesStarted[currentType]++;
    446   assert (trueNumberRows);
    447   workMultiplier[currentType]+=static_cast<double>(numberRowsEtc)/static_cast<double>(trueNumberRows);
     446  assert(trueNumberRows);
     447  workMultiplier[currentType] += static_cast< double >(numberRowsEtc) / static_cast< double >(trueNumberRows);
    448448}
    449449void instrument_add(int count)
    450450{
    451   assert (currentType>=0);
    452   currentWork+=count;
     451  assert(currentType >= 0);
     452  currentWork += count;
    453453  timesCalled[currentType]++;
    454   if (count<MAX_COUNT-1)
     454  if (count < MAX_COUNT - 1)
    455455    currentCounts[count]++;
    456456  else
    457     currentCounts[MAX_COUNT-1]++;
    458   assert(count/fractionDivider>=0&&count/fractionDivider<MAX_FRACTION);
    459   currentCountsFraction[count/fractionDivider]++;
     457    currentCounts[MAX_COUNT - 1]++;
     458  assert(count / fractionDivider >= 0 && count / fractionDivider < MAX_FRACTION);
     459  currentCountsFraction[count / fractionDivider]++;
    460460}
    461 void instrument_do(const char * type,double count)
     461void instrument_do(const char *type, double count)
    462462{
    463   int iType=-1;
    464   for (int i=0;i<numberTypes;i++) {
    465     if (!strcmp(types[i],type)) {
    466       iType=i;
     463  int iType = -1;
     464  for (int i = 0; i < numberTypes; i++) {
     465    if (!strcmp(types[i], type)) {
     466      iType = i;
    467467      break;
    468468    }
    469469  }
    470   if (iType==-1) {
    471     assert (numberTypes<MAX_TYPES);
    472     iType=numberTypes;
    473     types[numberTypes++]=strdup(type);
     470  if (iType == -1) {
     471    assert(numberTypes < MAX_TYPES);
     472    iType = numberTypes;
     473    types[numberTypes++] = strdup(type);
    474474  }
    475475  timesStarted[iType]++;
    476   otherWork[iType]+=count;
     476  otherWork[iType] += count;
    477477}
    478478void instrument_end()
    479479{
    480   work[currentType]+=currentWork;
    481   currentType=-1;
     480  work[currentType] += currentWork;
     481  currentType = -1;
    482482}
    483483void instrument_end_and_adjust(double factor)
    484484{
    485   work[currentType]+=currentWork*factor;
    486   currentType=-1;
     485  work[currentType] += currentWork * factor;
     486  currentType = -1;
    487487}
    488488void instrument_print()
    489489{
    490   for (int iType=0;iType<numberTypes;iType++) {
     490  for (int iType = 0; iType < numberTypes; iType++) {
    491491    currentCounts = &counts[iType][0];
    492492    currentCountsFraction = &countsFraction[iType][0];
    493493    if (!otherWork[iType]) {
    494494      printf("%s started %d times, used %d times, work %g (average length %.1f) multiplier %g\n",
    495              types[iType],timesStarted[iType],timesCalled[iType],
    496              work[iType],work[iType]/(timesCalled[iType]+1.0e-100),workMultiplier[iType]/(timesStarted[iType]+1.0e-100));
    497       int n=0;
    498       for (int i=0;i<MAX_COUNT-1;i++) {
    499         if (currentCounts[i]) {
    500           if (n==5) {
    501             n=0;
    502             printf("\n");
    503           }
    504           n++;
    505           printf("(%d els,%.0f times) ",i,currentCounts[i]);
    506         }
    507       }
    508       if (currentCounts[MAX_COUNT-1]) {
    509         if (n==5) {
    510           n=0;
    511           printf("\n");
    512         }
    513         n++;
    514         printf("(>=%d els,%.0f times) ",MAX_COUNT-1,currentCounts[MAX_COUNT-1]);
     495        types[iType], timesStarted[iType], timesCalled[iType],
     496        work[iType], work[iType] / (timesCalled[iType] + 1.0e-100), workMultiplier[iType] / (timesStarted[iType] + 1.0e-100));
     497      int n = 0;
     498      for (int i = 0; i < MAX_COUNT - 1; i++) {
     499        if (currentCounts[i]) {
     500          if (n == 5) {
     501            n = 0;
     502            printf("\n");
     503          }
     504          n++;
     505          printf("(%d els,%.0f times) ", i, currentCounts[i]);
     506        }
     507      }
     508      if (currentCounts[MAX_COUNT - 1]) {
     509        if (n == 5) {
     510          n = 0;
     511          printf("\n");
     512        }
     513        n++;
     514        printf("(>=%d els,%.0f times) ", MAX_COUNT - 1, currentCounts[MAX_COUNT - 1]);
    515515      }
    516516      printf("\n");
    517517      int largestFraction;
    518       int nBig=0;
    519       for (largestFraction=MAX_FRACTION-1;largestFraction>=10;largestFraction--) {
    520         double count = currentCountsFraction[largestFraction];
    521         if (count&&largestFraction>10)
    522           nBig++;
    523         if (nBig>4)
    524           break;
    525       }
    526       int chunk=(largestFraction+5)/10;
    527       int lo=0;
    528       for (int iChunk=0;iChunk<largestFraction;iChunk+=chunk) {
    529         int hi=CoinMin(lo+chunk*fractionDivider,trueNumberRows);
    530         double sum=0.0;
    531         for (int i=iChunk;i<CoinMin(iChunk+chunk,MAX_FRACTION);i++)
    532           sum += currentCountsFraction[i];
    533         if (sum)
    534           printf("(%d-%d %.0f) ",lo,hi,sum);
    535         lo=hi;
    536       }
    537       for (int i=lo/fractionDivider;i<MAX_FRACTION;i++) {
    538         if (currentCountsFraction[i])
    539           printf("(%d %.0f) ",i*fractionDivider,currentCountsFraction[i]);
     518      int nBig = 0;
     519      for (largestFraction = MAX_FRACTION - 1; largestFraction >= 10; largestFraction--) {
     520        double count = currentCountsFraction[largestFraction];
     521        if (count && largestFraction > 10)
     522          nBig++;
     523        if (nBig > 4)
     524          break;
     525      }
     526      int chunk = (largestFraction + 5) / 10;
     527      int lo = 0;
     528      for (int iChunk = 0; iChunk < largestFraction; iChunk += chunk) {
     529        int hi = CoinMin(lo + chunk * fractionDivider, trueNumberRows);
     530        double sum = 0.0;
     531        for (int i = iChunk; i < CoinMin(iChunk + chunk, MAX_FRACTION); i++)
     532          sum += currentCountsFraction[i];
     533        if (sum)
     534          printf("(%d-%d %.0f) ", lo, hi, sum);
     535        lo = hi;
     536      }
     537      for (int i = lo / fractionDivider; i < MAX_FRACTION; i++) {
     538        if (currentCountsFraction[i])
     539          printf("(%d %.0f) ", i * fractionDivider, currentCountsFraction[i]);
    540540      }
    541541      printf("\n");
    542542    } else {
    543543      printf("%s started %d times, used %d times, work %g multiplier %g other work %g\n",
    544              types[iType],timesStarted[iType],timesCalled[iType],
    545              work[iType],workMultiplier[iType],otherWork[iType]);
     544        types[iType], timesStarted[iType], timesCalled[iType],
     545        work[iType], workMultiplier[iType], otherWork[iType]);
    546546    }
    547547    free(types[iType]);
     
    549549}
    550550#endif
    551 #if ABC_PARALLEL==2
     551#if ABC_PARALLEL == 2
    552552#ifndef FAKE_CILK
    553 int number_cilk_workers=0;
     553int number_cilk_workers = 0;
    554554#include <cilk/cilk_api.h>
    555555#endif
     
    558558AbcSimplex *
    559559ClpSimplex::dealWithAbc(int solveType, int startUp,
    560                         bool interrupt)
     560  bool interrupt)
    561561{
    562   bool keepAbc = startUp>=10;
    563   startUp=startUp%10;
    564   AbcSimplex * abcModel2 = NULL;
    565   if (!this->abcState()||!numberRows_||!numberColumns_) {
     562  bool keepAbc = startUp >= 10;
     563  startUp = startUp % 10;
     564  AbcSimplex *abcModel2 = NULL;
     565  if (!this->abcState() || !numberRows_ || !numberColumns_) {
    566566    //this->readBasis("aaa.bas");
    567567    if (!solveType) {
    568568      this->dual(0);
    569     } else if (solveType==1) {
    570       int ifValuesPass=startUp ? 1 : 0;
    571       if (startUp==3)
    572         ifValuesPass=2;
     569    } else if (solveType == 1) {
     570      int ifValuesPass = startUp ? 1 : 0;
     571      if (startUp == 3)
     572        ifValuesPass = 2;
    573573      this->primal(ifValuesPass);
    574574    }
    575575    //this->writeBasis("a.bas",true);
    576576  } else {
    577     abcModel2=new AbcSimplex(*this);
     577    abcModel2 = new AbcSimplex(*this);
    578578    if (interrupt)
    579579      currentAbcModel = abcModel2;
    580580    //if (abcSimplex_) {
    581581    // move factorization stuff
    582     abcModel2->factorization()->synchronize(this->factorization(),abcModel2);
     582    abcModel2->factorization()->synchronize(this->factorization(), abcModel2);
    583583    //}
    584584    //abcModel2->startPermanentArrays();
    585     int crashState=abcModel2->abcState()&(256+512+1024);
    586     abcModel2->setAbcState(CLP_ABC_WANTED|crashState|(abcModel2->abcState()&15));
    587     int ifValuesPass=startUp ? 1 : 0;
    588     if (startUp==3)
    589       ifValuesPass=2;
     585    int crashState = abcModel2->abcState() & (256 + 512 + 1024);
     586    abcModel2->setAbcState(CLP_ABC_WANTED | crashState | (abcModel2->abcState() & 15));
     587    int ifValuesPass = startUp ? 1 : 0;
     588    if (startUp == 3)
     589      ifValuesPass = 2;
    590590    // temp
    591     if (fabs(this->primalTolerance()-1.001e-6)<0.999e-9) {
    592       int type=1;
    593       double diff=this->primalTolerance()-1.001e-6;
    594       printf("Diff %g\n",diff);
    595       if (fabs(diff-1.0e-10)<1.0e-13)
    596         type=2;
    597       else if (fabs(diff-1.0e-11)<1.0e-13)
    598         type=3;
     591    if (fabs(this->primalTolerance() - 1.001e-6) < 0.999e-9) {
     592      int type = 1;
     593      double diff = this->primalTolerance() - 1.001e-6;
     594      printf("Diff %g\n", diff);
     595      if (fabs(diff - 1.0e-10) < 1.0e-13)
     596        type = 2;
     597      else if (fabs(diff - 1.0e-11) < 1.0e-13)
     598        type = 3;
    599599#if 0
    600600      ClpSimplex * thisModel = static_cast<ClpSimplexOther *> (this)->dualOfModel(1.0,1.0);
     
    619619#else
    620620      if (!solveType) {
    621         this->dual(0);
    622         abcModel2->setupDualValuesPass(this->dualRowSolution(),
    623                                        this->primalColumnSolution(),type);
    624       } else if (solveType==1) {
    625         ifValuesPass=1;
    626         abcModel2->setStateOfProblem(abcModel2->stateOfProblem() | VALUES_PASS);
    627         Idiot info(*abcModel2);
    628         info.setStrategy(512 | info.getStrategy());
    629         // Allow for scaling
    630         info.setStrategy(32 | info.getStrategy());
    631         info.setStartingWeight(1.0e3);
    632         info.setReduceIterations(6);
    633         info.crash(200, abcModel2->messageHandler(), abcModel2->messagesPointer(),false);
    634         //memcpy(abcModel2->primalColumnSolution(),this->primalColumnSolution(),
    635         //     this->numberColumns()*sizeof(double));
    636       }
    637 #endif
    638     }
    639     int numberCpu=this->abcState()&15;
    640     if (numberCpu==9) {
    641       numberCpu=1;
    642 #if ABC_PARALLEL==2
     621        this->dual(0);
     622        abcModel2->setupDualValuesPass(this->dualRowSolution(),
     623          this->primalColumnSolution(), type);
     624      } else if (solveType == 1) {
     625        ifValuesPass = 1;
     626        abcModel2->setStateOfProblem(abcModel2->stateOfProblem() | VALUES_PASS);
     627        Idiot info(*abcModel2);
     628        info.setStrategy(512 | info.getStrategy());
     629        // Allow for scaling
     630        info.setStrategy(32 | info.getStrategy());
     631        info.setStartingWeight(1.0e3);
     632        info.setReduceIterations(6);
     633        info.crash(200, abcModel2->messageHandler(), abcModel2->messagesPointer(), false);
     634        //memcpy(abcModel2->primalColumnSolution(),this->primalColumnSolution(),
     635        //     this->numberColumns()*sizeof(double));
     636      }
     637#endif
     638    }
     639    int numberCpu = this->abcState() & 15;
     640    if (numberCpu == 9) {
     641      numberCpu = 1;
     642#if ABC_PARALLEL == 2
    643643#ifndef FAKE_CILK
    644       if (number_cilk_workers>1)
    645       numberCpu=CoinMin(2*number_cilk_workers,8);
    646 #endif
    647 #endif
    648     } else if (numberCpu==10) {
     644      if (number_cilk_workers > 1)
     645        numberCpu = CoinMin(2 * number_cilk_workers, 8);
     646#endif
     647#endif
     648    } else if (numberCpu == 10) {
    649649      // maximum
    650       numberCpu=4;
    651     } else if (numberCpu==10) {
     650      numberCpu = 4;
     651    } else if (numberCpu == 10) {
    652652      // decide
    653       if (abcModel2->getNumElements()<5000)
    654         numberCpu=1;
    655 #if ABC_PARALLEL==2
     653      if (abcModel2->getNumElements() < 5000)
     654        numberCpu = 1;
     655#if ABC_PARALLEL == 2
    656656#ifndef FAKE_CILK
    657       else if (number_cilk_workers>1)
    658         numberCpu=CoinMin(2*number_cilk_workers,8);
     657      else if (number_cilk_workers > 1)
     658        numberCpu = CoinMin(2 * number_cilk_workers, 8);
    659659#endif
    660660#endif
    661661      else
    662         numberCpu=1;
     662        numberCpu = 1;
    663663    } else {
    664 #if ABC_PARALLEL==2
     664#if ABC_PARALLEL == 2
    665665#ifndef FAKE_CILK
    666666      char temp[3];
    667       sprintf(temp,"%d",numberCpu);
    668       __cilkrts_set_param("nworkers",temp);
    669       printf("setting cilk workers to %d\n",numberCpu);
    670       number_cilk_workers=numberCpu;
     667      sprintf(temp, "%d", numberCpu);
     668      __cilkrts_set_param("nworkers", temp);
     669      printf("setting cilk workers to %d\n", numberCpu);
     670      number_cilk_workers = numberCpu;
    671671
    672672#endif
     
    675675    char line[200];
    676676#if ABC_PARALLEL
    677 #if ABC_PARALLEL==2
     677#if ABC_PARALLEL == 2
    678678#ifndef FAKE_CILK
    679679    if (!number_cilk_workers) {
    680       number_cilk_workers=__cilkrts_get_nworkers();
    681       sprintf(line,"%d cilk workers",number_cilk_workers);
     680      number_cilk_workers = __cilkrts_get_nworkers();
     681      sprintf(line, "%d cilk workers", number_cilk_workers);
    682682      handler_->message(CLP_GENERAL, messages_)
    683         << line
    684         << CoinMessageEol;
    685     }
    686 #endif
    687 #endif
    688     abcModel2->setParallelMode(numberCpu-1);
     683        << line
     684        << CoinMessageEol;
     685    }
     686#endif
     687#endif
     688    abcModel2->setParallelMode(numberCpu - 1);
    689689#endif
    690690    //if (abcState()==3||abcState()==4) {
     
    694694    //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|131072*processTune);
    695695#if ABC_INSTRUMENT
    696     double startTimeCpu=CoinCpuTime();
    697     double startTimeElapsed=CoinGetTimeOfDay();
    698 #if ABC_INSTRUMENT>1
    699     memset(abcPricing,0,sizeof(abcPricing));
    700     memset(abcPricingDense,0,sizeof(abcPricing));
     696    double startTimeCpu = CoinCpuTime();
     697    double startTimeElapsed = CoinGetTimeOfDay();
     698#if ABC_INSTRUMENT > 1
     699    memset(abcPricing, 0, sizeof(abcPricing));
     700    memset(abcPricingDense, 0, sizeof(abcPricing));
    701701    instrument_initialize(abcModel2->numberRows());
    702702#endif
     
    704704    if (!solveType) {
    705705      abcModel2->ClpSimplex::doAbcDual();
    706     } else if (solveType==1) {
    707       int saveOptions=abcModel2->specialOptions();
     706    } else if (solveType == 1) {
     707      int saveOptions = abcModel2->specialOptions();
    708708      //if (startUp==2)
    709709      //abcModel2->setSpecialOptions(8192|saveOptions);
     
    712712    }
    713713#if ABC_INSTRUMENT
    714     if (solveType<2) {
    715       double timeCpu=CoinCpuTime()-startTimeCpu;
    716       double timeElapsed=CoinGetTimeOfDay()-startTimeElapsed;
    717       sprintf(line,"Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
    718               this->problemName().c_str(),this->numberRows(),this->numberColumns(),
    719               this->getNumElements(),
    720               timeCpu,timeElapsed,timeElapsed ? timeCpu/timeElapsed : 1.0,abcModel2->numberIterations());
     714    if (solveType < 2) {
     715      double timeCpu = CoinCpuTime() - startTimeCpu;
     716      double timeElapsed = CoinGetTimeOfDay() - startTimeElapsed;
     717      sprintf(line, "Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
     718        this->problemName().c_str(), this->numberRows(), this->numberColumns(),
     719        this->getNumElements(),
     720        timeCpu, timeElapsed, timeElapsed ? timeCpu / timeElapsed : 1.0, abcModel2->numberIterations());
    721721      handler_->message(CLP_GENERAL, messages_)
    722         << line
    723         << CoinMessageEol;
    724 #if ABC_INSTRUMENT>1
     722        << line
     723        << CoinMessageEol;
     724#if ABC_INSTRUMENT > 1
    725725      {
    726         int n;
    727         n=0;
    728         for (int i=0;i<20;i++)
    729           n+= abcPricing[i];
    730         printf("CCSparse pricing done %d times",n);
    731         int n2=0;
    732         for (int i=0;i<20;i++)
    733           n2+= abcPricingDense[i];
    734         if (n2)
    735           printf(" and dense pricing done %d times\n",n2);
    736         else
    737           printf("\n");
    738         n=0;
    739         printf ("CCS");
    740         for (int i=0;i<19;i++) {
    741           if (abcPricing[i]) {
    742             if (n==5) {
    743               n=0;
    744               printf("\nCCS");
    745             }
    746             n++;
    747             printf("(%d els,%d times) ",i,abcPricing[i]);
    748           }
    749         }
    750         if (abcPricing[19]) {
    751           if (n==5) {
    752             n=0;
    753             printf("\nCCS");
    754           }
    755           n++;
    756           printf("(>=19 els,%d times) ",abcPricing[19]);
    757         }
    758         if (n2) {
    759           printf ("CCD");
    760           for (int i=0;i<19;i++) {
    761             if (abcPricingDense[i]) {
    762               if (n==5) {
    763                 n=0;
    764                 printf("\nCCD");
    765               }
    766               n++;
    767               int k1=(numberRows_/16)*i;;
    768               int k2=CoinMin(numberRows_,k1+(numberRows_/16)-1);
    769               printf("(%d-%d els,%d times) ",k1,k2,abcPricingDense[i]);
    770             }
    771           }
    772         }
    773         printf("\n");
     726        int n;
     727        n = 0;
     728        for (int i = 0; i < 20; i++)
     729          n += abcPricing[i];
     730        printf("CCSparse pricing done %d times", n);
     731        int n2 = 0;
     732        for (int i = 0; i < 20; i++)
     733          n2 += abcPricingDense[i];
     734        if (n2)
     735          printf(" and dense pricing done %d times\n", n2);
     736        else
     737          printf("\n");
     738        n = 0;
     739        printf("CCS");
     740        for (int i = 0; i < 19; i++) {
     741          if (abcPricing[i]) {
     742            if (n == 5) {
     743              n = 0;
     744              printf("\nCCS");
     745            }
     746            n++;
     747            printf("(%d els,%d times) ", i, abcPricing[i]);
     748          }
     749        }
     750        if (abcPricing[19]) {
     751          if (n == 5) {
     752            n = 0;
     753            printf("\nCCS");
     754          }
     755          n++;
     756          printf("(>=19 els,%d times) ", abcPricing[19]);
     757        }
     758        if (n2) {
     759          printf("CCD");
     760          for (int i = 0; i < 19; i++) {
     761            if (abcPricingDense[i]) {
     762              if (n == 5) {
     763                n = 0;
     764                printf("\nCCD");
     765              }
     766              n++;
     767              int k1 = (numberRows_ / 16) * i;
     768              ;
     769              int k2 = CoinMin(numberRows_, k1 + (numberRows_ / 16) - 1);
     770              printf("(%d-%d els,%d times) ", k1, k2, abcPricingDense[i]);
     771            }
     772          }
     773        }
     774        printf("\n");
    774775      }
    775776      instrument_print();
     
    782783      double offset;
    783784      CoinMemcpyN(this->objectiveAsObject()->gradient(this,
    784                                                       this->primalColumnSolution(), offset, true),
    785                   numberColumns_, this->dualColumnSolution());
     785                    this->primalColumnSolution(), offset, true),
     786        numberColumns_, this->dualColumnSolution());
    786787      this->clpMatrix()->transposeTimes(-1.0,
    787                                         this->dualRowSolution(),
    788                                         this->dualColumnSolution());
     788        this->dualRowSolution(),
     789        this->dualColumnSolution());
    789790      memset(this->primalRowSolution(), 0, numberRows_ * sizeof(double));
    790       this->clpMatrix()->times(1.0, 
    791                                this->primalColumnSolution(),
    792                                this->primalRowSolution());
     791      this->clpMatrix()->times(1.0,
     792        this->primalColumnSolution(),
     793        this->primalRowSolution());
    793794      this->checkSolutionInternal();
    794       if (sumDualInfeasibilities_>100.0*dualTolerance_) {
    795 #if CBC_USEFUL_PRINTING>0
    796         printf("internal check on duals failed %d %g\n",
    797                numberDualInfeasibilities_,sumDualInfeasibilities_);
     795      if (sumDualInfeasibilities_ > 100.0 * dualTolerance_) {
     796#if CBC_USEFUL_PRINTING > 0
     797        printf("internal check on duals failed %d %g\n",
     798          numberDualInfeasibilities_, sumDualInfeasibilities_);
    798799#endif
    799800      } else {
    800         sumDualInfeasibilities_=0.0;
    801         numberDualInfeasibilities_=0;
    802       }
    803       if (sumPrimalInfeasibilities_>100.0*primalTolerance_) {
    804 #if CBC_USEFUL_PRINTING>0
    805         printf("internal check on primals failed %d %g\n",
    806                numberPrimalInfeasibilities_,sumPrimalInfeasibilities_);
     801        sumDualInfeasibilities_ = 0.0;
     802        numberDualInfeasibilities_ = 0;
     803      }
     804      if (sumPrimalInfeasibilities_ > 100.0 * primalTolerance_) {
     805#if CBC_USEFUL_PRINTING > 0
     806        printf("internal check on primals failed %d %g\n",
     807          numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
    807808#endif
    808809      } else {
    809         sumPrimalInfeasibilities_=0.0;
    810         numberPrimalInfeasibilities_=0;
    811       }
    812       problemStatus_=0;
     810        sumPrimalInfeasibilities_ = 0.0;
     811        numberPrimalInfeasibilities_ = 0;
     812      }
     813      problemStatus_ = 0;
    813814    }
    814815#endif
    815816    //ClpModel::stopPermanentArrays();
    816     this->setSpecialOptions(this->specialOptions()&~65536);
     817    this->setSpecialOptions(this->specialOptions() & ~65536);
    817818#if 0
    818819    this->writeBasis("a.bas",true);
     
    829830    if (!keepAbc) {
    830831      delete abcModel2;
    831       abcModel2=NULL;
     832      abcModel2 = NULL;
    832833    }
    833834  }
     
    845846    64 - do not use sprint even if problem looks good
    846847 */
    847 int
    848 ClpSimplex::initialSolve(ClpSolve & options)
     848int ClpSimplex::initialSolve(ClpSolve &options)
    849849{
    850      ClpSolve::SolveType method = options.getSolveType();
    851      //ClpSolve::SolveType originalMethod=method;
    852      ClpSolve::PresolveType presolve = options.getPresolveType();
    853      int saveMaxIterations = maximumIterations();
    854      int finalStatus = -1;
    855      int numberIterations = 0;
    856      double time1 = CoinCpuTime();
    857      double timeX = time1;
    858      double time2=0.0;
    859      ClpMatrixBase * saveMatrix = NULL;
    860      ClpObjective * savedObjective = NULL;
    861      int idiotOptions=0;
    862      if(options.getSpecialOption(6))
    863        idiotOptions=options.getExtraInfo(6)*32768;
     850  ClpSolve::SolveType method = options.getSolveType();
     851  //ClpSolve::SolveType originalMethod=method;
     852  ClpSolve::PresolveType presolve = options.getPresolveType();
     853  int saveMaxIterations = maximumIterations();
     854  int finalStatus = -1;
     855  int numberIterations = 0;
     856  double time1 = CoinCpuTime();
     857  double timeX = time1;
     858  double time2 = 0.0;
     859  ClpMatrixBase *saveMatrix = NULL;
     860  ClpObjective *savedObjective = NULL;
     861  int idiotOptions = 0;
     862  if (options.getSpecialOption(6))
     863    idiotOptions = options.getExtraInfo(6) * 32768;
    864864#ifdef CLP_USEFUL_PRINTOUT
    865      debugInt[0]=numberRows();
    866      debugInt[1]=numberColumns();
    867      debugInt[2]=matrix()->getNumElements();
    868 #endif
    869      if (!objective_ || !matrix_) {
    870           // totally empty
    871           handler_->message(CLP_EMPTY_PROBLEM, messages_)
    872                     << 0
    873                     << 0
    874                     << 0
    875                     << CoinMessageEol;
    876           return -1;
    877      } else if (!numberRows_ || !numberColumns_ || !getNumElements()) {
    878           presolve = ClpSolve::presolveOff;
    879      }
    880      if (objective_->type() >= 2 && optimizationDirection_ == 0) {
    881           // pretend linear
    882           savedObjective = objective_;
    883           // make up objective
    884           double * obj = new double[numberColumns_];
    885           for (int i = 0; i < numberColumns_; i++) {
    886                double l = fabs(columnLower_[i]);
    887                double u = fabs(columnUpper_[i]);
    888                obj[i] = 0.0;
    889                if (CoinMin(l, u) < 1.0e20) {
    890                     if (l < u)
    891                          obj[i] = 1.0 + randomNumberGenerator_.randomDouble() * 1.0e-2;
    892                     else
    893                          obj[i] = -1.0 - randomNumberGenerator_.randomDouble() * 1.0e-2;
    894                }
    895           }
    896           objective_ = new ClpLinearObjective(obj, numberColumns_);
    897           delete [] obj;
    898      }
    899      ClpSimplex * model2 = this;
    900      bool interrupt = (options.getSpecialOption(2) == 0);
    901      CoinSighandler_t saveSignal = static_cast<CoinSighandler_t> (0);
    902      if (interrupt) {
    903           currentModel = model2;
    904           // register signal handler
    905           saveSignal = signal(SIGINT, signal_handler);
    906      }
    907      // If no status array - set up basis
    908      if (!status_)
    909           allSlackBasis();
    910      ClpPresolve * pinfo = new ClpPresolve();
    911      pinfo->setSubstitution(options.substitution());
    912      int presolveOptions = options.presolveActions();
    913      bool presolveToFile = (presolveOptions & 0x40000000) != 0;
    914      presolveOptions &= ~0x40000000;
    915      if ((presolveOptions & 0xffffff) != 0)
    916           pinfo->setPresolveActions(presolveOptions);
    917      // switch off singletons to slacks
    918      //pinfo->setDoSingletonColumn(false); // done by bits
    919      int printOptions = options.getSpecialOption(5);
    920      if ((printOptions & 1) != 0)
    921           pinfo->statistics();
    922      double timePresolve = 0.0;
    923      double timeIdiot = 0.0;
    924      double timeCore = 0.0;
    925      eventHandler()->event(ClpEventHandler::presolveStart);
    926      int savePerturbation = perturbation_;
    927      int saveScaling = scalingFlag_;
     865  debugInt[0] = numberRows();
     866  debugInt[1] = numberColumns();
     867  debugInt[2] = matrix()->getNumElements();
     868#endif
     869  if (!objective_ || !matrix_) {
     870    // totally empty
     871    handler_->message(CLP_EMPTY_PROBLEM, messages_)
     872      << 0
     873      << 0
     874      << 0
     875      << CoinMessageEol;
     876    return -1;
     877  } else if (!numberRows_ || !numberColumns_ || !getNumElements()) {
     878    presolve = ClpSolve::presolveOff;
     879  }
     880  if (objective_->type() >= 2 && optimizationDirection_ == 0) {
     881    // pretend linear
     882    savedObjective = objective_;
     883    // make up objective
     884    double *obj = new double[numberColumns_];
     885    for (int i = 0; i < numberColumns_; i++) {
     886      double l = fabs(columnLower_[i]);
     887      double u = fabs(columnUpper_[i]);
     888      obj[i] = 0.0;
     889      if (CoinMin(l, u) < 1.0e20) {
     890        if (l < u)
     891          obj[i] = 1.0 + randomNumberGenerator_.randomDouble() * 1.0e-2;
     892        else
     893          obj[i] = -1.0 - randomNumberGenerator_.randomDouble() * 1.0e-2;
     894      }
     895    }
     896    objective_ = new ClpLinearObjective(obj, numberColumns_);
     897    delete[] obj;
     898  }
     899  ClpSimplex *model2 = this;
     900  bool interrupt = (options.getSpecialOption(2) == 0);
     901  CoinSighandler_t saveSignal = static_cast< CoinSighandler_t >(0);
     902  if (interrupt) {
     903    currentModel = model2;
     904    // register signal handler
     905    saveSignal = signal(SIGINT, signal_handler);
     906  }
     907  // If no status array - set up basis
     908  if (!status_)
     909    allSlackBasis();
     910  ClpPresolve *pinfo = new ClpPresolve();
     911  pinfo->setSubstitution(options.substitution());
     912  int presolveOptions = options.presolveActions();
     913  bool presolveToFile = (presolveOptions & 0x40000000) != 0;
     914  presolveOptions &= ~0x40000000;
     915  if ((presolveOptions & 0xffffff) != 0)
     916    pinfo->setPresolveActions(presolveOptions);
     917  // switch off singletons to slacks
     918  //pinfo->setDoSingletonColumn(false); // done by bits
     919  int printOptions = options.getSpecialOption(5);
     920  if ((printOptions & 1) != 0)
     921    pinfo->statistics();
     922  double timePresolve = 0.0;
     923  double timeIdiot = 0.0;
     924  double timeCore = 0.0;
     925  eventHandler()->event(ClpEventHandler::presolveStart);
     926  int savePerturbation = perturbation_;
     927  int saveScaling = scalingFlag_;
    928928#ifndef SLIM_CLP
    929929#ifndef NO_RTTI
    930      if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
    931           // network - switch off stuff
    932           presolve = ClpSolve::presolveOff;
    933      }
     930  if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) {
     931    // network - switch off stuff
     932    presolve = ClpSolve::presolveOff;
     933  }
    934934#else
    935      if (matrix_->type() == 11) {
    936           // network - switch off stuff
    937           presolve = ClpSolve::presolveOff;
    938      }
    939 #endif
    940 #endif
    941 #ifndef CLPSOLVE_ACTIONS 
     935  if (matrix_->type() == 11) {
     936    // network - switch off stuff
     937    presolve = ClpSolve::presolveOff;
     938  }
     939#endif
     940#endif
     941#ifndef CLPSOLVE_ACTIONS
    942942#define CLPSOLVE_ACTIONS 2
    943943#endif
    944 #if CLPSOLVE_ACTIONS 
    945      bool wasAutomatic = (method==ClpSolve::automatic);
    946 #endif
    947      if (presolve != ClpSolve::presolveOff) {
    948           bool costedSlacks = false;
    949 #if CLP_INHERIT_MODE>1
    950           int numberPasses = 20;
     944#if CLPSOLVE_ACTIONS
     945  bool wasAutomatic = (method == ClpSolve::automatic);
     946#endif
     947  if (presolve != ClpSolve::presolveOff) {
     948    bool costedSlacks = false;
     949#if CLP_INHERIT_MODE > 1
     950    int numberPasses = 20;
    951951#else
    952           int numberPasses = 5;
    953 #endif
    954           if (presolve == ClpSolve::presolveNumber) {
    955                numberPasses = options.getPresolvePasses();
    956                presolve = ClpSolve::presolveOn;
    957           } else if (presolve == ClpSolve::presolveNumberCost) {
    958                numberPasses = options.getPresolvePasses();
    959                presolve = ClpSolve::presolveOn;
    960                costedSlacks = true;
    961                // switch on singletons to slacks
    962                pinfo->setDoSingletonColumn(true);
    963                // gub stuff for testing
    964                //pinfo->setDoGubrow(true);
    965           }
     952    int numberPasses = 5;
     953#endif
     954    if (presolve == ClpSolve::presolveNumber) {
     955      numberPasses = options.getPresolvePasses();
     956      presolve = ClpSolve::presolveOn;
     957    } else if (presolve == ClpSolve::presolveNumberCost) {
     958      numberPasses = options.getPresolvePasses();
     959      presolve = ClpSolve::presolveOn;
     960      costedSlacks = true;
     961      // switch on singletons to slacks
     962      pinfo->setDoSingletonColumn(true);
     963      // gub stuff for testing
     964      //pinfo->setDoGubrow(true);
     965    }
    966966#ifndef CLP_NO_STD
    967           if (presolveToFile) {
    968                // PreSolve to file - not fully tested
    969                printf("Presolving to file - presolve.save\n");
    970                pinfo->presolvedModelToFile(*this, "presolve.save", dblParam_[ClpPresolveTolerance],
    971                                           false, numberPasses);
    972                model2 = this;
    973           } else {
    974 #endif
    975                model2 = pinfo->presolvedModel(*this, dblParam_[ClpPresolveTolerance],
    976                                              false, numberPasses, true, costedSlacks);
     967    if (presolveToFile) {
     968      // PreSolve to file - not fully tested
     969      printf("Presolving to file - presolve.save\n");
     970      pinfo->presolvedModelToFile(*this, "presolve.save", dblParam_[ClpPresolveTolerance],
     971        false, numberPasses);
     972      model2 = this;
     973    } else {
     974#endif
     975      model2 = pinfo->presolvedModel(*this, dblParam_[ClpPresolveTolerance],
     976        false, numberPasses, true, costedSlacks);
    977977#ifndef CLP_NO_STD
    978           }
    979 #endif
    980           time2 = CoinCpuTime();
    981           timePresolve = time2 - timeX;
    982           handler_->message(CLP_INTERVAL_TIMING, messages_)
    983                     << "Presolve" << timePresolve << time2 - time1
    984                     << CoinMessageEol;
    985           timeX = time2;
    986           if (!model2) {
    987                handler_->message(CLP_INFEASIBLE, messages_)
    988                          << CoinMessageEol;
    989                model2 = this;
    990                eventHandler()->event(ClpEventHandler::presolveInfeasible);
    991                problemStatus_ = pinfo->presolveStatus();
    992                if (options.infeasibleReturn() || (moreSpecialOptions_ & 1) != 0) {
    993                 delete pinfo;
    994                     return -1;
    995                }
    996                presolve = ClpSolve::presolveOff;
    997           } else {
     978    }
     979#endif
     980    time2 = CoinCpuTime();
     981    timePresolve = time2 - timeX;
     982    handler_->message(CLP_INTERVAL_TIMING, messages_)
     983      << "Presolve" << timePresolve << time2 - time1
     984      << CoinMessageEol;
     985    timeX = time2;
     986    if (!model2) {
     987      handler_->message(CLP_INFEASIBLE, messages_)
     988        << CoinMessageEol;
     989      model2 = this;
     990      eventHandler()->event(ClpEventHandler::presolveInfeasible);
     991      problemStatus_ = pinfo->presolveStatus();
     992      if (options.infeasibleReturn() || (moreSpecialOptions_ & 1) != 0) {
     993        delete pinfo;
     994        return -1;
     995      }
     996      presolve = ClpSolve::presolveOff;
     997    } else {
    998998#if 0 //def ABC_INHERIT
    999999            {
     
    10041004            }
    10051005#else
    1006             //ClpModel::stopPermanentArrays();
    1007             //setSpecialOptions(specialOptions()&~65536);
    1008             // try setting tolerances up
     1006      //ClpModel::stopPermanentArrays();
     1007      //setSpecialOptions(specialOptions()&~65536);
     1008      // try setting tolerances up
    10091009#if CLPSOLVE_ACTIONS
    1010             bool changeTolerances=wasAutomatic;
    1011 #if CLPSOLVE_ACTIONS>1
    1012             changeTolerances=true;
    1013 #endif
    1014             if (changeTolerances && model2 != this) {
     1010      bool changeTolerances = wasAutomatic;
     1011#if CLPSOLVE_ACTIONS > 1
     1012      changeTolerances = true;
     1013#endif
     1014      if (changeTolerances && model2 != this) {
    10151015#define CLP_NEW_TOLERANCE 1.0e-6
    1016               if (model2->primalTolerance()==1.0e-7&&model2->dualTolerance()==1.0e-7) {
    1017                 model2->setPrimalTolerance(CLP_NEW_TOLERANCE);
    1018                 model2->setDualTolerance(CLP_NEW_TOLERANCE);
    1019               }
    1020             }
    1021 #endif
    1022 #endif
    1023               model2->eventHandler()->setSimplex(model2);
    1024               int rcode=model2->eventHandler()->event(ClpEventHandler::presolveSize);
    1025               // see if too big or small
    1026               if (rcode==2) {
    1027                   delete model2;
    1028                  delete pinfo;
    1029                   return -2;
    1030               } else if (rcode==3) {
    1031                   delete model2;
    1032                  delete pinfo;
    1033                   return -3;
    1034               }
     1016        if (model2->primalTolerance() == 1.0e-7 && model2->dualTolerance() == 1.0e-7) {
     1017          model2->setPrimalTolerance(CLP_NEW_TOLERANCE);
     1018          model2->setDualTolerance(CLP_NEW_TOLERANCE);
     1019        }
     1020      }
     1021#endif
     1022#endif
     1023      model2->eventHandler()->setSimplex(model2);
     1024      int rcode = model2->eventHandler()->event(ClpEventHandler::presolveSize);
     1025      // see if too big or small
     1026      if (rcode == 2) {
     1027        delete model2;
     1028        delete pinfo;
     1029        return -2;
     1030      } else if (rcode == 3) {
     1031        delete model2;
     1032        delete pinfo;
     1033        return -3;
     1034      }
     1035    }
     1036    model2->setMoreSpecialOptions(model2->moreSpecialOptions() & (~1024));
     1037    model2->eventHandler()->setSimplex(model2);
     1038    // We may be better off using original (but if dual leave because of bounds)
     1039    if (presolve != ClpSolve::presolveOff && numberRows_ < 1.01 * model2->numberRows_ && numberColumns_ < 1.01 * model2->numberColumns_
     1040      && model2 != this) {
     1041      if (method != ClpSolve::useDual || (numberRows_ == model2->numberRows_ && numberColumns_ == model2->numberColumns_)) {
     1042        delete model2;
     1043        model2 = this;
     1044        presolve = ClpSolve::presolveOff;
     1045      }
     1046    }
     1047  }
     1048#ifdef CLP_USEFUL_PRINTOUT
     1049  debugInt[3] = model2->numberRows();
     1050  debugInt[4] = model2->numberColumns();
     1051  debugInt[5] = model2->matrix()->getNumElements();
     1052  // analyze
     1053  {
     1054    double time1 = CoinCpuTime();
     1055    int numberColumns = model2->numberColumns();
     1056    const double *columnLower = model2->columnLower();
     1057    const double *columnUpper = model2->columnUpper();
     1058    int numberRows = model2->numberRows();
     1059    const double *rowLower = model2->rowLower();
     1060    const double *rowUpper = model2->rowUpper();
     1061    const double *objective = model2->objective();
     1062    CoinPackedMatrix *matrix = model2->matrix();
     1063    CoinBigIndex numberElements = matrix->getNumElements();
     1064    const int *columnLength = matrix->getVectorLengths();
     1065    //const CoinBigIndex * columnStart = matrix->getVectorStarts();
     1066    const double *elementByColumn = matrix->getElements();
     1067    const int *row = matrix->getIndices();
     1068    int *rowCount = new int[numberRows];
     1069    memset(rowCount, 0, numberRows * sizeof(int));
     1070    int n = CoinMax(2 * numberRows, numberElements);
     1071    n = CoinMax(2 * numberColumns, n);
     1072    double *check = new double[n];
     1073    memcpy(check, elementByColumn, numberElements * sizeof(double));
     1074    for (int i = 0; i < numberElements; i++) {
     1075      check[i] = fabs(check[i]);
     1076      rowCount[row[i]]++;
     1077    }
     1078    int largestIndex = 0;
     1079    for (int i = 0; i < numberColumns; i++) {
     1080      largestIndex = CoinMax(largestIndex, columnLength[i]);
     1081    }
     1082    debugInt[12] = largestIndex;
     1083    largestIndex = 0;
     1084    for (int i = 0; i < numberRows; i++) {
     1085      largestIndex = CoinMax(largestIndex, rowCount[i]);
     1086    }
     1087    n = numberElements;
     1088    delete[] rowCount;
     1089    debugInt[11] = largestIndex;
     1090    std::sort(check, check + n);
     1091    debugDouble[4] = check[0];
     1092    debugDouble[5] = check[n - 1];
     1093    int nAtOne = 0;
     1094    int nDifferent = 0;
     1095    double last = -COIN_DBL_MAX;
     1096    for (int i = 0; i < n; i++) {
     1097      if (fabs(last - check[i]) > 1.0e-12) {
     1098        nDifferent++;
     1099        last = check[i];
     1100      }
     1101      if (check[i] == 1.0)
     1102        nAtOne++;
     1103    }
     1104    debugInt[10] = nDifferent;
     1105    debugInt[15] = nAtOne;
     1106    int nInf = 0;
     1107    int nZero = 0;
     1108    n = 0;
     1109    for (int i = 0; i < numberRows; i++) {
     1110      double value = fabs(rowLower[i]);
     1111      if (!value)
     1112        nZero++;
     1113      else if (value != COIN_DBL_MAX)
     1114        check[n++] = value;
     1115      else
     1116        nInf++;
     1117    }
     1118    for (int i = 0; i < numberRows; i++) {
     1119      double value = fabs(rowUpper[i]);
     1120      if (!value)
     1121        nZero++;
     1122      else if (value != COIN_DBL_MAX)
     1123        check[n++] = value;
     1124      else
     1125        nInf++;
     1126    }
     1127    debugInt[16] = nInf;
     1128    debugInt[20] = nZero;
     1129    if (n) {
     1130      std::sort(check, check + n);
     1131      debugDouble[0] = check[0];
     1132      debugDouble[1] = check[n - 1];
     1133    } else {
     1134      debugDouble[0] = 0.0;
     1135      debugDouble[1] = 0.0;
     1136    }
     1137    nAtOne = 0;
     1138    nDifferent = 0;
     1139    last = -COIN_DBL_MAX;
     1140    for (int i = 0; i < n; i++) {
     1141      if (fabs(last - check[i]) > 1.0e-12) {
     1142        nDifferent++;
     1143        last = check[i];
     1144      }
     1145      if (check[i] == 1.0)
     1146        nAtOne++;
     1147    }
     1148    debugInt[8] = nDifferent;
     1149    debugInt[13] = nAtOne;
     1150    nZero = 0;
     1151    n = 0;
     1152    for (int i = 0; i < numberColumns; i++) {
     1153      double value = fabs(objective[i]);
     1154      if (value)
     1155        check[n++] = value;
     1156      else
     1157        nZero++;
     1158    }
     1159    debugInt[21] = nZero;
     1160    if (n) {
     1161      std::sort(check, check + n);
     1162      debugDouble[2] = check[0];
     1163      debugDouble[3] = check[n - 1];
     1164    } else {
     1165      debugDouble[2] = 0.0;
     1166      debugDouble[3] = 0.0;
     1167    }
     1168    nAtOne = 0;
     1169    nDifferent = 0;
     1170    last = -COIN_DBL_MAX;
     1171    for (int i = 0; i < n; i++) {
     1172      if (fabs(last - check[i]) > 1.0e-12) {
     1173        nDifferent++;
     1174        last = check[i];
     1175      }
     1176      if (check[i] == 1.0)
     1177        nAtOne++;
     1178    }
     1179    debugInt[9] = nDifferent;
     1180    debugInt[14] = nAtOne;
     1181    nInf = 0;
     1182    nZero = 0;
     1183    n = 0;
     1184    for (int i = 0; i < numberColumns; i++) {
     1185      double value = fabs(columnLower[i]);
     1186      if (!value)
     1187        nZero++;
     1188      else if (value != COIN_DBL_MAX)
     1189        check[n++] = value;
     1190      else
     1191        nInf++;
     1192    }
     1193    for (int i = 0; i < numberColumns; i++) {
     1194      double value = fabs(columnUpper[i]);
     1195      if (!value)
     1196        nZero++;
     1197      else if (value != COIN_DBL_MAX)
     1198        check[n++] = value;
     1199      else
     1200        nInf++;
     1201    }
     1202    debugInt[17] = nInf;
     1203    double smallestColBound;
     1204    double largestColBound;
     1205    if (n) {
     1206      std::sort(check, check + n);
     1207      smallestColBound = check[0];
     1208      largestColBound = check[n - 1];
     1209    } else {
     1210      smallestColBound = 0.0;
     1211      largestColBound = 0.0;
     1212    }
     1213    nAtOne = 0;
     1214    nDifferent = 0;
     1215    last = -COIN_DBL_MAX;
     1216    for (int i = 0; i < n; i++) {
     1217      if (fabs(last - check[i]) > 1.0e-12) {
     1218        nDifferent++;
     1219        last = check[i];
     1220      }
     1221      if (check[i] == 1.0)
     1222        nAtOne++;
     1223    }
     1224    //debugInt[8]=nDifferent;
     1225    //debugInt[13]=nAtOne;
     1226    printf("BENCHMARK_STATS rhs %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
     1227      debugInt[8], debugDouble[0], debugDouble[1], debugInt[13], debugInt[16], debugInt[20]);
     1228    printf("BENCHMARK_STATS col %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
     1229      nDifferent, smallestColBound, largestColBound, nAtOne, nInf, nZero);
     1230    printf("BENCHMARK_STATS els %d different - %g -> %g (%d at one) - longest r,c %d,%d\n",
     1231      debugInt[10], debugDouble[4], debugDouble[5], debugInt[15],
     1232      debugInt[11], debugInt[12]);
     1233    printf("BENCHMARK_STATS obj %d different - %g -> %g (%d at one, %d zero) - time %g\n",
     1234      debugInt[9], debugDouble[2], debugDouble[3], debugInt[14], debugInt[21],
     1235      CoinCpuTime() - time1);
     1236    delete[] check;
     1237  }
     1238#endif
     1239  if (interrupt)
     1240    currentModel = model2;
     1241  int saveMoreOptions = moreSpecialOptions_;
     1242  // For below >0 overrides
     1243  // 0 means no, -1 means maybe
     1244  int doIdiot = 0;
     1245  int doCrash = 0;
     1246  int doSprint = 0;
     1247  int doSlp = 0;
     1248  int primalStartup = 1;
     1249  model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve);
     1250#if CLP_POOL_MATRIX
     1251  if (vectorMode() >= 10) {
     1252    ClpPoolMatrix *poolMatrix = new ClpPoolMatrix(*model2->matrix());
     1253    char output[80];
     1254    int numberDifferent = poolMatrix->getNumDifferentElements();
     1255    if (numberDifferent > 0) {
     1256      sprintf(output, "Pool matrix has %d different values",
     1257        numberDifferent);
     1258      model2->replaceMatrix(poolMatrix, true);
     1259    } else {
     1260      delete poolMatrix;
     1261      sprintf(output, "Pool matrix has more than %d different values - no good",
     1262        -numberDifferent);
     1263    }
     1264    handler_->message(CLP_GENERAL, messages_) << output
     1265                                              << CoinMessageEol;
     1266  }
     1267#endif
     1268  int tryItSave = 0;
     1269#if CLPSOLVE_ACTIONS
     1270  if (method == ClpSolve::automatic)
     1271    model2->moreSpecialOptions_ |= 8192; // stop switch over
     1272#endif
     1273  // switch to primal from automatic if just one cost entry
     1274  if (method == ClpSolve::automatic && model2->numberColumns() > 5000
     1275#ifndef CLPSOLVE_ACTIONS
     1276    && (specialOptions_ & 1024) != 0
     1277#endif
     1278  ) {
     1279    // look at original model for objective
     1280    int numberColumns = model2->numberColumns();
     1281    int numberRows = model2->numberRows();
     1282    int numberColumnsOrig = this->numberColumns();
     1283    const double *obj = this->objective();
     1284    int nNon = 0;
     1285    bool allOnes = true;
     1286    for (int i = 0; i < numberColumnsOrig; i++) {
     1287      if (obj[i]) {
     1288        nNon++;
     1289        if (fabs(obj[i]) != 1.0)
     1290          allOnes = false;
     1291      }
     1292    }
     1293    if (nNon <= 1 || allOnes || (options.getExtraInfo(1) > 0 && options.getSpecialOption(1) == 2)) {
     1294#ifdef COIN_DEVELOP
     1295      printf("Forcing primal\n");
     1296#endif
     1297      method = ClpSolve::usePrimal;
     1298#ifndef CLPSOLVE_ACTIONS
     1299      tryItSave = (numberRows > 200 && numberColumns > 2000 && (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0)) ? 3 : 0;
     1300#else
     1301      if (numberRows > 200 && numberColumns > 2000 && numberColumns > 2 * numberRows) {
     1302        tryItSave = 3;
     1303      } else {
     1304        // If rhs also rubbish then maybe
     1305        int numberRowsOrig = this->numberRows();
     1306        const double *rowLower = this->rowLower();
     1307        const double *rowUpper = this->rowUpper();
     1308        double last = COIN_DBL_MAX;
     1309        int nDifferent = 0;
     1310        for (int i = 0; i < numberRowsOrig; i++) {
     1311          double value = fabs(rowLower[i]);
     1312          if (value && value != COIN_DBL_MAX) {
     1313            if (value != last) {
     1314              nDifferent++;
     1315              last = value;
     1316            }
    10351317          }
    1036           model2->setMoreSpecialOptions(model2->moreSpecialOptions()&(~1024));
    1037           model2->eventHandler()->setSimplex(model2);
    1038           // We may be better off using original (but if dual leave because of bounds)
    1039           if (presolve != ClpSolve::presolveOff &&
    1040                     numberRows_ < 1.01 * model2->numberRows_ && numberColumns_ < 1.01 * model2->numberColumns_
    1041                     && model2 != this) {
    1042                if(method != ClpSolve::useDual ||
    1043                          (numberRows_ == model2->numberRows_ && numberColumns_ == model2->numberColumns_)) {
    1044                     delete model2;
    1045                     model2 = this;
    1046                     presolve = ClpSolve::presolveOff;
    1047                }
     1318          value = fabs(rowUpper[i]);
     1319          if (value && value != COIN_DBL_MAX) {
     1320            if (value != last) {
     1321              nDifferent++;
     1322              last = value;
     1323            }
    10481324          }
    1049      }
    1050 #ifdef CLP_USEFUL_PRINTOUT
    1051      debugInt[3]=model2->numberRows();
    1052      debugInt[4]=model2->numberColumns();
    1053      debugInt[5]=model2->matrix()->getNumElements();
    1054      // analyze
    1055      {
    1056        double time1 =CoinCpuTime();
    1057        int numberColumns = model2->numberColumns();
    1058        const double * columnLower = model2->columnLower();
    1059        const double * columnUpper = model2->columnUpper();
    1060        int numberRows = model2->numberRows();
    1061        const double * rowLower = model2->rowLower();
    1062        const double * rowUpper = model2->rowUpper();
    1063        const double * objective = model2->objective();
    1064        CoinPackedMatrix * matrix = model2->matrix();
    1065        CoinBigIndex numberElements = matrix->getNumElements();
    1066        const int * columnLength = matrix->getVectorLengths();
    1067        //const CoinBigIndex * columnStart = matrix->getVectorStarts();
    1068        const double * elementByColumn = matrix->getElements();
    1069        const int * row = matrix->getIndices();
    1070        int * rowCount = new int [numberRows];
    1071        memset(rowCount,0,numberRows*sizeof(int));
    1072        int n=CoinMax (2*numberRows,numberElements);
    1073        n= CoinMax(2*numberColumns,n);
    1074        double * check = new double[n];
    1075        memcpy(check,elementByColumn,numberElements*sizeof(double));
    1076        for (int i=0;i<numberElements;i++) {
    1077          check[i]=fabs(check[i]);
    1078          rowCount[row[i]]++;
    1079        }
    1080        int largestIndex=0;
    1081        for (int i=0;i<numberColumns;i++) {
    1082          largestIndex=CoinMax(largestIndex,columnLength[i]);
    1083        }
    1084        debugInt[12]=largestIndex;
    1085        largestIndex=0;
    1086        for (int i=0;i<numberRows;i++) {
    1087          largestIndex=CoinMax(largestIndex,rowCount[i]);
    1088        }
    1089        n=numberElements;
    1090        delete [] rowCount;
    1091        debugInt[11]=largestIndex;
    1092        std::sort(check,check+n);
    1093        debugDouble[4]=check[0];
    1094        debugDouble[5]=check[n-1];
    1095        int nAtOne=0;
    1096        int nDifferent=0;
    1097        double last=-COIN_DBL_MAX;
    1098        for (int i=0;i<n;i++) {
    1099          if (fabs(last-check[i])>1.0e-12) {
    1100            nDifferent++;
    1101            last=check[i];
    1102          }
    1103          if (check[i]==1.0)
    1104            nAtOne++;
    1105        }
    1106        debugInt[10]=nDifferent;
    1107        debugInt[15]=nAtOne;
    1108        int nInf=0;
    1109        int nZero=0;
    1110        n=0;
    1111        for (int i=0;i<numberRows;i++) {
    1112          double value=fabs(rowLower[i]);
    1113          if (!value)
    1114            nZero++;
    1115          else if (value!=COIN_DBL_MAX)
    1116            check[n++]=value;
    1117          else
    1118            nInf++;
    1119        }
    1120        for (int i=0;i<numberRows;i++) {
    1121          double value=fabs(rowUpper[i]);
    1122          if (!value)
    1123            nZero++;
    1124          else if (value!=COIN_DBL_MAX)
    1125            check[n++]=value;
    1126          else
    1127            nInf++;
    1128        }
    1129        debugInt[16]=nInf;
    1130        debugInt[20]=nZero;
    1131        if (n) {
    1132          std::sort(check,check+n);
    1133          debugDouble[0]=check[0];
    1134          debugDouble[1]=check[n-1];
    1135        } else {
    1136          debugDouble[0]=0.0;
    1137          debugDouble[1]=0.0;
    1138        }
    1139        nAtOne=0;
    1140        nDifferent=0;
    1141        last=-COIN_DBL_MAX;
    1142        for (int i=0;i<n;i++) {
    1143          if (fabs(last-check[i])>1.0e-12) {
    1144            nDifferent++;
    1145            last=check[i];
    1146          }
    1147          if (check[i]==1.0)
    1148            nAtOne++;
    1149        }
    1150        debugInt[8]=nDifferent;
    1151        debugInt[13]=nAtOne;
    1152        nZero=0;
    1153        n=0;
    1154        for (int i=0;i<numberColumns;i++) {
    1155          double value=fabs(objective[i]);
    1156          if (value)
    1157            check[n++]=value;
    1158          else
    1159            nZero++;
    1160        }
    1161        debugInt[21]=nZero;
    1162        if (n) {
    1163          std::sort(check,check+n);
    1164          debugDouble[2]=check[0];
    1165          debugDouble[3]=check[n-1];
    1166        } else {
    1167          debugDouble[2]=0.0;
    1168          debugDouble[3]=0.0;
    1169        }
    1170        nAtOne=0;
    1171        nDifferent=0;
    1172        last=-COIN_DBL_MAX;
    1173        for (int i=0;i<n;i++) {
    1174          if (fabs(last-check[i])>1.0e-12) {
    1175            nDifferent++;
    1176            last=check[i];
    1177          }
    1178          if (check[i]==1.0)
    1179            nAtOne++;
    1180        }
    1181        debugInt[9]=nDifferent;
    1182        debugInt[14]=nAtOne;
    1183        nInf=0;
    1184        nZero=0;
    1185        n=0;
    1186        for (int i=0;i<numberColumns;i++) {
    1187          double value=fabs(columnLower[i]);
    1188          if (!value)
    1189            nZero++;
    1190          else if (value!=COIN_DBL_MAX)
    1191            check[n++]=value;
    1192          else
    1193            nInf++;
    1194        }
    1195        for (int i=0;i<numberColumns;i++) {
    1196          double value=fabs(columnUpper[i]);
    1197          if (!value)
    1198            nZero++;
    1199          else if (value!=COIN_DBL_MAX)
    1200            check[n++]=value;
    1201          else
    1202            nInf++;
    1203        }
    1204        debugInt[17]=nInf;
    1205        double smallestColBound;
    1206        double largestColBound;
    1207        if (n) {
    1208          std::sort(check,check+n);
    1209          smallestColBound=check[0];
    1210          largestColBound=check[n-1];
    1211        } else {
    1212          smallestColBound=0.0;
    1213          largestColBound=0.0;
    1214        }
    1215        nAtOne=0;
    1216        nDifferent=0;
    1217        last=-COIN_DBL_MAX;
    1218        for (int i=0;i<n;i++) {
    1219          if (fabs(last-check[i])>1.0e-12) {
    1220            nDifferent++;
    1221            last=check[i];
    1222          }
    1223          if (check[i]==1.0)
    1224            nAtOne++;
    1225        }
    1226        //debugInt[8]=nDifferent;
    1227        //debugInt[13]=nAtOne;
    1228        printf("BENCHMARK_STATS rhs %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
    1229               debugInt[8],debugDouble[0],debugDouble[1],debugInt[13],debugInt[16],debugInt[20]);
    1230        printf("BENCHMARK_STATS col %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
    1231               nDifferent,smallestColBound,largestColBound,nAtOne,nInf,nZero);
    1232        printf("BENCHMARK_STATS els %d different - %g -> %g (%d at one) - longest r,c %d,%d\n",
    1233               debugInt[10],debugDouble[4],debugDouble[5],debugInt[15],
    1234               debugInt[11],debugInt[12]);
    1235        printf("BENCHMARK_STATS obj %d different - %g -> %g (%d at one, %d zero) - time %g\n",
    1236               debugInt[9],debugDouble[2],debugDouble[3],debugInt[14],debugInt[21],
    1237               CoinCpuTime()-time1);
    1238        delete [] check;
    1239      }
    1240 #endif
    1241      if (interrupt)
    1242           currentModel = model2;
    1243      int saveMoreOptions = moreSpecialOptions_;
    1244      // For below >0 overrides
    1245      // 0 means no, -1 means maybe
    1246      int doIdiot = 0;
    1247      int doCrash = 0;
    1248      int doSprint = 0;
    1249      int doSlp = 0;
    1250      int primalStartup = 1;
    1251      model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve);
    1252 #if CLP_POOL_MATRIX
    1253      if (vectorMode()>=10) {
    1254        ClpPoolMatrix * poolMatrix = new ClpPoolMatrix(*model2->matrix());
    1255        char output[80];
    1256        int numberDifferent=poolMatrix->getNumDifferentElements();
    1257        if (numberDifferent>0) {
    1258          sprintf(output,"Pool matrix has %d different values",
    1259                  numberDifferent);
    1260          model2->replaceMatrix(poolMatrix,true);
    1261        } else {
    1262          delete poolMatrix;
    1263          sprintf(output,"Pool matrix has more than %d different values - no good",
    1264                  -numberDifferent);
    1265        }
    1266        handler_->message(CLP_GENERAL, messages_) << output
    1267          << CoinMessageEol;
    1268      }
    1269 #endif
    1270      int tryItSave = 0;
    1271 #if CLPSOLVE_ACTIONS
    1272      if (method == ClpSolve::automatic )
    1273        model2->moreSpecialOptions_ |= 8192; // stop switch over
    1274 #endif
    1275      // switch to primal from automatic if just one cost entry
    1276      if (method == ClpSolve::automatic && model2->numberColumns() > 5000
    1277 #ifndef CLPSOLVE_ACTIONS
    1278          && (specialOptions_ & 1024) != 0
    1279 #endif
    1280          ) {
    1281           // look at original model for objective
    1282           int numberColumns = model2->numberColumns();
    1283           int numberRows = model2->numberRows();
    1284           int numberColumnsOrig = this->numberColumns();
    1285           const double * obj = this->objective();
    1286           int nNon = 0;
    1287           bool allOnes=true;
    1288           for (int i = 0; i < numberColumnsOrig; i++) {
    1289             if (obj[i]) {
    1290               nNon++;
    1291               if (fabs(obj[i])!=1.0)
    1292                 allOnes=false;
    1293             }
    1294           }
    1295           if (nNon <= 1 || allOnes ||
    1296               (options.getExtraInfo(1) > 0 && options.getSpecialOption(1)==2)) {
    1297 #ifdef COIN_DEVELOP
    1298                printf("Forcing primal\n");
    1299 #endif
    1300                method = ClpSolve::usePrimal;
    1301 #ifndef CLPSOLVE_ACTIONS
    1302                tryItSave = (numberRows > 200 && numberColumns > 2000 &&
    1303                             (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0)) ? 3 : 0;
     1325        }
     1326        if (nDifferent < 2)
     1327          tryItSave = 1;
     1328      }
     1329#endif
     1330    }
     1331  }
     1332  if (method != ClpSolve::useDual && method != ClpSolve::useBarrier
     1333    && method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe
     1334    && method != ClpSolve::useBarrierNoCross) {
     1335    switch (options.getSpecialOption(1)) {
     1336    case 0:
     1337      doIdiot = -1;
     1338      doCrash = -1;
     1339      doSprint = -1;
     1340      break;
     1341    case 1:
     1342      doIdiot = 0;
     1343      doCrash = 1;
     1344      if (options.getExtraInfo(1) > 0)
     1345        doCrash = options.getExtraInfo(1);
     1346      doSprint = 0;
     1347      break;
     1348    case 2:
     1349      doIdiot = 1;
     1350      if (options.getExtraInfo(1) > 0)
     1351        doIdiot = options.getExtraInfo(1);
     1352      doCrash = 0;
     1353      doSprint = 0;
     1354      break;
     1355    case 3:
     1356      doIdiot = 0;
     1357      doCrash = 0;
     1358      doSprint = 1;
     1359      break;
     1360    case 4:
     1361      doIdiot = 0;
     1362      doCrash = 0;
     1363      doSprint = 0;
     1364      break;
     1365    case 5:
     1366      doIdiot = 0;
     1367      doCrash = -1;
     1368      doSprint = -1;
     1369      break;
     1370    case 6:
     1371      doIdiot = -1;
     1372      doCrash = -1;
     1373      doSprint = 0;
     1374      break;
     1375    case 7:
     1376      doIdiot = -1;
     1377      doCrash = 0;
     1378      doSprint = -1;
     1379      break;
     1380    case 8:
     1381      doIdiot = -1;
     1382      doCrash = 0;
     1383      doSprint = 0;
     1384      break;
     1385    case 9:
     1386      doIdiot = 0;
     1387      doCrash = 0;
     1388      doSprint = -1;
     1389      break;
     1390    case 10:
     1391      doIdiot = 0;
     1392      doCrash = 0;
     1393      doSprint = 0;
     1394      if (options.getExtraInfo(1))
     1395        doSlp = options.getExtraInfo(1);
     1396      break;
     1397    case 11:
     1398      doIdiot = 0;
     1399      doCrash = 0;
     1400      doSprint = 0;
     1401      primalStartup = 0;
     1402      break;
     1403    default:
     1404      abort();
     1405    }
     1406  } else if (method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe) {
     1407    // Dual
     1408    switch (options.getSpecialOption(0)) {
     1409    case 0:
     1410      doIdiot = 0;
     1411      doCrash = 0;
     1412      doSprint = 0;
     1413      break;
     1414    case 1:
     1415      doIdiot = 0;
     1416      doCrash = 1;
     1417      if (options.getExtraInfo(0) > 0)
     1418        doCrash = options.getExtraInfo(0);
     1419      doSprint = 0;
     1420      break;
     1421    case 2:
     1422      doIdiot = -1;
     1423      if (options.getExtraInfo(0) > 0)
     1424        doIdiot = options.getExtraInfo(0);
     1425      doCrash = 0;
     1426      doSprint = 0;
     1427      break;
     1428    default:
     1429      abort();
     1430    }
     1431  } else {
     1432    // decomposition
     1433  }
     1434#ifndef NO_RTTI
     1435  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objectiveAsObject()));
    13041436#else
    1305                if (numberRows > 200 && numberColumns > 2000 &&
    1306                    numberColumns > 2 * numberRows) {
    1307                  tryItSave= 3;
    1308                } else {
    1309                  // If rhs also rubbish then maybe
    1310                  int numberRowsOrig = this->numberRows();
    1311                  const double * rowLower = this->rowLower();
    1312                  const double * rowUpper = this->rowUpper();
    1313                  double last=COIN_DBL_MAX;
    1314                  int nDifferent=0;
    1315                  for (int i=0;i<numberRowsOrig;i++) {
    1316                    double value=fabs(rowLower[i]);
    1317                    if (value&&value!=COIN_DBL_MAX) {
    1318                      if (value!=last) {
    1319                        nDifferent++;
    1320                        last=value;
    1321                      }
    1322                    }
    1323                    value=fabs(rowUpper[i]);
    1324                    if (value&&value!=COIN_DBL_MAX) {
    1325                      if (value!=last) {
    1326                        nDifferent++;
    1327                        last=value;
    1328                      }
    1329                    }
    1330                  }
    1331                  if (nDifferent<2)
    1332                    tryItSave=1;
    1333                }
    1334 #endif
    1335           }
    1336      }
    1337      if (method != ClpSolve::useDual && method != ClpSolve::useBarrier
    1338          && method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe
    1339                && method != ClpSolve::useBarrierNoCross) {
    1340           switch (options.getSpecialOption(1)) {
    1341           case 0:
    1342                doIdiot = -1;
    1343                doCrash = -1;
    1344                doSprint = -1;
    1345                break;
    1346           case 1:
    1347                doIdiot = 0;
    1348                doCrash = 1;
    1349                if (options.getExtraInfo(1) > 0)
    1350                     doCrash = options.getExtraInfo(1);
    1351                doSprint = 0;
    1352                break;
    1353           case 2:
    1354                doIdiot = 1;
    1355                if (options.getExtraInfo(1) > 0)
    1356                     doIdiot = options.getExtraInfo(1);
    1357                doCrash = 0;
    1358                doSprint = 0;
    1359                break;
    1360           case 3:
    1361                doIdiot = 0;
    1362                doCrash = 0;
    1363                doSprint = 1;
    1364                break;
    1365           case 4:
    1366                doIdiot = 0;
    1367                doCrash = 0;
    1368                doSprint = 0;
    1369                break;
    1370           case 5:
    1371                doIdiot = 0;
    1372                doCrash = -1;
    1373                doSprint = -1;
    1374                break;
    1375           case 6:
    1376                doIdiot = -1;
    1377                doCrash = -1;
    1378                doSprint = 0;
    1379                break;
    1380           case 7:
    1381                doIdiot = -1;
    1382                doCrash = 0;
    1383                doSprint = -1;
    1384                break;
    1385           case 8:
    1386                doIdiot = -1;
    1387                doCrash = 0;
    1388                doSprint = 0;
    1389                break;
    1390           case 9:
    1391                doIdiot = 0;
    1392                doCrash = 0;
    1393                doSprint = -1;
    1394                break;
    1395           case 10:
    1396                doIdiot = 0;
    1397                doCrash = 0;
    1398                doSprint = 0;
    1399                if (options.getExtraInfo(1) )
    1400                     doSlp = options.getExtraInfo(1);
    1401                break;
    1402           case 11:
    1403                doIdiot = 0;
    1404                doCrash = 0;
    1405                doSprint = 0;
    1406                primalStartup = 0;
    1407                break;
    1408           default:
    1409                abort();
    1410           }
    1411      } else if (method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe) {
    1412           // Dual
    1413           switch (options.getSpecialOption(0)) {
    1414           case 0:
    1415                doIdiot = 0;
    1416                doCrash = 0;
    1417                doSprint = 0;
    1418                break;
    1419           case 1:
    1420                doIdiot = 0;
    1421                doCrash = 1;
    1422                if (options.getExtraInfo(0) > 0)
    1423                     doCrash = options.getExtraInfo(0);
    1424                doSprint = 0;
    1425                break;
    1426           case 2:
    1427                doIdiot = -1;
    1428                if (options.getExtraInfo(0) > 0)
    1429                     doIdiot = options.getExtraInfo(0);
    1430                doCrash = 0;
    1431                doSprint = 0;
    1432                break;
    1433           default:
    1434                abort();
    1435           }
    1436      } else {
    1437        // decomposition
    1438      }
    1439 #ifndef NO_RTTI
    1440      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objectiveAsObject()));
    1441 #else
    1442      ClpQuadraticObjective * quadraticObj = NULL;
    1443      if (objective_->type() == 2)
    1444           quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    1445 #endif
    1446      // If quadratic then primal or barrier or slp
    1447      if (quadraticObj) {
    1448           doSprint = 0;
    1449           //doIdiot = 0;
    1450           // off
    1451           if (method == ClpSolve::useBarrier)
    1452                method = ClpSolve::useBarrierNoCross;
    1453           else if (method != ClpSolve::useBarrierNoCross)
    1454                method = ClpSolve::usePrimal;
    1455      }
     1437  ClpQuadraticObjective *quadraticObj = NULL;
     1438  if (objective_->type() == 2)
     1439    quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
     1440#endif
     1441  // If quadratic then primal or barrier or slp
     1442  if (quadraticObj) {
     1443    doSprint = 0;
     1444    //doIdiot = 0;
     1445    // off
     1446    if (method == ClpSolve::useBarrier)
     1447      method = ClpSolve::useBarrierNoCross;
     1448    else if (method != ClpSolve::useBarrierNoCross)
     1449      method = ClpSolve::usePrimal;
     1450  }
    14561451#ifdef COIN_HAS_VOL
    1457      // Save number of idiot
    1458      int saveDoIdiot = doIdiot;
    1459 #endif
    1460      // Just do this number of passes in Sprint
    1461      int maxSprintPass = 100;
    1462      // See if worth trying +- one matrix
    1463      bool plusMinus = false;
    1464      CoinBigIndex numberElements = model2->getNumElements();
     1452  // Save number of idiot
     1453  int saveDoIdiot = doIdiot;
     1454#endif
     1455  // Just do this number of passes in Sprint
     1456  int maxSprintPass = 100;
     1457  // See if worth trying +- one matrix
     1458  bool plusMinus = false;
     1459  CoinBigIndex numberElements = model2->getNumElements();
    14651460#ifndef SLIM_CLP
    14661461#ifndef NO_RTTI
    1467      if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
    1468           // network - switch off stuff
     1462  if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) {
     1463    // network - switch off stuff
     1464    doIdiot = 0;
     1465    if (doSprint < 0)
     1466      doSprint = 0;
     1467  }
     1468#else
     1469  if (matrix_->type() == 11) {
     1470    // network - switch off stuff
     1471    doIdiot = 0;
     1472    //doSprint=0;
     1473  }
     1474#endif
     1475#endif
     1476  int numberColumns = model2->numberColumns();
     1477  int numberRows = model2->numberRows();
     1478  // If not all slack basis - switch off all except sprint
     1479  int numberRowsBasic = 0;
     1480  int iRow;
     1481  for (iRow = 0; iRow < numberRows; iRow++)
     1482    if (model2->getRowStatus(iRow) == basic)
     1483      numberRowsBasic++;
     1484  if (numberRowsBasic < numberRows && objective_->type() < 2) {
     1485    doIdiot = 0;
     1486    doCrash = 0;
     1487    //doSprint=0;
     1488  }
     1489  if (options.getSpecialOption(3) == 0) {
     1490    if (numberElements > 100000)
     1491      plusMinus = true;
     1492    if (numberElements > 10000 && (doIdiot || doSprint))
     1493      plusMinus = true;
     1494  } else if ((specialOptions_ & 1024) != 0) {
     1495    plusMinus = true;
     1496  }
     1497#ifndef SLIM_CLP
     1498  // Statistics (+1,-1, other) - used to decide on strategy if not +-1
     1499  CoinBigIndex statistics[3] = { -1, 0, 0 };
     1500  if (plusMinus) {
     1501    saveMatrix = model2->clpMatrix();
     1502#ifndef NO_RTTI
     1503    ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(saveMatrix);
     1504#else
     1505    ClpPackedMatrix *clpMatrix = NULL;
     1506    if (saveMatrix->type() == 1)
     1507      clpMatrix = static_cast< ClpPackedMatrix * >(saveMatrix);
     1508#endif
     1509    if (clpMatrix) {
     1510      ClpPlusMinusOneMatrix *newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
     1511      if (newMatrix->getIndices()) {
     1512        // CHECKME This test of specialOptions and the one above
     1513        // don't seem compatible.
     1514#ifndef ABC_INHERIT
     1515        if ((specialOptions_ & 1024) == 0) {
     1516          model2->replaceMatrix(newMatrix);
     1517        } else {
     1518#endif
     1519          // in integer (or abc) - just use for sprint/idiot
     1520          saveMatrix = NULL;
     1521          delete newMatrix;
     1522#ifndef ABC_INHERIT
     1523        }
     1524#endif
     1525      } else {
     1526        handler_->message(CLP_MATRIX_CHANGE, messages_)
     1527          << "+- 1"
     1528          << CoinMessageEol;
     1529        CoinMemcpyN(newMatrix->startPositive(), 3, statistics);
     1530        saveMatrix = NULL;
     1531        plusMinus = false;
     1532        delete newMatrix;
     1533      }
     1534    } else {
     1535      saveMatrix = NULL;
     1536      plusMinus = false;
     1537    }
     1538  }
     1539#endif
     1540  if (this->factorizationFrequency() == 200) {
     1541    // User did not touch preset
     1542    model2->defaultFactorizationFrequency();
     1543  } else if (model2 != this) {
     1544    // make sure model2 has correct value
     1545    model2->setFactorizationFrequency(this->factorizationFrequency());
     1546  }
     1547  if (method == ClpSolve::automatic) {
     1548    if (doSprint == 0 && doIdiot == 0) {
     1549      // off
     1550      method = ClpSolve::useDual;
     1551    } else {
     1552      // only do primal if sprint or idiot
     1553      if (doSprint > 0) {
     1554        method = ClpSolve::usePrimalorSprint;
     1555      } else if (doIdiot > 0) {
     1556        method = ClpSolve::usePrimal;
     1557      } else {
     1558        if (numberElements < 500000) {
     1559          // Small problem
     1560          if (numberRows * 10 > numberColumns || numberColumns < 6000
     1561            || (numberRows * 20 > numberColumns && !plusMinus))
     1562            doSprint = 0; // switch off sprint
     1563        } else {
     1564          // larger problem
     1565          if (numberRows * 8 > numberColumns)
     1566            doSprint = 0; // switch off sprint
     1567        }
     1568        // switch off idiot or sprint if any free variable
     1569        // switch off sprint if very few with costs
     1570        // or great variation in cost
     1571        int iColumn;
     1572        const double *columnLower = model2->columnLower();
     1573        const double *columnUpper = model2->columnUpper();
     1574        const double *objective = model2->objective();
     1575        int nObj = 0;
     1576        int nFree = 0;
     1577        double smallestObj = COIN_DBL_MAX;
     1578        double largestObj = 0.0;
     1579        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1580          if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
     1581            nFree++;
     1582          } else if (objective[iColumn]) {
     1583            nObj++;
     1584            smallestObj = CoinMin(smallestObj, objective[iColumn]);
     1585            largestObj = CoinMax(largestObj, objective[iColumn]);
     1586          }
     1587        }
     1588        if (nObj * 10 < numberColumns || smallestObj * 10.0 < largestObj)
     1589          doSprint = 0;
     1590        if (nFree)
    14691591          doIdiot = 0;
    1470           if (doSprint < 0)
    1471                doSprint = 0;
    1472      }
     1592        if (nFree * 10 > numberRows)
     1593          doSprint = 0;
     1594        int nPasses = 0;
     1595        // look at rhs
     1596        int iRow;
     1597        double largest = 0.0;
     1598        double smallest = 1.0e30;
     1599        double largestGap = 0.0;
     1600        int numberNotE = 0;
     1601        bool notInteger = false;
     1602        for (iRow = 0; iRow < numberRows; iRow++) {
     1603          double value1 = model2->rowLower_[iRow];
     1604          if (value1 && value1 > -1.0e31) {
     1605            largest = CoinMax(largest, fabs(value1));
     1606            smallest = CoinMin(smallest, fabs(value1));
     1607            if (fabs(value1 - floor(value1 + 0.5)) > 1.0e-8) {
     1608              notInteger = true;
     1609              break;
     1610            }
     1611          }
     1612          double value2 = model2->rowUpper_[iRow];
     1613          if (value2 && value2 < 1.0e31) {
     1614            largest = CoinMax(largest, fabs(value2));
     1615            smallest = CoinMin(smallest, fabs(value2));
     1616            if (fabs(value2 - floor(value2 + 0.5)) > 1.0e-8) {
     1617              notInteger = true;
     1618              break;
     1619            }
     1620          }
     1621          // CHECKME This next bit can't be right...
     1622          if (value2 > value1) {
     1623            numberNotE++;
     1624            //if (value2 > 1.0e31 || value1 < -1.0e31)
     1625            //   largestGap = COIN_DBL_MAX;
     1626            //else
     1627            //   largestGap = value2 - value1;
     1628          }
     1629        }
     1630        int tryIt = (numberRows > 200 && numberColumns > 2000 && (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0))) ? 3 : 0;
     1631        tryItSave = tryIt;
     1632        if (numberRows < 1000 && numberColumns < 3000)
     1633          tryIt = 0;
     1634        if (notInteger)
     1635          tryIt = 0;
     1636        if (largest / smallest > 10 || (largest / smallest > 2.0 && largest > 50))
     1637          tryIt = 0;
     1638        if (tryIt) {
     1639          if (largest / smallest > 2.0) {
     1640            nPasses = 10 + numberColumns / 100000;
     1641            nPasses = CoinMin(nPasses, 50);
     1642            nPasses = CoinMax(nPasses, 15);
     1643            if (numberRows > 20000 && nPasses > 5) {
     1644              // Might as well go for it
     1645              nPasses = CoinMax(nPasses, 71);
     1646            } else if (numberRows > 2000 && nPasses > 5) {
     1647              nPasses = CoinMax(nPasses, 50);
     1648            } else if (numberElements < 3 * numberColumns) {
     1649              nPasses = CoinMin(nPasses, 10); // probably not worh it
     1650            }
     1651          } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) {
     1652            nPasses = 10 + numberColumns / 1000;
     1653            nPasses = CoinMin(nPasses, 100);
     1654            nPasses = CoinMax(nPasses, 30);
     1655            if (numberRows > 25000) {
     1656              // Might as well go for it
     1657              nPasses = CoinMax(nPasses, 71);
     1658            }
     1659            if (!largestGap)
     1660              nPasses *= 2;
     1661          } else {
     1662            nPasses = 10 + numberColumns / 1000;
     1663            nPasses = CoinMax(nPasses, 100);
     1664            if (!largestGap)
     1665              nPasses *= 2;
     1666            nPasses = CoinMin(nPasses, 200);
     1667          }
     1668        }
     1669        //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
     1670        //     numberRows,numberColumns,plusMinus ? 'Y' : 'N',
     1671        //     tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
     1672        //exit(0);
     1673        if (!tryIt || nPasses <= 5)
     1674          doIdiot = 0;
     1675#if CLPSOLVE_ACTIONS
     1676        if (doIdiot && doSprint < 0 && wasAutomatic && 20 * model2->numberRows() > model2->numberColumns())
     1677          doSprint = 0; // switch off sprint
     1678#endif
     1679        if (doSprint) {
     1680          method = ClpSolve::usePrimalorSprint;
     1681        } else if (doIdiot) {
     1682          method = ClpSolve::usePrimal;
     1683        } else {
     1684          method = ClpSolve::useDual;
     1685        }
     1686      }
     1687    }
     1688  }
     1689  if (method == ClpSolve::tryBenders) {
     1690    // Now build model
     1691    int lengthNames = model2->lengthNames();
     1692    model2->setLengthNames(0);
     1693    CoinModel *build = model2->createCoinModel();
     1694    model2->setLengthNames(lengthNames);
     1695    CoinStructuredModel benders;
     1696    build->convertMatrix();
     1697    int numberBlocks = options.independentOption(0);
     1698    benders.setMessageHandler(handler_);
     1699    numberBlocks = benders.decompose(*build, 2, numberBlocks, NULL);
     1700    delete build;
     1701    //exit(0);
     1702    if (numberBlocks) {
     1703      options.setIndependentOption(1, 1); // don't do final clean up
     1704      model2->solveBenders(&benders, options);
     1705      //move solution
     1706      method = ClpSolve::notImplemented;
     1707      time2 = CoinCpuTime();
     1708      timeCore = time2 - timeX;
     1709      handler_->message(CLP_INTERVAL_TIMING, messages_)
     1710        << "Crossover" << timeCore << time2 - time1
     1711        << CoinMessageEol;
     1712      timeX = time2;
     1713    } else {
     1714      printf("No structure\n");
     1715      method = ClpSolve::useDual;
     1716    }
     1717  } else if (method == ClpSolve::tryDantzigWolfe) {
     1718    abort();
     1719  }
     1720  if (method == ClpSolve::usePrimalorSprint) {
     1721    if (doSprint < 0) {
     1722      if (numberElements < 500000) {
     1723        // Small problem
     1724        if (numberRows * 10 > numberColumns || numberColumns < 6000
     1725          || (numberRows * 20 > numberColumns && !plusMinus))
     1726          method = ClpSolve::usePrimal; // switch off sprint
     1727      } else {
     1728        // larger problem
     1729        if (numberRows * 8 > numberColumns)
     1730          method = ClpSolve::usePrimal; // switch off sprint
     1731        // but make lightweight
     1732        if (numberRows * 10 > numberColumns || numberColumns < 6000
     1733          || (numberRows * 20 > numberColumns && !plusMinus))
     1734          maxSprintPass = 10;
     1735      }
     1736    } else if (doSprint == 0) {
     1737      method = ClpSolve::usePrimal; // switch off sprint
     1738    }
     1739  }
     1740  if (method == ClpSolve::useDual) {
     1741#ifdef CLP_USEFUL_PRINTOUT
     1742    debugInt[6] = 1;
     1743#endif
     1744    double *saveLower = NULL;
     1745    double *saveUpper = NULL;
     1746    if (presolve == ClpSolve::presolveOn) {
     1747      int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0);
     1748      if (numberInfeasibilities) {
     1749        handler_->message(CLP_INFEASIBLE, messages_)
     1750          << CoinMessageEol;
     1751        delete model2;
     1752        model2 = this;
     1753        presolve = ClpSolve::presolveOff;
     1754      }
     1755    } else if (numberRows_ + numberColumns_ > 5000) {
     1756      // do anyway
     1757      saveLower = new double[numberRows_ + numberColumns_];
     1758      CoinMemcpyN(model2->columnLower(), numberColumns_, saveLower);
     1759      CoinMemcpyN(model2->rowLower(), numberRows_, saveLower + numberColumns_);
     1760      saveUpper = new double[numberRows_ + numberColumns_];
     1761      CoinMemcpyN(model2->columnUpper(), numberColumns_, saveUpper);
     1762      CoinMemcpyN(model2->rowUpper(), numberRows_, saveUpper + numberColumns_);
     1763      int numberInfeasibilities = model2->tightenPrimalBounds();
     1764      if (numberInfeasibilities) {
     1765        handler_->message(CLP_INFEASIBLE, messages_)
     1766          << CoinMessageEol;
     1767        CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
     1768        CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
     1769        delete[] saveLower;
     1770        saveLower = NULL;
     1771        CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
     1772        CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
     1773        delete[] saveUpper;
     1774        saveUpper = NULL;
     1775        // return if wanted
     1776        if (options.infeasibleReturn() || (moreSpecialOptions_ & 1) != 0)
     1777          return -1;
     1778      }
     1779    }
     1780#ifndef COIN_HAS_VOL
     1781    // switch off idiot and volume for now
     1782    doIdiot = 0;
     1783#endif
     1784    // pick up number passes
     1785    int nPasses = 0;
     1786    int numberNotE = 0;
     1787#ifndef SLIM_CLP
     1788    if ((doIdiot < 0 && plusMinus) || doIdiot > 0) {
     1789      // See if candidate for idiot
     1790      nPasses = 0;
     1791      Idiot info(*model2);
     1792      info.setStrategy(idiotOptions | info.getStrategy());
     1793      // Get average number of elements per column
     1794      double ratio = static_cast< double >(numberElements) / static_cast< double >(numberColumns);
     1795      // look at rhs
     1796      int iRow;
     1797      double largest = 0.0;
     1798      double smallest = 1.0e30;
     1799      for (iRow = 0; iRow < numberRows; iRow++) {
     1800        double value1 = model2->rowLower_[iRow];
     1801        if (value1 && value1 > -1.0e31) {
     1802          largest = CoinMax(largest, fabs(value1));
     1803          smallest = CoinMin(smallest, fabs(value1));
     1804        }
     1805        double value2 = model2->rowUpper_[iRow];
     1806        if (value2 && value2 < 1.0e31) {
     1807          largest = CoinMax(largest, fabs(value2));
     1808          smallest = CoinMin(smallest, fabs(value2));
     1809        }
     1810        if (value2 > value1) {
     1811          numberNotE++;
     1812        }
     1813      }
     1814      if (doIdiot < 0) {
     1815        if (numberRows > 200 && numberColumns > 5000 && ratio >= 3.0 && largest / smallest < 1.1 && !numberNotE) {
     1816          nPasses = 71;
     1817        }
     1818      }
     1819      if (doIdiot > 0) {
     1820        nPasses = CoinMax(nPasses, doIdiot);
     1821        if (nPasses > 70) {
     1822          info.setStartingWeight(1.0e3);
     1823          info.setDropEnoughFeasibility(0.01);
     1824        }
     1825      }
     1826      if (nPasses > 20) {
     1827#ifdef COIN_HAS_VOL
     1828        int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
     1829        if (!returnCode) {
     1830          time2 = CoinCpuTime();
     1831          timeIdiot = time2 - timeX;
     1832          handler_->message(CLP_INTERVAL_TIMING, messages_)
     1833            << "Idiot Crash" << timeIdiot << time2 - time1
     1834            << CoinMessageEol;
     1835          timeX = time2;
     1836        } else {
     1837          nPasses = 0;
     1838        }
    14731839#else
    1474      if (matrix_->type() == 11) {
    1475           // network - switch off stuff
    1476           doIdiot = 0;
    1477           //doSprint=0;
    1478      }
    1479 #endif
    1480 #endif
    1481      int numberColumns = model2->numberColumns();
    1482      int numberRows = model2->numberRows();
    1483      // If not all slack basis - switch off all except sprint
    1484      int numberRowsBasic = 0;
    1485      int iRow;
    1486      for (iRow = 0; iRow < numberRows; iRow++)
    1487           if (model2->getRowStatus(iRow) == basic)
    1488                numberRowsBasic++;
    1489      if (numberRowsBasic < numberRows && objective_->type()<2) {
    1490           doIdiot = 0;
    1491           doCrash = 0;
    1492           //doSprint=0;
    1493      }
    1494      if (options.getSpecialOption(3) == 0) {
    1495           if(numberElements > 100000)
    1496                plusMinus = true;
    1497           if(numberElements > 10000 && (doIdiot || doSprint))
    1498                plusMinus = true;
    1499      } else if ((specialOptions_ & 1024) != 0) {
    1500           plusMinus = true;
    1501      }
     1840        nPasses = 0;
     1841#endif
     1842      } else {
     1843        nPasses = 0;
     1844      }
     1845    }
     1846#endif
     1847    if (doCrash) {
     1848#ifdef ABC_INHERIT
     1849      if (!model2->abcState()) {
     1850#endif
     1851        switch (doCrash) {
     1852          // standard
     1853        case 1:
     1854          model2->crash(1000, 1);
     1855          break;
     1856          // As in paper by Solow and Halim (approx)
     1857        case 2:
     1858        case 3:
     1859          model2->crash(model2->dualBound(), 0);
     1860          break;
     1861          // Just put free in basis
     1862        case 4:
     1863          model2->crash(0.0, 3);
     1864          break;
     1865        }
     1866#ifdef ABC_INHERIT
     1867      } else if (doCrash >= 0) {
     1868        model2->setAbcState(model2->abcState() | 256 * doCrash);
     1869      }
     1870#endif
     1871    }
     1872    if (!nPasses) {
     1873      int saveOptions = model2->specialOptions();
     1874      if (model2->numberRows() > 100)
     1875        model2->setSpecialOptions(saveOptions | 64); // go as far as possible
     1876      //int numberRows = model2->numberRows();
     1877      //int numberColumns = model2->numberColumns();
     1878      if (dynamic_cast< ClpPackedMatrix * >(matrix_)) {
     1879        // See if original wanted vector
     1880        ClpPackedMatrix *clpMatrixO = dynamic_cast< ClpPackedMatrix * >(matrix_);
     1881        ClpMatrixBase *matrix = model2->clpMatrix();
     1882        if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
     1883          ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     1884          clpMatrix->makeSpecialColumnCopy();
     1885          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
     1886          model2->dual(0);
     1887          clpMatrix->releaseSpecialColumnCopy();
     1888        } else {
     1889#ifndef ABC_INHERIT
     1890          model2->dual(0);
     1891#else
     1892          model2->dealWithAbc(0, 0, interrupt);
     1893#endif
     1894        }
     1895      } else {
     1896        model2->dual(0);
     1897      }
     1898    } else if (!numberNotE && 0) {
     1899      // E so we can do in another way
     1900      double *pi = model2->dualRowSolution();
     1901      int i;
     1902      int numberColumns = model2->numberColumns();
     1903      int numberRows = model2->numberRows();
     1904      double *saveObj = new double[numberColumns];
     1905      CoinMemcpyN(model2->objective(), numberColumns, saveObj);
     1906      CoinMemcpyN(model2->objective(),
     1907        numberColumns, model2->dualColumnSolution());
     1908      model2->clpMatrix()->transposeTimes(-1.0, pi, model2->dualColumnSolution());
     1909      CoinMemcpyN(model2->dualColumnSolution(),
     1910        numberColumns, model2->objective());
     1911      const double *rowsol = model2->primalRowSolution();
     1912      double offset = 0.0;
     1913      for (i = 0; i < numberRows; i++) {
     1914        offset += pi[i] * rowsol[i];
     1915      }
     1916      double value2;
     1917      model2->getDblParam(ClpObjOffset, value2);
     1918      //printf("Offset %g %g\n",offset,value2);
     1919      model2->setDblParam(ClpObjOffset, value2 - offset);
     1920      model2->setPerturbation(51);
     1921      //model2->setRowObjective(pi);
     1922      // zero out pi
     1923      //memset(pi,0,numberRows*sizeof(double));
     1924      // Could put some in basis - only partially tested
     1925      model2->allSlackBasis();
     1926      //model2->factorization()->maximumPivots(200);
     1927      //model2->setLogLevel(63);
     1928      // solve
     1929      model2->dual(0);
     1930      model2->setDblParam(ClpObjOffset, value2);
     1931      CoinMemcpyN(saveObj, numberColumns, model2->objective());
     1932      // zero out pi
     1933      //memset(pi,0,numberRows*sizeof(double));
     1934      //model2->setRowObjective(pi);
     1935      delete[] saveObj;
     1936      //model2->dual(0);
     1937      model2->setPerturbation(50);
     1938      model2->primal();
     1939    } else {
     1940      // solve
     1941      model2->setPerturbation(100);
     1942      model2->dual(2);
     1943      model2->setPerturbation(50);
     1944      model2->dual(0);
     1945    }
     1946    if (saveLower) {
     1947      CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
     1948      CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
     1949      delete[] saveLower;
     1950      saveLower = NULL;
     1951      CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
     1952      CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
     1953      delete[] saveUpper;
     1954      saveUpper = NULL;
     1955    }
     1956    time2 = CoinCpuTime();
     1957    timeCore = time2 - timeX;
     1958    handler_->message(CLP_INTERVAL_TIMING, messages_)
     1959      << "Dual" << timeCore << time2 - time1
     1960      << CoinMessageEol;
     1961    timeX = time2;
     1962  } else if (method == ClpSolve::usePrimal) {
     1963#ifdef CLP_USEFUL_PRINTOUT
     1964    debugInt[6] = 2;
     1965#endif
    15021966#ifndef SLIM_CLP
    1503      // Statistics (+1,-1, other) - used to decide on strategy if not +-1
    1504      CoinBigIndex statistics[3] = { -1, 0, 0};
    1505      if (plusMinus) {
    1506           saveMatrix = model2->clpMatrix();
    1507 #ifndef NO_RTTI
    1508           ClpPackedMatrix* clpMatrix =
    1509                dynamic_cast< ClpPackedMatrix*>(saveMatrix);
    1510 #else
    1511           ClpPackedMatrix* clpMatrix = NULL;
    1512           if (saveMatrix->type() == 1)
    1513                clpMatrix =
    1514                     static_cast< ClpPackedMatrix*>(saveMatrix);
    1515 #endif
    1516           if (clpMatrix) {
    1517                ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
    1518                if (newMatrix->getIndices()) {
    1519                   // CHECKME This test of specialOptions and the one above
    1520                   // don't seem compatible.
    1521 #ifndef ABC_INHERIT
    1522                     if ((specialOptions_ & 1024) == 0) {
    1523                          model2->replaceMatrix(newMatrix);
    1524                     } else {
    1525 #endif
    1526                          // in integer (or abc) - just use for sprint/idiot
    1527                          saveMatrix = NULL;
    1528                          delete newMatrix;
    1529 #ifndef ABC_INHERIT
    1530                     }
    1531 #endif
    1532                } else {
    1533                     handler_->message(CLP_MATRIX_CHANGE, messages_)
    1534                               << "+- 1"
    1535                               << CoinMessageEol;
    1536                     CoinMemcpyN(newMatrix->startPositive(), 3, statistics);
    1537                     saveMatrix = NULL;
    1538                     plusMinus = false;
    1539                     delete newMatrix;
    1540                }
     1967    if (doIdiot) {
     1968      int nPasses = 0;
     1969      Idiot info(*model2);
     1970      info.setStrategy(idiotOptions | info.getStrategy());
     1971      // Get average number of elements per column
     1972      double ratio = static_cast< double >(numberElements) / static_cast< double >(numberColumns);
     1973      // look at rhs
     1974      int iRow;
     1975      double largest = 0.0;
     1976      double smallest = 1.0e30;
     1977      double largestGap = 0.0;
     1978      int numberNotE = 0;
     1979      for (iRow = 0; iRow < numberRows; iRow++) {
     1980        double value1 = model2->rowLower_[iRow];
     1981        if (value1 && value1 > -1.0e31) {
     1982          largest = CoinMax(largest, fabs(value1));
     1983          smallest = CoinMin(smallest, fabs(value1));
     1984        }
     1985        double value2 = model2->rowUpper_[iRow];
     1986        if (value2 && value2 < 1.0e31) {
     1987          largest = CoinMax(largest, fabs(value2));
     1988          smallest = CoinMin(smallest, fabs(value2));
     1989        }
     1990        if (value2 > value1) {
     1991          numberNotE++;
     1992          if (value2 > 1.0e31 || value1 < -1.0e31)
     1993            largestGap = COIN_DBL_MAX;
     1994          else
     1995            largestGap = value2 - value1;
     1996        }
     1997      }
     1998      bool increaseSprint = plusMinus;
     1999      if ((specialOptions_ & 1024) != 0)
     2000        increaseSprint = false;
     2001      if (!plusMinus) {
     2002        // If 90% +- 1 then go for sprint
     2003        if (statistics[0] >= 0 && 10 * statistics[2] < statistics[0] + statistics[1])
     2004          increaseSprint = true;
     2005      }
     2006      int tryIt = tryItSave;
     2007      if (numberRows < 1000 && numberColumns < 3000)
     2008        tryIt = 0;
     2009      if (tryIt) {
     2010        if (increaseSprint) {
     2011          info.setStartingWeight(1.0e3);
     2012          info.setReduceIterations(6);
     2013          // also be more lenient on infeasibilities
     2014          info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
     2015          info.setDropEnoughWeighted(-2.0);
     2016          if (largest / smallest > 2.0) {
     2017            nPasses = 10 + numberColumns / 100000;
     2018            nPasses = CoinMin(nPasses, 50);
     2019            nPasses = CoinMax(nPasses, 15);
     2020            if (numberRows > 20000 && nPasses > 5) {
     2021              // Might as well go for it
     2022              nPasses = CoinMax(nPasses, 71);
     2023            } else if (numberRows > 2000 && nPasses > 5) {
     2024              nPasses = CoinMax(nPasses, 50);
     2025            } else if (numberElements < 3 * numberColumns) {
     2026              nPasses = CoinMin(nPasses, 10); // probably not worh it
     2027              if (doIdiot < 0)
     2028                info.setLightweight(1); // say lightweight idiot
     2029            } else {
     2030              if (doIdiot < 0)
     2031                info.setLightweight(1); // say lightweight idiot
     2032            }
     2033          } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) {
     2034            nPasses = 10 + numberColumns / 1000;
     2035            nPasses = CoinMin(nPasses, 100);
     2036            nPasses = CoinMax(nPasses, 30);
     2037            if (numberRows > 25000) {
     2038              // Might as well go for it
     2039              nPasses = CoinMax(nPasses, 71);
     2040            }
     2041            if (!largestGap)
     2042              nPasses *= 2;
    15412043          } else {
    1542                saveMatrix = NULL;
    1543                plusMinus = false;
     2044            nPasses = 10 + numberColumns / 1000;
     2045            nPasses = CoinMin(nPasses, 200);
     2046            nPasses = CoinMax(nPasses, 100);
     2047            info.setStartingWeight(1.0e-1);
     2048            info.setReduceIterations(6);
     2049            if (!largestGap && nPasses <= 50)
     2050              nPasses *= 2;
     2051            //info.setFeasibilityTolerance(1.0e-7);
    15442052          }
    1545      }
    1546 #endif
    1547      if (this->factorizationFrequency() == 200) {
    1548           // User did not touch preset
    1549           model2->defaultFactorizationFrequency();
    1550      } else if (model2 != this) {
    1551           // make sure model2 has correct value
    1552           model2->setFactorizationFrequency(this->factorizationFrequency());
    1553      }
    1554      if (method == ClpSolve::automatic) {
    1555           if (doSprint == 0 && doIdiot == 0) {
    1556                // off
    1557                method = ClpSolve::useDual;
     2053          // If few passes - don't bother
     2054          if (nPasses <= 5 && !plusMinus)
     2055            nPasses = 0;
     2056        } else {
     2057          if (doIdiot < 0)
     2058            info.setLightweight(1); // say lightweight idiot
     2059          if (largest / smallest > 1.01 || numberNotE || statistics[2] > statistics[0] + statistics[1]) {
     2060            if (numberRows > 25000 || numberColumns > 5 * numberRows) {
     2061              nPasses = 50;
     2062            } else if (numberColumns > 4 * numberRows) {
     2063              nPasses = 20;
     2064            } else {
     2065              nPasses = 5;
     2066            }
    15582067          } else {
    1559                // only do primal if sprint or idiot
    1560                if (doSprint > 0) {
    1561                     method = ClpSolve::usePrimalorSprint;
    1562                } else if (doIdiot > 0) {
    1563                     method = ClpSolve::usePrimal;
    1564                } else {
    1565                     if (numberElements < 500000) {
    1566                          // Small problem
    1567                          if(numberRows * 10 > numberColumns || numberColumns < 6000
    1568                                    || (numberRows * 20 > numberColumns && !plusMinus))
    1569                               doSprint = 0; // switch off sprint
    1570                     } else {
    1571                          // larger problem
    1572                          if(numberRows * 8 > numberColumns)
    1573                               doSprint = 0; // switch off sprint
    1574                     }
    1575                     // switch off idiot or sprint if any free variable
    1576                     // switch off sprint if very few with costs
    1577                     // or great variation in cost
    1578                     int iColumn;
    1579                     const double * columnLower = model2->columnLower();
    1580                     const double * columnUpper = model2->columnUpper();
    1581                     const double * objective = model2->objective();
    1582                     int nObj=0;
    1583                     int nFree=0;
    1584                     double smallestObj=COIN_DBL_MAX;
    1585                     double largestObj=0.0;
    1586                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    1587                          if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
    1588                            nFree++;
    1589                          } else if (objective[iColumn]) {
    1590                            nObj++;
    1591                            smallestObj=CoinMin(smallestObj,objective[iColumn]);
    1592                            largestObj=CoinMax(largestObj,objective[iColumn]);
    1593                          }
    1594                     }
    1595                     if (nObj*10<numberColumns ||
    1596                         smallestObj*10.0<largestObj)
    1597                       doSprint=0;
    1598                     if (nFree)
    1599                       doIdiot=0;
    1600                     if (nFree*10>numberRows)
    1601                       doSprint=0;
    1602                     int nPasses = 0;
    1603                     // look at rhs
    1604                     int iRow;
    1605                     double largest = 0.0;
    1606                     double smallest = 1.0e30;
    1607                     double largestGap = 0.0;
    1608                     int numberNotE = 0;
    1609                     bool notInteger = false;
    1610                     for (iRow = 0; iRow < numberRows; iRow++) {
    1611                          double value1 = model2->rowLower_[iRow];
    1612                          if (value1 && value1 > -1.0e31) {
    1613                               largest = CoinMax(largest, fabs(value1));
    1614                               smallest = CoinMin(smallest, fabs(value1));
    1615                               if (fabs(value1 - floor(value1 + 0.5)) > 1.0e-8) {
    1616                                    notInteger = true;
    1617                                    break;
    1618                               }
    1619                          }
    1620                          double value2 = model2->rowUpper_[iRow];
    1621                          if (value2 && value2 < 1.0e31) {
    1622                               largest = CoinMax(largest, fabs(value2));
    1623                               smallest = CoinMin(smallest, fabs(value2));
    1624                               if (fabs(value2 - floor(value2 + 0.5)) > 1.0e-8) {
    1625                                    notInteger = true;
    1626                                    break;
    1627                               }
    1628                          }
    1629                          // CHECKME This next bit can't be right...
    1630                          if (value2 > value1) {
    1631                               numberNotE++;
    1632                               //if (value2 > 1.0e31 || value1 < -1.0e31)
    1633                               //   largestGap = COIN_DBL_MAX;
    1634                               //else
    1635                               //   largestGap = value2 - value1;
    1636                          }
    1637                     }
    1638                     int tryIt = (numberRows > 200 && numberColumns > 2000 &&
    1639                                  (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0))) ? 3 : 0;
    1640                     tryItSave = tryIt;
    1641                     if (numberRows < 1000 && numberColumns < 3000)
    1642                          tryIt = 0;
    1643                     if (notInteger)
    1644                          tryIt = 0;
    1645                     if (largest / smallest > 10 || (largest / smallest > 2.0 && largest > 50))
    1646                          tryIt = 0;
    1647                     if (tryIt) {
    1648                          if (largest / smallest > 2.0) {
    1649                               nPasses = 10 + numberColumns / 100000;
    1650                               nPasses = CoinMin(nPasses, 50);
    1651                               nPasses = CoinMax(nPasses, 15);
    1652                               if (numberRows > 20000 && nPasses > 5) {
    1653                                    // Might as well go for it
    1654                                    nPasses = CoinMax(nPasses, 71);
    1655                               } else if (numberRows > 2000 && nPasses > 5) {
    1656                                    nPasses = CoinMax(nPasses, 50);
    1657                               } else if (numberElements < 3 * numberColumns) {
    1658                                    nPasses = CoinMin(nPasses, 10); // probably not worh it
    1659                               }
    1660                          } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) {
    1661                               nPasses = 10 + numberColumns / 1000;
    1662                               nPasses = CoinMin(nPasses, 100);
    1663                               nPasses = CoinMax(nPasses, 30);
    1664                               if (numberRows > 25000) {
    1665                                    // Might as well go for it
    1666                                    nPasses = CoinMax(nPasses, 71);
    1667                               }
    1668                               if (!largestGap)
    1669                                    nPasses *= 2;
    1670                          } else {
    1671                               nPasses = 10 + numberColumns / 1000;
    1672                               nPasses = CoinMax(nPasses, 100);
    1673                               if (!largestGap)
    1674                                    nPasses *= 2;
    1675                               nPasses = CoinMin(nPasses, 200);
    1676                          }
    1677                     }
    1678                     //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
    1679                     //     numberRows,numberColumns,plusMinus ? 'Y' : 'N',
    1680                     //     tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
    1681                     //exit(0);
    1682                     if (!tryIt || nPasses <= 5)
    1683                          doIdiot = 0;
    1684 #if CLPSOLVE_ACTIONS
    1685                     if (doIdiot&&doSprint<0&&wasAutomatic &&
    1686                         20*model2->numberRows()>model2->numberColumns())
    1687                       doSprint=0; // switch off sprint
    1688 #endif
    1689                     if (doSprint) {
    1690                          method = ClpSolve::usePrimalorSprint;
    1691                     } else if (doIdiot) {
    1692                          method = ClpSolve::usePrimal;
    1693                     } else {
    1694                          method = ClpSolve::useDual;
    1695                     }
    1696                }
     2068            if (numberRows > 25000 || numberColumns > 5 * numberRows) {
     2069              nPasses = 50;
     2070              info.setLightweight(0); // say not lightweight idiot
     2071            } else if (numberColumns > 4 * numberRows) {
     2072              nPasses = 20;
     2073            } else {
     2074              nPasses = 15;
     2075            }
    16972076          }
    1698      }
    1699      if (method == ClpSolve::tryBenders) {
    1700        // Now build model
    1701        int lengthNames=model2->lengthNames();
    1702        model2->setLengthNames(0);
    1703        CoinModel * build = model2->createCoinModel();
    1704        model2->setLengthNames(lengthNames);
    1705        CoinStructuredModel benders;
    1706        build->convertMatrix();
    1707        int numberBlocks = options.independentOption(0);
    1708        benders.setMessageHandler(handler_);
    1709        numberBlocks=benders.decompose(*build,2,numberBlocks,NULL);
    1710        delete build;
    1711        //exit(0);
    1712        if (numberBlocks) {
    1713          options.setIndependentOption(1,1); // don't do final clean up
    1714          model2->solveBenders(&benders,options);
    1715          //move solution
    1716          method=ClpSolve::notImplemented;
    1717          time2 = CoinCpuTime();
    1718          timeCore = time2 - timeX;
    1719          handler_->message(CLP_INTERVAL_TIMING, messages_)
    1720            << "Crossover" << timeCore << time2 - time1
    1721            << CoinMessageEol;
    1722          timeX = time2;
    1723        } else {
    1724          printf("No structure\n");
    1725          method=ClpSolve::useDual;
    1726        }
    1727      } else if (method == ClpSolve::tryDantzigWolfe) {
    1728          abort();
    1729      }
    1730      if (method == ClpSolve::usePrimalorSprint) {
    1731           if (doSprint < 0) {
    1732                if (numberElements < 500000) {
    1733                     // Small problem
    1734                     if(numberRows * 10 > numberColumns || numberColumns < 6000
    1735                               || (numberRows * 20 > numberColumns && !plusMinus))
    1736                          method = ClpSolve::usePrimal; // switch off sprint
    1737                } else {
    1738                     // larger problem
    1739                     if(numberRows * 8 > numberColumns)
    1740                          method = ClpSolve::usePrimal; // switch off sprint
    1741                     // but make lightweight
    1742                     if(numberRows * 10 > numberColumns || numberColumns < 6000
    1743                               || (numberRows * 20 > numberColumns && !plusMinus))
    1744                          maxSprintPass = 10;
    1745                }
    1746           } else if (doSprint == 0) {
    1747                method = ClpSolve::usePrimal; // switch off sprint
     2077          if (ratio < 3.0) {
     2078            nPasses = static_cast< int >(ratio * static_cast< double >(nPasses) / 4.0); // probably not worth it
     2079          } else {
     2080            nPasses = CoinMax(nPasses, 5);
    17482081          }
    1749      }
    1750      if (method == ClpSolve::useDual) {
     2082          if (numberRows > 25000 && nPasses > 5) {
     2083            // Might as well go for it
     2084            nPasses = CoinMax(nPasses, 71);
     2085          } else if (increaseSprint) {
     2086            nPasses *= 2;
     2087            nPasses = CoinMin(nPasses, 71);
     2088          } else if (nPasses == 5 && ratio > 5.0) {
     2089            nPasses = static_cast< int >(static_cast< double >(nPasses) * (ratio / 5.0)); // increase if lots of elements per column
     2090          }
     2091          if (nPasses <= 5 && !plusMinus)
     2092            nPasses = 0;
     2093          //info.setStartingWeight(1.0e-1);
     2094        }
     2095        if (tryIt == 1) {
     2096          idiotOptions |= 262144;
     2097          info.setStrategy(idiotOptions | info.getStrategy());
     2098          //model2->setSpecialOptions(model2->specialOptions()
     2099          //                    |8388608);
     2100        }
     2101      }
     2102      if (doIdiot > 0) {
     2103        // pick up number passes
     2104        nPasses = options.getExtraInfo(1) % 1000000;
     2105#ifdef COIN_HAS_VOL
     2106        int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
     2107        nPasses = 0;
     2108        if (!returnCode) {
     2109          time2 = CoinCpuTime();
     2110          timeIdiot = time2 - timeX;
     2111          handler_->message(CLP_INTERVAL_TIMING, messages_)
     2112            << "Idiot Crash" << timeIdiot << time2 - time1
     2113            << CoinMessageEol;
     2114          timeX = time2;
     2115        }
     2116#endif
    17512117#ifdef CLP_USEFUL_PRINTOUT
    1752        debugInt[6]=1;
    1753 #endif
    1754           double * saveLower = NULL;
    1755           double * saveUpper = NULL;
    1756           if (presolve == ClpSolve::presolveOn) {
    1757                int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0);
    1758                if (numberInfeasibilities) {
    1759                     handler_->message(CLP_INFEASIBLE, messages_)
    1760                               << CoinMessageEol;
    1761                     delete model2;
    1762                     model2 = this;
    1763                     presolve = ClpSolve::presolveOff;
    1764                }
    1765           } else if (numberRows_ + numberColumns_ > 5000) {
    1766                // do anyway
    1767                saveLower = new double[numberRows_+numberColumns_];
    1768                CoinMemcpyN(model2->columnLower(), numberColumns_, saveLower);
    1769                CoinMemcpyN(model2->rowLower(), numberRows_, saveLower + numberColumns_);
    1770                saveUpper = new double[numberRows_+numberColumns_];
    1771                CoinMemcpyN(model2->columnUpper(), numberColumns_, saveUpper);
    1772                CoinMemcpyN(model2->rowUpper(), numberRows_, saveUpper + numberColumns_);
    1773                int numberInfeasibilities = model2->tightenPrimalBounds();
    1774                if (numberInfeasibilities) {
    1775                     handler_->message(CLP_INFEASIBLE, messages_)
    1776                               << CoinMessageEol;
    1777                     CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
    1778                     CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
    1779                     delete [] saveLower;
    1780                     saveLower = NULL;
    1781                     CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
    1782                     CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
    1783                     delete [] saveUpper;
    1784                     saveUpper = NULL;
    1785                     // return if wanted
    1786                     if (options.infeasibleReturn() ||
    1787                         (moreSpecialOptions_ & 1) != 0)
    1788                       return -1;
    1789                }
     2118        debugInt[6] = 4;
     2119        debugInt[7] = nPasses;
     2120#endif
     2121        if (nPasses > 70) {
     2122          info.setStartingWeight(1.0e3);
     2123          info.setReduceIterations(6);
     2124          //if (nPasses > 200)
     2125          //info.setFeasibilityTolerance(1.0e-9);
     2126          //if (nPasses > 1900)
     2127          //info.setWeightFactor(0.93);
     2128          if (nPasses > 900) {
     2129            double reductions = nPasses / 6.0;
     2130            if (nPasses < 5000) {
     2131              reductions /= 12.0;
     2132            } else {
     2133              reductions /= 13.0;
     2134              info.setStartingWeight(1.0e4);
     2135            }
     2136            double ratio = 1.0 / std::pow(10.0, (1.0 / reductions));
     2137            printf("%d passes reduction factor %g\n", nPasses, ratio);
     2138            info.setWeightFactor(ratio);
     2139          } else if (nPasses > 500) {
     2140            info.setWeightFactor(0.7);
     2141          } else if (nPasses > 200) {
     2142            info.setWeightFactor(0.5);
    17902143          }
    1791 #ifndef COIN_HAS_VOL
    1792           // switch off idiot and volume for now
    1793           doIdiot = 0;
    1794 #endif
    1795           // pick up number passes
    1796           int nPasses = 0;
    1797           int numberNotE = 0;
    1798 #ifndef SLIM_CLP
    1799           if ((doIdiot < 0 && plusMinus) || doIdiot > 0) {
    1800                // See if candidate for idiot
    1801                nPasses = 0;
    1802                Idiot info(*model2);
    1803                info.setStrategy(idiotOptions | info.getStrategy());
    1804                // Get average number of elements per column
    1805                double ratio  = static_cast<double> (numberElements) / static_cast<double> (numberColumns);
    1806                // look at rhs
    1807                int iRow;
    1808                double largest = 0.0;
    1809                double smallest = 1.0e30;
    1810                for (iRow = 0; iRow < numberRows; iRow++) {
    1811                     double value1 = model2->rowLower_[iRow];
    1812                     if (value1 && value1 > -1.0e31) {
    1813                          largest = CoinMax(largest, fabs(value1));
    1814                          smallest = CoinMin(smallest, fabs(value1));
    1815                     }
    1816                     double value2 = model2->rowUpper_[iRow];
    1817                     if (value2 && value2 < 1.0e31) {
    1818                          largest = CoinMax(largest, fabs(value2));
    1819                          smallest = CoinMin(smallest, fabs(value2));
    1820                     }
    1821                     if (value2 > value1) {
    1822                          numberNotE++;
    1823                     }
    1824                }
    1825                if (doIdiot < 0) {
    1826                     if (numberRows > 200 && numberColumns > 5000 && ratio >= 3.0 &&
    1827                               largest / smallest < 1.1 && !numberNotE) {
    1828                          nPasses = 71;
    1829                     }
    1830                }
    1831                if (doIdiot > 0) {
    1832                     nPasses = CoinMax(nPasses, doIdiot);
    1833                     if (nPasses > 70) {
    1834                          info.setStartingWeight(1.0e3);
    1835                          info.setDropEnoughFeasibility(0.01);
    1836                     }
    1837                }
    1838                if (nPasses > 20) {
    1839 #ifdef COIN_HAS_VOL
    1840                     int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
    1841                     if (!returnCode) {
    1842                          time2 = CoinCpuTime();
    1843                          timeIdiot = time2 - timeX;
    1844                          handler_->message(CLP_INTERVAL_TIMING, messages_)
    1845                                    << "Idiot Crash" << timeIdiot << time2 - time1
    1846                                    << CoinMessageEol;
    1847                          timeX = time2;
    1848                     } else {
    1849                          nPasses = 0;
    1850                     }
    1851 #else
    1852                     nPasses = 0;
    1853 #endif
    1854                } else {
    1855                     nPasses = 0;
    1856                }
     2144          if (maximumIterations() < nPasses) {
     2145            printf("Presuming maximumIterations is just for Idiot\n");
     2146            nPasses = maximumIterations();
     2147            setMaximumIterations(COIN_INT_MAX);
     2148            model2->setMaximumIterations(COIN_INT_MAX);
    18572149          }
    1858 #endif
    1859           if (doCrash) {
    1860 #ifdef ABC_INHERIT
    1861             if (!model2->abcState()) {
    1862 #endif
    1863                switch(doCrash) {
    1864                     // standard
    1865                case 1:
    1866                     model2->crash(1000, 1);
    1867                     break;
    1868                     // As in paper by Solow and Halim (approx)
    1869                case 2:
    1870                case 3:
    1871                     model2->crash(model2->dualBound(), 0);
    1872                     break;
    1873                     // Just put free in basis
    1874                case 4:
    1875                     model2->crash(0.0, 3);
    1876                     break;
    1877                }
    1878 #ifdef ABC_INHERIT
    1879             } else if (doCrash>=0) {
    1880                model2->setAbcState(model2->abcState()|256*doCrash);
    1881             }
    1882 #endif
     2150          if (nPasses >= 10000 && nPasses < 100000) {
     2151            int k = nPasses % 100;
     2152            nPasses /= 200;
     2153            info.setReduceIterations(3);
     2154            if (k)
     2155              info.setStartingWeight(1.0e2);
    18832156          }
    1884           if (!nPasses) {
    1885                int saveOptions = model2->specialOptions();
    1886                if (model2->numberRows() > 100)
    1887                     model2->setSpecialOptions(saveOptions | 64); // go as far as possible
    1888                //int numberRows = model2->numberRows();
    1889                //int numberColumns = model2->numberColumns();
    1890                if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
    1891                     // See if original wanted vector
    1892                     ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
    1893                     ClpMatrixBase * matrix = model2->clpMatrix();
    1894                     if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
    1895                          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    1896                          clpMatrix->makeSpecialColumnCopy();
    1897                          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
    1898                          model2->dual(0);
    1899                          clpMatrix->releaseSpecialColumnCopy();
    1900                     } else {
    1901 #ifndef ABC_INHERIT
    1902                       model2->dual(0);
    1903 #else
    1904                       model2->dealWithAbc(0,0,interrupt);
    1905 #endif
    1906                     }
    1907                } else {
    1908                     model2->dual(0);
    1909                }
    1910           } else if (!numberNotE && 0) {
    1911                // E so we can do in another way
    1912                double * pi = model2->dualRowSolution();
    1913                int i;
    1914                int numberColumns = model2->numberColumns();
    1915                int numberRows = model2->numberRows();
    1916                double * saveObj = new double[numberColumns];
    1917                CoinMemcpyN(model2->objective(), numberColumns, saveObj);
    1918                CoinMemcpyN(model2->objective(),
    1919                            numberColumns, model2->dualColumnSolution());
    1920                model2->clpMatrix()->transposeTimes(-1.0, pi, model2->dualColumnSolution());
    1921                CoinMemcpyN(model2->dualColumnSolution(),
    1922                            numberColumns, model2->objective());
    1923                const double * rowsol = model2->primalRowSolution();
    1924                double offset = 0.0;
    1925                for (i = 0; i < numberRows; i++) {
    1926                     offset += pi[i] * rowsol[i];
    1927                }
    1928                double value2;
    1929                model2->getDblParam(ClpObjOffset, value2);
    1930                //printf("Offset %g %g\n",offset,value2);
    1931                model2->setDblParam(ClpObjOffset, value2 - offset);
    1932                model2->setPerturbation(51);
    1933                //model2->setRowObjective(pi);
    1934                // zero out pi
    1935                //memset(pi,0,numberRows*sizeof(double));
    1936                // Could put some in basis - only partially tested
    1937                model2->allSlackBasis();
    1938                //model2->factorization()->maximumPivots(200);
    1939                //model2->setLogLevel(63);
    1940                // solve
    1941                model2->dual(0);
    1942                model2->setDblParam(ClpObjOffset, value2);
    1943                CoinMemcpyN(saveObj, numberColumns, model2->objective());
    1944                // zero out pi
    1945                //memset(pi,0,numberRows*sizeof(double));
    1946                //model2->setRowObjective(pi);
    1947                delete [] saveObj;
    1948                //model2->dual(0);
    1949                model2->setPerturbation(50);
    1950                model2->primal();
    1951           } else {
    1952                // solve
    1953                model2->setPerturbation(100);
    1954                model2->dual(2);
    1955                model2->setPerturbation(50);
    1956                model2->dual(0);
    1957           }
    1958           if (saveLower) {
    1959                CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
    1960                CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
    1961                delete [] saveLower;
    1962                saveLower = NULL;
    1963                CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
    1964                CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
    1965                delete [] saveUpper;
    1966                saveUpper = NULL;
    1967           }
    1968           time2 = CoinCpuTime();
    1969           timeCore = time2 - timeX;
    1970           handler_->message(CLP_INTERVAL_TIMING, messages_)
    1971                     << "Dual" << timeCore << time2 - time1
    1972                     << CoinMessageEol;
    1973           timeX = time2;
    1974      } else if (method == ClpSolve::usePrimal) {
    1975 #ifdef CLP_USEFUL_PRINTOUT
    1976        debugInt[6]=2;
    1977 #endif
    1978 #ifndef SLIM_CLP
    1979           if (doIdiot) {
    1980                int nPasses = 0;
    1981                Idiot info(*model2);
    1982                info.setStrategy(idiotOptions | info.getStrategy());
    1983                // Get average number of elements per column
    1984                double ratio  = static_cast<double> (numberElements) / static_cast<double> (numberColumns);
    1985                // look at rhs
    1986                int iRow;
    1987                double largest = 0.0;
    1988                double smallest = 1.0e30;
    1989                double largestGap = 0.0;
    1990                int numberNotE = 0;
    1991                for (iRow = 0; iRow < numberRows; iRow++) {
    1992                     double value1 = model2->rowLower_[iRow];
    1993                     if (value1 && value1 > -1.0e31) {
    1994                          largest = CoinMax(largest, fabs(value1));
    1995                          smallest = CoinMin(smallest, fabs(value1));
    1996                     }
    1997                     double value2 = model2->rowUpper_[iRow];
    1998                     if (value2 && value2 < 1.0e31) {
    1999                          largest = CoinMax(largest, fabs(value2));
    2000                          smallest = CoinMin(smallest, fabs(value2));
    2001                     }
    2002                     if (value2 > value1) {
    2003                          numberNotE++;
    2004                          if (value2 > 1.0e31 || value1 < -1.0e31)
    2005                               largestGap = COIN_DBL_MAX;
    2006                          else
    2007                               largestGap = value2 - value1;
    2008                     }
    2009                }
    2010                bool increaseSprint = plusMinus;
    2011                if ((specialOptions_ & 1024) != 0)
    2012                     increaseSprint = false;
    2013                if (!plusMinus) {
    2014                     // If 90% +- 1 then go for sprint
    2015                     if (statistics[0] >= 0 && 10 * statistics[2] < statistics[0] + statistics[1])
    2016                          increaseSprint = true;
    2017                }
    2018                int tryIt = tryItSave;
    2019                if (numberRows < 1000 && numberColumns < 3000)
    2020                     tryIt = 0;
    2021                if (tryIt) {
    2022                     if (increaseSprint) {
    2023                          info.setStartingWeight(1.0e3);
    2024                          info.setReduceIterations(6);
    2025                          // also be more lenient on infeasibilities
    2026                          info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
    2027                          info.setDropEnoughWeighted(-2.0);
    2028                          if (largest / smallest > 2.0) {
    2029                               nPasses = 10 + numberColumns / 100000;
    2030                               nPasses = CoinMin(nPasses, 50);
    2031                               nPasses = CoinMax(nPasses, 15);
    2032                               if (numberRows > 20000 && nPasses > 5) {
    2033                                    // Might as well go for it
    2034                                    nPasses = CoinMax(nPasses, 71);
    2035                               } else if (numberRows > 2000 && nPasses > 5) {
    2036                                    nPasses = CoinMax(nPasses, 50);
    2037                               } else if (numberElements < 3 * numberColumns) {
    2038                                    nPasses = CoinMin(nPasses, 10); // probably not worh it
    2039                                    if (doIdiot < 0)
    2040                                         info.setLightweight(1); // say lightweight idiot
    2041                               } else {
    2042                                    if (doIdiot < 0)
    2043                                         info.setLightweight(1); // say lightweight idiot
    2044                               }
    2045                          } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) {
    2046                               nPasses = 10 + numberColumns / 1000;
    2047                               nPasses = CoinMin(nPasses, 100);
    2048                               nPasses = CoinMax(nPasses, 30);
    2049                               if (numberRows > 25000) {
    2050                                    // Might as well go for it
    2051                                    nPasses = CoinMax(nPasses, 71);
    2052                               }
    2053                               if (!largestGap)
    2054                                    nPasses *= 2;
    2055                          } else {
    2056                               nPasses = 10 + numberColumns / 1000;
    2057                               nPasses = CoinMin(nPasses, 200);
    2058                               nPasses = CoinMax(nPasses, 100);
    2059                               info.setStartingWeight(1.0e-1);
    2060                               info.setReduceIterations(6);
    2061                               if (!largestGap && nPasses <= 50)
    2062                                    nPasses *= 2;
    2063                               //info.setFeasibilityTolerance(1.0e-7);
    2064                          }
    2065                          // If few passes - don't bother
    2066                          if (nPasses <= 5 && !plusMinus)
    2067                               nPasses = 0;
    2068                     } else {
    2069                          if (doIdiot < 0)
    2070                               info.setLightweight(1); // say lightweight idiot
    2071                          if (largest / smallest > 1.01 || numberNotE || statistics[2] > statistics[0] + statistics[1]) {
    2072                               if (numberRows > 25000 || numberColumns > 5 * numberRows) {
    2073                                    nPasses = 50;
    2074                               } else if (numberColumns > 4 * numberRows) {
    2075                                    nPasses = 20;
    2076                               } else {
    2077                                    nPasses = 5;
    2078                               }
    2079                          } else {
    2080                               if (numberRows > 25000 || numberColumns > 5 * numberRows) {
    2081                                    nPasses = 50;
    2082                                    info.setLightweight(0); // say not lightweight idiot
    2083                               } else if (numberColumns > 4 * numberRows) {
    2084                                    nPasses = 20;
    2085                               } else {
    2086                                    nPasses = 15;
    2087                               }
    2088                          }
    2089                          if (ratio < 3.0) {
    2090                               nPasses = static_cast<int> (ratio * static_cast<double> (nPasses) / 4.0); // probably not worth it
    2091                          } else {
    2092                               nPasses = CoinMax(nPasses, 5);
    2093                          }
    2094                          if (numberRows > 25000 && nPasses > 5) {
    2095                               // Might as well go for it
    2096                               nPasses = CoinMax(nPasses, 71);
    2097                          } else if (increaseSprint) {
    2098                               nPasses *= 2;
    2099                               nPasses = CoinMin(nPasses, 71);
    2100                          } else if (nPasses == 5 && ratio > 5.0) {
    2101                               nPasses = static_cast<int> (static_cast<double> (nPasses) * (ratio / 5.0)); // increase if lots of elements per column
    2102                          }
    2103                          if (nPasses <= 5 && !plusMinus)
    2104                               nPasses = 0;
    2105                          //info.setStartingWeight(1.0e-1);
    2106                     }
    2107                     if (tryIt==1) {
    2108                       idiotOptions |= 262144;
    2109                       info.setStrategy(idiotOptions | info.getStrategy());
    2110                       //model2->setSpecialOptions(model2->specialOptions()
    2111                       //                        |8388608);
    2112                     }
    2113                }
    2114                if (doIdiot > 0) {
    2115                     // pick up number passes
    2116                     nPasses = options.getExtraInfo(1) % 1000000;
    2117 #ifdef COIN_HAS_VOL
    2118                     int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
    2119                     nPasses=0;
    2120                     if (!returnCode) {
    2121                       time2 = CoinCpuTime();
    2122                       timeIdiot = time2 - timeX;
    2123                       handler_->message(CLP_INTERVAL_TIMING, messages_)
    2124                         << "Idiot Crash" << timeIdiot << time2 - time1
    2125                         << CoinMessageEol;
    2126                       timeX = time2;
    2127                     }
    2128 #endif
    2129 #ifdef CLP_USEFUL_PRINTOUT
    2130                     debugInt[6]=4;
    2131                     debugInt[7]=nPasses;
    2132 #endif
    2133                     if (nPasses > 70) {
    2134                          info.setStartingWeight(1.0e3);
    2135                          info.setReduceIterations(6);
    2136                          //if (nPasses > 200)
    2137                          //info.setFeasibilityTolerance(1.0e-9);
    2138                          //if (nPasses > 1900)
    2139                          //info.setWeightFactor(0.93);
    2140                          if (nPasses > 900) {
    2141                            double reductions=nPasses/6.0;
    2142                            if (nPasses<5000) {
    2143                              reductions /= 12.0;
    2144                            } else {
    2145                              reductions /= 13.0;
    2146                              info.setStartingWeight(1.0e4);
    2147                            }
    2148                            double ratio=1.0/std::pow(10.0,(1.0/reductions));
    2149                            printf("%d passes reduction factor %g\n",nPasses,ratio);
    2150                            info.setWeightFactor(ratio);
    2151                          } else if (nPasses > 500) {
    2152                            info.setWeightFactor(0.7);
    2153                          } else if (nPasses > 200) {
    2154                            info.setWeightFactor(0.5);
    2155                          }
    2156                          if (maximumIterations()<nPasses) {
    2157                            printf("Presuming maximumIterations is just for Idiot\n");
    2158                            nPasses=maximumIterations();
    2159                            setMaximumIterations(COIN_INT_MAX);
    2160                            model2->setMaximumIterations(COIN_INT_MAX);
    2161                          }
    2162                          if (nPasses >= 10000&&nPasses<100000) {
    2163                               int k = nPasses % 100;
    2164                               nPasses /= 200;
    2165                               info.setReduceIterations(3);
    2166                               if (k)
    2167                                    info.setStartingWeight(1.0e2);
    2168                          }
    2169                          // also be more lenient on infeasibilities
    2170                          info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
    2171                          info.setDropEnoughWeighted(-2.0);
    2172                     } else if (nPasses >= 50) {
    2173                          info.setStartingWeight(1.0e3);
    2174                          //info.setReduceIterations(6);
    2175                     }
    2176                     // For experimenting
    2177                     if (nPasses < 70 && (nPasses % 10) > 0 && (nPasses % 10) < 4) {
    2178                          info.setStartingWeight(1.0e3);
    2179                          info.setLightweight(nPasses % 10); // special testing
     2157          // also be more lenient on infeasibilities
     2158          info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
     2159          info.setDropEnoughWeighted(-2.0);
     2160        } else if (nPasses >= 50) {
     2161          info.setStartingWeight(1.0e3);
     2162          //info.setReduceIterations(6);
     2163        }
     2164        // For experimenting
     2165        if (nPasses < 70 && (nPasses % 10) > 0 && (nPasses % 10) < 4) {
     2166          info.setStartingWeight(1.0e3);
     2167          info.setLightweight(nPasses % 10); // special testing
    21802168#ifdef COIN_DEVELOP
    2181                          printf("warning - odd lightweight %d\n", nPasses % 10);
    2182                          //info.setReduceIterations(6);
    2183 #endif
    2184                     }
    2185                }
    2186                if (options.getExtraInfo(1) > 1000000)
    2187                     nPasses += 1000000;
    2188                if (nPasses) {
    2189                     doCrash = 0;
     2169          printf("warning - odd lightweight %d\n", nPasses % 10);
     2170          //info.setReduceIterations(6);
     2171#endif
     2172        }
     2173      }
     2174      if (options.getExtraInfo(1) > 1000000)
     2175        nPasses += 1000000;
     2176      if (nPasses) {
     2177        doCrash = 0;
    21902178#if 0
    21912179                    double * solution = model2->primalColumnSolution();
     
    22132201                    delete [] saveUpper;
    22142202#else
    2215                     // Allow for crossover
    2216                     //#define LACI_TRY
     2203        // Allow for crossover
     2204        //#define LACI_TRY
    22172205#ifndef LACI_TRY
    2218                     //if (doIdiot>0)
     2206        //if (doIdiot>0)
    22192207#if 0 //def ABC_INHERIT
    2220                     if (!model2->abcState()) 
    2221 #endif
    2222                     info.setStrategy(512 | info.getStrategy());
    2223 #endif
    2224                     // Allow for scaling
    2225                     info.setStrategy(32 | info.getStrategy());
    2226                     int saveScalingFlag=model2->scalingFlag();
    2227                     bool linearObjective = objective_->type()<2;
    2228                     if (!linearObjective)
    2229                       model2->scaling(0);
     2208                    if (!model2->abcState())
     2209#endif
     2210        info.setStrategy(512 | info.getStrategy());
     2211#endif
     2212        // Allow for scaling
     2213        info.setStrategy(32 | info.getStrategy());
     2214        int saveScalingFlag = model2->scalingFlag();
     2215        bool linearObjective = objective_->type() < 2;
     2216        if (!linearObjective)
     2217          model2->scaling(0);
    22302218#define CLP_DUAL_IDIOT
    22312219#ifdef CLP_DUAL_IDIOT
    2232                     bool doubleIdiot=false;
    2233                     if (nPasses==99&&linearObjective)
    2234                       doubleIdiot=true;
    2235                     if (doubleIdiot) {
    2236                       ClpSimplex * dualModel2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(1.0,1.0);
    2237                       if (dualModel2) {
    2238                         printf("Dual of model has %d rows and %d columns\n",
    2239                                dualModel2->numberRows(), dualModel2->numberColumns());
    2240                         dualModel2->setOptimizationDirection(1.0);
    2241                         Idiot infoDual(info);
    2242                         infoDual.setModel(dualModel2);
     2220        bool doubleIdiot = false;
     2221        if (nPasses == 99 && linearObjective)
     2222          doubleIdiot = true;
     2223        if (doubleIdiot) {
     2224          ClpSimplex *dualModel2 = static_cast< ClpSimplexOther * >(model2)->dualOfModel(1.0, 1.0);
     2225          if (dualModel2) {
     2226            printf("Dual of model has %d rows and %d columns\n",
     2227              dualModel2->numberRows(), dualModel2->numberColumns());
     2228            dualModel2->setOptimizationDirection(1.0);
     2229            Idiot infoDual(info);
     2230            infoDual.setModel(dualModel2);
    22432231#if 0
    22442232                        info.setStrategy(512 | info.getStrategy());
     
    22482236                        info.setReduceIterations(6);
    22492237#endif
    2250                         info.crash(nPasses, model2->messageHandler(),
    2251                                    model2->messagesPointer(),false);
    2252                         infoDual.crash(nPasses, model2->messageHandler(),
    2253                                    model2->messagesPointer(),false);
    2254                         // two copies of solutions
    2255                         ClpSimplex temp(*model2);
     2238            info.crash(nPasses, model2->messageHandler(),
     2239              model2->messagesPointer(), false);
     2240            infoDual.crash(nPasses, model2->messageHandler(),
     2241              model2->messagesPointer(), false);
     2242            // two copies of solutions
     2243            ClpSimplex temp(*model2);
    22562244#if 0
    22572245                        static_cast<ClpSimplexOther *> (&temp)->restoreFromDual(dualModel2);
    22582246#else
    2259                         // move duals and just copy primal
    2260                         memcpy(temp.dualRowSolution(),dualModel2->primalColumnSolution(),
    2261                                numberRows*sizeof(double));
    2262                         memcpy(temp.primalColumnSolution(),model2->primalColumnSolution(),
    2263                                numberColumns*sizeof(double));
    2264 #endif
    2265                         delete dualModel2;
    2266                         int numberRows=model2->numberRows();
    2267                         int numberColumns=model2->numberColumns();
    2268                         ClpSimplex * tempModel[2];
    2269                         tempModel[0]=model2;
    2270                         tempModel[1]=&temp;
    2271                         const double * primalColumn[2];
    2272                         const double * dualRow[2];
    2273                         double * dualColumn[2];
    2274                         double * primalRow[2];
    2275                         for (int i=0;i<2;i++) {
    2276                           primalColumn[i]=tempModel[i]->primalColumnSolution();
    2277                           dualRow[i]=tempModel[i]->dualRowSolution();
    2278                           dualColumn[i]=tempModel[i]->dualColumnSolution();
    2279                           primalRow[i]=tempModel[i]->primalRowSolution();
    2280                           memcpy(dualColumn[i],model2->objective(),
    2281                                  numberColumns*sizeof(double));
    2282                           memset(primalRow[i], 0, numberRows * sizeof(double));
    2283                           tempModel[i]->clpMatrix()->transposeTimes(-1.0,
    2284                                                               dualRow[i],
    2285                                                               dualColumn[i]);
    2286                           tempModel[i]->clpMatrix()->times(1.0,
    2287                                                      primalColumn[i],
    2288                                                      primalRow[i]);
    2289                           tempModel[i]->checkSolutionInternal();
    2290                           printf("model %d - dual inf %g primal inf %g\n",
    2291                                  i,tempModel[i]->sumDualInfeasibilities(),
    2292                                 tempModel[i]->sumPrimalInfeasibilities());
    2293                         }
    2294                         printf("What now\n");
    2295                       } else {
    2296                         doubleIdiot=false;
    2297                       }
    2298                     }
    2299                     if (!doubleIdiot)
    2300                       info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),
    2301                                 (objective_->type() <2));
     2247            // move duals and just copy primal
     2248            memcpy(temp.dualRowSolution(), dualModel2->primalColumnSolution(),
     2249              numberRows * sizeof(double));
     2250            memcpy(temp.primalColumnSolution(), model2->primalColumnSolution(),
     2251              numberColumns * sizeof(double));
     2252#endif
     2253            delete dualModel2;
     2254            int numberRows = model2->numberRows();
     2255            int numberColumns = model2->numberColumns();
     2256            ClpSimplex *tempModel[2];
     2257            tempModel[0] = model2;
     2258            tempModel[1] = &temp;
     2259            const double *primalColumn[2];
     2260            const double *dualRow[2];
     2261            double *dualColumn[2];
     2262            double *primalRow[2];
     2263            for (int i = 0; i < 2; i++) {
     2264              primalColumn[i] = tempModel[i]->primalColumnSolution();
     2265              dualRow[i] = tempModel[i]->dualRowSolution();
     2266              dualColumn[i] = tempModel[i]->dualColumnSolution();
     2267              primalRow[i] = tempModel[i]->primalRowSolution();
     2268              memcpy(dualColumn[i], model2->objective(),
     2269                numberColumns * sizeof(double));
     2270              memset(primalRow[i], 0, numberRows * sizeof(double));
     2271              tempModel[i]->clpMatrix()->transposeTimes(-1.0,
     2272                dualRow[i],
     2273                dualColumn[i]);
     2274              tempModel[i]->clpMatrix()->times(1.0,
     2275                primalColumn[i],
     2276                primalRow[i]);
     2277              tempModel[i]->checkSolutionInternal();
     2278              printf("model %d - dual inf %g primal inf %g\n",
     2279                i, tempModel[i]->sumDualInfeasibilities(),
     2280                tempModel[i]->sumPrimalInfeasibilities());
     2281            }
     2282            printf("What now\n");
     2283          } else {
     2284            doubleIdiot = false;
     2285          }
     2286        }
     2287        if (!doubleIdiot)
     2288          info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),
     2289            (objective_->type() < 2));
    23022290#else
    2303                     info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),(objective_->type() <2));
    2304 #endif
    2305                       model2->scaling(saveScalingFlag);
    2306 #endif
    2307                     time2 = CoinCpuTime();
    2308                     timeIdiot = time2 - timeX;
    2309                     handler_->message(CLP_INTERVAL_TIMING, messages_)
    2310                               << "Idiot Crash" << timeIdiot << time2 - time1
    2311                               << CoinMessageEol;
    2312                     timeX = time2;
    2313                     if (nPasses>100000&&nPasses<100500) {
    2314                       // make sure no status left
    2315                       model2->createStatus();
    2316                       // solve
    2317                       if (model2->factorizationFrequency() == 200) {
    2318                         // User did not touch preset
    2319                         model2->defaultFactorizationFrequency();
    2320                       }
    2321                       //int numberRows = model2->numberRows();
    2322                       int numberColumns = model2->numberColumns();
    2323                       // save duals
    2324                       //double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows);
    2325                       // for moment this only works on nug etc (i.e. all ==)
    2326                       // needs row objective
    2327                       double * saveObj = CoinCopyOfArray(model2->objective(),numberColumns);
    2328                       double * pi = model2->dualRowSolution();
    2329                       model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective());
    2330                       // just primal values pass
    2331                       double saveScale = model2->objectiveScale();
    2332                       model2->setObjectiveScale(1.0e-3);
    2333                       model2->primal(2);
    2334                       model2->writeMps("xx.mps");
    2335                       double * solution = model2->primalColumnSolution();
    2336                       double * upper = model2->columnUpper();
    2337                       for (int i=0;i<numberColumns;i++) {
    2338                         if (solution[i]<100.0)
    2339                           upper[i]=1000.0;
    2340                       }
    2341                       model2->setProblemStatus(-1);
    2342                       model2->setObjectiveScale(saveScale);
     2291        info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(), (objective_->type() < 2));
     2292#endif
     2293        model2->scaling(saveScalingFlag);
     2294#endif
     2295        time2 = CoinCpuTime();
     2296        timeIdiot = time2 - timeX;
     2297        handler_->message(CLP_INTERVAL_TIMING, messages_)
     2298          << "Idiot Crash" << timeIdiot << time2 - time1
     2299          << CoinMessageEol;
     2300        timeX = time2;
     2301        if (nPasses > 100000 && nPasses < 100500) {
     2302          // make sure no status left
     2303          model2->createStatus();
     2304          // solve
     2305          if (model2->factorizationFrequency() == 200) {
     2306            // User did not touch preset
     2307            model2->defaultFactorizationFrequency();
     2308          }
     2309          //int numberRows = model2->numberRows();
     2310          int numberColumns = model2->numberColumns();
     2311          // save duals
     2312          //double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows);
     2313          // for moment this only works on nug etc (i.e. all ==)
     2314          // needs row objective
     2315          double *saveObj = CoinCopyOfArray(model2->objective(), numberColumns);
     2316          double *pi = model2->dualRowSolution();
     2317          model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective());
     2318          // just primal values pass
     2319          double saveScale = model2->objectiveScale();
     2320          model2->setObjectiveScale(1.0e-3);
     2321          model2->primal(2);
     2322          model2->writeMps("xx.mps");
     2323          double *solution = model2->primalColumnSolution();
     2324          double *upper = model2->columnUpper();
     2325          for (int i = 0; i < numberColumns; i++) {
     2326            if (solution[i] < 100.0)
     2327              upper[i] = 1000.0;
     2328          }
     2329          model2->setProblemStatus(-1);
     2330          model2->setObjectiveScale(saveScale);
    23432331#ifdef ABC_INHERIT
    2344                       AbcSimplex * abcModel2=new AbcSimplex(*model2);
    2345                       if (interrupt)
    2346                         currentAbcModel = abcModel2;
    2347                       if (abcSimplex_) {
    2348                         // move factorization stuff
    2349                         abcModel2->factorization()->synchronize(model2->factorization(),abcModel2);
    2350                       }
    2351                       abcModel2->startPermanentArrays();
    2352                       abcModel2->setAbcState(CLP_ABC_WANTED);
     2332          AbcSimplex *abcModel2 = new AbcSimplex(*model2);
     2333          if (interrupt)
     2334            currentAbcModel = abcModel2;
     2335          if (abcSimplex_) {
     2336            // move factorization stuff
     2337            abcModel2->factorization()->synchronize(model2->factorization(), abcModel2);
     2338          }
     2339          abcModel2->startPermanentArrays();
     2340          abcModel2->setAbcState(CLP_ABC_WANTED);
    23532341#if ABC_PARALLEL
    2354                       int parallelMode=1;
    2355                       printf("Parallel mode %d\n",parallelMode);
    2356                       abcModel2->setParallelMode(parallelMode);
    2357 #endif
    2358                       //if (processTune>0&&processTune<8)
    2359                       //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|65536*processTune);
    2360                       abcModel2->doAbcDual();
    2361                       abcModel2->moveStatusToClp(model2);
    2362                       //ClpModel::stopPermanentArrays();
    2363                       model2->setSpecialOptions(model2->specialOptions()&~65536);
    2364                       //model2->dual();
    2365                       //model2->setNumberIterations(abcModel2->numberIterations()+model2->numberIterations());
    2366                       delete abcModel2;
    2367 #endif
    2368                       memcpy(model2->objective(),saveObj,numberColumns*sizeof(double));
    2369                       //delete [] saveDuals;
    2370                       delete [] saveObj;
    2371                       model2->dual(2);
    2372                     } // end dubious idiot
    2373                }
    2374           }
    2375 #endif
    2376           // ?
    2377           if (doCrash) {
    2378                switch(doCrash) {
    2379                     // standard
    2380                case 1:
    2381                     model2->crash(1000, 1);
    2382                     break;
    2383                     // As in paper by Solow and Halim (approx)
    2384                case 2:
    2385                     model2->crash(model2->dualBound(), 0);
    2386                     break;
    2387                     // My take on it
    2388                case 3:
    2389                     model2->crash(model2->dualBound(), -1);
    2390                     break;
    2391                     // Just put free in basis
    2392                case 4:
    2393                     model2->crash(0.0, 3);
    2394                     break;
    2395                }
    2396           }
     2342          int parallelMode = 1;
     2343          printf("Parallel mode %d\n", parallelMode);
     2344          abcModel2->setParallelMode(parallelMode);
     2345#endif
     2346          //if (processTune>0&&processTune<8)
     2347          //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|65536*processTune);
     2348          abcModel2->doAbcDual();
     2349          abcModel2->moveStatusToClp(model2);
     2350          //ClpModel::stopPermanentArrays();
     2351          model2->setSpecialOptions(model2->specialOptions() & ~65536);
     2352          //model2->dual();
     2353          //model2->setNumberIterations(abcModel2->numberIterations()+model2->numberIterations());
     2354          delete abcModel2;
     2355#endif
     2356          memcpy(model2->objective(), saveObj, numberColumns * sizeof(double));
     2357          //delete [] saveDuals;
     2358          delete[] saveObj;
     2359          model2->dual(2);
     2360        } // end dubious idiot
     2361      }
     2362    }
     2363#endif
     2364    // ?
     2365    if (doCrash) {
     2366      switch (doCrash) {
     2367        // standard
     2368      case 1:
     2369        model2->crash(1000, 1);
     2370        break;
     2371        // As in paper by Solow and Halim (approx)
     2372      case 2:
     2373        model2->crash(model2->dualBound(), 0);
     2374        break;
     2375        // My take on it
     2376      case 3:
     2377        model2->crash(model2->dualBound(), -1);
     2378        break;
     2379        // Just put free in basis
     2380      case 4:
     2381        model2->crash(0.0, 3);
     2382        break;
     2383      }
     2384    }
    23972385#ifndef SLIM_CLP
    2398           if (doSlp && objective_->type() == 2) {
    2399                model2->nonlinearSLP(doSlp, 1.0e-5);
    2400           }
     2386    if (doSlp && objective_->type() == 2) {
     2387      model2->nonlinearSLP(doSlp, 1.0e-5);
     2388    }
    24012389#endif
    24022390#ifndef LACI_TRY
    2403           if (options.getSpecialOption(1) != 2 ||
    2404                     options.getExtraInfo(1) < 1000000) {
    2405                if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
    2406                     // See if original wanted vector
    2407                     ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
    2408                     ClpMatrixBase * matrix = model2->clpMatrix();
    2409                     if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
    2410                          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    2411                          clpMatrix->makeSpecialColumnCopy();
    2412                          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
    2413                          model2->primal(primalStartup);
    2414                          clpMatrix->releaseSpecialColumnCopy();
    2415                     } else {
     2391    if (options.getSpecialOption(1) != 2 || options.getExtraInfo(1) < 1000000) {
     2392      if (dynamic_cast< ClpPackedMatrix * >(matrix_)) {
     2393        // See if original wanted vector
     2394        ClpPackedMatrix *clpMatrixO = dynamic_cast< ClpPackedMatrix * >(matrix_);
     2395        ClpMatrixBase *matrix = model2->clpMatrix();
     2396        if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
     2397          ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     2398          clpMatrix->makeSpecialColumnCopy();
     2399          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
     2400          model2->primal(primalStartup);
     2401          clpMatrix->releaseSpecialColumnCopy();
     2402        } else {
    24162403#ifndef ABC_INHERIT
    2417                         model2->primal(primalStartup);
     2404          model2->primal(primalStartup);
    24182405#else
    2419                         model2->dealWithAbc(1,primalStartup,interrupt);
    2420 #endif
    2421                     }
    2422                } else {
     2406          model2->dealWithAbc(1, primalStartup, interrupt);
     2407#endif
     2408        }
     2409      } else {
    24232410#ifndef ABC_INHERIT
    2424                     model2->primal(primalStartup);
     2411        model2->primal(primalStartup);
    24252412#else
    2426                     model2->dealWithAbc(1,primalStartup,interrupt);
    2427 #endif
    2428                }
    2429           }
    2430 #endif
    2431           time2 = CoinCpuTime();
    2432           timeCore = time2 - timeX;
    2433           handler_->message(CLP_INTERVAL_TIMING, messages_)
    2434                     << "Primal" << timeCore << time2 - time1
    2435                     << CoinMessageEol;
    2436           timeX = time2;
    2437      } else if (method == ClpSolve::usePrimalorSprint) {
    2438           // Sprint
    2439           /*
     2413        model2->dealWithAbc(1, primalStartup, interrupt);
     2414#endif
     2415      }
     2416    }
     2417#endif
     2418    time2 = CoinCpuTime();
     2419    timeCore = time2 - timeX;
     2420    handler_->message(CLP_INTERVAL_TIMING, messages_)
     2421      << "Primal" << timeCore << time2 - time1
     2422      << CoinMessageEol;
     2423    timeX = time2;
     2424  } else if (method == ClpSolve::usePrimalorSprint) {
     2425    // Sprint
     2426    /*
    24402427            This driver implements what I called Sprint when I introduced the idea
    24412428            many years ago.  Cplex calls it "sifting" which I think is just as silly.
     
    24492436          */
    24502437
    2451           /* The idea works best if you can get feasible easily.  To make it
     2438    /* The idea works best if you can get feasible easily.  To make it
    24522439             more general we can add in costed slacks */
    24532440
    2454           int originalNumberColumns = model2->numberColumns();
    2455           int numberRows = model2->numberRows();
    2456           ClpSimplex * originalModel2 = model2;
     2441    int originalNumberColumns = model2->numberColumns();
     2442    int numberRows = model2->numberRows();
     2443    ClpSimplex *originalModel2 = model2;
    24572444
    2458           // We will need arrays to choose variables.  These are too big but ..
    2459           double * weight = new double [numberRows+originalNumberColumns];
    2460           int * sort = new int [numberRows+originalNumberColumns];
    2461           int numberSort = 0;
    2462           // We are going to add slacks to get feasible.
    2463           // initial list will just be artificials
    2464           int iColumn;
    2465           const double * columnLower = model2->columnLower();
    2466           const double * columnUpper = model2->columnUpper();
    2467           double * columnSolution = model2->primalColumnSolution();
     2445    // We will need arrays to choose variables.  These are too big but ..
     2446    double *weight = new double[numberRows + originalNumberColumns];
     2447    int *sort = new int[numberRows + originalNumberColumns];
     2448    int numberSort = 0;
     2449    // We are going to add slacks to get feasible.
     2450    // initial list will just be artificials
     2451    int iColumn;
     2452    const double *columnLower = model2->columnLower();
     2453    const double *columnUpper = model2->columnUpper();
     2454    double *columnSolution = model2->primalColumnSolution();
    24682455
    2469           // See if we have costed slacks
    2470           int * negSlack = new int[numberRows+numberColumns];
    2471           int * posSlack = new int[numberRows+numberColumns];
    2472           int * nextNegSlack = negSlack+numberRows;
    2473           int * nextPosSlack = posSlack+numberRows;
    2474           int iRow;
    2475           for (iRow = 0; iRow < numberRows+numberColumns; iRow++) {
    2476                negSlack[iRow] = -1;
    2477                posSlack[iRow] = -1;
     2456    // See if we have costed slacks
     2457    int *negSlack = new int[numberRows + numberColumns];
     2458    int *posSlack = new int[numberRows + numberColumns];
     2459    int *nextNegSlack = negSlack + numberRows;
     2460    int *nextPosSlack = posSlack + numberRows;
     2461    int iRow;
     2462    for (iRow = 0; iRow < numberRows + numberColumns; iRow++) {
     2463      negSlack[iRow] = -1;
     2464      posSlack[iRow] = -1;
     2465    }
     2466    const double *element = model2->matrix()->getElements();
     2467    const int *row = model2->matrix()->getIndices();
     2468    const CoinBigIndex *columnStart = model2->matrix()->getVectorStarts();
     2469    const int *columnLength = model2->matrix()->getVectorLengths();
     2470    //bool allSlack = (numberRowsBasic==numberRows);
     2471    for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
     2472      if (!columnSolution[iColumn] || fabs(columnSolution[iColumn]) > 1.0e20) {
     2473        double value = 0.0;
     2474        if (columnLower[iColumn] > 0.0)
     2475          value = columnLower[iColumn];
     2476        else if (columnUpper[iColumn] < 0.0)
     2477          value = columnUpper[iColumn];
     2478        columnSolution[iColumn] = value;
     2479      }
     2480      if (columnLength[iColumn] == 1) {
     2481        int jRow = row[columnStart[iColumn]];
     2482        if (!columnLower[iColumn]) {
     2483          if (element[columnStart[iColumn]] > 0.0) {
     2484            if (posSlack[jRow] < 0) {
     2485              posSlack[jRow] = iColumn;
     2486            } else {
     2487              int jColumn = posSlack[jRow];
     2488              while (nextPosSlack[jColumn] >= 0)
     2489                jColumn = nextPosSlack[jColumn];
     2490              nextPosSlack[jColumn] = iColumn;
     2491            }
     2492          } else if (element[columnStart[iColumn]] < 0.0) {
     2493            if (negSlack[jRow] < 0) {
     2494              negSlack[jRow] = iColumn;
     2495            } else {
     2496              int jColumn = negSlack[jRow];
     2497              while (nextNegSlack[jColumn] >= 0)
     2498                jColumn = nextNegSlack[jColumn];
     2499              nextNegSlack[jColumn] = iColumn;
     2500            }
    24782501          }
    2479           const double * element = model2->matrix()->getElements();
    2480           const int * row = model2->matrix()->getIndices();
    2481           const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
    2482           const int * columnLength = model2->matrix()->getVectorLengths();
    2483           //bool allSlack = (numberRowsBasic==numberRows);
    2484           for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
    2485                if (!columnSolution[iColumn] || fabs(columnSolution[iColumn]) > 1.0e20) {
    2486                     double value = 0.0;
    2487                     if (columnLower[iColumn] > 0.0)
    2488                          value = columnLower[iColumn];
    2489                     else if (columnUpper[iColumn] < 0.0)
    2490                          value = columnUpper[iColumn];
    2491                     columnSolution[iColumn] = value;
    2492                }
    2493                if (columnLength[iColumn] == 1) {
    2494                     int jRow = row[columnStart[iColumn]];
    2495                     if (!columnLower[iColumn]) {
    2496                       if (element[columnStart[iColumn]] > 0.0) {
    2497                         if(posSlack[jRow] < 0) {
    2498                           posSlack[jRow] = iColumn;
    2499                         } else {
    2500                           int jColumn=posSlack[jRow];
    2501                           while (nextPosSlack[jColumn]>=0)
    2502                             jColumn=nextPosSlack[jColumn];
    2503                           nextPosSlack[jColumn]=iColumn;
    2504                         }
    2505                       } else if (element[columnStart[iColumn]] < 0.0) {
    2506                         if(negSlack[jRow] < 0) {
    2507                           negSlack[jRow] = iColumn;
    2508                         } else {
    2509                           int jColumn=negSlack[jRow];
    2510                           while (nextNegSlack[jColumn]>=0)
    2511                             jColumn=nextNegSlack[jColumn];
    2512                           nextNegSlack[jColumn]=iColumn;
    2513                         }
    2514                       }
    2515                     } else if (!columnUpper[iColumn]&& false) {// out for testing
    2516                          if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0)
    2517                               posSlack[jRow] = iColumn;
    2518                          else if (element[columnStart[iColumn]] > 0.0 && negSlack[jRow] < 0)
    2519                               negSlack[jRow] = iColumn;
    2520                     }
    2521                }
     2502        } else if (!columnUpper[iColumn] && false) { // out for testing
     2503          if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0)
     2504            posSlack[jRow] = iColumn;
     2505          else if (element[columnStart[iColumn]] > 0.0 && negSlack[jRow] < 0)
     2506            negSlack[jRow] = iColumn;
     2507        }
     2508      }
     2509    }
     2510    // now see what that does to row solution
     2511    const double *objective = model2->objective();
     2512    double *rowSolution = model2->primalRowSolution();
     2513    CoinZeroN(rowSolution, numberRows);
     2514    model2->clpMatrix()->times(1.0, columnSolution, rowSolution);
     2515    // See if we can adjust using costed slacks
     2516    double penalty = CoinMax(1.0e5, CoinMin(infeasibilityCost_ * 0.01, 1.0e10)) * optimizationDirection_;
     2517    const double *lower = model2->rowLower();
     2518    const double *upper = model2->rowUpper();
     2519    for (iRow = 0; iRow < numberRows; iRow++) {
     2520      if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
     2521        int jColumn = posSlack[iRow];
     2522        if (jColumn >= 0) {
     2523          // sort if more than one
     2524          int nPos = 1;
     2525          sort[0] = jColumn;
     2526          weight[0] = objective[jColumn];
     2527          while (nextPosSlack[jColumn] >= 0) {
     2528            jColumn = nextPosSlack[jColumn];
     2529            sort[nPos] = jColumn;
     2530            weight[nPos++] = objective[jColumn];
    25222531          }
    2523           // now see what that does to row solution
    2524           const double * objective = model2->objective();
    2525           double * rowSolution = model2->primalRowSolution();
    2526           CoinZeroN (rowSolution, numberRows);
    2527           model2->clpMatrix()->times(1.0, columnSolution, rowSolution);
    2528           // See if we can adjust using costed slacks
    2529           double penalty = CoinMax(1.0e5, CoinMin(infeasibilityCost_ * 0.01, 1.0e10)) * optimizationDirection_;
    2530           const double * lower = model2->rowLower();
    2531           const double * upper = model2->rowUpper();
    2532           for (iRow = 0; iRow < numberRows; iRow++) {
    2533                if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
    2534                     int jColumn = posSlack[iRow];
    2535                     if (jColumn >= 0) {
    2536                       // sort if more than one
    2537                       int nPos=1;
    2538                       sort[0]=jColumn;
    2539                       weight[0]=objective[jColumn];
    2540                       while (nextPosSlack[jColumn]>=0) {
    2541                         jColumn=nextPosSlack[jColumn];
    2542                         sort[nPos]=jColumn;
    2543                         weight[nPos++]=objective[jColumn];
    2544                       }
    2545                       if (nPos>1) {
    2546                         CoinSort_2(weight,weight+nPos,sort);
    2547                         for (int i=0;i<nPos;i++) {
    2548                           double difference = lower[iRow] - rowSolution[iRow];
    2549                           jColumn=sort[i];
    2550                           double elementValue = element[columnStart[jColumn]];
    2551                           assert (elementValue > 0.0);
    2552                           double value=columnSolution[jColumn];
    2553                           double movement = CoinMin(difference / elementValue, columnUpper[jColumn]-value);
    2554                           columnSolution[jColumn] += movement;
    2555                           rowSolution[iRow] += movement * elementValue;
    2556                         }
    2557                         continue;
    2558                       }
    2559                         if (jColumn<0||columnSolution[jColumn])
    2560                               continue;
    2561                          double difference = lower[iRow] - rowSolution[iRow];
    2562                          double elementValue = element[columnStart[jColumn]];
    2563                          if (elementValue > 0.0) {
    2564                               double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
    2565                               columnSolution[jColumn] = movement;
    2566                               rowSolution[iRow] += movement * elementValue;
    2567                          } else {
    2568                               double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
    2569                               columnSolution[jColumn] = movement;
    2570                               rowSolution[iRow] += movement * elementValue;
    2571                          }
    2572                     }
    2573                } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
    2574                     int jColumn = negSlack[iRow];
    2575                     if (jColumn >= 0) {
    2576                       // sort if more than one
    2577                       int nNeg=1;
    2578                       sort[0]=jColumn;
    2579                       weight[0]=objective[jColumn];
    2580                       while (nextNegSlack[jColumn]>=0) {
    2581                         jColumn=nextNegSlack[jColumn];
    2582                         sort[nNeg]=jColumn;
    2583                         weight[nNeg++]=objective[jColumn];
    2584                       }
    2585                       if (nNeg>1) {
    2586                         CoinSort_2(weight,weight+nNeg,sort);
    2587                         for (int i=0;i<nNeg;i++) {
    2588                           double difference = rowSolution[iRow]-upper[iRow];
    2589                           jColumn=sort[i];
    2590                           double elementValue = element[columnStart[jColumn]];
    2591                           assert (elementValue < 0.0);
    2592                           double value=columnSolution[jColumn];
    2593                           double movement = CoinMin(difference / -elementValue, columnUpper[jColumn]-value);
    2594                           columnSolution[jColumn] += movement;
    2595                           rowSolution[iRow] += movement * elementValue;
    2596                         }
    2597                         continue;
    2598                       }
    2599                         if (jColumn<0||columnSolution[jColumn])
    2600                               continue;
    2601                          double difference = upper[iRow] - rowSolution[iRow];
    2602                          double elementValue = element[columnStart[jColumn]];
    2603                          if (elementValue < 0.0) {
    2604                               double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
    2605                               columnSolution[jColumn] = movement;
    2606                               rowSolution[iRow] += movement * elementValue;
    2607                          } else {
    2608                               double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
    2609                               columnSolution[jColumn] = movement;
    2610                               rowSolution[iRow] += movement * elementValue;
    2611                          }
    2612                     }
    2613                }
     2532          if (nPos > 1) {
     2533            CoinSort_2(weight, weight + nPos, sort);
     2534            for (int i = 0; i < nPos; i++) {
     2535              double difference = lower[iRow] - rowSolution[iRow];
     2536              jColumn = sort[i];
     2537              double elementValue = element[columnStart[jColumn]];
     2538              assert(elementValue > 0.0);
     2539              double value = columnSolution[jColumn];
     2540              double movement = CoinMin(difference / elementValue, columnUpper[jColumn] - value);
     2541              columnSolution[jColumn] += movement;
     2542              rowSolution[iRow] += movement * elementValue;
     2543            }
     2544            continue;
    26142545          }
    2615           delete [] negSlack;
    2616           delete [] posSlack;
    2617           int nRow = numberRows;
    2618           bool network = false;
    2619           if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
    2620                network = true;
    2621                nRow *= 2;
     2546          if (jColumn < 0 || columnSolution[jColumn])
     2547            continue;
     2548          double difference = lower[iRow] - rowSolution[iRow];
     2549          double elementValue = element[columnStart[jColumn]];
     2550          if (elementValue > 0.0) {
     2551            double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
     2552            columnSolution[jColumn] = movement;
     2553            rowSolution[iRow] += movement * elementValue;
     2554          } else {
     2555            double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
     2556            columnSolution[jColumn] = movement;
     2557            rowSolution[iRow] += movement * elementValue;
    26222558          }
    2623           CoinBigIndex * addStarts = new CoinBigIndex [nRow+1];
    2624           int * addRow = new int[nRow];
    2625           double * addElement = new double[nRow];
    2626           addStarts[0] = 0;
    2627           int numberArtificials = 0;
    2628           int numberAdd = 0;
    2629           double * addCost = new double [numberRows];
    2630           for (iRow = 0; iRow < numberRows; iRow++) {
    2631                if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
    2632                     addRow[numberAdd] = iRow;
    2633                     addElement[numberAdd++] = 1.0;
    2634                     if (network) {
    2635                          addRow[numberAdd] = numberRows;
    2636                          addElement[numberAdd++] = -1.0;
    2637                     }
    2638                     addCost[numberArtificials] = penalty;
    2639                     numberArtificials++;
    2640                     addStarts[numberArtificials] = numberAdd;
    2641                } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
    2642                     addRow[numberAdd] = iRow;
    2643                     addElement[numberAdd++] = -1.0;
    2644                     if (network) {
    2645                          addRow[numberAdd] = numberRows;
    2646                          addElement[numberAdd++] = 1.0;
    2647                     }
    2648                     addCost[numberArtificials] = penalty;
    2649                     numberArtificials++;
    2650                     addStarts[numberArtificials] = numberAdd;
    2651                }
     2559        }
     2560      } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
     2561        int jColumn = negSlack[iRow];
     2562        if (jColumn >= 0) {
     2563          // sort if more than one
     2564          int nNeg = 1;
     2565          sort[0] = jColumn;
     2566          weight[0] = objective[jColumn];
     2567          while (nextNegSlack[jColumn] >= 0) {
     2568            jColumn = nextNegSlack[jColumn];
     2569            sort[nNeg] = jColumn;
     2570            weight[nNeg++] = objective[jColumn];
    26522571          }
    2653           if (numberArtificials) {
    2654                // need copy so as not to disturb original
    2655                model2 = new ClpSimplex(*model2);
    2656                if (network) {
    2657                     // network - add a null row
    2658                     model2->addRow(0, NULL, NULL, -COIN_DBL_MAX, COIN_DBL_MAX);
    2659                     numberRows++;
    2660                }
    2661                model2->addColumns(numberArtificials, NULL, NULL, addCost,
    2662                                   addStarts, addRow, addElement);
     2572          if (nNeg > 1) {
     2573            CoinSort_2(weight, weight + nNeg, sort);
     2574            for (int i = 0; i < nNeg; i++) {
     2575              double difference = rowSolution[iRow] - upper[iRow];
     2576              jColumn = sort[i];
     2577              double elementValue = element[columnStart[jColumn]];
     2578              assert(elementValue < 0.0);
     2579              double value = columnSolution[jColumn];
     2580              double movement = CoinMin(difference / -elementValue, columnUpper[jColumn] - value);
     2581              columnSolution[jColumn] += movement;
     2582              rowSolution[iRow] += movement * elementValue;
     2583            }
     2584            continue;
    26632585          }
    2664           delete [] addStarts;
    2665           delete [] addRow;
    2666           delete [] addElement;
    2667           delete [] addCost;
    2668           // look at rhs to see if to perturb
    2669           double largest = 0.0;
    2670           double smallest = 1.0e30;
    2671           for (iRow = 0; iRow < numberRows; iRow++) {
    2672                double value;
    2673                value = fabs(model2->rowLower_[iRow]);
    2674                if (value && value < 1.0e30) {
    2675                     largest = CoinMax(largest, value);
    2676                     smallest = CoinMin(smallest, value);
    2677                }
    2678                value = fabs(model2->rowUpper_[iRow]);
    2679                if (value && value < 1.0e30) {
    2680                     largest = CoinMax(largest, value);
    2681                     smallest = CoinMin(smallest, value);
    2682                }
     2586          if (jColumn < 0 || columnSolution[jColumn])
     2587            continue;
     2588          double difference = upper[iRow] - rowSolution[iRow];
     2589          double elementValue = element[columnStart[jColumn]];
     2590          if (elementValue < 0.0) {
     2591            double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
     2592            columnSolution[jColumn] = movement;
     2593            rowSolution[iRow] += movement * elementValue;
     2594          } else {
     2595            double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
     2596            columnSolution[jColumn] = movement;
     2597            rowSolution[iRow] += movement * elementValue;
    26832598          }
    2684           double * saveLower = NULL;
    2685           double * saveUpper = NULL;
    2686           if (largest < -2.01 * smallest) { // switch off for now
    2687                // perturb - so switch off standard
    2688                model2->setPerturbation(100);
    2689                saveLower = new double[numberRows];
    2690                CoinMemcpyN(model2->rowLower_, numberRows, saveLower);
    2691                saveUpper = new double[numberRows];
    2692                CoinMemcpyN(model2->rowUpper_, numberRows, saveUpper);
    2693                double * lower = model2->rowLower();
    2694                double * upper = model2->rowUpper();
    2695                for (iRow = 0; iRow < numberRows; iRow++) {
    2696                     double lowerValue = lower[iRow], upperValue = upper[iRow];
    2697                     double value = randomNumberGenerator_.randomDouble();
    2698                     if (upperValue > lowerValue + primalTolerance_) {
    2699                          if (lowerValue > -1.0e20 && lowerValue)
    2700                               lowerValue -= value * 1.0e-4 * fabs(lowerValue);
    2701                          if (upperValue < 1.0e20 && upperValue)
    2702                               upperValue += value * 1.0e-4 * fabs(upperValue);
    2703                     } else if (upperValue > 0.0) {
    2704                          upperValue -= value * 1.0e-4 * fabs(lowerValue);
    2705                          lowerValue -= value * 1.0e-4 * fabs(lowerValue);
    2706                     } else if (upperValue < 0.0) {
    2707                          upperValue += value * 1.0e-4 * fabs(lowerValue);
    2708                          lowerValue += value * 1.0e-4 * fabs(lowerValue);
    2709                     } else {
    2710                     }
    2711                     lower[iRow] = lowerValue;
    2712                     upper[iRow] = upperValue;
    2713                }
    2714           }
    2715           int i;
    2716           // Just do this number of passes in Sprint
    2717           if (doSprint > 0)
    2718                maxSprintPass = options.getExtraInfo(1);
    2719           // but if big use to get ratio
    2720           double ratio = 3;
     2599        }
     2600      }
     2601    }
     2602    delete[] negSlack;
     2603    delete[] posSlack;
     2604    int nRow = numberRows;
     2605    bool network = false;
     2606    if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) {
     2607      network = true;
     2608      nRow *= 2;
     2609    }
     2610    CoinBigIndex *addStarts = new CoinBigIndex[nRow + 1];
     2611    int *addRow = new int[nRow];
     2612    double *addElement = new double[nRow];
     2613    addStarts[0] = 0;
     2614    int numberArtificials = 0;
     2615    int numberAdd = 0;
     2616    double *addCost = new double[numberRows];
     2617    for (iRow = 0; iRow < numberRows; iRow++) {
     2618      if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
     2619        addRow[numberAdd] = iRow;
     2620        addElement[numberAdd++] = 1.0;
     2621        if (network) {
     2622          addRow[numberAdd] = numberRows;
     2623          addElement[numberAdd++] = -1.0;
     2624        }
     2625        addCost[numberArtificials] = penalty;
     2626        numberArtificials++;
     2627        addStarts[numberArtificials] = numberAdd;
     2628      } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
     2629        addRow[numberAdd] = iRow;
     2630        addElement[numberAdd++] = -1.0;
     2631        if (network) {
     2632          addRow[numberAdd] = numberRows;
     2633          addElement[numberAdd++] = 1.0;
     2634        }
     2635        addCost[numberArtificials] = penalty;
     2636        numberArtificials++;
     2637        addStarts[numberArtificials] = numberAdd;
     2638      }
     2639    }
     2640    if (numberArtificials) {
     2641      // need copy so as not to disturb original
     2642      model2 = new ClpSimplex(*model2);
     2643      if (network) {
     2644        // network - add a null row
     2645        model2->addRow(0, NULL, NULL, -COIN_DBL_MAX, COIN_DBL_MAX);
     2646        numberRows++;
     2647      }
     2648      model2->addColumns(numberArtificials, NULL, NULL, addCost,
     2649        addStarts, addRow, addElement);
     2650    }
     2651    delete[] addStarts;
     2652    delete[] addRow;
     2653    delete[] addElement;
     2654    delete[] addCost;
     2655    // look at rhs to see if to perturb
     2656    double largest = 0.0;
     2657    double smallest = 1.0e30;
     2658    for (iRow = 0; iRow < numberRows; iRow++) {
     2659      double value;
     2660      value = fabs(model2->rowLower_[iRow]);
     2661      if (value && value < 1.0e30) {
     2662        largest = CoinMax(largest, value);
     2663        smallest = CoinMin(smallest, value);
     2664      }
     2665      value = fabs(model2->rowUpper_[iRow]);
     2666      if (value && value < 1.0e30) {
     2667        largest = CoinMax(largest, value);
     2668        smallest = CoinMin(smallest, value);
     2669      }
     2670    }
     2671    double *saveLower = NULL;
     2672    double *saveUpper = NULL;
     2673    if (largest < -2.01 * smallest) { // switch off for now
     2674      // perturb - so switch off standard
     2675      model2->setPerturbation(100);
     2676      saveLower = new double[numberRows];
     2677      CoinMemcpyN(model2->rowLower_, numberRows, saveLower);
     2678      saveUpper = new double[numberRows];
     2679      CoinMemcpyN(model2->rowUpper_, numberRows, saveUpper);
     2680      double *lower = model2->rowLower();
     2681      double *upper = model2->rowUpper();
     2682      for (iRow = 0; iRow < numberRows; iRow++) {
     2683        double lowerValue = lower[iRow], upperValue = upper[iRow];
     2684        double value = randomNumberGenerator_.randomDouble();
     2685        if (upperValue > lowerValue + primalTolerance_) {
     2686          if (lowerValue > -1.0e20 && lowerValue)
     2687            lowerValue -= value * 1.0e-4 * fabs(lowerValue);
     2688          if (upperValue < 1.0e20 && upperValue)
     2689            upperValue += value * 1.0e-4 * fabs(upperValue);
     2690        } else if (upperValue > 0.0) {
     2691          upperValue -= value * 1.0e-4 * fabs(lowerValue);
     2692          lowerValue -= value * 1.0e-4 * fabs(lowerValue);
     2693        } else if (upperValue < 0.0) {
     2694          upperValue += value * 1.0e-4 * fabs(lowerValue);
     2695          lowerValue += value * 1.0e-4 * fabs(lowerValue);
     2696        } else {
     2697        }
     2698        lower[iRow] = lowerValue;
     2699        upper[iRow] = upperValue;
     2700      }
     2701    }
     2702    int i;
     2703    // Just do this number of passes in Sprint
     2704    if (doSprint > 0)
     2705      maxSprintPass = options.getExtraInfo(1);
     2706    // but if big use to get ratio
     2707    double ratio = 3;
    27212708#ifdef CLP_USEFUL_PRINTOUT
    2722           debugInt[6]=3;
    2723           debugInt[7]=maxSprintPass;
    2724 #endif
    2725           if (maxSprintPass > 1000) {
    2726                ratio = static_cast<double> (maxSprintPass) * 0.0001;
    2727                ratio = CoinMax(ratio, 1.1);
    2728                maxSprintPass = maxSprintPass % 1000;
     2709    debugInt[6] = 3;
     2710    debugInt[7] = maxSprintPass;
     2711#endif
     2712    if (maxSprintPass > 1000) {
     2713      ratio = static_cast< double >(maxSprintPass) * 0.0001;
     2714      ratio = CoinMax(ratio, 1.1);
     2715      maxSprintPass = maxSprintPass % 1000;
    27292716#ifdef COIN_DEVELOP
    2730                printf("%d passes wanted with ratio of %g\n", maxSprintPass, ratio);
    2731 #endif
    2732           }
    2733           // Just take this number of columns in small problem
    2734           int smallNumberColumns = static_cast<int> (CoinMin(ratio * numberRows, static_cast<double> (numberColumns)));
    2735           smallNumberColumns = CoinMax(smallNumberColumns, 3000);
    2736           smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
    2737           int saveSmallNumber = smallNumberColumns;
    2738           bool emergencyMode=false;
    2739           //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
    2740           //smallNumberColumns = CoinMax(smallNumberColumns,3000);
    2741           //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
    2742           // redo as may have changed
    2743           columnLower = model2->columnLower();
    2744           columnUpper = model2->columnUpper();
    2745           columnSolution = model2->primalColumnSolution();
    2746           // Set up initial list
    2747           numberSort = 0;
    2748           if (numberArtificials) {
    2749                numberSort = numberArtificials;
    2750                for (i = 0; i < numberSort; i++)
    2751                     sort[i] = i + originalNumberColumns;
    2752           }
    2753           // put in free
    2754           // maybe a solution there already
    2755           for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
    2756                if (model2->getColumnStatus(iColumn) == basic||
    2757                    (columnLower[iColumn]<-1.0e30&&columnUpper[iColumn]>1.0e30))
    2758                     sort[numberSort++] = iColumn;
    2759           }
    2760           for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
    2761                if (model2->getColumnStatus(iColumn) != basic) {
    2762                     if (columnSolution[iColumn] > columnLower[iColumn] &&
    2763                               columnSolution[iColumn] < columnUpper[iColumn] &&
    2764                               columnSolution[iColumn])
    2765                          sort[numberSort++] = iColumn;
    2766                }
    2767           }
    2768           numberSort = CoinMin(numberSort, smallNumberColumns);
     2717      printf("%d passes wanted with ratio of %g\n", maxSprintPass, ratio);
     2718#endif
     2719    }
     2720    // Just take this number of columns in small problem
     2721    int smallNumberColumns = static_cast< int >(CoinMin(ratio * numberRows, static_cast< double >(numberColumns)));
     2722    smallNumberColumns = CoinMax(smallNumberColumns, 3000);
     2723    smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
     2724    int saveSmallNumber = smallNumberColumns;
     2725    bool emergencyMode = false;
     2726    //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
     2727    //smallNumberColumns = CoinMax(smallNumberColumns,3000);
     2728    //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
     2729    // redo as may have changed
     2730    columnLower = model2->columnLower();
     2731    columnUpper = model2->columnUpper();
     2732    columnSolution = model2->primalColumnSolution();
     2733    // Set up initial list
     2734    numberSort = 0;
     2735    if (numberArtificials) {
     2736      numberSort = numberArtificials;
     2737      for (i = 0; i < numberSort; i++)
     2738        sort[i] = i + originalNumberColumns;
     2739    }
     2740    // put in free
     2741    // maybe a solution there already
     2742    for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
     2743      if (model2->getColumnStatus(iColumn) == basic || (columnLower[iColumn] < -1.0e30 && columnUpper[iColumn] > 1.0e30))
     2744        sort[numberSort++] = iColumn;
     2745    }
     2746    for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
     2747      if (model2->getColumnStatus(iColumn) != basic) {
     2748        if (columnSolution[iColumn] > columnLower[iColumn] && columnSolution[iColumn] < columnUpper[iColumn] && columnSolution[iColumn])
     2749          sort[numberSort++] = iColumn;
     2750      }
     2751    }
     2752    numberSort = CoinMin(numberSort, smallNumberColumns);
    27692753
    2770           int numberColumns = model2->numberColumns();
    2771           double * fullSolution = model2->primalColumnSolution();
     2754    int numberColumns = model2->numberColumns();
     2755    double *fullSolution = model2->primalColumnSolution();
    27722756
     2757    int iPass;
     2758    double lastObjective[] = { 1.0e31, 1.0e31 };
     2759    // It will be safe to allow dense
     2760    model2->setInitialDenseFactorization(true);
    27732761
    2774           int iPass;
    2775           double lastObjective[] = {1.0e31,1.0e31};
    2776           // It will be safe to allow dense
    2777           model2->setInitialDenseFactorization(true);
    2778 
    2779           // We will be using all rows
    2780           int * whichRows = new int [numberRows];
    2781           for (iRow = 0; iRow < numberRows; iRow++)
    2782                whichRows[iRow] = iRow;
    2783           double originalOffset;
    2784           model2->getDblParam(ClpObjOffset, originalOffset);
    2785           int totalIterations = 0;
    2786           double lastSumArtificials = COIN_DBL_MAX;
    2787           int originalMaxSprintPass = maxSprintPass;
    2788           maxSprintPass = 20; // so we do that many if infeasible
    2789           for (iPass = 0; iPass < maxSprintPass; iPass++) {
    2790                //printf("Bug until submodel new version\n");
    2791                //CoinSort_2(sort,sort+numberSort,weight);
    2792                // Create small problem
    2793                ClpSimplex small(model2, numberRows, whichRows, numberSort, sort);
    2794                small.setPerturbation(model2->perturbation());
    2795                small.setInfeasibilityCost(model2->infeasibilityCost());
    2796                if (model2->factorizationFrequency() == 200) {
    2797                     // User did not touch preset
    2798                     small.defaultFactorizationFrequency();
    2799                }
    2800                // now see what variables left out do to row solution
    2801                double * rowSolution = model2->primalRowSolution();
    2802                double * sumFixed = new double[numberRows];
    2803                CoinZeroN (sumFixed, numberRows);
    2804                int iRow, iColumn;
    2805                // zero out ones in small problem
    2806                for (iColumn = 0; iColumn < numberSort; iColumn++) {
    2807                     int kColumn = sort[iColumn];
    2808                     fullSolution[kColumn] = 0.0;
    2809                }
    2810                // Get objective offset
    2811                const double * objective = model2->objective();
    2812                double offset = 0.0;
    2813                for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
    2814                     offset += fullSolution[iColumn] * objective[iColumn];
     2762    // We will be using all rows
     2763    int *whichRows = new int[numberRows];
     2764    for (iRow = 0; iRow < numberRows; iRow++)
     2765      whichRows[iRow] = iRow;
     2766    double originalOffset;
     2767    model2->getDblParam(ClpObjOffset, originalOffset);
     2768    int totalIterations = 0;
     2769    double lastSumArtificials = COIN_DBL_MAX;
     2770    int originalMaxSprintPass = maxSprintPass;
     2771    maxSprintPass = 20; // so we do that many if infeasible
     2772    for (iPass = 0; iPass < maxSprintPass; iPass++) {
     2773      //printf("Bug until submodel new version\n");
     2774      //CoinSort_2(sort,sort+numberSort,weight);
     2775      // Create small problem
     2776      ClpSimplex small(model2, numberRows, whichRows, numberSort, sort);
     2777      small.setPerturbation(model2->perturbation());
     2778      small.setInfeasibilityCost(model2->infeasibilityCost());
     2779      if (model2->factorizationFrequency() == 200) {
     2780        // User did not touch preset
     2781        small.defaultFactorizationFrequency();
     2782      }
     2783      // now see what variables left out do to row solution
     2784      double *rowSolution = model2->primalRowSolution();
     2785      double *sumFixed = new double[numberRows];
     2786      CoinZeroN(sumFixed, numberRows);
     2787      int iRow, iColumn;
     2788      // zero out ones in small problem
     2789      for (iColumn = 0; iColumn < numberSort; iColumn++) {
     2790        int kColumn = sort[iColumn];
     2791        fullSolution[kColumn] = 0.0;
     2792      }
     2793      // Get objective offset
     2794      const double *objective = model2->objective();
     2795      double offset = 0.0;
     2796      for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
     2797        offset += fullSolution[iColumn] * objective[iColumn];
    28152798#if 0
    28162799               // Set artificials to zero if first time close to zero
     
    28222805               }
    28232806#endif
    2824                small.setDblParam(ClpObjOffset, originalOffset - offset);
    2825                model2->clpMatrix()->times(1.0, fullSolution, sumFixed);
     2807      small.setDblParam(ClpObjOffset, originalOffset - offset);
     2808      model2->clpMatrix()->times(1.0, fullSolution, sumFixed);
    28262809
    2827                double * lower = small.rowLower();
    2828                double * upper = small.rowUpper();
    2829                for (iRow = 0; iRow < numberRows; iRow++) {
    2830                     if (lower[iRow] > -1.0e50)
    2831                          lower[iRow] -= sumFixed[iRow];
    2832                     if (upper[iRow] < 1.0e50)
    2833                          upper[iRow] -= sumFixed[iRow];
    2834                     rowSolution[iRow] -= sumFixed[iRow];
    2835                }
    2836                delete [] sumFixed;
    2837                // Solve
    2838                if (interrupt)
    2839                     currentModel = &small;
    2840                small.defaultFactorizationFrequency();
    2841                if (emergencyMode) {
    2842                 // not much happening so big model
    2843                  int options=small.moreSpecialOptions();
    2844                  small.setMoreSpecialOptions(options|1048576);
    2845                 small.setMaximumIterations(1000400);
    2846                 small.setPerturbation(100);
    2847                }
    2848                if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
    2849                     // See if original wanted vector
    2850                     ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
    2851                     ClpMatrixBase * matrix = small.clpMatrix();
    2852                     if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
    2853                          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    2854                          clpMatrix->makeSpecialColumnCopy();
    2855                          small.primal(1);
    2856                          clpMatrix->releaseSpecialColumnCopy();
    2857                     } else {
     2810      double *lower = small.rowLower();
     2811      double *upper = small.rowUpper();
     2812      for (iRow = 0; iRow < numberRows; iRow++) {
     2813        if (lower[iRow] > -1.0e50)
     2814          lower[iRow] -= sumFixed[iRow];
     2815        if (upper[iRow] < 1.0e50)
     2816          upper[iRow] -= sumFixed[iRow];
     2817        rowSolution[iRow] -= sumFixed[iRow];
     2818      }
     2819      delete[] sumFixed;
     2820      // Solve
     2821      if (interrupt)
     2822        currentModel = &small;
     2823      small.defaultFactorizationFrequency();
     2824      if (emergencyMode) {
     2825        // not much happening so big model
     2826        int options = small.moreSpecialOptions();
     2827        small.setMoreSpecialOptions(options | 1048576);
     2828        small.setMaximumIterations(1000400);
     2829        small.setPerturbation(100);
     2830      }
     2831      if (dynamic_cast< ClpPackedMatrix * >(matrix_)) {
     2832        // See if original wanted vector
     2833        ClpPackedMatrix *clpMatrixO = dynamic_cast< ClpPackedMatrix * >(matrix_);
     2834        ClpMatrixBase *matrix = small.clpMatrix();
     2835        if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
     2836          ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     2837          clpMatrix->makeSpecialColumnCopy();
     2838          small.primal(1);
     2839          clpMatrix->releaseSpecialColumnCopy();
     2840        } else {
    28582841#if 1
    28592842#ifdef ABC_INHERIT
    2860                       //small.writeMps("try.mps");
    2861                       if (iPass||!numberArtificials)
    2862                          small.dealWithAbc(1,1);
    2863                        else
    2864                          small.dealWithAbc(0,0);
     2843          //small.writeMps("try.mps");
     2844          if (iPass || !numberArtificials)
     2845            small.dealWithAbc(1, 1);
     2846          else
     2847            small.dealWithAbc(0, 0);
    28652848#else
    2866                       if (iPass||!numberArtificials)
    2867                          small.primal(1);
    2868                       else
    2869                          small.dual(0);
    2870 #endif
    2871                       if (emergencyMode) {
    2872                         if (small.problemStatus()==3)
    2873                           small.setProblemStatus(0);
    2874                         smallNumberColumns=saveSmallNumber;
    2875                         emergencyMode=false;
    2876                         double * temp = new double [numberRows];
    2877                         memset(temp,0,numberRows*sizeof(double));
    2878                         double * solution = small.primalColumnSolution();
    2879                         small.matrix()->times(solution,temp);
    2880                         double sumInf=0.0;
    2881                         double * lower = small.rowLower();
    2882                         double * upper = small.rowUpper();
    2883                         for (int iRow=0;iRow<numberRows;iRow++) {
    2884                           if (temp[iRow]>upper[iRow])
    2885                             sumInf += temp[iRow]-upper[iRow];
    2886                           else if (temp[iRow]<lower[iRow])
    2887                             sumInf += lower[iRow]-temp[iRow];
    2888                         }
    2889                         printf("row inf %g\n",sumInf);
    2890                         sumInf=0.0;
    2891                         lower = small.columnLower();
    2892                         upper = small.columnUpper();
    2893                         for (int iColumn=0;iColumn<small.numberColumns();iColumn++) {
    2894                           if (solution[iColumn]>upper[iColumn])
    2895                             sumInf += solution[iColumn]-upper[iColumn];
    2896                           else if (solution[iColumn]<lower[iColumn])
    2897                             sumInf += lower[iColumn]-solution[iColumn];
    2898                         }
    2899                         printf("column inf %g\n",sumInf);
    2900                         delete [] temp;
    2901                       }
    2902                       if (small.problemStatus()) {
    2903                         int numberIterations=small.numberIterations();
    2904                         small.dual(0);
    2905                         small.setNumberIterations(small.numberIterations()+numberIterations);
    2906                       }
     2849          if (iPass || !numberArtificials)
     2850            small.primal(1);
     2851          else
     2852            small.dual(0);
     2853#endif
     2854          if (emergencyMode) {
     2855            if (small.problemStatus() == 3)
     2856              small.setProblemStatus(0);
     2857            smallNumberColumns = saveSmallNumber;
     2858            emergencyMode = false;
     2859            double *temp = new double[numberRows];
     2860            memset(temp, 0, numberRows * sizeof(double));
     2861            double *solution = small.primalColumnSolution();
     2862            small.matrix()->times(solution, temp);
     2863            double sumInf = 0.0;
     2864            double *lower = small.rowLower();
     2865            double *upper = small.rowUpper();
     2866            for (int iRow = 0; iRow < numberRows; iRow++) {
     2867              if (temp[iRow] > upper[iRow])
     2868                sumInf += temp[iRow] - upper[iRow];
     2869              else if (temp[iRow] < lower[iRow])
     2870                sumInf += lower[iRow] - temp[iRow];
     2871            }
     2872            printf("row inf %g\n", sumInf);
     2873            sumInf = 0.0;
     2874            lower = small.columnLower();
     2875            upper = small.columnUpper();
     2876            for (int iColumn = 0; iColumn < small.numberColumns(); iColumn++) {
     2877              if (solution[iColumn] > upper[iColumn])
     2878                sumInf += solution[iColumn] - upper[iColumn];
     2879              else if (solution[iColumn] < lower[iColumn])
     2880                sumInf += lower[iColumn] - solution[iColumn];
     2881            }
     2882            printf("column inf %g\n", sumInf);
     2883            delete[] temp;
     2884          }
     2885          if (small.problemStatus()) {
     2886            int numberIterations = small.numberIterations();
     2887            small.dual(0);
     2888            small.setNumberIterations(small.numberIterations() + numberIterations);
     2889          }
    29072890#else
    2908                          int numberColumns = small.numberColumns();
    2909                          int numberRows = small.numberRows();
    2910                          // Use dual region
    2911                          double * rhs = small.dualRowSolution();
    2912                          int * whichRow = new int[3*numberRows];
    2913                          int * whichColumn = new int[2*numberColumns];
    2914                          int nBound;
    2915                          ClpSimplex * small2 = ((ClpSimplexOther *) (&small))->crunch(rhs, whichRow, whichColumn,
    2916                                                nBound, false, false);
    2917                          if (small2) {
     2891          int numberColumns = small.numberColumns();
     2892          int numberRows = small.numberRows();
     2893          // Use dual region
     2894          double *rhs = small.dualRowSolution();
     2895          int *whichRow = new int[3 * numberRows];
     2896          int *whichColumn = new int[2 * numberColumns];
     2897          int nBound;
     2898          ClpSimplex *small2 = ((ClpSimplexOther *)(&small))->crunch(rhs, whichRow, whichColumn, nBound, false, false);
     2899          if (small2) {
    29182900#ifdef ABC_INHERIT
    2919                               small2->dealWithAbc(1,1);
     2901            small2->dealWithAbc(1, 1);
    29202902#else
    2921                               small.primal(1);
    2922 #endif
    2923                               if (small2->problemStatus() == 0) {
    2924                                    small.setProblemStatus(0);
    2925                                    ((ClpSimplexOther *) (&small))->afterCrunch(*small2, whichRow, whichColumn, nBound);
    2926                               } else {
     2903            small.primal(1);
     2904#endif
     2905            if (small2->problemStatus() == 0) {
     2906              small.setProblemStatus(0);
     2907              ((ClpSimplexOther *)(&small))->afterCrunch(*small2, whichRow, whichColumn, nBound);
     2908            } else {
    29272909#ifdef ABC_INHERIT
    2928                                    small2->dealWithAbc(1,1);
     2910              small2->dealWithAbc(1, 1);
    29292911#else
    2930                                    small.primal(1);
    2931 #endif
    2932                                    if (small2->problemStatus())
    2933                                         small.primal(1);
    2934                               }
    2935                               delete small2;
    2936                          } else {
    2937                               small.primal(1);
    2938                          }
    2939                          delete [] whichRow;
    2940                          delete [] whichColumn;
    2941 #endif
    2942                     }
    2943                } else {
    2944                     small.primal(1);
    2945                }
    2946                int smallIterations=small.numberIterations();
    2947                totalIterations += smallIterations;
    2948                if (2*smallIterations<CoinMin(numberRows,1000) && iPass) {
    2949                  int oldNumber=smallNumberColumns;
    2950                  if (smallIterations<100)
    2951                    smallNumberColumns *= 1.2;
    2952                  else
    2953                    smallNumberColumns *= 1.1;
    2954                  smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
    2955                  if (smallIterations<200 && smallNumberColumns>2*saveSmallNumber) {
    2956                    // try kicking it
    2957                    emergencyMode=true;
    2958                    smallNumberColumns=numberColumns;
    2959                  }
    2960                  //              smallNumberColumns = CoinMin(smallNumberColumns, 3*saveSmallNumber);
    2961                  char line[100];
    2962                  sprintf(line,"sample size increased from %d to %d",
    2963                          oldNumber,smallNumberColumns); 
    2964                  handler_->message(CLP_GENERAL, messages_)
    2965                    << line
    2966                    << CoinMessageEol;
    2967                }
    2968                // move solution back
    2969                const double * solution = small.primalColumnSolution();
    2970                for (iColumn = 0; iColumn < numberSort; iColumn++) {
    2971                     int kColumn = sort[iColumn];
    2972                     model2->setColumnStatus(kColumn, small.getColumnStatus(iColumn));
    2973                     fullSolution[kColumn] = solution[iColumn];
    2974                }
    2975                for (iRow = 0; iRow < numberRows; iRow++)
    2976                     model2->setRowStatus(iRow, small.getRowStatus(iRow));
    2977                CoinMemcpyN(small.primalRowSolution(),
    2978                            numberRows, model2->primalRowSolution());
    2979                double sumArtificials = 0.0;
    2980                for (i = 0; i < numberArtificials; i++)
    2981                     sumArtificials += fullSolution[i + originalNumberColumns];
    2982                if (sumArtificials && iPass > 5 && sumArtificials >= lastSumArtificials) {
    2983                     // increase costs
    2984                     double * cost = model2->objective() + originalNumberColumns;
    2985                     double newCost = CoinMin(1.0e10, cost[0] * 1.5);
    2986                     for (i = 0; i < numberArtificials; i++)
    2987                          cost[i] = newCost;
    2988                }
    2989                lastSumArtificials = sumArtificials;
    2990                // get reduced cost for large problem
    2991                double * djs = model2->dualColumnSolution();
    2992                CoinMemcpyN(model2->objective(), numberColumns, djs);
    2993                model2->clpMatrix()->transposeTimes(-1.0, small.dualRowSolution(), djs);
    2994                int numberNegative = 0;
    2995                double sumNegative = 0.0;
    2996                // now massage weight so all basic in plus good djs
    2997                // first count and do basic
    2998                numberSort = 0;
    2999                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    3000                     double dj = djs[iColumn] * optimizationDirection_;
    3001                     double value = fullSolution[iColumn];
    3002                     if (model2->getColumnStatus(iColumn) == ClpSimplex::basic) {
    3003                          sort[numberSort++] = iColumn;
    3004                     } else if (dj < -dualTolerance_ && value < columnUpper[iColumn]) {
    3005                          numberNegative++;
    3006                          sumNegative -= dj;
    3007                     } else if (dj > dualTolerance_ && value > columnLower[iColumn]) {
    3008                          numberNegative++;
    3009                          sumNegative += dj;
    3010                     }
    3011                }
    3012                handler_->message(CLP_SPRINT, messages_)
    3013                          << iPass + 1 << small.numberIterations() << small.objectiveValue() << sumNegative
    3014                          << numberNegative
    3015                          << CoinMessageEol;
    3016                if (sumArtificials < 1.0e-8 && originalMaxSprintPass >= 0) {
    3017                     maxSprintPass = iPass + originalMaxSprintPass;
    3018                     originalMaxSprintPass = -1;
    3019                }
    3020                if (iPass > 20)
    3021                     sumArtificials = 0.0;
    3022                if ((small.objectiveValue()*optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 15 && sumArtificials < 1.0e-8 && maxSprintPass<200) ||
    3023                          (!small.numberIterations() && iPass) ||
    3024                          iPass == maxSprintPass - 1 || small.status() == 3) {
     2912              small.primal(1);
     2913#endif
     2914              if (small2->problemStatus())
     2915                small.primal(1);
     2916            }
     2917            delete small2;
     2918          } else {
     2919            small.primal(1);
     2920          }
     2921          delete[] whichRow;
     2922          delete[] whichColumn;
     2923#endif
     2924        }
     2925      } else {
     2926        small.primal(1);
     2927      }
     2928      int smallIterations = small.numberIterations();
     2929      totalIterations += smallIterations;
     2930      if (2 * smallIterations < CoinMin(numberRows, 1000) && iPass) {
     2931        int oldNumber = smallNumberColumns;
     2932        if (smallIterations < 100)
     2933          smallNumberColumns *= 1.2;
     2934        else
     2935          smallNumberColumns *= 1.1;
     2936        smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
     2937        if (smallIterations < 200 && smallNumberColumns > 2 * saveSmallNumber) {
     2938          // try kicking it
     2939          emergencyMode = true;
     2940          smallNumberColumns = numberColumns;
     2941        }
     2942        //               smallNumberColumns = CoinMin(smallNumberColumns, 3*saveSmallNumber);
     2943        char line[100];
     2944        sprintf(line, "sample size increased from %d to %d",
     2945          oldNumber, smallNumberColumns);
     2946        handler_->message(CLP_GENERAL, messages_)
     2947          << line
     2948          << CoinMessageEol;
     2949      }
     2950      // move solution back
     2951      const double *solution = small.primalColumnSolution();
     2952      for (iColumn = 0; iColumn < numberSort; iColumn++) {
     2953        int kColumn = sort[iColumn];
     2954        model2->setColumnStatus(kColumn, small.getColumnStatus(iColumn));
     2955        fullSolution[kColumn] = solution[iColumn];
     2956      }
     2957      for (iRow = 0; iRow < numberRows; iRow++)
     2958        model2->setRowStatus(iRow, small.getRowStatus(iRow));
     2959      CoinMemcpyN(small.primalRowSolution(),
     2960        numberRows, model2->primalRowSolution());
     2961      double sumArtificials = 0.0;
     2962      for (i = 0; i < numberArtificials; i++)
     2963        sumArtificials += fullSolution[i + originalNumberColumns];
     2964      if (sumArtificials && iPass > 5 && sumArtificials >= lastSumArtificials) {
     2965        // increase costs
     2966        double *cost = model2->objective() + originalNumberColumns;
     2967        double newCost = CoinMin(1.0e10, cost[0] * 1.5);
     2968        for (i = 0; i < numberArtificials; i++)
     2969          cost[i] = newCost;
     2970      }
     2971      lastSumArtificials = sumArtificials;
     2972      // get reduced cost for large problem
     2973      double *djs = model2->dualColumnSolution();
     2974      CoinMemcpyN(model2->objective(), numberColumns, djs);
     2975      model2->clpMatrix()->transposeTimes(-1.0, small.dualRowSolution(), djs);
     2976      int numberNegative = 0;
     2977      double sumNegative = 0.0;
     2978      // now massage weight so all basic in plus good djs
     2979      // first count and do basic
     2980      numberSort = 0;
     2981      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2982        double dj = djs[iColumn] * optimizationDirection_;
     2983        double value = fullSolution[iColumn];
     2984        if (model2->getColumnStatus(iColumn) == ClpSimplex::basic) {
     2985          sort[numberSort++] = iColumn;
     2986        } else if (dj < -dualTolerance_ && value < columnUpper[iColumn]) {
     2987          numberNegative++;
     2988          sumNegative -= dj;
     2989        } else if (dj > dualTolerance_ && value > columnLower[iColumn]) {
     2990          numberNegative++;
     2991          sumNegative += dj;
     2992        }
     2993      }
     2994      handler_->message(CLP_SPRINT, messages_)
     2995        << iPass + 1 << small.numberIterations() << small.objectiveValue() << sumNegative
     2996        << numberNegative
     2997        << CoinMessageEol;
     2998      if (sumArtificials < 1.0e-8 && originalMaxSprintPass >= 0) {
     2999        maxSprintPass = iPass + originalMaxSprintPass;
     3000        originalMaxSprintPass = -1;
     3001      }
     3002      if (iPass > 20)
     3003        sumArtificials = 0.0;
     3004      if ((small.objectiveValue() * optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 15 && sumArtificials < 1.0e-8 && maxSprintPass < 200) || (!small.numberIterations() && iPass) || iPass == maxSprintPass - 1 || small.status() == 3) {
    30253005
    3026                     break; // finished
    3027                } else {
    3028                     lastObjective[1] = lastObjective[0];
    3029                     lastObjective[0] = small.objectiveValue() * optimizationDirection_;
    3030                     double tolerance;
    3031                     double averageNegDj = sumNegative / static_cast<double> (numberNegative + 1);
    3032                     if (numberNegative + numberSort > smallNumberColumns && false)
    3033                          tolerance = -dualTolerance_;
    3034                     else
    3035                          tolerance = 10.0 * averageNegDj;
    3036                     if (emergencyMode)
    3037                       tolerance=1.0e100;
    3038                     int saveN = numberSort;
    3039                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    3040                          double dj = djs[iColumn] * optimizationDirection_;
    3041                          double value = fullSolution[iColumn];
    3042                          if (model2->getColumnStatus(iColumn) != ClpSimplex::basic) {
    3043                               if (dj < -dualTolerance_ && value < columnUpper[iColumn])
    3044                                    dj = dj;
    3045                               else if (dj > dualTolerance_ && value > columnLower[iColumn])
    3046                                    dj = -dj;
    3047                               else if (columnUpper[iColumn] > columnLower[iColumn])
    3048                                    dj = fabs(dj);
    3049                               else
    3050                                    dj = 1.0e50;
    3051                               if (dj < tolerance) {
    3052                                    weight[numberSort] = dj;
    3053                                    sort[numberSort++] = iColumn;
    3054                               }
    3055                          }
    3056                     }
    3057                     // sort
    3058                     CoinSort_2(weight + saveN, weight + numberSort, sort + saveN);
    3059                     if (numberSort<smallNumberColumns)
    3060                       printf("using %d columns not %d\n",numberSort,smallNumberColumns);
    3061                     numberSort = CoinMin(smallNumberColumns, numberSort);
    3062                     // try singletons
    3063                     char * markX = new char[numberColumns];
    3064                     memset(markX,0,numberColumns);
    3065                     for (int i=0;i<numberSort;i++)
    3066                       markX[sort[i]]=1;
    3067                     int n=numberSort;
    3068                     for (int i=0;i<numberColumns;i++) {
    3069                       if (columnLength[i]==1&&!markX[i])
    3070                         sort[numberSort++]=i;
    3071                     }
    3072                     if (n<numberSort)
    3073                       printf("%d slacks added\n",numberSort-n);
    3074                     delete [] markX;
    3075                }
     3006        break; // finished
     3007      } else {
     3008        lastObjective[1] = lastObjective[0];
     3009        lastObjective[0] = small.objectiveValue() * optimizationDirection_;
     3010        double tolerance;
     3011        double averageNegDj = sumNegative / static_cast< double >(numberNegative + 1);
     3012        if (numberNegative + numberSort > smallNumberColumns && false)
     3013          tolerance = -dualTolerance_;
     3014        else
     3015          tolerance = 10.0 * averageNegDj;
     3016        if (emergencyMode)
     3017          tolerance = 1.0e100;
     3018        int saveN = numberSort;
     3019        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     3020          double dj = djs[iColumn] * optimizationDirection_;
     3021          double value = fullSolution[iColumn];
     3022          if (model2->getColumnStatus(iColumn) != ClpSimplex::basic) {
     3023            if (dj < -dualTolerance_ && value < columnUpper[iColumn])
     3024              dj = dj;
     3025            else if (dj > dualTolerance_ && value > columnLower[iColumn])
     3026              dj = -dj;
     3027            else if (columnUpper[iColumn] > columnLower[iColumn])
     3028              dj = fabs(dj);
     3029            else
     3030              dj = 1.0e50;
     3031            if (dj < tolerance) {
     3032              weight[numberSort] = dj;
     3033              sort[numberSort++] = iColumn;
     3034            }
    30763035          }
    3077           if (interrupt)
    3078                currentModel = model2;
    3079           for (i = 0; i < numberArtificials; i++)
    3080                sort[i] = i + originalNumberColumns;
    3081           model2->deleteColumns(numberArtificials, sort);
    3082           if (network) {
    3083                int iRow = numberRows - 1;
    3084                model2->deleteRows(1, &iRow);
    3085           }
    3086           delete [] weight;
    3087           delete [] sort;
    3088           delete [] whichRows;
    3089           if (saveLower) {
    3090                // unperturb and clean
    3091                for (iRow = 0; iRow < numberRows; iRow++) {
    3092                     model2->rowLower_[iRow] = saveLower[iRow];
    3093                     model2->rowUpper_[iRow] = saveUpper[iRow];
    3094                }
    3095                delete [] saveLower;
    3096                delete [] saveUpper;
    3097           }
     3036        }
     3037        // sort
     3038        CoinSort_2(weight + saveN, weight + numberSort, sort + saveN);
     3039        if (numberSort < smallNumberColumns)
     3040          printf("using %d columns not %d\n", numberSort, smallNumberColumns);
     3041        numberSort = CoinMin(smallNumberColumns, numberSort);
     3042        // try singletons
     3043        char *markX = new char[numberColumns];
     3044        memset(markX, 0, numberColumns);
     3045        for (int i = 0; i < numberSort; i++)
     3046          markX[sort[i]] = 1;
     3047        int n = numberSort;
     3048        for (int i = 0; i < numberColumns; i++) {
     3049          if (columnLength[i] == 1 && !markX[i])
     3050            sort[numberSort++] = i;
     3051        }
     3052        if (n < numberSort)
     3053          printf("%d slacks added\n", numberSort - n);
     3054        delete[] markX;
     3055      }
     3056    }
     3057    if (interrupt)
     3058      currentModel = model2;
     3059    for (i = 0; i < numberArtificials; i++)
     3060      sort[i] = i + originalNumberColumns;
     3061    model2->deleteColumns(numberArtificials, sort);
     3062    if (network) {
     3063      int iRow = numberRows - 1;
     3064      model2->deleteRows(1, &iRow);
     3065    }
     3066    delete[] weight;
     3067    delete[] sort;
     3068    delete[] whichRows;
     3069    if (saveLower) {
     3070      // unperturb and clean
     3071      for (iRow = 0; iRow < numberRows; iRow++) {
     3072        model2->rowLower_[iRow] = saveLower[iRow];
     3073&n