source: trunk/Clp/src/ClpSolve.cpp @ 2439

Last change on this file since 2439 was 2439, checked in by stefan, 3 months ago

merge r2435,r2437,r2438 from stable/1.17 to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 280.1 KB
Line 
1/* $Id: ClpSolve.cpp 2439 2019-03-28 02:10:04Z stefan $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6// This file has higher level solve functions
7
8#include "CoinPragma.hpp"
9#include "ClpConfig.h"
10
11// check already here if COIN_HAS_GLPK is defined, since we do not want to get confused by a COIN_HAS_GLPK in config_coinutils.h
12#if defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD) || defined(COIN_HAS_GLPK)
13#define UFL_BARRIER
14#endif
15
16#include <math.h>
17#ifdef _MSC_VER
18#include <windows.h> // for Sleep()
19#ifdef small
20#undef small
21#endif
22#else
23#include <unistd.h> // for usleep()
24#endif
25
26#include "CoinHelperFunctions.hpp"
27#include "ClpHelperFunctions.hpp"
28#include "CoinSort.hpp"
29#include "ClpFactorization.hpp"
30#include "ClpSimplex.hpp"
31#include "ClpSimplexOther.hpp"
32#include "ClpSimplexDual.hpp"
33#ifndef SLIM_CLP
34#include "ClpQuadraticObjective.hpp"
35#include "ClpInterior.hpp"
36#include "ClpCholeskyDense.hpp"
37#include "ClpCholeskyBase.hpp"
38#include "ClpPlusMinusOneMatrix.hpp"
39#include "ClpNetworkMatrix.hpp"
40#endif
41#include "ClpEventHandler.hpp"
42#include "ClpLinearObjective.hpp"
43#include "ClpSolve.hpp"
44#include "ClpPackedMatrix.hpp"
45#include "ClpMessage.hpp"
46#include "CoinTime.hpp"
47#if CLP_HAS_ABC
48#include "CoinAbcCommon.hpp"
49#endif
50#ifdef ABC_INHERIT
51#include "AbcSimplex.hpp"
52#include "AbcSimplexFactorization.hpp"
53#endif
54#include "CoinStructuredModel.hpp"
55double zz_slack_value = 0.0;
56#ifdef CLP_USEFUL_PRINTOUT
57double debugDouble[10];
58int debugInt[24];
59#endif
60
61#include "ClpPresolve.hpp"
62#ifndef SLIM_CLP
63#include "Idiot.hpp"
64#ifdef COIN_HAS_WSMP
65#include "ClpCholeskyWssmp.hpp"
66#include "ClpCholeskyWssmpKKT.hpp"
67#endif
68#include "ClpCholeskyUfl.hpp"
69#ifdef TAUCS_BARRIER
70#include "ClpCholeskyTaucs.hpp"
71#endif
72#ifdef PARDISO_BARRIER
73#include "ClpCholeskyPardiso.hpp"
74#endif
75#include "ClpCholeskyMumps.hpp"
76#ifdef COIN_HAS_VOL
77#include "VolVolume.hpp"
78#include "CoinHelperFunctions.hpp"
79#include "CoinPackedMatrix.hpp"
80#include "CoinMpsIO.hpp"
81//#############################################################################
82
83class lpHook : public VOL_user_hooks {
84private:
85  lpHook(const lpHook &);
86  lpHook &operator=(const lpHook &);
87
88private:
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_;
99
100  /// The problem matrix in a row ordered form
101  CoinPackedMatrix rowMatrix_;
102  /// The problem matrix in a column ordered form
103  CoinPackedMatrix colMatrix_;
104
105public:
106  lpHook(double *clb, double *cub, double *obj,
107    double *rhs, char *sense, const CoinPackedMatrix &mat);
108  virtual ~lpHook();
109
110public:
111  // for all hooks: return value of -1 means that volume should quit
112  /** compute reduced costs
113         @param u (IN) the dual variables
114         @param rc (OUT) the reduced cost with respect to the dual values
115     */
116  virtual int compute_rc(const VOL_dvector &u, VOL_dvector &rc);
117
118  /** Solve the subproblem for the subgradient step.
119         @param dual (IN) the dual variables
120         @param rc (IN) the reduced cost with respect to the dual values
121         @param lcost (OUT) the lagrangean cost with respect to the dual values
122         @param x (OUT) the primal result of solving the subproblem
123         @param v (OUT) b-Ax for the relaxed constraints
124         @param pcost (OUT) the primal objective value of <code>x</code>
125     */
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
130         an integer solution
131         @param x (IN) the primal vector
132         @param heur_val (OUT) the value of the integer solution (return
133         <code>DBL_MAX</code> here if no feas sol was found
134     */
135  virtual int heuristics(const VOL_problem &p,
136    const VOL_dvector &x, double &heur_val)
137  {
138    return 0;
139  }
140};
141
142//#############################################################################
143
144lpHook::lpHook(double *clb, double *cub, double *obj,
145  double *rhs, char *sense,
146  const CoinPackedMatrix &mat)
147{
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);
156}
157
158//-----------------------------------------------------------------------------
159
160lpHook::~lpHook()
161{
162}
163
164//#############################################################################
165
166int lpHook::compute_rc(const VOL_dvector &u, VOL_dvector &rc)
167{
168  rowMatrix_.transposeTimes(u.v, rc.v);
169  const int psize = rowMatrix_.getNumCols();
170
171  for (int i = 0; i < psize; ++i)
172    rc[i] = objcoeffs_[i] - rc[i];
173  return 0;
174}
175
176//-----------------------------------------------------------------------------
177
178int lpHook::solve_subproblem(const VOL_dvector &dual, const VOL_dvector &rc,
179  double &lcost, VOL_dvector &x, VOL_dvector &v,
180  double &pcost)
181{
182  int i;
183  const int psize = x.size();
184  const int dsize = v.size();
185
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];
189
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];
196
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];
201
202  // compute the lagrangean primal objective
203  pcost = 0;
204  for (i = 0; i < psize; ++i)
205    pcost += x[i] * objcoeffs_[i];
206
207  return 0;
208}
209
210//#############################################################################
211/** A quick inlined function to convert from lb/ub style constraint
212    definition to sense/rhs/range style */
213inline void
214convertBoundToSense(const double lower, const double upper,
215  char &sense, double &right,
216  double &range)
217{
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  }
241}
242
243static int
244solveWithVolume(ClpSimplex *model, int numberPasses, int doIdiot)
245{
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];
265
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);
319
320  volprob.solve(myHook, false /* no warmstart */);
321
322  if (saveLower != model->columnLower()) {
323    delete[] saveLower;
324    delete[] saveUpper;
325  }
326  //------------- extract the solution ---------------------------
327
328  //printf("Best lagrangean value: %f\n", volprob.value);
329
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  }
346
347  //printf("Average primal constraint violation: %f\n", avg/dsize);
348
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;
355}
356#endif
357static ClpInterior *currentModel2 = NULL;
358#endif
359//#############################################################################
360// Allow for interrupts
361// But is this threadsafe ? (so switched off by option)
362
363#include "CoinSignal.hpp"
364static ClpSimplex *currentModel = NULL;
365#ifdef ABC_INHERIT
366static AbcSimplex *currentAbcModel = NULL;
367#endif
368
369extern "C" {
370static void
371#if defined(_MSC_VER)
372  __cdecl
373#endif // _MSC_VER
374  signal_handler(int /*whichSignal*/)
375{
376  if (currentModel != NULL)
377    currentModel->setMaximumIterations(0); // stop at next iterations
378#ifdef ABC_INHERIT
379  if (currentAbcModel != NULL)
380    currentAbcModel->setMaximumIterations(0); // stop at next iterations
381#endif
382#ifndef SLIM_CLP
383  if (currentModel2 != NULL)
384    currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
385#endif
386  return;
387}
388}
389#if ABC_INSTRUMENT > 1
390int abcPricing[20];
391int abcPricingDense[20];
392static int trueNumberRows;
393static int numberTypes;
394#define MAX_TYPES 25
395#define MAX_COUNT 20
396#define MAX_FRACTION 101
397static char *types[MAX_TYPES];
398static double counts[MAX_TYPES][MAX_COUNT];
399static double countsFraction[MAX_TYPES][MAX_FRACTION];
400static double *currentCounts;
401static double *currentCountsFraction;
402static int currentType;
403static double workMultiplier[MAX_TYPES];
404static double work[MAX_TYPES];
405static double currentWork;
406static double otherWork[MAX_TYPES];
407static int timesCalled[MAX_TYPES];
408static int timesStarted[MAX_TYPES];
409static int fractionDivider;
410void instrument_initialize(int numberRows)
411{
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);
425}
426void instrument_start(const char *type, int numberRowsEtc)
427{
428  if (currentType >= 0)
429    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;
435      break;
436    }
437  }
438  if (currentType == -1) {
439    assert(numberTypes < MAX_TYPES);
440    currentType = numberTypes;
441    types[numberTypes++] = strdup(type);
442  }
443  currentCounts = &counts[currentType][0];
444  currentCountsFraction = &countsFraction[currentType][0];
445  timesStarted[currentType]++;
446  assert(trueNumberRows);
447  workMultiplier[currentType] += static_cast< double >(numberRowsEtc) / static_cast< double >(trueNumberRows);
448}
449void instrument_add(int count)
450{
451  assert(currentType >= 0);
452  currentWork += count;
453  timesCalled[currentType]++;
454  if (count < MAX_COUNT - 1)
455    currentCounts[count]++;
456  else
457    currentCounts[MAX_COUNT - 1]++;
458  assert(count / fractionDivider >= 0 && count / fractionDivider < MAX_FRACTION);
459  currentCountsFraction[count / fractionDivider]++;
460}
461void instrument_do(const char *type, double count)
462{
463  int iType = -1;
464  for (int i = 0; i < numberTypes; i++) {
465    if (!strcmp(types[i], type)) {
466      iType = i;
467      break;
468    }
469  }
470  if (iType == -1) {
471    assert(numberTypes < MAX_TYPES);
472    iType = numberTypes;
473    types[numberTypes++] = strdup(type);
474  }
475  timesStarted[iType]++;
476  otherWork[iType] += count;
477}
478void instrument_end()
479{
480  work[currentType] += currentWork;
481  currentType = -1;
482}
483void instrument_end_and_adjust(double factor)
484{
485  work[currentType] += currentWork * factor;
486  currentType = -1;
487}
488void instrument_print()
489{
490  for (int iType = 0; iType < numberTypes; iType++) {
491    currentCounts = &counts[iType][0];
492    currentCountsFraction = &countsFraction[iType][0];
493    if (!otherWork[iType]) {
494      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]);
515      }
516      printf("\n");
517      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]);
540      }
541      printf("\n");
542    } else {
543      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]);
546    }
547    free(types[iType]);
548  }
549}
550#endif
551#if ABC_PARALLEL == 2
552#ifndef FAKE_CILK
553int number_cilk_workers = 0;
554#include <cilk/cilk_api.h>
555#endif
556#endif
557#ifdef ABC_INHERIT
558AbcSimplex *
559ClpSimplex::dealWithAbc(int solveType, int startUp,
560  bool interrupt)
561{
562  bool keepAbc = startUp >= 10;
563  startUp = startUp % 10;
564  AbcSimplex *abcModel2 = NULL;
565  if (!this->abcState() || !numberRows_ || !numberColumns_) {
566    //this->readBasis("aaa.bas");
567    if (!solveType) {
568      this->dual(0);
569    } else if (solveType == 1) {
570      int ifValuesPass = startUp ? 1 : 0;
571      if (startUp == 3)
572        ifValuesPass = 2;
573      this->primal(ifValuesPass);
574    }
575    //this->writeBasis("a.bas",true);
576  } else {
577    abcModel2 = new AbcSimplex(*this);
578    if (interrupt)
579      currentAbcModel = abcModel2;
580    //if (abcSimplex_) {
581    // move factorization stuff
582    abcModel2->factorization()->synchronize(this->factorization(), abcModel2);
583    //}
584    //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;
590    // 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;
599#if 0
600      ClpSimplex * thisModel = static_cast<ClpSimplexOther *> (this)->dualOfModel(1.0,1.0);
601      if (thisModel) {
602        printf("Dual of model has %d rows and %d columns\n",
603               thisModel->numberRows(), thisModel->numberColumns());
604        thisModel->setOptimizationDirection(1.0);
605        Idiot info(*thisModel);
606        info.setStrategy(512 | info.getStrategy());
607        // Allow for scaling
608        info.setStrategy(32 | info.getStrategy());
609        info.setStartingWeight(1.0e3);
610        info.setReduceIterations(6);
611        info.crash(50, this->messageHandler(), this->messagesPointer(),false);
612        // make sure later that primal solution in correct place
613        // and has correct sign
614        abcModel2->setupDualValuesPass(thisModel->primalColumnSolution(),
615                                       thisModel->dualRowSolution(),type);
616        //thisModel->dual();
617        delete thisModel;
618      }
619#else
620      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
643#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) {
649      // maximum
650      numberCpu = 4;
651    } else if (numberCpu == 10) {
652      // decide
653      if (abcModel2->getNumElements() < 5000)
654        numberCpu = 1;
655#if ABC_PARALLEL == 2
656#ifndef FAKE_CILK
657      else if (number_cilk_workers > 1)
658        numberCpu = CoinMin(2 * number_cilk_workers, 8);
659#endif
660#endif
661      else
662        numberCpu = 1;
663    } else {
664#if ABC_PARALLEL == 2
665#ifndef FAKE_CILK
666      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;
671
672#endif
673#endif
674    }
675    char line[200];
676#if ABC_PARALLEL
677#if ABC_PARALLEL == 2
678#ifndef FAKE_CILK
679    if (!number_cilk_workers) {
680      number_cilk_workers = __cilkrts_get_nworkers();
681      sprintf(line, "%d cilk workers", number_cilk_workers);
682      handler_->message(CLP_GENERAL, messages_)
683        << line
684        << CoinMessageEol;
685    }
686#endif
687#endif
688    abcModel2->setParallelMode(numberCpu - 1);
689#endif
690    //if (abcState()==3||abcState()==4) {
691    //abcModel2->setMoreSpecialOptions((131072*2)|abcModel2->moreSpecialOptions());
692    //}
693    //if (processTune>0&&processTune<8)
694    //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|131072*processTune);
695#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));
701    instrument_initialize(abcModel2->numberRows());
702#endif
703#endif
704    if (!solveType) {
705      abcModel2->ClpSimplex::doAbcDual();
706    } else if (solveType == 1) {
707      int saveOptions = abcModel2->specialOptions();
708      //if (startUp==2)
709      //abcModel2->setSpecialOptions(8192|saveOptions);
710      abcModel2->ClpSimplex::doAbcPrimal(ifValuesPass);
711      abcModel2->setSpecialOptions(saveOptions);
712    }
713#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());
721      handler_->message(CLP_GENERAL, messages_)
722        << line
723        << CoinMessageEol;
724#if ABC_INSTRUMENT > 1
725      {
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");
775      }
776      instrument_print();
777#endif
778    }
779#endif
780    abcModel2->moveStatusToClp(this);
781#if 1
782    if (!problemStatus_) {
783      double offset;
784      CoinMemcpyN(this->objectiveAsObject()->gradient(this,
785                    this->primalColumnSolution(), offset, true),
786        numberColumns_, this->dualColumnSolution());
787      this->clpMatrix()->transposeTimes(-1.0,
788        this->dualRowSolution(),
789        this->dualColumnSolution());
790      memset(this->primalRowSolution(), 0, numberRows_ * sizeof(double));
791      this->clpMatrix()->times(1.0,
792        this->primalColumnSolution(),
793        this->primalRowSolution());
794      this->checkSolutionInternal();
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_);
799#endif
800      } else {
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_);
808#endif
809      } else {
810        sumPrimalInfeasibilities_ = 0.0;
811        numberPrimalInfeasibilities_ = 0;
812      }
813      problemStatus_ = 0;
814    }
815#endif
816    //ClpModel::stopPermanentArrays();
817    this->setSpecialOptions(this->specialOptions() & ~65536);
818#if 0
819    this->writeBasis("a.bas",true);
820    for (int i=0;i<this->numberRows();i++)
821      printf("%d %g\n",i,this->dualRowSolution()[i]);
822    this->dual();
823    this->writeBasis("b.bas",true);
824    for (int i=0;i<this->numberRows();i++)
825      printf("%d %g\n",i,this->dualRowSolution()[i]);
826#endif
827    // switch off initialSolve flag
828    moreSpecialOptions_ &= ~16384;
829    //this->setNumberIterations(abcModel2->numberIterations()+this->numberIterations());
830    if (!keepAbc) {
831      delete abcModel2;
832      abcModel2 = NULL;
833    }
834  }
835  return abcModel2;
836}
837#endif
838/** General solve algorithm which can do presolve
839    special options (bits)
840    1 - do not perturb
841    2 - do not scale
842    4 - use crash (default allslack in dual, idiot in primal)
843    8 - all slack basis in primal
844    16 - switch off interrupt handling
845    32 - do not try and make plus minus one matrix
846    64 - do not use sprint even if problem looks good
847 */
848int ClpSimplex::initialSolve(ClpSolve &options)
849{
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;
864#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_;
928#ifndef SLIM_CLP
929#ifndef NO_RTTI
930  if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) {
931    // network - switch off stuff
932    presolve = ClpSolve::presolveOff;
933  }
934#else
935  if (matrix_->type() == 11) {
936    // network - switch off stuff
937    presolve = ClpSolve::presolveOff;
938  }
939#endif
940#endif
941#ifndef CLPSOLVE_ACTIONS
942#define CLPSOLVE_ACTIONS 2
943#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;
951#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    }
966#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);
977#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 {
998#if 0 //def ABC_INHERIT
999            {
1000              AbcSimplex * abcModel2=new AbcSimplex(*model2);
1001              delete model2;
1002              model2=abcModel2;
1003              pinfo->setPresolvedModel(model2);
1004            }
1005#else
1006      //ClpModel::stopPermanentArrays();
1007      //setSpecialOptions(specialOptions()&~65536);
1008      // try setting tolerances up
1009#if CLPSOLVE_ACTIONS
1010      bool changeTolerances = wasAutomatic;
1011#if CLPSOLVE_ACTIONS > 1
1012      changeTolerances = true;
1013#endif
1014      if (changeTolerances && model2 != this) {
1015#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      }
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            }
1317          }
1318          value = fabs(rowUpper[i]);
1319          if (value && value != COIN_DBL_MAX) {
1320            if (value != last) {
1321              nDifferent++;
1322              last = value;
1323            }
1324          }
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()));
1436#else
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  }
1451#ifdef COIN_HAS_VOL
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();
1460#ifndef SLIM_CLP
1461#ifndef NO_RTTI
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)
1591          doIdiot = 0;
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        }
1839#else
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
1966#ifndef SLIM_CLP
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;
2043          } else {
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);
2052          }
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            }
2067          } else {
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            }
2076          }
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);
2081          }
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
2117#ifdef CLP_USEFUL_PRINTOUT
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);
2143          }
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);
2149          }
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);
2156          }
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
2168#ifdef COIN_DEVELOP
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;
2178#if 0
2179                    double * solution = model2->primalColumnSolution();
2180                    int iColumn;
2181                    double * saveLower = new double[numberColumns];
2182                    CoinMemcpyN(model2->columnLower(), numberColumns, saveLower);
2183                    double * saveUpper = new double[numberColumns];
2184                    CoinMemcpyN(model2->columnUpper(), numberColumns, saveUpper);
2185                    printf("doing tighten before idiot\n");
2186                    model2->tightenPrimalBounds();
2187                    // Move solution
2188                    double * columnLower = model2->columnLower();
2189                    double * columnUpper = model2->columnUpper();
2190                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2191                         if (columnLower[iColumn] > 0.0)
2192                              solution[iColumn] = columnLower[iColumn];
2193                         else if (columnUpper[iColumn] < 0.0)
2194                              solution[iColumn] = columnUpper[iColumn];
2195                         else
2196                              solution[iColumn] = 0.0;
2197                    }
2198                    CoinMemcpyN(saveLower, numberColumns, columnLower);
2199                    CoinMemcpyN(saveUpper, numberColumns, columnUpper);
2200                    delete [] saveLower;
2201                    delete [] saveUpper;
2202#else
2203        // Allow for crossover
2204        //#define LACI_TRY
2205#ifndef LACI_TRY
2206        //if (doIdiot>0)
2207#if 0 //def ABC_INHERIT
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);
2218#define CLP_DUAL_IDIOT
2219#ifdef CLP_DUAL_IDIOT
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);
2231#if 0
2232                        info.setStrategy(512 | info.getStrategy());
2233                        // Allow for scaling
2234                        info.setStrategy(32 | info.getStrategy());
2235                        info.setStartingWeight(1.0e3);
2236                        info.setReduceIterations(6);
2237#endif
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);
2244#if 0
2245                        static_cast<ClpSimplexOther *> (&temp)->restoreFromDual(dualModel2);
2246#else
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));
2290#else
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);
2331#ifdef ABC_INHERIT
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);
2341#if ABC_PARALLEL
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    }
2385#ifndef SLIM_CLP
2386    if (doSlp && objective_->type() == 2) {
2387      model2->nonlinearSLP(doSlp, 1.0e-5);
2388    }
2389#endif
2390#ifndef LACI_TRY
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 {
2403#ifndef ABC_INHERIT
2404          model2->primal(primalStartup);
2405#else
2406          model2->dealWithAbc(1, primalStartup, interrupt);
2407#endif
2408        }
2409      } else {
2410#ifndef ABC_INHERIT
2411        model2->primal(primalStartup);
2412#else
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    /*
2427            This driver implements what I called Sprint when I introduced the idea
2428            many years ago.  Cplex calls it "sifting" which I think is just as silly.
2429            When I thought of this trivial idea
2430            it reminded me of an LP code of the 60's called sprint which after
2431            every factorization took a subset of the matrix into memory (all
2432            64K words!) and then iterated very fast on that subset.  On the
2433            problems of those days it did not work very well, but it worked very
2434            well on aircrew scheduling problems where there were very large numbers
2435            of columns all with the same flavor.
2436          */
2437
2438    /* The idea works best if you can get feasible easily.  To make it
2439             more general we can add in costed slacks */
2440
2441    int originalNumberColumns = model2->numberColumns();
2442    int numberRows = model2->numberRows();
2443    ClpSimplex *originalModel2 = model2;
2444
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();
2455
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            }
2501          }
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];
2531          }
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;
2545          }
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;
2558          }
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];
2571          }
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;
2585          }
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;
2598          }
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;
2708#ifdef CLP_USEFUL_PRINTOUT
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;
2716#ifdef COIN_DEVELOP
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);
2753
2754    int numberColumns = model2->numberColumns();
2755    double *fullSolution = model2->primalColumnSolution();
2756
2757    int iPass;
2758    double lastObjective[] = { 1.0e31, 1.0e31 };
2759    // It will be safe to allow dense
2760    model2->setInitialDenseFactorization(true);
2761
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];
2798#if 0
2799               // Set artificials to zero if first time close to zero
2800               for (iColumn = originalNumberColumns; iColumn < numberColumns; iColumn++) {
2801                 if (fullSolution[iColumn]<primalTolerance_&&objective[iColumn]==penalty) {
2802                   model2->objective()[iColumn]=2.0*penalty;
2803                   fullSolution[iColumn]=0.0;
2804                 }
2805               }
2806#endif
2807      small.setDblParam(ClpObjOffset, originalOffset - offset);
2808      model2->clpMatrix()->times(1.0, fullSolution, sumFixed);
2809
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 {
2841#if 1
2842#ifdef ABC_INHERIT
2843          //small.writeMps("try.mps");
2844          if (iPass || !numberArtificials)
2845            small.dealWithAbc(1, 1);
2846          else
2847            small.dealWithAbc(0, 0);
2848#else
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          }
2890#else
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) {
2900#ifdef ABC_INHERIT
2901            small2->dealWithAbc(1, 1);
2902#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 {
2909#ifdef ABC_INHERIT
2910              small2->dealWithAbc(1, 1);
2911#else
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) {
3005
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            }
3035          }
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        model2->rowUpper_[iRow] = saveUpper[iRow];
3074      }
3075      delete[] saveLower;
3076      delete[] saveUpper;
3077    }
3078#ifdef ABC_INHERIT
3079    model2->dealWithAbc(1, 1);
3080#else
3081    model2->primal(1);
3082#endif
3083    model2->setPerturbation(savePerturbation);
3084    if (model2 != originalModel2) {
3085      originalModel2->moveInfo(*model2);
3086      delete model2;
3087      model2 = originalModel2;
3088    }
3089    time2 = CoinCpuTime();
3090    timeCore = time2 - timeX;
3091    handler_->message(CLP_INTERVAL_TIMING, messages_)
3092      << "Sprint" << timeCore << time2 - time1
3093      << CoinMessageEol;
3094    timeX = time2;
3095    model2->setNumberIterations(model2->numberIterations() + totalIterations);
3096  } else if (method == ClpSolve::useBarrier || method == ClpSolve::useBarrierNoCross) {
3097    if (presolve == ClpSolve::presolveOn) {
3098      int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0);
3099      if (numberInfeasibilities) {
3100        handler_->message(CLP_INFEASIBLE, messages_)
3101          << CoinMessageEol;
3102        delete model2;
3103        model2 = this;
3104        presolve = ClpSolve::presolveOff;
3105      }
3106    }
3107#ifndef SLIM_CLP
3108    //printf("***** experimental pretty crude barrier\n");
3109    //#define SAVEIT 2
3110#ifndef SAVEIT
3111#define BORROW
3112#endif
3113#ifdef BORROW
3114    ClpInterior barrier;
3115    barrier.borrowModel(*model2);
3116#else
3117    ClpInterior barrier(*model2);
3118#endif
3119    barrier.eventHandler()->setSimplex(NULL);
3120    if (interrupt)
3121      currentModel2 = &barrier;
3122    if (barrier.numberRows() + barrier.numberColumns() > 10000)
3123      barrier.setMaximumBarrierIterations(1000);
3124    int barrierOptions = options.getSpecialOption(4);
3125    int aggressiveGamma = 0;
3126    bool presolveInCrossover = false;
3127    bool scale = false;
3128    bool doKKT = false;
3129    bool forceFixing = false;
3130    int speed = 0;
3131    if (barrierOptions & 16) {
3132      barrierOptions &= ~16;
3133      doKKT = true;
3134    }
3135    if (barrierOptions & (32 + 64 + 128)) {
3136      aggressiveGamma = (barrierOptions & (32 + 64 + 128)) >> 5;
3137      barrierOptions &= ~(32 + 64 + 128);
3138    }
3139    if (barrierOptions & 256) {
3140      barrierOptions &= ~256;
3141      presolveInCrossover = true;
3142    }
3143    if (barrierOptions & 512) {
3144      barrierOptions &= ~512;
3145      forceFixing = true;
3146    }
3147    if (barrierOptions & 1024) {
3148      barrierOptions &= ~1024;
3149      barrier.setProjectionTolerance(1.0e-9);
3150    }
3151    if (barrierOptions & (2048 | 4096)) {
3152      speed = (barrierOptions & (2048 | 4096)) >> 11;
3153      barrierOptions &= ~(2048 | 4096);
3154    }
3155    if (barrierOptions & 8) {
3156      barrierOptions &= ~8;
3157      scale = true;
3158    }
3159    // If quadratic force KKT
3160    if (quadraticObj) {
3161      doKKT = true;
3162    }
3163    switch (barrierOptions) {
3164    case 0:
3165    default:
3166      if (!doKKT) {
3167        ClpCholeskyBase *cholesky = new ClpCholeskyBase(options.getExtraInfo(1));
3168        cholesky->setIntegerParameter(0, speed);
3169        barrier.setCholesky(cholesky);
3170      } else {
3171        ClpCholeskyBase *cholesky = new ClpCholeskyBase();
3172        cholesky->setKKT(true);
3173        barrier.setCholesky(cholesky);
3174      }
3175      break;
3176    case 1:
3177      if (!doKKT) {
3178        ClpCholeskyDense *cholesky = new ClpCholeskyDense();
3179        barrier.setCholesky(cholesky);
3180      } else {
3181        ClpCholeskyDense *cholesky = new ClpCholeskyDense();
3182        cholesky->setKKT(true);
3183        barrier.setCholesky(cholesky);
3184      }
3185      break;
3186#ifdef COIN_HAS_WSMP
3187    case 2: {
3188      if (!doKKT) {
3189        ClpCholeskyWssmp *cholesky = new ClpCholeskyWssmp(CoinMax(100, model2->numberRows() / 10));
3190        barrier.setCholesky(cholesky);
3191      } else {
3192        ClpCholeskyWssmpKKT *cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3193        barrier.setCholesky(cholesky);
3194      }
3195    } break;
3196    case 3:
3197      if (!doKKT) {
3198        ClpCholeskyWssmp *cholesky = new ClpCholeskyWssmp();
3199        barrier.setCholesky(cholesky);
3200      } else {
3201        ClpCholeskyWssmpKKT *cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3202        barrier.setCholesky(cholesky);
3203      }
3204      break;
3205#endif
3206#ifdef UFL_BARRIER
3207    case 4:
3208      if (!doKKT) {
3209        ClpCholeskyUfl *cholesky = new ClpCholeskyUfl(options.getExtraInfo(1));
3210        barrier.setCholesky(cholesky);
3211      } else {
3212        ClpCholeskyUfl *cholesky = new ClpCholeskyUfl();
3213        cholesky->setKKT(true);
3214        barrier.setCholesky(cholesky);
3215      }
3216      break;
3217#endif
3218#ifdef TAUCS_BARRIER
3219    case 5: {
3220      ClpCholeskyTaucs *cholesky = new ClpCholeskyTaucs();
3221      barrier.setCholesky(cholesky);
3222      assert(!doKKT);
3223    } break;
3224#endif
3225#ifdef COIN_HAS_MUMPS
3226    case 6: {
3227      if (!doKKT) {
3228        int logLevel = 0;
3229        if (barrier.logLevel() > 3) {
3230          logLevel = (barrier.logLevel() > 4) ? 2 : 1;
3231        }
3232        ClpCholeskyMumps *cholesky = new ClpCholeskyMumps(-1, logLevel);
3233        barrier.setCholesky(cholesky);
3234      } else {
3235        ClpCholeskyMumps *cholesky = new ClpCholeskyMumps();
3236        cholesky->setKKT(true);
3237        barrier.setCholesky(cholesky);
3238      }
3239    } break;
3240#endif
3241#ifdef PARDISO_BARRIER
3242    case 7: {
3243      ClpCholeskyPardiso *cholesky = new ClpCholeskyPardiso();
3244      barrier.setCholesky(cholesky);
3245      assert(!doKKT);
3246    } break;
3247#endif
3248    }
3249    int numberRows = model2->numberRows();
3250    int numberColumns = model2->numberColumns();
3251    int saveMaxIts = model2->maximumIterations();
3252    if (saveMaxIts < 1000) {
3253      barrier.setMaximumBarrierIterations(saveMaxIts);
3254      model2->setMaximumIterations(10000000);
3255    }
3256#ifndef SAVEIT
3257    //barrier.setDiagonalPerturbation(1.0e-25);
3258    if (aggressiveGamma) {
3259      switch (aggressiveGamma) {
3260      case 1:
3261        barrier.setGamma(1.0e-5);
3262        barrier.setDelta(1.0e-5);
3263        break;
3264      case 2:
3265        barrier.setGamma(1.0e-7);
3266        break;
3267      case 3:
3268        barrier.setDelta(1.0e-5);
3269        break;
3270      case 4:
3271        barrier.setGamma(1.0e-3);
3272        barrier.setDelta(1.0e-3);
3273        break;
3274      case 5:
3275        barrier.setGamma(1.0e-3);
3276        break;
3277      case 6:
3278        barrier.setDelta(1.0e-3);
3279        break;
3280      }
3281    }
3282    if (scale)
3283      barrier.scaling(1);
3284    else
3285      barrier.scaling(0);
3286    barrier.primalDual();
3287#elif SAVEIT == 1
3288    barrier.primalDual();
3289#else
3290    model2->restoreModel("xx.save");
3291    // move solutions
3292    CoinMemcpyN(model2->primalRowSolution(),
3293      numberRows, barrier.primalRowSolution());
3294    CoinMemcpyN(model2->dualRowSolution(),
3295      numberRows, barrier.dualRowSolution());
3296    CoinMemcpyN(model2->primalColumnSolution(),
3297      numberColumns, barrier.primalColumnSolution());
3298    CoinMemcpyN(model2->dualColumnSolution(),
3299      numberColumns, barrier.dualColumnSolution());
3300#endif
3301    time2 = CoinCpuTime();
3302    timeCore = time2 - timeX;
3303    handler_->message(CLP_INTERVAL_TIMING, messages_)
3304      << "Barrier" << timeCore << time2 - time1
3305      << CoinMessageEol;
3306    timeX = time2;
3307    int maxIts = barrier.maximumBarrierIterations();
3308    int barrierStatus = barrier.status();
3309    double gap = barrier.complementarityGap();
3310    // get which variables are fixed
3311    double *saveLower = NULL;
3312    double *saveUpper = NULL;
3313    ClpPresolve pinfo2;
3314    ClpSimplex *saveModel2 = NULL;
3315    bool extraPresolve = false;
3316    int numberFixed = barrier.numberFixed();
3317    if (numberFixed) {
3318      int numberRows = barrier.numberRows();
3319      int numberColumns = barrier.numberColumns();
3320      int numberTotal = numberRows + numberColumns;
3321      saveLower = new double[numberTotal];
3322      saveUpper = new double[numberTotal];
3323      CoinMemcpyN(barrier.columnLower(), numberColumns, saveLower);
3324      CoinMemcpyN(barrier.rowLower(), numberRows, saveLower + numberColumns);
3325      CoinMemcpyN(barrier.columnUpper(), numberColumns, saveUpper);
3326      CoinMemcpyN(barrier.rowUpper(), numberRows, saveUpper + numberColumns);
3327    }
3328    if (((numberFixed * 20 > barrier.numberRows() && numberFixed > 5000) || forceFixing) && presolveInCrossover) {
3329      // may as well do presolve
3330      if (!forceFixing) {
3331        barrier.fixFixed();
3332      } else {
3333        // Fix
3334        int n = barrier.numberColumns();
3335        double *lower = barrier.columnLower();
3336        double *upper = barrier.columnUpper();
3337        double *solution = barrier.primalColumnSolution();
3338        int nFix = 0;
3339        for (int i = 0; i < n; i++) {
3340          if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3341            double value = solution[i];
3342            if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3343              solution[i] = lower[i];
3344              upper[i] = lower[i];
3345              nFix++;
3346            } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3347              solution[i] = upper[i];
3348              lower[i] = upper[i];
3349              nFix++;
3350            }
3351          }
3352        }
3353#ifdef CLP_INVESTIGATE
3354        printf("%d columns fixed\n", nFix);
3355#endif
3356        int nr = barrier.numberRows();
3357        lower = barrier.rowLower();
3358        upper = barrier.rowUpper();
3359        solution = barrier.primalRowSolution();
3360        nFix = 0;
3361        for (int i = 0; i < nr; i++) {
3362          if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3363            double value = solution[i];
3364            if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3365              solution[i] = lower[i];
3366              upper[i] = lower[i];
3367              nFix++;
3368            } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3369              solution[i] = upper[i];
3370              lower[i] = upper[i];
3371              nFix++;
3372            }
3373          }
3374        }
3375#ifdef CLP_INVESTIGATE
3376        printf("%d row slacks fixed\n", nFix);
3377#endif
3378      }
3379      saveModel2 = model2;
3380      extraPresolve = true;
3381    } else if (numberFixed) {
3382      // Set fixed to bounds (may have restored earlier solution)
3383      if (!forceFixing) {
3384        barrier.fixFixed(false);
3385      } else {
3386        // Fix
3387        int n = barrier.numberColumns();
3388        double *lower = barrier.columnLower();
3389        double *upper = barrier.columnUpper();
3390        double *solution = barrier.primalColumnSolution();
3391        int nFix = 0;
3392        for (int i = 0; i < n; i++) {
3393          if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3394            double value = solution[i];
3395            if (value < lower[i] + 1.0e-8 && value - lower[i] < upper[i] - value) {
3396              solution[i] = lower[i];
3397              upper[i] = lower[i];
3398              nFix++;
3399            } else if (value > upper[i] - 1.0e-8 && value - lower[i] > upper[i] - value) {
3400              solution[i] = upper[i];
3401              lower[i] = upper[i];
3402              nFix++;
3403            } else {
3404              //printf("fixcol %d %g <= %g <= %g\n",
3405              //     i,lower[i],solution[i],upper[i]);
3406            }
3407          }
3408        }
3409#ifdef CLP_INVESTIGATE
3410        printf("%d columns fixed\n", nFix);
3411#endif
3412        int nr = barrier.numberRows();
3413        lower = barrier.rowLower();
3414        upper = barrier.rowUpper();
3415        solution = barrier.primalRowSolution();
3416        nFix = 0;
3417        for (int i = 0; i < nr; i++) {
3418          if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3419            double value = solution[i];
3420            if (value < lower[i] + 1.0e-5 && value - lower[i] < upper[i] - value) {
3421              solution[i] = lower[i];
3422              upper[i] = lower[i];
3423              nFix++;
3424            } else if (value > upper[i] - 1.0e-5 && value - lower[i] > upper[i] - value) {
3425              solution[i] = upper[i];
3426              lower[i] = upper[i];
3427              nFix++;
3428            } else {
3429              //printf("fixrow %d %g <= %g <= %g\n",
3430              //     i,lower[i],solution[i],upper[i]);
3431            }
3432          }
3433        }
3434#ifdef CLP_INVESTIGATE
3435        printf("%d row slacks fixed\n", nFix);
3436#endif
3437      }
3438    }
3439#ifdef BORROW
3440    int saveNumberIterations = barrier.numberIterations();
3441    barrier.returnModel(*model2);
3442    double *rowPrimal = new double[numberRows];
3443    double *columnPrimal = new double[numberColumns];
3444    double *rowDual = new double[numberRows];
3445    double *columnDual = new double[numberColumns];
3446    // move solutions other way
3447    CoinMemcpyN(model2->primalRowSolution(),
3448      numberRows, rowPrimal);
3449    CoinMemcpyN(model2->dualRowSolution(),
3450      numberRows, rowDual);
3451    CoinMemcpyN(model2->primalColumnSolution(),
3452      numberColumns, columnPrimal);
3453    CoinMemcpyN(model2->dualColumnSolution(),
3454      numberColumns, columnDual);
3455#else
3456    double *rowPrimal = barrier.primalRowSolution();
3457    double *columnPrimal = barrier.primalColumnSolution();
3458    double *rowDual = barrier.dualRowSolution();
3459    double *columnDual = barrier.dualColumnSolution();
3460    // move solutions
3461    CoinMemcpyN(rowPrimal,
3462      numberRows, model2->primalRowSolution());
3463    CoinMemcpyN(rowDual,
3464      numberRows, model2->dualRowSolution());
3465    CoinMemcpyN(columnPrimal,
3466      numberColumns, model2->primalColumnSolution());
3467    CoinMemcpyN(columnDual,
3468      numberColumns, model2->dualColumnSolution());
3469#endif
3470    if (saveModel2) {
3471      // do presolve
3472      model2 = pinfo2.presolvedModel(*model2, dblParam_[ClpPresolveTolerance],
3473        false, 5, true);
3474      if (!model2) {
3475        model2 = saveModel2;
3476        saveModel2 = NULL;
3477        int numberRows = model2->numberRows();
3478        int numberColumns = model2->numberColumns();
3479        CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3480        CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3481        delete[] saveLower;
3482        CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3483        CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3484        delete[] saveUpper;
3485        saveLower = NULL;
3486        saveUpper = NULL;
3487      }
3488    }
3489    if (method == ClpSolve::useBarrier || barrierStatus < 0) {
3490      if (maxIts && barrierStatus < 4 && !quadraticObj) {
3491        //printf("***** crossover - needs more thought on difficult models\n");
3492#if SAVEIT == 1
3493        model2->ClpSimplex::saveModel("xx.save");
3494#endif
3495        // make sure no status left
3496        model2->createStatus();
3497        // solve
3498        if (!forceFixing)
3499          model2->setPerturbation(100);
3500        if (model2->factorizationFrequency() == 200) {
3501          // User did not touch preset
3502          model2->defaultFactorizationFrequency();
3503        }
3504#if 1 //ndef ABC_INHERIT //#if 1 \
3505  // throw some into basis
3506        if (!forceFixing) {
3507          int numberRows = model2->numberRows();
3508          int numberColumns = model2->numberColumns();
3509          double *dsort = new double[numberColumns];
3510          int *sort = new int[numberColumns];
3511          int n = 0;
3512          const double *columnLower = model2->columnLower();
3513          const double *columnUpper = model2->columnUpper();
3514          double *primalSolution = model2->primalColumnSolution();
3515          const double *dualSolution = model2->dualColumnSolution();
3516          double tolerance = 10.0 * primalTolerance_;
3517          int i;
3518          for (i = 0; i < numberRows; i++)
3519            model2->setRowStatus(i, superBasic);
3520          for (i = 0; i < numberColumns; i++) {
3521            double distance = CoinMin(columnUpper[i] - primalSolution[i],
3522              primalSolution[i] - columnLower[i]);
3523            if (distance > tolerance) {
3524              if (fabs(dualSolution[i]) < 1.0e-5)
3525                distance *= 100.0;
3526              dsort[n] = -distance;
3527              sort[n++] = i;
3528              model2->setStatus(i, superBasic);
3529            } else if (distance > primalTolerance_) {
3530              model2->setStatus(i, superBasic);
3531            } else if (primalSolution[i] <= columnLower[i] + primalTolerance_) {
3532              model2->setStatus(i, atLowerBound);
3533              primalSolution[i] = columnLower[i];
3534            } else {
3535              model2->setStatus(i, atUpperBound);
3536              primalSolution[i] = columnUpper[i];
3537            }
3538          }
3539          CoinSort_2(dsort, dsort + n, sort);
3540          n = CoinMin(numberRows, n);
3541          for (i = 0; i < n; i++) {
3542            int iColumn = sort[i];
3543            model2->setStatus(iColumn, basic);
3544          }
3545          delete[] sort;
3546          delete[] dsort;
3547          // model2->allSlackBasis();
3548          if (gap < 1.0e-3 * static_cast< double >(numberRows + numberColumns)) {
3549            if (saveUpper) {
3550              int numberRows = model2->numberRows();
3551              int numberColumns = model2->numberColumns();
3552              CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3553              CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3554              CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3555              CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3556              delete[] saveLower;
3557              delete[] saveUpper;
3558              saveLower = NULL;
3559              saveUpper = NULL;
3560            }
3561            //int numberRows = model2->numberRows();
3562            //int numberColumns = model2->numberColumns();
3563#ifdef ABC_INHERIT
3564            model2->checkSolution(0);
3565            printf("%d primal infeasibilities summing to %g\n",
3566              model2->numberPrimalInfeasibilities(),
3567              model2->sumPrimalInfeasibilities());
3568            model2->dealWithAbc(1, 1);
3569          }
3570        }
3571#else
3572            // just primal values pass
3573            double saveScale = model2->objectiveScale();
3574            model2->setObjectiveScale(1.0e-3);
3575            model2->primal(2);
3576            model2->setObjectiveScale(saveScale);
3577            // save primal solution and copy back dual
3578            CoinMemcpyN(model2->primalRowSolution(),
3579              numberRows, rowPrimal);
3580            CoinMemcpyN(rowDual,
3581              numberRows, model2->dualRowSolution());
3582            CoinMemcpyN(model2->primalColumnSolution(),
3583              numberColumns, columnPrimal);
3584            CoinMemcpyN(columnDual,
3585              numberColumns, model2->dualColumnSolution());
3586            //model2->primal(1);
3587            // clean up reduced costs and flag variables
3588            {
3589              double *dj = model2->dualColumnSolution();
3590              double *cost = model2->objective();
3591              double *saveCost = new double[numberColumns];
3592              CoinMemcpyN(cost, numberColumns, saveCost);
3593              double *saveLower = new double[numberColumns];
3594              double *lower = model2->columnLower();
3595              CoinMemcpyN(lower, numberColumns, saveLower);
3596              double *saveUpper = new double[numberColumns];
3597              double *upper = model2->columnUpper();
3598              CoinMemcpyN(upper, numberColumns, saveUpper);
3599              int i;
3600              double tolerance = 10.0 * dualTolerance_;
3601              for (i = 0; i < numberColumns; i++) {
3602                if (model2->getStatus(i) == basic) {
3603                  dj[i] = 0.0;
3604                } else if (model2->getStatus(i) == atLowerBound) {
3605                  if (optimizationDirection_ * dj[i] < tolerance) {
3606                    if (optimizationDirection_ * dj[i] < 0.0) {
3607                      //if (dj[i]<-1.0e-3)
3608                      //printf("bad dj at lb %d %g\n",i,dj[i]);
3609                      cost[i] -= dj[i];
3610                      dj[i] = 0.0;
3611                    }
3612                  } else {
3613                    upper[i] = lower[i];
3614                  }
3615                } else if (model2->getStatus(i) == atUpperBound) {
3616                  if (optimizationDirection_ * dj[i] > tolerance) {
3617                    if (optimizationDirection_ * dj[i] > 0.0) {
3618                      //if (dj[i]>1.0e-3)
3619                      //printf("bad dj at ub %d %g\n",i,dj[i]);
3620                      cost[i] -= dj[i];
3621                      dj[i] = 0.0;
3622                    }
3623                  } else {
3624                    lower[i] = upper[i];
3625                  }
3626                }
3627              }
3628              // just dual values pass
3629              //model2->setLogLevel(63);
3630              //model2->setFactorizationFrequency(1);
3631              if (!gap)
3632                model2->dual(2);
3633              CoinMemcpyN(saveCost, numberColumns, cost);
3634              delete[] saveCost;
3635              CoinMemcpyN(saveLower, numberColumns, lower);
3636              delete[] saveLower;
3637              CoinMemcpyN(saveUpper, numberColumns, upper);
3638              delete[] saveUpper;
3639            }
3640          }
3641          // and finish
3642          // move solutions
3643          CoinMemcpyN(rowPrimal,
3644            numberRows, model2->primalRowSolution());
3645          CoinMemcpyN(columnPrimal,
3646            numberColumns, model2->primalColumnSolution());
3647        }
3648        double saveScale = model2->objectiveScale();
3649        model2->setObjectiveScale(1.0e-3);
3650        model2->primal(2);
3651        model2->setObjectiveScale(saveScale);
3652        model2->primal(1);
3653#endif
3654#else
3655        // just primal
3656#ifdef ABC_INHERIT
3657        model2->checkSolution(0);
3658        printf("%d primal infeasibilities summing to %g\n",
3659          model2->numberPrimalInfeasibilities(),
3660          model2->sumPrimalInfeasibilities());
3661        model2->dealWithAbc(1, 1);
3662#else
3663        model2->primal(1);
3664#endif
3665        //model2->primal(1);
3666#endif
3667      } else if (barrierStatus == 4) {
3668        // memory problems
3669        model2->setPerturbation(savePerturbation);
3670        model2->createStatus();
3671        model2->dual();
3672      } else if (maxIts && quadraticObj) {
3673        // make sure no status left
3674        model2->createStatus();
3675        // solve
3676        model2->setPerturbation(100);
3677        model2->reducedGradient(1);
3678      }
3679    }
3680
3681    //model2->setMaximumIterations(saveMaxIts);
3682#ifdef BORROW
3683    model2->setNumberIterations(model2->numberIterations() + saveNumberIterations);
3684    delete[] rowPrimal;
3685    delete[] columnPrimal;
3686    delete[] rowDual;
3687    delete[] columnDual;
3688#endif
3689    if (extraPresolve) {
3690      pinfo2.postsolve(true);
3691      delete model2;
3692      model2 = saveModel2;
3693    }
3694    if (saveUpper) {
3695      if (!forceFixing) {
3696        int numberRows = model2->numberRows();
3697        int numberColumns = model2->numberColumns();
3698        CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3699        CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3700        CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3701        CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3702      }
3703      delete[] saveLower;
3704      delete[] saveUpper;
3705      saveLower = NULL;
3706      saveUpper = NULL;
3707      if (method != ClpSolve::useBarrierNoCross)
3708        model2->primal(1);
3709    }
3710    model2->setPerturbation(savePerturbation);
3711    time2 = CoinCpuTime();
3712    timeCore = time2 - timeX;
3713    handler_->message(CLP_INTERVAL_TIMING, messages_)
3714      << "Crossover" << timeCore << time2 - time1
3715      << CoinMessageEol;
3716    timeX = time2;
3717#else
3718    abort();
3719#endif
3720  } else if (method == ClpSolve::notImplemented) {
3721    printf("done decomposition\n");
3722  } else {
3723    assert(method != ClpSolve::automatic); // later
3724    time2 = 0.0;
3725  }
3726  if (saveMatrix) {
3727    if (model2 == this) {
3728      // delete and replace
3729      delete model2->clpMatrix();
3730      model2->replaceMatrix(saveMatrix);
3731    } else {
3732      delete saveMatrix;
3733    }
3734  }
3735  numberIterations = model2->numberIterations();
3736  finalStatus = model2->status();
3737  int finalSecondaryStatus = model2->secondaryStatus();
3738  if (presolve == ClpSolve::presolveOn) {
3739    int saveLevel = logLevel();
3740    if ((specialOptions_ & 1024) == 0)
3741      setLogLevel(CoinMin(1, saveLevel));
3742    else
3743      setLogLevel(CoinMin(0, saveLevel));
3744    pinfo->postsolve(true);
3745    numberIterations_ = 0;
3746    delete pinfo;
3747    pinfo = NULL;
3748    factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
3749    time2 = CoinCpuTime();
3750    timePresolve += time2 - timeX;
3751    handler_->message(CLP_INTERVAL_TIMING, messages_)
3752      << "Postsolve" << time2 - timeX << time2 - time1
3753      << CoinMessageEol;
3754    timeX = time2;
3755    if (!presolveToFile) {
3756#if 1 //ndef ABC_INHERIT
3757      delete model2;
3758#else
3759      if (model2->abcSimplex())
3760        delete model2->abcSimplex();
3761      else
3762        delete model2;
3763#endif
3764    }
3765    if (interrupt)
3766      currentModel = this;
3767    // checkSolution(); already done by postSolve
3768    setLogLevel(saveLevel);
3769    int oldStatus = problemStatus_;
3770    setProblemStatus(finalStatus);
3771    setSecondaryStatus(finalSecondaryStatus);
3772    /* Code modified so rcode -1 as normal, >0 clean up, 0 say optimal */
3773    int rcode = eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
3774    //#define TREAT_AS_OPTIMAL_TOLERANCE 1.0e-4
3775#ifdef TREAT_AS_OPTIMAL_TOLERANCE
3776    if (rcode == -1 && sumPrimalInfeasibilities_ < TREAT_AS_OPTIMAL_TOLERANCE && sumDualInfeasibilities_ < TREAT_AS_OPTIMAL_TOLERANCE)
3777      rcode = 0;
3778#endif
3779    if (finalStatus != 3 && rcode < 0 && (finalStatus || oldStatus == -1)) {
3780      double sumPrimal = sumPrimalInfeasibilities_;
3781      double sumDual = sumDualInfeasibilities_;
3782      if (sumDual > 1.0e-6 && sumPrimal > 1.0e-6)
3783        moreSpecialOptions_ &= ~2; // be safe and do final solve
3784      // ignore some parts of solution
3785      if (finalStatus == 1) {
3786        // infeasible
3787        sumDual = 0.0;
3788      } else if (finalStatus == 2) {
3789        sumPrimal = 0.0;
3790      }
3791      int savePerturbation = perturbation();
3792      if (savePerturbation == 50)
3793        setPerturbation(51); // small
3794      if (!finalStatus || finalStatus == 2 || (moreSpecialOptions_ & 2) == 0 || fabs(sumDual) + fabs(sumPrimal) < 1.0e-3) {
3795        if (finalStatus == 2) {
3796          if (sumDual > 1.0e-4) {
3797            // unbounded - get feasible first
3798            double save = optimizationDirection_;
3799            optimizationDirection_ = 0.0;
3800            primal(1);
3801            optimizationDirection_ = save;
3802          }
3803          primal(1);
3804        } else if (finalStatus == 1) {
3805          dual();
3806        } else {
3807          if ((moreSpecialOptions_ & 65536) == 0) {
3808            if (numberRows_ < 10000 || true)
3809              setPerturbation(100); // probably better to perturb after n its
3810            else if (savePerturbation < 100)
3811              setPerturbation(51); // probably better to perturb after n its
3812          }
3813#ifndef ABC_INHERIT
3814          // use method thought suitable
3815          int numberSuperBasic = 0;
3816          for (int i = 0; i < numberColumns_; i++) {
3817            if (getColumnStatus(i) == superBasic)
3818              numberSuperBasic++;
3819          }
3820          if (sumDual > 1000.0 * sumPrimal || numberSuperBasic) {
3821            primal(1);
3822          } else if (sumPrimal > 1000.0 * sumDual) {
3823            dual();
3824          } else {
3825            if (method != ClpSolve::useDual)
3826              primal(1);
3827            else
3828              dual();
3829          }
3830#else
3831          dealWithAbc(1, 2, interrupt);
3832#endif
3833        }
3834      } else {
3835        // just set status
3836        problemStatus_ = finalStatus;
3837      }
3838      setPerturbation(savePerturbation);
3839      numberIterations += numberIterations_;
3840      numberIterations_ = numberIterations;
3841      finalStatus = status();
3842      time2 = CoinCpuTime();
3843      handler_->message(CLP_INTERVAL_TIMING, messages_)
3844        << "Cleanup" << time2 - timeX << time2 - time1
3845        << CoinMessageEol;
3846      timeX = time2;
3847    } else if (rcode > 0) { // was >= 0
3848#ifdef ABC_INHERIT
3849      dealWithAbc(1, 2, true);
3850#else
3851      primal(1);
3852#endif
3853    } else {
3854      secondaryStatus_ = finalSecondaryStatus;
3855    }
3856  } else if (model2 != this) {
3857    // not presolved - but different model used (sprint probably)
3858    CoinMemcpyN(model2->primalRowSolution(),
3859      numberRows_, this->primalRowSolution());
3860    CoinMemcpyN(model2->dualRowSolution(),
3861      numberRows_, this->dualRowSolution());
3862    CoinMemcpyN(model2->primalColumnSolution(),
3863      numberColumns_, this->primalColumnSolution());
3864    CoinMemcpyN(model2->dualColumnSolution(),
3865      numberColumns_, this->dualColumnSolution());
3866    CoinMemcpyN(model2->statusArray(),
3867      numberColumns_ + numberRows_, this->statusArray());
3868    objectiveValue_ = model2->objectiveValue_;
3869    numberIterations_ = model2->numberIterations_;
3870    problemStatus_ = model2->problemStatus_;
3871    secondaryStatus_ = model2->secondaryStatus_;
3872    delete model2;
3873  }
3874  if (method != ClpSolve::useBarrierNoCross && method != ClpSolve::useBarrier)
3875    setMaximumIterations(saveMaxIterations);
3876  std::string statusMessage[] = { "Unknown", "Optimal", "PrimalInfeasible", "DualInfeasible", "Stopped",
3877    "Errors", "User stopped" };
3878  assert(finalStatus >= -1 && finalStatus <= 5);
3879  numberIterations_ = numberIterations;
3880  handler_->message(CLP_TIMING, messages_)
3881    << statusMessage[finalStatus + 1] << objectiveValue() << numberIterations << time2 - time1;
3882  handler_->printing(presolve == ClpSolve::presolveOn)
3883    << timePresolve;
3884  handler_->printing(timeIdiot != 0.0)
3885    << timeIdiot;
3886  handler_->message() << CoinMessageEol;
3887  if (interrupt)
3888    signal(SIGINT, saveSignal);
3889  perturbation_ = savePerturbation;
3890  scalingFlag_ = saveScaling;
3891  // If faking objective - put back correct one
3892  if (savedObjective) {
3893    delete objective_;
3894    objective_ = savedObjective;
3895  }
3896  if (options.getSpecialOption(1) == 2 && options.getExtraInfo(1) > 1000000) {
3897    ClpObjective *savedObjective = objective_;
3898    // make up zero objective
3899    double *obj = new double[numberColumns_];
3900    for (int i = 0; i < numberColumns_; i++)
3901      obj[i] = 0.0;
3902    objective_ = new ClpLinearObjective(obj, numberColumns_);
3903    delete[] obj;
3904    primal(1);
3905    delete objective_;
3906    objective_ = savedObjective;
3907    finalStatus = status();
3908  }
3909  eventHandler()->event(ClpEventHandler::presolveEnd);
3910  delete pinfo;
3911  moreSpecialOptions_ = saveMoreOptions;
3912#ifdef CLP_USEFUL_PRINTOUT
3913  debugInt[23] = numberIterations_;
3914#endif
3915  return finalStatus;
3916}
3917// General solve
3918int ClpSimplex::initialSolve()
3919{
3920  // Default so use dual
3921  ClpSolve options;
3922  return initialSolve(options);
3923}
3924// General dual solve
3925int ClpSimplex::initialDualSolve()
3926{
3927  ClpSolve options;
3928  // Use dual
3929  options.setSolveType(ClpSolve::useDual);
3930  return initialSolve(options);
3931}
3932// General primal solve
3933int ClpSimplex::initialPrimalSolve()
3934{
3935  ClpSolve options;
3936  // Use primal
3937  options.setSolveType(ClpSolve::usePrimal);
3938  return initialSolve(options);
3939}
3940// barrier solve, not to be followed by crossover
3941int ClpSimplex::initialBarrierNoCrossSolve()
3942{
3943  ClpSolve options;
3944  // Use primal
3945  options.setSolveType(ClpSolve::useBarrierNoCross);
3946  return initialSolve(options);
3947}
3948
3949// General barrier solve
3950int ClpSimplex::initialBarrierSolve()
3951{
3952  ClpSolve options;
3953  // Use primal
3954  options.setSolveType(ClpSolve::useBarrier);
3955  return initialSolve(options);
3956}
3957
3958// Default constructor
3959ClpSolve::ClpSolve()
3960{
3961  method_ = automatic;
3962  presolveType_ = presolveOn;
3963  numberPasses_ = 5;
3964  int i;
3965  for (i = 0; i < 7; i++)
3966    options_[i] = 0;
3967  // say no +-1 matrix
3968  options_[3] = 1;
3969  for (i = 0; i < 7; i++)
3970    extraInfo_[i] = -1;
3971  independentOptions_[0] = 0;
3972  // But switch off slacks
3973  independentOptions_[1] = 512;
3974  // Substitute up to 3
3975  independentOptions_[2] = 3;
3976}
3977// Constructor when you really know what you are doing
3978ClpSolve::ClpSolve(SolveType method, PresolveType presolveType,
3979  int numberPasses, int options[6],
3980  int extraInfo[6], int independentOptions[3])
3981{
3982  method_ = method;
3983  presolveType_ = presolveType;
3984  numberPasses_ = numberPasses;
3985  int i;
3986  for (i = 0; i < 6; i++)
3987    options_[i] = options[i];
3988  options_[6] = 0;
3989  for (i = 0; i < 6; i++)
3990    extraInfo_[i] = extraInfo[i];
3991  extraInfo_[6] = 0;
3992  for (i = 0; i < 3; i++)
3993    independentOptions_[i] = independentOptions[i];
3994}
3995
3996// Copy constructor.
3997ClpSolve::ClpSolve(const ClpSolve &rhs)
3998{
3999  method_ = rhs.method_;
4000  presolveType_ = rhs.presolveType_;
4001  numberPasses_ = rhs.numberPasses_;
4002  int i;
4003  for (i = 0; i < 7; i++)
4004    options_[i] = rhs.options_[i];
4005  for (i = 0; i < 7; i++)
4006    extraInfo_[i] = rhs.extraInfo_[i];
4007  for (i = 0; i < 3; i++)
4008    independentOptions_[i] = rhs.independentOptions_[i];
4009}
4010// Assignment operator. This copies the data
4011ClpSolve &
4012ClpSolve::operator=(const ClpSolve &rhs)
4013{
4014  if (this != &rhs) {
4015    method_ = rhs.method_;
4016    presolveType_ = rhs.presolveType_;
4017    numberPasses_ = rhs.numberPasses_;
4018    int i;
4019    for (i = 0; i < 7; i++)
4020      options_[i] = rhs.options_[i];
4021    for (i = 0; i < 7; i++)
4022      extraInfo_[i] = rhs.extraInfo_[i];
4023    for (i = 0; i < 3; i++)
4024      independentOptions_[i] = rhs.independentOptions_[i];
4025  }
4026  return *this;
4027}
4028// Destructor
4029ClpSolve::~ClpSolve()
4030{
4031}
4032// See header file for details
4033void ClpSolve::setSpecialOption(int which, int value, int extraInfo)
4034{
4035  options_[which] = value;
4036  extraInfo_[which] = extraInfo;
4037}
4038int ClpSolve::getSpecialOption(int which) const
4039{
4040  return options_[which];
4041}
4042
4043// Solve types
4044void ClpSolve::setSolveType(SolveType method, int /*extraInfo*/)
4045{
4046  method_ = method;
4047}
4048
4049ClpSolve::SolveType
4050ClpSolve::getSolveType()
4051{
4052  return method_;
4053}
4054
4055// Presolve types
4056void ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
4057{
4058  presolveType_ = amount;
4059  numberPasses_ = extraInfo;
4060}
4061ClpSolve::PresolveType
4062ClpSolve::getPresolveType()
4063{
4064  return presolveType_;
4065}
4066// Extra info for idiot (or sprint)
4067int ClpSolve::getExtraInfo(int which) const
4068{
4069  return extraInfo_[which];
4070}
4071int ClpSolve::getPresolvePasses() const
4072{
4073  return numberPasses_;
4074}
4075/* Say to return at once if infeasible,
4076   default is to solve */
4077void ClpSolve::setInfeasibleReturn(bool trueFalse)
4078{
4079  independentOptions_[0] = trueFalse ? 1 : 0;
4080}
4081#include <string>
4082// Generates code for above constructor
4083void ClpSolve::generateCpp(FILE *fp)
4084{
4085  std::string solveType[] = {
4086    "ClpSolve::useDual",
4087    "ClpSolve::usePrimal",
4088    "ClpSolve::usePrimalorSprint",
4089    "ClpSolve::useBarrier",
4090    "ClpSolve::useBarrierNoCross",
4091    "ClpSolve::automatic",
4092    "ClpSolve::notImplemented"
4093  };
4094  std::string presolveType[] = {
4095    "ClpSolve::presolveOn",
4096    "ClpSolve::presolveOff",
4097    "ClpSolve::presolveNumber",
4098    "ClpSolve::presolveNumberCost"
4099  };
4100  fprintf(fp, "3  ClpSolve::SolveType method = %s;\n", solveType[method_].c_str());
4101  fprintf(fp, "3  ClpSolve::PresolveType presolveType = %s;\n",
4102    presolveType[presolveType_].c_str());
4103  fprintf(fp, "3  int numberPasses = %d;\n", numberPasses_);
4104  fprintf(fp, "3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
4105    options_[0], options_[1], options_[2],
4106    options_[3], options_[4], options_[5]);
4107  fprintf(fp, "3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
4108    extraInfo_[0], extraInfo_[1], extraInfo_[2],
4109    extraInfo_[3], extraInfo_[4], extraInfo_[5]);
4110  fprintf(fp, "3  int independentOptions[] = {%d,%d,%d};\n",
4111    independentOptions_[0], independentOptions_[1], independentOptions_[2]);
4112  fprintf(fp, "3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
4113  fprintf(fp, "3                    options,extraInfo,independentOptions);\n");
4114}
4115//#############################################################################
4116#include "ClpNonLinearCost.hpp"
4117
4118ClpSimplexProgress::ClpSimplexProgress()
4119{
4120  int i;
4121  for (i = 0; i < CLP_PROGRESS; i++) {
4122    objective_[i] = COIN_DBL_MAX * 1.0e-50;
4123    infeasibility_[i] = -1.0; // set to an impossible value
4124    realInfeasibility_[i] = COIN_DBL_MAX * 1.0e-50;
4125    numberInfeasibilities_[i] = -1;
4126    iterationNumber_[i] = -1;
4127  }
4128#ifdef CLP_PROGRESS_WEIGHT
4129  for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4130    objectiveWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4131    infeasibilityWeight_[i] = -1.0; // set to an impossible value
4132    realInfeasibilityWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4133    numberInfeasibilitiesWeight_[i] = -1;
4134    iterationNumberWeight_[i] = -1;
4135  }
4136  drop_ = 0.0;
4137  best_ = 0.0;
4138#endif
4139  initialWeight_ = 0.0;
4140  for (i = 0; i < CLP_CYCLE; i++) {
4141    //obj_[i]=COIN_DBL_MAX*1.0e-50;
4142    in_[i] = -1;
4143    out_[i] = -1;
4144    way_[i] = 0;
4145  }
4146  numberTimes_ = 0;
4147  numberBadTimes_ = 0;
4148  numberReallyBadTimes_ = 0;
4149  numberTimesFlagged_ = 0;
4150  model_ = NULL;
4151  oddState_ = 0;
4152}
4153
4154//-----------------------------------------------------------------------------
4155
4156ClpSimplexProgress::~ClpSimplexProgress()
4157{
4158}
4159// Copy constructor.
4160ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
4161{
4162  int i;
4163  for (i = 0; i < CLP_PROGRESS; i++) {
4164    objective_[i] = rhs.objective_[i];
4165    infeasibility_[i] = rhs.infeasibility_[i];
4166    realInfeasibility_[i] = rhs.realInfeasibility_[i];
4167    numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
4168    iterationNumber_[i] = rhs.iterationNumber_[i];
4169  }
4170#ifdef CLP_PROGRESS_WEIGHT
4171  for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4172    objectiveWeight_[i] = rhs.objectiveWeight_[i];
4173    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4174    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4175    numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4176    iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4177  }
4178  drop_ = rhs.drop_;
4179  best_ = rhs.best_;
4180#endif
4181  initialWeight_ = rhs.initialWeight_;
4182  for (i = 0; i < CLP_CYCLE; i++) {
4183    //obj_[i]=rhs.obj_[i];
4184    in_[i] = rhs.in_[i];
4185    out_[i] = rhs.out_[i];
4186    way_[i] = rhs.way_[i];
4187  }
4188  numberTimes_ = rhs.numberTimes_;
4189  numberBadTimes_ = rhs.numberBadTimes_;
4190  numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4191  numberTimesFlagged_ = rhs.numberTimesFlagged_;
4192  model_ = rhs.model_;
4193  oddState_ = rhs.oddState_;
4194}
4195// Copy constructor.from model
4196ClpSimplexProgress::ClpSimplexProgress(ClpSimplex *model)
4197{
4198  model_ = model;
4199  reset();
4200  initialWeight_ = 0.0;
4201}
4202// Fill from model
4203void ClpSimplexProgress::fillFromModel(ClpSimplex *model)
4204{
4205  model_ = model;
4206  reset();
4207  initialWeight_ = 0.0;
4208}
4209// Assignment operator. This copies the data
4210ClpSimplexProgress &
4211ClpSimplexProgress::operator=(const ClpSimplexProgress &rhs)
4212{
4213  if (this != &rhs) {
4214    int i;
4215    for (i = 0; i < CLP_PROGRESS; i++) {
4216      objective_[i] = rhs.objective_[i];
4217      infeasibility_[i] = rhs.infeasibility_[i];
4218      realInfeasibility_[i] = rhs.realInfeasibility_[i];
4219      numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
4220      iterationNumber_[i] = rhs.iterationNumber_[i];
4221    }
4222#ifdef CLP_PROGRESS_WEIGHT
4223    for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4224      objectiveWeight_[i] = rhs.objectiveWeight_[i];
4225      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4226      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4227      numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4228      iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4229    }
4230    drop_ = rhs.drop_;
4231    best_ = rhs.best_;
4232#endif
4233    initialWeight_ = rhs.initialWeight_;
4234    for (i = 0; i < CLP_CYCLE; i++) {
4235      //obj_[i]=rhs.obj_[i];
4236      in_[i] = rhs.in_[i];
4237      out_[i] = rhs.out_[i];
4238      way_[i] = rhs.way_[i];
4239    }
4240    numberTimes_ = rhs.numberTimes_;
4241    numberBadTimes_ = rhs.numberBadTimes_;
4242    numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4243    numberTimesFlagged_ = rhs.numberTimesFlagged_;
4244    model_ = rhs.model_;
4245    oddState_ = rhs.oddState_;
4246  }
4247  return *this;
4248}
4249// Seems to be something odd about exact comparison of doubles on linux
4250static bool equalDouble(double value1, double value2)
4251{
4252
4253  union {
4254    double d;
4255    int i[2];
4256  } v1, v2;
4257  v1.d = value1;
4258  v2.d = value2;
4259  if (sizeof(int) * 2 == sizeof(double))
4260    return (v1.i[0] == v2.i[0] && v1.i[1] == v2.i[1]);
4261  else
4262    return (v1.i[0] == v2.i[0]);
4263}
4264int ClpSimplexProgress::looping()
4265{
4266  if (!model_)
4267    return -1;
4268  double objective;
4269  if (model_->algorithm() < 0) {
4270    objective = model_->rawObjectiveValue();
4271    objective -= model_->bestPossibleImprovement();
4272  } else {
4273    objective = model_->nonLinearCost()->feasibleReportCost();
4274  }
4275  double infeasibility;
4276  double realInfeasibility = 0.0;
4277  int numberInfeasibilities;
4278  int iterationNumber = model_->numberIterations();
4279  //numberTimesFlagged_ = 0;
4280  if (model_->algorithm() < 0) {
4281    // dual
4282    infeasibility = model_->sumPrimalInfeasibilities();
4283    numberInfeasibilities = model_->numberPrimalInfeasibilities();
4284  } else {
4285    //primal
4286    infeasibility = model_->sumDualInfeasibilities();
4287    realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
4288    numberInfeasibilities = model_->numberDualInfeasibilities();
4289  }
4290  int i;
4291  int numberMatched = 0;
4292  int matched = 0;
4293  int nsame = 0;
4294  for (i = 0; i < CLP_PROGRESS; i++) {
4295    bool matchedOnObjective = equalDouble(objective, objective_[i]);
4296    bool matchedOnInfeasibility = equalDouble(infeasibility, infeasibility_[i]);
4297    bool matchedOnInfeasibilities = (numberInfeasibilities == numberInfeasibilities_[i]);
4298
4299    if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
4300      matched |= (1 << i);
4301      // Check not same iteration
4302      if (iterationNumber != iterationNumber_[i]) {
4303        numberMatched++;
4304        // here mainly to get over compiler bug?
4305        if (model_->messageHandler()->logLevel() > 10)
4306          printf("%d %d %d %d %d loop check\n", i, numberMatched,
4307            matchedOnObjective, matchedOnInfeasibility,
4308            matchedOnInfeasibilities);
4309      } else {
4310        // stuck but code should notice
4311        nsame++;
4312      }
4313    }
4314    if (i) {
4315      objective_[i - 1] = objective_[i];
4316      infeasibility_[i - 1] = infeasibility_[i];
4317      realInfeasibility_[i - 1] = realInfeasibility_[i];
4318      numberInfeasibilities_[i - 1] = numberInfeasibilities_[i];
4319      iterationNumber_[i - 1] = iterationNumber_[i];
4320    }
4321  }
4322  objective_[CLP_PROGRESS - 1] = objective;
4323  infeasibility_[CLP_PROGRESS - 1] = infeasibility;
4324  realInfeasibility_[CLP_PROGRESS - 1] = realInfeasibility;
4325  numberInfeasibilities_[CLP_PROGRESS - 1] = numberInfeasibilities;
4326  iterationNumber_[CLP_PROGRESS - 1] = iterationNumber;
4327  if (nsame == CLP_PROGRESS)
4328    numberMatched = CLP_PROGRESS; // really stuck
4329  if (model_->progressFlag())
4330    numberMatched = 0;
4331  numberTimes_++;
4332  if (numberTimes_ < 10)
4333    numberMatched = 0;
4334  // skip if just last time as may be checking something
4335  if (matched == (1 << (CLP_PROGRESS - 1)))
4336    numberMatched = 0;
4337  if (numberMatched && model_->clpMatrix()->type() < 15) {
4338    model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
4339      << numberMatched
4340      << matched
4341      << numberTimes_
4342      << CoinMessageEol;
4343    numberBadTimes_++;
4344    if (numberBadTimes_ < 10) {
4345      // make factorize every iteration
4346      model_->forceFactorization(1);
4347      if (numberBadTimes_ < 2) {
4348        startCheck(); // clear other loop check
4349        if (model_->algorithm() < 0) {
4350          // dual - change tolerance
4351          model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
4352          // if infeasible increase dual bound
4353          if (model_->dualBound() < 1.0e17) {
4354            model_->setDualBound(model_->dualBound() * 1.1);
4355            static_cast< ClpSimplexDual * >(model_)->resetFakeBounds(0);
4356          }
4357        } else {
4358          // primal - change tolerance
4359          if (numberBadTimes_ > 3)
4360            model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
4361          // if infeasible increase infeasibility cost
4362          if (model_->nonLinearCost()->numberInfeasibilities() && model_->infeasibilityCost() < 1.0e17) {
4363            model_->setInfeasibilityCost(model_->infeasibilityCost() * 1.1);
4364          }
4365        }
4366      } else {
4367        // flag
4368        int iSequence;
4369        if (model_->algorithm() < 0) {
4370          // dual
4371          if (model_->dualBound() > 1.0e14)
4372            model_->setDualBound(1.0e14);
4373          iSequence = in_[CLP_CYCLE - 1];
4374        } else {
4375          // primal
4376          //if (model_->infeasibilityCost() > 1.0e14)
4377          //   model_->setInfeasibilityCost(1.0e14);
4378          iSequence = out_[CLP_CYCLE - 1];
4379        }
4380        if (iSequence >= 0) {
4381          char x = model_->isColumn(iSequence) ? 'C' : 'R';
4382          if (model_->messageHandler()->logLevel() >= 63)
4383            model_->messageHandler()->message(CLP_SIMPLEX_FLAG, model_->messages())
4384              << x << model_->sequenceWithin(iSequence)
4385              << CoinMessageEol;
4386          // if Gub then needs to be sequenceIn_
4387          int save = model_->sequenceIn();
4388          model_->setSequenceIn(iSequence);
4389          model_->setFlagged(iSequence);
4390          model_->setSequenceIn(save);
4391          //printf("flagging %d from loop\n",iSequence);
4392          startCheck();
4393        } else {
4394          // Give up
4395          if (model_->messageHandler()->logLevel() >= 63)
4396            printf("***** All flagged?\n");
4397          return 4;
4398        }
4399        // reset
4400        numberBadTimes_ = 2;
4401      }
4402      return -2;
4403    } else {
4404      // look at solution and maybe declare victory
4405      if (infeasibility < 1.0e-4) {
4406        return 0;
4407      } else {
4408        model_->messageHandler()->message(CLP_LOOP, model_->messages())
4409          << CoinMessageEol;
4410#ifndef NDEBUG
4411        printf("debug loop ClpSimplex A\n");
4412        abort();
4413#endif
4414        return 3;
4415      }
4416    }
4417  }
4418  return -1;
4419}
4420// Resets as much as possible
4421void ClpSimplexProgress::reset()
4422{
4423  int i;
4424  for (i = 0; i < CLP_PROGRESS; i++) {
4425    if (model_->algorithm() >= 0)
4426      objective_[i] = COIN_DBL_MAX * 1.0e-50;
4427    else
4428      objective_[i] = -COIN_DBL_MAX * 1.0e-50;
4429    infeasibility_[i] = -1.0; // set to an impossible value
4430    realInfeasibility_[i] = COIN_DBL_MAX * 1.0e-50;
4431    numberInfeasibilities_[i] = -1;
4432    iterationNumber_[i] = -1;
4433  }
4434#ifdef CLP_PROGRESS_WEIGHT
4435  for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4436    objectiveWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4437    infeasibilityWeight_[i] = -1.0; // set to an impossible value
4438    realInfeasibilityWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4439    numberInfeasibilitiesWeight_[i] = -1;
4440    iterationNumberWeight_[i] = -1;
4441  }
4442  drop_ = 0.0;
4443  best_ = 0.0;
4444#endif
4445  for (i = 0; i < CLP_CYCLE; i++) {
4446    //obj_[i]=COIN_DBL_MAX*1.0e-50;
4447    in_[i] = -1;
4448    out_[i] = -1;
4449    way_[i] = 0;
4450  }
4451  numberTimes_ = 0;
4452  numberBadTimes_ = 0;
4453  numberReallyBadTimes_ = 0;
4454  numberTimesFlagged_ = 0;
4455  oddState_ = 0;
4456}
4457// Returns previous objective (if -1) - current if (0)
4458double
4459ClpSimplexProgress::lastObjective(int back) const
4460{
4461  return objective_[CLP_PROGRESS - 1 - back];
4462}
4463// Returns previous infeasibility (if -1) - current if (0)
4464double
4465ClpSimplexProgress::lastInfeasibility(int back) const
4466{
4467  return realInfeasibility_[CLP_PROGRESS - 1 - back];
4468}
4469// Sets real primal infeasibility
4470void ClpSimplexProgress::setInfeasibility(double value)
4471{
4472  for (int i = 1; i < CLP_PROGRESS; i++)
4473    realInfeasibility_[i - 1] = realInfeasibility_[i];
4474  realInfeasibility_[CLP_PROGRESS - 1] = value;
4475}
4476// Returns number of primal infeasibilities (if -1) - current if (0)
4477int ClpSimplexProgress::numberInfeasibilities(int back) const
4478{
4479  return numberInfeasibilities_[CLP_PROGRESS - 1 - back];
4480}
4481// Modify objective e.g. if dual infeasible in dual
4482void ClpSimplexProgress::modifyObjective(double value)
4483{
4484  objective_[CLP_PROGRESS - 1] = value;
4485}
4486// Returns previous iteration number (if -1) - current if (0)
4487int ClpSimplexProgress::lastIterationNumber(int back) const
4488{
4489  return iterationNumber_[CLP_PROGRESS - 1 - back];
4490}
4491// clears iteration numbers (to switch off panic)
4492void ClpSimplexProgress::clearIterationNumbers()
4493{
4494  for (int i = 0; i < CLP_PROGRESS; i++)
4495    iterationNumber_[i] = -1;
4496}
4497// Start check at beginning of whileIterating
4498void ClpSimplexProgress::startCheck()
4499{
4500  int i;
4501  for (i = 0; i < CLP_CYCLE; i++) {
4502    //obj_[i]=COIN_DBL_MAX*1.0e-50;
4503    in_[i] = -1;
4504    out_[i] = -1;
4505    way_[i] = 0;
4506  }
4507}
4508// Returns cycle length in whileIterating
4509int ClpSimplexProgress::cycle(int in, int out, int wayIn, int wayOut)
4510{
4511  int i;
4512#if 0
4513     if (model_->numberIterations() > 206571) {
4514          printf("in %d out %d\n", in, out);
4515          for (i = 0; i < CLP_CYCLE; i++)
4516               printf("cy %d in %d out %d\n", i, in_[i], out_[i]);
4517     }
4518#endif
4519  int matched = 0;
4520  // first see if in matches any out
4521  for (i = 1; i < CLP_CYCLE; i++) {
4522    if (in == out_[i]) {
4523      // even if flip then suspicious
4524      matched = -1;
4525      break;
4526    }
4527  }
4528#if 0
4529     if (!matched || in_[0] < 0) {
4530          // can't be cycle
4531          for (i = 0; i < CLP_CYCLE - 1; i++) {
4532               //obj_[i]=obj_[i+1];
4533               in_[i] = in_[i+1];
4534               out_[i] = out_[i+1];
4535               way_[i] = way_[i+1];
4536          }
4537     } else {
4538          // possible cycle
4539          matched = 0;
4540          for (i = 0; i < CLP_CYCLE - 1; i++) {
4541               int k;
4542               char wayThis = way_[i];
4543               int inThis = in_[i];
4544               int outThis = out_[i];
4545               //double objThis = obj_[i];
4546               for(k = i + 1; k < CLP_CYCLE; k++) {
4547                    if (inThis == in_[k] && outThis == out_[k] && wayThis == way_[k]) {
4548                         int distance = k - i;
4549                         if (k + distance < CLP_CYCLE) {
4550                              // See if repeats
4551                              int j = k + distance;
4552                              if (inThis == in_[j] && outThis == out_[j] && wayThis == way_[j]) {
4553                                   matched = distance;
4554                                   break;
4555                              }
4556                         } else {
4557                              matched = distance;
4558                              break;
4559                         }
4560                    }
4561               }
4562               //obj_[i]=obj_[i+1];
4563               in_[i] = in_[i+1];
4564               out_[i] = out_[i+1];
4565               way_[i] = way_[i+1];
4566          }
4567     }
4568#else
4569  if (matched && in_[0] >= 0) {
4570    // possible cycle - only check [0] against all
4571    matched = 0;
4572    int nMatched = 0;
4573    char way0 = way_[0];
4574    int in0 = in_[0];
4575    int out0 = out_[0];
4576    //double obj0 = obj_[i];
4577    for (int k = 1; k < CLP_CYCLE - 4; k++) {
4578      if (in0 == in_[k] && out0 == out_[k] && way0 == way_[k]) {
4579        nMatched++;
4580        // See if repeats
4581        int end = CLP_CYCLE - k;
4582        int j;
4583        for (j = 1; j < end; j++) {
4584          if (in_[j + k] != in_[j] || out_[j + k] != out_[j] || way_[j + k] != way_[j])
4585            break;
4586        }
4587        if (j == end) {
4588          matched = k;
4589          break;
4590        }
4591      }
4592    }
4593    // If three times then that is too much even if not regular
4594    if (matched <= 0 && nMatched > 1)
4595      matched = 100;
4596  }
4597  for (i = 0; i < CLP_CYCLE - 1; i++) {
4598    //obj_[i]=obj_[i+1];
4599    in_[i] = in_[i + 1];
4600    out_[i] = out_[i + 1];
4601    way_[i] = way_[i + 1];
4602  }
4603#endif
4604  int way = 1 - wayIn + 4 * (1 - wayOut);
4605  //obj_[i]=model_->objectiveValue();
4606  in_[CLP_CYCLE - 1] = in;
4607  out_[CLP_CYCLE - 1] = out;
4608  way_[CLP_CYCLE - 1] = static_cast< char >(way);
4609  return matched;
4610}
4611#include "CoinStructuredModel.hpp"
4612// Solve using structure of model and maybe in parallel
4613int ClpSimplex::solve(CoinStructuredModel *model)
4614{
4615  // analyze structure
4616  int numberRowBlocks = model->numberRowBlocks();
4617  int numberColumnBlocks = model->numberColumnBlocks();
4618  int numberElementBlocks = model->numberElementBlocks();
4619  if (numberElementBlocks == 1) {
4620    loadProblem(*model, false);
4621    return dual();
4622  }
4623  // For now just get top level structure
4624  CoinModelBlockInfo *blockInfo = new CoinModelBlockInfo[numberElementBlocks];
4625  for (int i = 0; i < numberElementBlocks; i++) {
4626    CoinStructuredModel *subModel = dynamic_cast< CoinStructuredModel * >(model->block(i));
4627    CoinModel *thisBlock;
4628    if (subModel) {
4629      thisBlock = subModel->coinModelBlock(blockInfo[i]);
4630      model->setCoinModel(thisBlock, i);
4631    } else {
4632      thisBlock = dynamic_cast< CoinModel * >(model->block(i));
4633      assert(thisBlock);
4634      // just fill in info
4635      CoinModelBlockInfo info = CoinModelBlockInfo();
4636      int whatsSet = thisBlock->whatIsSet();
4637      info.matrix = static_cast< char >(((whatsSet & 1) != 0) ? 1 : 0);
4638      info.rhs = static_cast< char >(((whatsSet & 2) != 0) ? 1 : 0);
4639      info.rowName = static_cast< char >(((whatsSet & 4) != 0) ? 1 : 0);
4640      info.integer = static_cast< char >(((whatsSet & 32) != 0) ? 1 : 0);
4641      info.bounds = static_cast< char >(((whatsSet & 8) != 0) ? 1 : 0);
4642      info.columnName = static_cast< char >(((whatsSet & 16) != 0) ? 1 : 0);
4643      // Which block
4644      int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4645      info.rowBlock = iRowBlock;
4646      int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4647      info.columnBlock = iColumnBlock;
4648      blockInfo[i] = info;
4649    }
4650  }
4651  int *rowCounts = new int[numberRowBlocks];
4652  CoinZeroN(rowCounts, numberRowBlocks);
4653  int *columnCounts = new int[numberColumnBlocks + 1];
4654  CoinZeroN(columnCounts, numberColumnBlocks);
4655  int decomposeType = 0;
4656  for (int i = 0; i < numberElementBlocks; i++) {
4657    int iRowBlock = blockInfo[i].rowBlock;
4658    int iColumnBlock = blockInfo[i].columnBlock;
4659    rowCounts[iRowBlock]++;
4660    columnCounts[iColumnBlock]++;
4661  }
4662  if (numberRowBlocks == numberColumnBlocks || numberRowBlocks == numberColumnBlocks + 1) {
4663    // could be Dantzig-Wolfe
4664    int numberG1 = 0;
4665    for (int i = 0; i < numberRowBlocks; i++) {
4666      if (rowCounts[i] > 1)
4667        numberG1++;
4668    }
4669    bool masterColumns = (numberColumnBlocks == numberRowBlocks);
4670    if ((masterColumns && numberElementBlocks == 2 * numberRowBlocks - 1)
4671      || (!masterColumns && numberElementBlocks == 2 * numberRowBlocks)) {
4672      if (numberG1 < 2)
4673        decomposeType = 1;
4674    }
4675  }
4676  if (!decomposeType && (numberRowBlocks == numberColumnBlocks || numberRowBlocks == numberColumnBlocks - 1)) {
4677    // could be Benders
4678    int numberG1 = 0;
4679    for (int i = 0; i < numberColumnBlocks; i++) {
4680      if (columnCounts[i] > 1)
4681        numberG1++;
4682    }
4683    bool masterRows = (numberColumnBlocks == numberRowBlocks);
4684    if ((masterRows && numberElementBlocks == 2 * numberColumnBlocks - 1)
4685      || (!masterRows && numberElementBlocks == 2 * numberColumnBlocks)) {
4686      if (numberG1 < 2)
4687        decomposeType = 2;
4688    }
4689  }
4690  delete[] rowCounts;
4691  delete[] columnCounts;
4692  delete[] blockInfo;
4693  // decide what to do
4694  ClpSolve options;
4695  options.setIndependentOption(2, 100);
4696  switch (decomposeType) {
4697    // No good
4698  case 0:
4699    loadProblem(*model, false);
4700    return dual();
4701    // DW
4702  case 1:
4703    return solveDW(model, options);
4704    // Benders
4705  case 2:
4706    return solveBenders(model, options);
4707  }
4708  return 0; // to stop compiler warning
4709}
4710/* This loads a model from a CoinStructuredModel object - returns number of errors.
4711   If originalOrder then keep to order stored in blocks,
4712   otherwise first column/rows correspond to first block - etc.
4713   If keepSolution true and size is same as current then
4714   keeps current status and solution
4715*/
4716int ClpSimplex::loadProblem(CoinStructuredModel &coinModel,
4717  bool originalOrder,
4718  bool keepSolution)
4719{
4720  unsigned char *status = NULL;
4721  double *psol = NULL;
4722  double *dsol = NULL;
4723  int numberRows = coinModel.numberRows();
4724  int numberColumns = coinModel.numberColumns();
4725  int numberRowBlocks = coinModel.numberRowBlocks();
4726  int numberColumnBlocks = coinModel.numberColumnBlocks();
4727  int numberElementBlocks = coinModel.numberElementBlocks();
4728  if (status_ && numberRows_ && numberRows_ == numberRows && numberColumns_ == numberColumns && keepSolution) {
4729    status = new unsigned char[numberRows_ + numberColumns_];
4730    CoinMemcpyN(status_, numberRows_ + numberColumns_, status);
4731    psol = new double[numberRows_ + numberColumns_];
4732    CoinMemcpyN(columnActivity_, numberColumns_, psol);
4733    CoinMemcpyN(rowActivity_, numberRows_, psol + numberColumns_);
4734    dsol = new double[numberRows_ + numberColumns_];
4735    CoinMemcpyN(reducedCost_, numberColumns_, dsol);
4736    CoinMemcpyN(dual_, numberRows_, dsol + numberColumns_);
4737  }
4738  int returnCode = 0;
4739  double *rowLower = new double[numberRows];
4740  double *rowUpper = new double[numberRows];
4741  double *columnLower = new double[numberColumns];
4742  double *columnUpper = new double[numberColumns];
4743  double *objective = new double[numberColumns];
4744  int *integerType = new int[numberColumns];
4745  CoinBigIndex numberElements = 0;
4746  // Bases for blocks
4747  int *rowBase = new int[numberRowBlocks];
4748  CoinFillN(rowBase, numberRowBlocks, -1);
4749  // And row to put it
4750  int *whichRow = new int[numberRows + numberRowBlocks];
4751  int *columnBase = new int[numberColumnBlocks];
4752  CoinFillN(columnBase, numberColumnBlocks, -1);
4753  // And column to put it
4754  int *whichColumn = new int[numberColumns + numberColumnBlocks];
4755  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4756    CoinModel *block = coinModel.coinBlock(iBlock);
4757    numberElements += block->numberElements();
4758    //and set up elements etc
4759    double *associated = block->associatedArray();
4760    // If strings then do copies
4761    if (block->stringsExist())
4762      returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
4763        objective, integerType, associated);
4764    const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4765    int iRowBlock = info.rowBlock;
4766    int iColumnBlock = info.columnBlock;
4767    if (rowBase[iRowBlock] < 0) {
4768      rowBase[iRowBlock] = block->numberRows();
4769      // Save block number
4770      whichRow[numberRows + iRowBlock] = iBlock;
4771    } else {
4772      assert(rowBase[iRowBlock] == block->numberRows());
4773    }
4774    if (columnBase[iColumnBlock] < 0) {
4775      columnBase[iColumnBlock] = block->numberColumns();
4776      // Save block number
4777      whichColumn[numberColumns + iColumnBlock] = iBlock;
4778    } else {
4779      assert(columnBase[iColumnBlock] == block->numberColumns());
4780    }
4781  }
4782  // Fill arrays with defaults
4783  CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
4784  CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
4785  CoinFillN(columnLower, numberColumns, 0.0);
4786  CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
4787  CoinFillN(objective, numberColumns, 0.0);
4788  CoinFillN(integerType, numberColumns, 0);
4789  int n = 0;
4790  for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4791    int k = rowBase[iBlock];
4792    rowBase[iBlock] = n;
4793    assert(k >= 0);
4794    // block number
4795    int jBlock = whichRow[numberRows + iBlock];
4796    if (originalOrder) {
4797      memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
4798    } else {
4799      CoinIotaN(whichRow + n, k, n);
4800    }
4801    n += k;
4802  }
4803  assert(n == numberRows);
4804  n = 0;
4805  for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4806    int k = columnBase[iBlock];
4807    columnBase[iBlock] = n;
4808    assert(k >= 0);
4809    if (k) {
4810      // block number
4811      int jBlock = whichColumn[numberColumns + iBlock];
4812      if (originalOrder) {
4813        memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
4814          k * sizeof(int));
4815      } else {
4816        CoinIotaN(whichColumn + n, k, n);
4817      }
4818      n += k;
4819    }
4820  }
4821  assert(n == numberColumns);
4822  bool gotIntegers = false;
4823  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4824    CoinModel *block = coinModel.coinBlock(iBlock);
4825    const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4826    int iRowBlock = info.rowBlock;
4827    int iRowBase = rowBase[iRowBlock];
4828    int iColumnBlock = info.columnBlock;
4829    int iColumnBase = columnBase[iColumnBlock];
4830    if (info.rhs) {
4831      int nRows = block->numberRows();
4832      const double *lower = block->rowLowerArray();
4833      const double *upper = block->rowUpperArray();
4834      for (int i = 0; i < nRows; i++) {
4835        int put = whichRow[i + iRowBase];
4836        rowLower[put] = lower[i];
4837        rowUpper[put] = upper[i];
4838      }
4839    }
4840    if (info.bounds) {
4841      int nColumns = block->numberColumns();
4842      const double *lower = block->columnLowerArray();
4843      const double *upper = block->columnUpperArray();
4844      const double *obj = block->objectiveArray();
4845      for (int i = 0; i < nColumns; i++) {
4846        int put = whichColumn[i + iColumnBase];
4847        columnLower[put] = lower[i];
4848        columnUpper[put] = upper[i];
4849        objective[put] = obj[i];
4850      }
4851    }
4852    if (info.integer) {
4853      gotIntegers = true;
4854      int nColumns = block->numberColumns();
4855      const int *type = block->integerTypeArray();
4856      for (int i = 0; i < nColumns; i++) {
4857        int put = whichColumn[i + iColumnBase];
4858        integerType[put] = type[i];
4859      }
4860    }
4861  }
4862  gutsOfLoadModel(numberRows, numberColumns,
4863    columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
4864  delete[] rowLower;
4865  delete[] rowUpper;
4866  delete[] columnLower;
4867  delete[] columnUpper;
4868  delete[] objective;
4869  // Do integers if wanted
4870  if (gotIntegers) {
4871    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4872      if (integerType[iColumn])
4873        setInteger(iColumn);
4874    }
4875  }
4876  delete[] integerType;
4877  setObjectiveOffset(coinModel.objectiveOffset());
4878  // Space for elements
4879  int *row = new int[numberElements];
4880  int *column = new int[numberElements];
4881  double *element = new double[numberElements];
4882  numberElements = 0;
4883  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4884    CoinModel *block = coinModel.coinBlock(iBlock);
4885    const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4886    int iRowBlock = info.rowBlock;
4887    int iRowBase = rowBase[iRowBlock];
4888    int iColumnBlock = info.columnBlock;
4889    int iColumnBase = columnBase[iColumnBlock];
4890    if (info.rowName) {
4891      int numberItems = block->rowNames()->numberItems();
4892      assert(block->numberRows() >= numberItems);
4893      if (numberItems) {
4894        const char *const *rowNames = block->rowNames()->names();
4895        for (int i = 0; i < numberItems; i++) {
4896          int put = whichRow[i + iRowBase];
4897          std::string name = rowNames[i];
4898          setRowName(put, name);
4899        }
4900      }
4901    }
4902    if (info.columnName) {
4903      int numberItems = block->columnNames()->numberItems();
4904      assert(block->numberColumns() >= numberItems);
4905      if (numberItems) {
4906        const char *const *columnNames = block->columnNames()->names();
4907        for (int i = 0; i < numberItems; i++) {
4908          int put = whichColumn[i + iColumnBase];
4909          std::string name = columnNames[i];
4910          setColumnName(put, name);
4911        }
4912      }
4913    }
4914    if (info.matrix) {
4915      CoinPackedMatrix matrix2;
4916      const CoinPackedMatrix *matrix = block->packedMatrix();
4917      if (!matrix) {
4918        double *associated = block->associatedArray();
4919        block->createPackedMatrix(matrix2, associated);
4920        matrix = &matrix2;
4921      }
4922      // get matrix data pointers
4923      const int *row2 = matrix->getIndices();
4924      const CoinBigIndex *columnStart = matrix->getVectorStarts();
4925      const double *elementByColumn = matrix->getElements();
4926      const int *columnLength = matrix->getVectorLengths();
4927      int n = matrix->getNumCols();
4928      assert(matrix->isColOrdered());
4929      for (int iColumn = 0; iColumn < n; iColumn++) {
4930        CoinBigIndex j;
4931        int jColumn = whichColumn[iColumn + iColumnBase];
4932        for (j = columnStart[iColumn];
4933             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4934          row[numberElements] = whichRow[row2[j] + iRowBase];
4935          column[numberElements] = jColumn;
4936          element[numberElements++] = elementByColumn[j];
4937        }
4938      }
4939    }
4940  }
4941  delete[] whichRow;
4942  delete[] whichColumn;
4943  delete[] rowBase;
4944  delete[] columnBase;
4945  CoinPackedMatrix *matrix = new CoinPackedMatrix(true, row, column, element, numberElements);
4946  matrix_ = new ClpPackedMatrix(matrix);
4947  matrix_->setDimensions(numberRows, numberColumns);
4948  delete[] row;
4949  delete[] column;
4950  delete[] element;
4951  createStatus();
4952  if (status) {
4953    // copy back
4954    CoinMemcpyN(status, numberRows_ + numberColumns_, status_);
4955    CoinMemcpyN(psol, numberColumns_, columnActivity_);
4956    CoinMemcpyN(psol + numberColumns_, numberRows_, rowActivity_);
4957    CoinMemcpyN(dsol, numberColumns_, reducedCost_);
4958    CoinMemcpyN(dsol + numberColumns_, numberRows_, dual_);
4959    delete[] status;
4960    delete[] psol;
4961    delete[] dsol;
4962  }
4963  optimizationDirection_ = coinModel.optimizationDirection();
4964  return returnCode;
4965}
4966/*  If input negative scales objective so maximum <= -value
4967    and returns scale factor used.  If positive unscales and also
4968    redoes dual stuff
4969*/
4970double
4971ClpSimplex::scaleObjective(double value)
4972{
4973  double *obj = objective();
4974  double largest = 0.0;
4975  if (value < 0.0) {
4976    value = -value;
4977    for (int i = 0; i < numberColumns_; i++) {
4978      largest = CoinMax(largest, fabs(obj[i]));
4979    }
4980    if (largest > value) {
4981      double scaleFactor = value / largest;
4982      for (int i = 0; i < numberColumns_; i++) {
4983        obj[i] *= scaleFactor;
4984        reducedCost_[i] *= scaleFactor;
4985      }
4986      for (int i = 0; i < numberRows_; i++) {
4987        dual_[i] *= scaleFactor;
4988      }
4989      largest /= value;
4990    } else {
4991      // no need
4992      largest = 1.0;
4993    }
4994  } else {
4995    // at end
4996    if (value != 1.0) {
4997      for (int i = 0; i < numberColumns_; i++) {
4998        obj[i] *= value;
4999        reducedCost_[i] *= value;
5000      }
5001      for (int i = 0; i < numberRows_; i++) {
5002        dual_[i] *= value;
5003      }
5004      computeObjectiveValue();
5005    }
5006  }
5007  return largest;
5008}
5009#if defined(ABC_INHERIT) || defined(CBC_THREAD) || defined(THREADS_IN_ANALYZE)
5010void *clp_parallelManager(void *stuff)
5011{
5012  CoinPthreadStuff *driver = reinterpret_cast< CoinPthreadStuff * >(stuff);
5013  int whichThread = driver->whichThread();
5014  CoinThreadInfo *threadInfo = driver->threadInfoPointer(whichThread);
5015  threadInfo->status = -1;
5016  int *which = threadInfo->stuff;
5017#ifdef PTHREAD_BARRIER_SERIAL_THREAD
5018  pthread_barrier_wait(driver->barrierPointer());
5019#endif
5020#if 0
5021  int status=-1;
5022  while (status!=100)
5023    status=timedWait(driver,1000,2);
5024  pthread_cond_signal(driver->conditionPointer(1));
5025  pthread_mutex_unlock(driver->mutexPointer(1,whichThread));
5026#endif
5027  // so now mutex_ is locked
5028  int whichLocked = 0;
5029  while (true) {
5030    pthread_mutex_t *mutexPointer = driver->mutexPointer(whichLocked, whichThread);
5031    // wait
5032    //printf("Child waiting for %d - status %d %d %d\n",
5033    //     whichLocked,lockedX[0],lockedX[1],lockedX[2]);
5034#ifdef DETAIL_THREAD
5035    printf("thread %d about to lock mutex %d\n", whichThread, whichLocked);
5036#endif
5037    pthread_mutex_lock(mutexPointer);
5038    whichLocked++;
5039    if (whichLocked == 3)
5040      whichLocked = 0;
5041    int unLock = whichLocked + 1;
5042    if (unLock == 3)
5043      unLock = 0;
5044    //printf("child pointer %x status %d\n",threadInfo,threadInfo->status);
5045    assert(threadInfo->status >= 0);
5046    if (threadInfo->status == 1000)
5047      pthread_exit(NULL);
5048    int type = threadInfo->status;
5049    int &returnCode = which[0];
5050    int iPass = which[1];
5051    ClpSimplex *clpSimplex = reinterpret_cast< ClpSimplex * >(threadInfo->extraInfo);
5052    //CoinIndexedVector * array;
5053    //double dummy;
5054    switch (type) {
5055      // dummy
5056    case 0:
5057      break;
5058    case 1:
5059      if (!clpSimplex->problemStatus() || !iPass)
5060        returnCode = clpSimplex->dual();
5061      else
5062        returnCode = clpSimplex->primal();
5063      break;
5064    case 100:
5065      // initialization
5066      break;
5067    }
5068    threadInfo->status = -1;
5069#ifdef DETAIL_THREAD
5070    printf("thread %d about to unlock mutex %d\n", whichThread, unLock);
5071#endif
5072    pthread_mutex_unlock(driver->mutexPointer(unLock, whichThread));
5073  }
5074}
5075#endif
5076// Solve using Dantzig-Wolfe decomposition and maybe in parallel
5077int ClpSimplex::solveDW(CoinStructuredModel *model, ClpSolve &options)
5078{
5079  double time1 = CoinCpuTime();
5080  int numberColumns = model->numberColumns();
5081  int numberRowBlocks = model->numberRowBlocks();
5082  int numberColumnBlocks = model->numberColumnBlocks();
5083  int numberElementBlocks = model->numberElementBlocks();
5084  // We already have top level structure
5085  CoinModelBlockInfo *blockInfo = new CoinModelBlockInfo[numberElementBlocks];
5086  for (int i = 0; i < numberElementBlocks; i++) {
5087    CoinModel *thisBlock = model->coinBlock(i);
5088    assert(thisBlock);
5089    // just fill in info
5090    CoinModelBlockInfo info = CoinModelBlockInfo();
5091    int whatsSet = thisBlock->whatIsSet();
5092    info.matrix = static_cast< char >(((whatsSet & 1) != 0) ? 1 : 0);
5093    info.rhs = static_cast< char >(((whatsSet & 2) != 0) ? 1 : 0);
5094    info.rowName = static_cast< char >(((whatsSet & 4) != 0) ? 1 : 0);
5095    info.integer = static_cast< char >(((whatsSet & 32) != 0) ? 1 : 0);
5096    info.bounds = static_cast< char >(((whatsSet & 8) != 0) ? 1 : 0);
5097    info.columnName = static_cast< char >(((whatsSet & 16) != 0) ? 1 : 0);
5098    // Which block
5099    int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
5100    info.rowBlock = iRowBlock;
5101    int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
5102    info.columnBlock = iColumnBlock;
5103    blockInfo[i] = info;
5104  }
5105  // make up problems
5106  int numberBlocks = numberRowBlocks - 1;
5107  // Find master rows and columns
5108  int *rowCounts = new int[numberRowBlocks];
5109  CoinZeroN(rowCounts, numberRowBlocks);
5110  int *columnCounts = new int[numberColumnBlocks + 1];
5111  CoinZeroN(columnCounts, numberColumnBlocks);
5112  int iBlock;
5113  for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5114    int iRowBlock = blockInfo[iBlock].rowBlock;
5115    rowCounts[iRowBlock]++;
5116    int iColumnBlock = blockInfo[iBlock].columnBlock;
5117    columnCounts[iColumnBlock]++;
5118  }
5119  int *whichBlock = new int[numberElementBlocks];
5120  int masterRowBlock = -1;
5121  for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5122    if (rowCounts[iBlock] > 1) {
5123      if (masterRowBlock == -1) {
5124        masterRowBlock = iBlock;
5125      } else {
5126        // Can't decode
5127        masterRowBlock = -2;
5128        break;
5129      }
5130    }
5131  }
5132  int masterColumnBlock = -1;
5133  int kBlock = 0;
5134  for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5135    int count = columnCounts[iBlock];
5136    columnCounts[iBlock] = kBlock;
5137    kBlock += count;
5138  }
5139  for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5140    int iColumnBlock = blockInfo[iBlock].columnBlock;
5141    whichBlock[columnCounts[iColumnBlock]] = iBlock;
5142    columnCounts[iColumnBlock]++;
5143  }
5144  for (iBlock = numberColumnBlocks - 1; iBlock >= 0; iBlock--)
5145    columnCounts[iBlock + 1] = columnCounts[iBlock];
5146  columnCounts[0] = 0;
5147  for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5148    int count = columnCounts[iBlock + 1] - columnCounts[iBlock];
5149    if (count == 1) {
5150      int kBlock = whichBlock[columnCounts[iBlock]];
5151      int iRowBlock = blockInfo[kBlock].rowBlock;
5152      if (iRowBlock == masterRowBlock) {
5153        if (masterColumnBlock == -1) {
5154          masterColumnBlock = iBlock;
5155        } else {
5156          // Can't decode
5157          masterColumnBlock = -2;
5158          break;
5159        }
5160      }
5161    }
5162  }
5163  if (masterRowBlock < 0 || masterColumnBlock == -2) {
5164    // What now
5165    abort();
5166  }
5167  delete[] rowCounts;
5168  // create all data
5169  const CoinPackedMatrix **top = new const CoinPackedMatrix *[numberColumnBlocks];
5170  ClpSimplex *sub = new ClpSimplex[numberBlocks];
5171  ClpSimplex master;
5172  // Set offset
5173  master.setObjectiveOffset(model->objectiveOffset());
5174  bool reducePrint = logLevel() == 7;
5175  if (reducePrint) {
5176    // special
5177    setLogLevel(1);
5178    master.setLogLevel(1);
5179  }
5180  kBlock = 0;
5181  int masterBlock = -1;
5182  for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5183    top[kBlock] = NULL;
5184    int start = columnCounts[iBlock];
5185    int end = columnCounts[iBlock + 1];
5186    assert(end - start <= 2);
5187    for (int j = start; j < end; j++) {
5188      int jBlock = whichBlock[j];
5189      int iRowBlock = blockInfo[jBlock].rowBlock;
5190      int iColumnBlock = blockInfo[jBlock].columnBlock;
5191      assert(iColumnBlock == iBlock);
5192      if (iColumnBlock != masterColumnBlock && iRowBlock == masterRowBlock) {
5193        // top matrix
5194        top[kBlock] = model->coinBlock(jBlock)->packedMatrix();
5195      } else {
5196        const CoinPackedMatrix *matrix
5197          = model->coinBlock(jBlock)->packedMatrix();
5198        // Get pointers to arrays
5199        const double *rowLower;
5200        const double *rowUpper;
5201        const double *columnLower;
5202        const double *columnUpper;
5203        const double *objective;
5204        model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
5205          columnLower, columnUpper, objective);
5206        if (iColumnBlock != masterColumnBlock) {
5207          // diagonal block
5208          sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
5209            objective, rowLower, rowUpper);
5210          if (true) {
5211            double *lower = sub[kBlock].columnLower();
5212            double *upper = sub[kBlock].columnUpper();
5213            int n = sub[kBlock].numberColumns();
5214            for (int i = 0; i < n; i++) {
5215              lower[i] = CoinMax(-1.0e8, lower[i]);
5216              upper[i] = CoinMin(1.0e8, upper[i]);
5217            }
5218          }
5219          if (optimizationDirection_ < 0.0) {
5220            double *obj = sub[kBlock].objective();
5221            int n = sub[kBlock].numberColumns();
5222            for (int i = 0; i < n; i++)
5223              obj[i] = -obj[i];
5224          }
5225          if (this->factorizationFrequency() == 200) {
5226            // User did not touch preset
5227            sub[kBlock].defaultFactorizationFrequency();
5228          } else {
5229            // make sure model has correct value
5230            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
5231          }
5232          sub[kBlock].setPerturbation(50);
5233          // Set columnCounts to be diagonal block index for cleanup
5234          columnCounts[kBlock] = jBlock;
5235        } else {
5236          // master
5237          masterBlock = jBlock;
5238          master.loadProblem(*matrix, columnLower, columnUpper,
5239            objective, rowLower, rowUpper);
5240          if (optimizationDirection_ < 0.0) {
5241            double *obj = master.objective();
5242            int n = master.numberColumns();
5243            for (int i = 0; i < n; i++)
5244              obj[i] = -obj[i];
5245          }
5246        }
5247      }
5248    }
5249    if (iBlock != masterColumnBlock)
5250      kBlock++;
5251  }
5252  delete[] whichBlock;
5253  delete[] blockInfo;
5254  // For now master must have been defined (does not have to have columns)
5255  assert(master.numberRows());
5256  assert(masterBlock >= 0);
5257  int numberMasterRows = master.numberRows();
5258  // Overkill in terms of space
5259  int spaceNeeded = CoinMax(numberBlocks * (numberMasterRows + 1),
5260    2 * numberMasterRows);
5261  int *rowAdd = new int[spaceNeeded];
5262  double *elementAdd = new double[spaceNeeded];
5263  spaceNeeded = numberBlocks;
5264  CoinBigIndex *columnAdd = new CoinBigIndex[spaceNeeded + 1];
5265  double *objective = new double[spaceNeeded];
5266  // Add in costed slacks
5267  int firstArtificial = master.numberColumns();
5268  int lastArtificial = firstArtificial;
5269  if (true) {
5270    const double *lower = master.rowLower();
5271    const double *upper = master.rowUpper();
5272    int kCol = 0;
5273    for (int iRow = 0; iRow < numberMasterRows; iRow++) {
5274      if (lower[iRow] > -1.0e10) {
5275        rowAdd[kCol] = iRow;
5276        elementAdd[kCol++] = 1.0;
5277      }
5278      if (upper[iRow] < 1.0e10) {
5279        rowAdd[kCol] = iRow;
5280        elementAdd[kCol++] = -1.0;
5281      }
5282    }
5283    if (kCol > spaceNeeded) {
5284      spaceNeeded = kCol;
5285      delete[] columnAdd;
5286      delete[] objective;
5287      columnAdd = new CoinBigIndex[spaceNeeded + 1];
5288      objective = new double[spaceNeeded];
5289    }
5290    for (int i = 0; i < kCol; i++) {
5291      columnAdd[i] = i;
5292      objective[i] = 1.0e13;
5293    }
5294    columnAdd[kCol] = kCol;
5295    master.addColumns(kCol, NULL, NULL, objective,
5296      columnAdd, rowAdd, elementAdd);
5297    lastArtificial = master.numberColumns();
5298  }
5299  int maxPass = options.independentOption(2);
5300  if (maxPass < 2)
5301    maxPass = 100;
5302  int iPass;
5303  double lastObjective = 1.0e31;
5304  // Create convexity rows for proposals
5305  int numberMasterColumns = master.numberColumns();
5306  master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
5307  if (this->factorizationFrequency() == 200) {
5308    // User did not touch preset
5309    master.defaultFactorizationFrequency();
5310  } else {
5311    // make sure model has correct value
5312    master.setFactorizationFrequency(this->factorizationFrequency());
5313  }
5314  master.setPerturbation(50);
5315  // Arrays to say which block and when created
5316  int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
5317  whichBlock = new int[maximumColumns];
5318  int *when = new int[maximumColumns];
5319  int numberColumnsGenerated = numberBlocks;
5320  // fill in rhs and add in artificials
5321  {
5322    double *rowLower = master.rowLower();
5323    double *rowUpper = master.rowUpper();
5324    int iBlock;
5325    columnAdd[0] = 0;
5326    for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
5327      int iRow = iBlock + numberMasterRows;
5328      ;
5329      rowLower[iRow] = 1.0;
5330      rowUpper[iRow] = 1.0;
5331      rowAdd[iBlock] = iRow;
5332      elementAdd[iBlock] = 1.0;
5333      objective[iBlock] = 1.0e13;
5334      columnAdd[iBlock + 1] = iBlock + 1;
5335      when[iBlock] = -1;
5336      whichBlock[iBlock] = iBlock;
5337    }
5338    master.addColumns(numberBlocks, NULL, NULL, objective,
5339      columnAdd, rowAdd, elementAdd);
5340  }
5341  char generalPrint[200];
5342  // and resize matrix to double check clp will be happy
5343  //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
5344  //                     numberMasterColumns+numberBlocks);
5345  sprintf(generalPrint, "Time to decompose %.2f seconds", CoinCpuTime() - time1);
5346  handler_->message(CLP_GENERAL, messages_)
5347    << generalPrint
5348    << CoinMessageEol;
5349#ifdef ABC_INHERIT
5350  //AbcSimplex abcMaster;
5351  //if (!this->abcState())
5352  //setAbcState(1);
5353  int numberCpu = CoinMin((this->abcState() & 15), 4);
5354  CoinPthreadStuff threadInfo(numberCpu, clp_parallelManager);
5355  master.setAbcState(this->abcState());
5356  //AbcSimplex * tempMaster=master.dealWithAbc(2,10,true);
5357  //abcMaster=*tempMaster;
5358  //delete tempMaster;
5359  //abcMaster.startThreads(numberCpu);
5360  //#define master abcMaster
5361#endif
5362  for (iPass = 0; iPass < maxPass; iPass++) {
5363    sprintf(generalPrint, "Start of pass %d", iPass);
5364    handler_->message(CLP_GENERAL, messages_)
5365      << generalPrint
5366      << CoinMessageEol;
5367    // Solve master - may be infeasible
5368    //master.scaling(0);
5369    if (0) {
5370      master.writeMps("yy.mps");
5371    }
5372    // Correct artificials
5373    double sumArtificials = 0.0;
5374    if (iPass) {
5375      double *upper = master.columnUpper();
5376      double *solution = master.primalColumnSolution();
5377      double *obj = master.objective();
5378      sumArtificials = 0.0;
5379      for (int i = firstArtificial; i < lastArtificial; i++) {
5380        sumArtificials += solution[i];
5381        //assert (solution[i]>-1.0e-2);
5382        if (solution[i] < 1.0e-6) {
5383#if 0
5384                         // Could take out
5385                         obj[i] = 0.0;
5386                         upper[i] = 0.0;
5387#else
5388          obj[i] = 1.0e7;
5389          upper[i] = 1.0e-1;
5390#endif
5391          solution[i] = 0.0;
5392          master.setColumnStatus(i, isFixed);
5393        } else {
5394          upper[i] = solution[i] + 1.0e-5 * (1.0 + solution[i]);
5395        }
5396      }
5397      sprintf(generalPrint, "Sum of artificials before solve is %g", sumArtificials);
5398      handler_->message(CLP_GENERAL, messages_)
5399        << generalPrint
5400        << CoinMessageEol;
5401    }
5402    // scale objective to be reasonable
5403    double scaleFactor = master.scaleObjective(-1.0e9);
5404    {
5405      double *dual = master.dualRowSolution();
5406      int n = master.numberRows();
5407      memset(dual, 0, n * sizeof(double));
5408      double *solution = master.primalColumnSolution();
5409      master.clpMatrix()->times(1.0, solution, dual);
5410      double sum = 0.0;
5411      double *lower = master.rowLower();
5412      double *upper = master.rowUpper();
5413      for (int iRow = 0; iRow < n; iRow++) {
5414        double value = dual[iRow];
5415        if (value > upper[iRow])
5416          sum += value - upper[iRow];
5417        else if (value < lower[iRow])
5418          sum -= value - lower[iRow];
5419      }
5420      printf("** suminf %g\n", sum);
5421      lower = master.columnLower();
5422      upper = master.columnUpper();
5423      n = master.numberColumns();
5424      for (int iColumn = 0; iColumn < n; iColumn++) {
5425        double value = solution[iColumn];
5426        if (value > upper[iColumn] + 1.0e-5)
5427          sum += value - upper[iColumn];
5428        else if (value < lower[iColumn] - 1.0e-5)
5429          sum -= value - lower[iColumn];
5430      }
5431      printf("** suminf %g\n", sum);
5432    }
5433#ifdef ABC_INHERIT
5434    master.dealWithAbc(1, 0, true);
5435#else
5436    master.primal();
5437#endif
5438    //master.primal(1);
5439    // Correct artificials
5440    sumArtificials = 0.0;
5441    {
5442      double *solution = master.primalColumnSolution();
5443      for (int i = firstArtificial; i < lastArtificial; i++) {
5444        sumArtificials += solution[i];
5445      }
5446      printf("** Sum of artificials after solve is %g\n", sumArtificials);
5447    }
5448    master.scaleObjective(scaleFactor);
5449    int problemStatus = master.status(); // do here as can change (delcols)
5450    if (problemStatus == 2 && master.numberColumns()) {
5451#ifdef ABC_INHERIT
5452      master.dealWithAbc(1, 1, true);
5453#else
5454      master