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

Last change on this file since 2385 was 2385, checked in by unxusr, 3 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 280.1 KB
Line 
1/* $Id: ClpSolve.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
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    if (interrupt)
3120      currentModel2 = &barrier;
3121    if (barrier.numberRows() + barrier.numberColumns() > 10000)
3122      barrier.setMaximumBarrierIterations(1000);
3123    int barrierOptions = options.getSpecialOption(4);
3124    int aggressiveGamma = 0;
3125    bool presolveInCrossover = false;
3126    bool scale = false;
3127    bool doKKT = false;
3128    bool forceFixing = false;
3129    int speed = 0;
3130    if (barrierOptions & 16) {
3131      barrierOptions &= ~16;
3132      doKKT = true;
3133    }
3134    if (barrierOptions & (32 + 64 + 128)) {
3135      aggressiveGamma = (barrierOptions & (32 + 64 + 128)) >> 5;
3136      barrierOptions &= ~(32 + 64 + 128);
3137    }
3138    if (barrierOptions & 256) {
3139      barrierOptions &= ~256;
3140      presolveInCrossover = true;
3141    }
3142    if (barrierOptions & 512) {
3143      barrierOptions &= ~512;
3144      forceFixing = true;
3145    }
3146    if (barrierOptions & 1024) {
3147      barrierOptions &= ~1024;
3148      barrier.setProjectionTolerance(1.0e-9);
3149    }
3150    if (barrierOptions & (2048 | 4096)) {
3151      speed = (barrierOptions & (2048 | 4096)) >> 11;
3152      barrierOptions &= ~(2048 | 4096);
3153    }
3154    if (barrierOptions & 8) {
3155      barrierOptions &= ~8;
3156      scale = true;
3157    }
3158    // If quadratic force KKT
3159    if (quadraticObj) {
3160      doKKT = true;
3161    }
3162    switch (barrierOptions) {
3163    case 0:
3164    default:
3165      if (!doKKT) {
3166        ClpCholeskyBase *cholesky = new ClpCholeskyBase(options.getExtraInfo(1));
3167        cholesky->setIntegerParameter(0, speed);
3168        barrier.setCholesky(cholesky);
3169      } else {
3170        ClpCholeskyBase *cholesky = new ClpCholeskyBase();
3171        cholesky->setKKT(true);
3172        barrier.setCholesky(cholesky);
3173      }
3174      break;
3175    case 1:
3176      if (!doKKT) {
3177        ClpCholeskyDense *cholesky = new ClpCholeskyDense();
3178        barrier.setCholesky(cholesky);
3179      } else {
3180        ClpCholeskyDense *cholesky = new ClpCholeskyDense();
3181        cholesky->setKKT(true);
3182        barrier.setCholesky(cholesky);
3183      }
3184      break;
3185#ifdef COIN_HAS_WSMP
3186    case 2: {
3187      if (!doKKT) {
3188        ClpCholeskyWssmp *cholesky = new ClpCholeskyWssmp(CoinMax(100, model2->numberRows() / 10));
3189        barrier.setCholesky(cholesky);
3190      } else {
3191        ClpCholeskyWssmpKKT *cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3192        barrier.setCholesky(cholesky);
3193      }
3194    } break;
3195    case 3:
3196      if (!doKKT) {
3197        ClpCholeskyWssmp *cholesky = new ClpCholeskyWssmp();
3198        barrier.setCholesky(cholesky);
3199      } else {
3200        ClpCholeskyWssmpKKT *cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3201        barrier.setCholesky(cholesky);
3202      }
3203      break;
3204#endif
3205#ifdef UFL_BARRIER
3206    case 4:
3207      if (!doKKT) {
3208        ClpCholeskyUfl *cholesky = new ClpCholeskyUfl(options.getExtraInfo(1));
3209        barrier.setCholesky(cholesky);
3210      } else {
3211        ClpCholeskyUfl *cholesky = new ClpCholeskyUfl();
3212        cholesky->setKKT(true);
3213        barrier.setCholesky(cholesky);
3214      }
3215      break;
3216#endif
3217#ifdef TAUCS_BARRIER
3218    case 5: {
3219      ClpCholeskyTaucs *cholesky = new ClpCholeskyTaucs();
3220      barrier.setCholesky(cholesky);
3221      assert(!doKKT);
3222    } break;
3223#endif
3224#ifdef COIN_HAS_MUMPS
3225    case 6: {
3226      if (!doKKT) {
3227        int logLevel = 0;
3228        if (barrier.logLevel() > 3) {
3229          logLevel = (barrier.logLevel() > 4) ? 2 : 1;
3230        }
3231        ClpCholeskyMumps *cholesky = new ClpCholeskyMumps(-1, logLevel);
3232        barrier.setCholesky(cholesky);
3233      } else {
3234        ClpCholeskyMumps *cholesky = new ClpCholeskyMumps();
3235        cholesky->setKKT(true);
3236        barrier.setCholesky(cholesky);
3237      }
3238    } break;
3239#endif
3240#ifdef PARDISO_BARRIER
3241    case 7: {
3242      ClpCholeskyPardiso *cholesky = new ClpCholeskyPardiso();
3243      barrier.setCholesky(cholesky);
3244      assert(!doKKT);
3245    } break;
3246#endif
3247    }
3248    int numberRows = model2->numberRows();
3249    int numberColumns = model2->numberColumns();
3250    int saveMaxIts = model2->maximumIterations();
3251    if (saveMaxIts < 1000) {
3252      barrier.setMaximumBarrierIterations(saveMaxIts);
3253      model2->setMaximumIterations(10000000);
3254    }
3255#ifndef SAVEIT
3256    //barrier.setDiagonalPerturbation(1.0e-25);
3257    if (aggressiveGamma) {
3258      switch (aggressiveGamma) {
3259      case 1:
3260        barrier.setGamma(1.0e-5);
3261        barrier.setDelta(1.0e-5);
3262        break;
3263      case 2:
3264        barrier.setGamma(1.0e-7);
3265        break;
3266      case 3:
3267        barrier.setDelta(1.0e-5);
3268        break;
3269      case 4:
3270        barrier.setGamma(1.0e-3);
3271        barrier.setDelta(1.0e-3);
3272        break;
3273      case 5:
3274        barrier.setGamma(1.0e-3);
3275        break;
3276      case 6:
3277        barrier.setDelta(1.0e-3);
3278        break;
3279      }
3280    }
3281    if (scale)
3282      barrier.scaling(1);
3283    else
3284      barrier.scaling(0);
3285    barrier.primalDual();
3286#elif SAVEIT == 1
3287    barrier.primalDual();
3288#else
3289    model2->restoreModel("xx.save");
3290    // move solutions
3291    CoinMemcpyN(model2->primalRowSolution(),
3292      numberRows, barrier.primalRowSolution());
3293    CoinMemcpyN(model2->dualRowSolution(),
3294      numberRows, barrier.dualRowSolution());
3295    CoinMemcpyN(model2->primalColumnSolution(),
3296      numberColumns, barrier.primalColumnSolution());
3297    CoinMemcpyN(model2->dualColumnSolution(),
3298      numberColumns, barrier.dualColumnSolution());
3299#endif
3300    time2 = CoinCpuTime();
3301    timeCore = time2 - timeX;
3302    handler_->message(CLP_INTERVAL_TIMING, messages_)
3303      << "Barrier" << timeCore << time2 - time1
3304      << CoinMessageEol;
3305    timeX = time2;
3306    int maxIts = barrier.maximumBarrierIterations();
3307    int barrierStatus = barrier.status();
3308    double gap = barrier.complementarityGap();
3309    // get which variables are fixed
3310    double *saveLower = NULL;
3311    double *saveUpper = NULL;
3312    ClpPresolve pinfo2;
3313    ClpSimplex *saveModel2 = NULL;
3314    bool extraPresolve = false;
3315    int numberFixed = barrier.numberFixed();
3316    if (numberFixed) {
3317      int numberRows = barrier.numberRows();
3318      int numberColumns = barrier.numberColumns();
3319      int numberTotal = numberRows + numberColumns;
3320      saveLower = new double[numberTotal];
3321      saveUpper = new double[numberTotal];
3322      CoinMemcpyN(barrier.columnLower(), numberColumns, saveLower);
3323      CoinMemcpyN(barrier.rowLower(), numberRows, saveLower + numberColumns);
3324      CoinMemcpyN(barrier.columnUpper(), numberColumns, saveUpper);
3325      CoinMemcpyN(barrier.rowUpper(), numberRows, saveUpper + numberColumns);
3326    }
3327    if (((numberFixed * 20 > barrier.numberRows() && numberFixed > 5000) || forceFixing) && presolveInCrossover) {
3328      // may as well do presolve
3329      if (!forceFixing) {
3330        barrier.fixFixed();
3331      } else {
3332        // Fix
3333        int n = barrier.numberColumns();
3334        double *lower = barrier.columnLower();
3335        double *upper = barrier.columnUpper();
3336        double *solution = barrier.primalColumnSolution();
3337        int nFix = 0;
3338        for (int i = 0; i < n; i++) {
3339          if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3340            double value = solution[i];
3341            if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3342              solution[i] = lower[i];
3343              upper[i] = lower[i];
3344              nFix++;
3345            } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3346              solution[i] = upper[i];
3347              lower[i] = upper[i];
3348              nFix++;
3349            }
3350          }
3351        }
3352#ifdef CLP_INVESTIGATE
3353        printf("%d columns fixed\n", nFix);
3354#endif
3355        int nr = barrier.numberRows();
3356        lower = barrier.rowLower();
3357        upper = barrier.rowUpper();
3358        solution = barrier.primalRowSolution();
3359        nFix = 0;
3360        for (int i = 0; i < nr; i++) {
3361          if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3362            double value = solution[i];
3363            if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3364              solution[i] = lower[i];
3365              upper[i] = lower[i];
3366              nFix++;
3367            } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3368              solution[i] = upper[i];
3369              lower[i] = upper[i];
3370              nFix++;
3371            }
3372          }
3373        }
3374#ifdef CLP_INVESTIGATE
3375        printf("%d row slacks fixed\n", nFix);
3376#endif
3377      }
3378      saveModel2 = model2;
3379      extraPresolve = true;
3380    } else if (numberFixed) {
3381      // Set fixed to bounds (may have restored earlier solution)
3382      if (!forceFixing) {
3383        barrier.fixFixed(false);
3384      } else {
3385        // Fix
3386        int n = barrier.numberColumns();
3387        double *lower = barrier.columnLower();
3388        double *upper = barrier.columnUpper();
3389        double *solution = barrier.primalColumnSolution();
3390        int nFix = 0;
3391        for (int i = 0; i < n; i++) {
3392          if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3393            double value = solution[i];
3394            if (value < lower[i] + 1.0e-8 && value - lower[i] < upper[i] - value) {
3395              solution[i] = lower[i];
3396              upper[i] = lower[i];
3397              nFix++;
3398            } else if (value > upper[i] - 1.0e-8 && value - lower[i] > upper[i] - value) {
3399              solution[i] = upper[i];
3400              lower[i] = upper[i];
3401              nFix++;
3402            } else {
3403              //printf("fixcol %d %g <= %g <= %g\n",
3404              //     i,lower[i],solution[i],upper[i]);
3405            }
3406          }
3407        }
3408#ifdef CLP_INVESTIGATE
3409        printf("%d columns fixed\n", nFix);
3410#endif
3411        int nr = barrier.numberRows();
3412        lower = barrier.rowLower();
3413        upper = barrier.rowUpper();
3414        solution = barrier.primalRowSolution();
3415        nFix = 0;
3416        for (int i = 0; i < nr; i++) {
3417          if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3418            double value = solution[i];
3419            if (value < lower[i] + 1.0e-5 && value - lower[i] < upper[i] - value) {
3420              solution[i] = lower[i];
3421              upper[i] = lower[i];
3422              nFix++;
3423            } else if (value > upper[i] - 1.0e-5 && value - lower[i] > upper[i] - value) {
3424              solution[i] = upper[i];
3425              lower[i] = upper[i];
3426              nFix++;
3427            } else {
3428              //printf("fixrow %d %g <= %g <= %g\n",
3429              //     i,lower[i],solution[i],upper[i]);
3430            }
3431          }
3432        }
3433#ifdef CLP_INVESTIGATE
3434        printf("%d row slacks fixed\n", nFix);
3435#endif
3436      }
3437    }
3438#ifdef BORROW
3439    int saveNumberIterations = barrier.numberIterations();
3440    barrier.returnModel(*model2);
3441    double *rowPrimal = new double[numberRows];
3442    double *columnPrimal = new double[numberColumns];
3443    double *rowDual = new double[numberRows];
3444    double *columnDual = new double[numberColumns];
3445    // move solutions other way
3446    CoinMemcpyN(model2->primalRowSolution(),
3447      numberRows, rowPrimal);
3448    CoinMemcpyN(model2->dualRowSolution(),
3449      numberRows, rowDual);
3450    CoinMemcpyN(model2->primalColumnSolution(),
3451      numberColumns, columnPrimal);
3452    CoinMemcpyN(model2->dualColumnSolution(),
3453      numberColumns, columnDual);
3454#else
3455    double *rowPrimal = barrier.primalRowSolution();
3456    double *columnPrimal = barrier.primalColumnSolution();
3457    double *rowDual = barrier.dualRowSolution();
3458    double *columnDual = barrier.dualColumnSolution();
3459    // move solutions
3460    CoinMemcpyN(rowPrimal,
3461      numberRows, model2->primalRowSolution());
3462    CoinMemcpyN(rowDual,
3463      numberRows, model2->dualRowSolution());
3464    CoinMemcpyN(columnPrimal,
3465      numberColumns, model2->primalColumnSolution());
3466    CoinMemcpyN(columnDual,
3467      numberColumns, model2->dualColumnSolution());
3468#endif
3469    if (saveModel2) {
3470      // do presolve
3471      model2 = pinfo2.presolvedModel(*model2, dblParam_[ClpPresolveTolerance],
3472        false, 5, true);
3473      if (!model2) {
3474        model2 = saveModel2;
3475        saveModel2 = NULL;
3476        int numberRows = model2->numberRows();
3477        int numberColumns = model2->numberColumns();
3478        CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3479        CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3480        delete[] saveLower;
3481        CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3482        CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3483        delete[] saveUpper;
3484        saveLower = NULL;
3485        saveUpper = NULL;
3486      }
3487    }
3488    if (method == ClpSolve::useBarrier || barrierStatus < 0) {
3489      if (maxIts && barrierStatus < 4 && !quadraticObj) {
3490        //printf("***** crossover - needs more thought on difficult models\n");
3491#if SAVEIT == 1
3492        model2->ClpSimplex::saveModel("xx.save");
3493#endif
3494        // make sure no status left
3495        model2->createStatus();
3496        // solve
3497        if (!forceFixing)
3498          model2->setPerturbation(100);
3499        if (model2->factorizationFrequency() == 200) {
3500          // User did not touch preset
3501          model2->defaultFactorizationFrequency();
3502        }
3503#if 1 //ndef ABC_INHERIT //#if 1 \
3504  // throw some into basis
3505        if (!forceFixing) {
3506          int numberRows = model2->numberRows();
3507          int numberColumns = model2->numberColumns();
3508          double *dsort = new double[numberColumns];
3509          int *sort = new int[numberColumns];
3510          int n = 0;
3511          const double *columnLower = model2->columnLower();
3512          const double *columnUpper = model2->columnUpper();
3513          double *primalSolution = model2->primalColumnSolution();
3514          const double *dualSolution = model2->dualColumnSolution();
3515          double tolerance = 10.0 * primalTolerance_;
3516          int i;
3517          for (i = 0; i < numberRows; i++)
3518            model2->setRowStatus(i, superBasic);
3519          for (i = 0; i < numberColumns; i++) {
3520            double distance = CoinMin(columnUpper[i] - primalSolution[i],
3521              primalSolution[i] - columnLower[i]);
3522            if (distance > tolerance) {
3523              if (fabs(dualSolution[i]) < 1.0e-5)
3524                distance *= 100.0;
3525              dsort[n] = -distance;
3526              sort[n++] = i;
3527              model2->setStatus(i, superBasic);
3528            } else if (distance > primalTolerance_) {
3529              model2->setStatus(i, superBasic);
3530            } else if (primalSolution[i] <= columnLower[i] + primalTolerance_) {
3531              model2->setStatus(i, atLowerBound);
3532              primalSolution[i] = columnLower[i];
3533            } else {
3534              model2->setStatus(i, atUpperBound);
3535              primalSolution[i] = columnUpper[i];
3536            }
3537          }
3538          CoinSort_2(dsort, dsort + n, sort);
3539          n = CoinMin(numberRows, n);
3540          for (i = 0; i < n; i++) {
3541            int iColumn = sort[i];
3542            model2->setStatus(iColumn, basic);
3543          }
3544          delete[] sort;
3545          delete[] dsort;
3546          // model2->allSlackBasis();
3547          if (gap < 1.0e-3 * static_cast< double >(numberRows + numberColumns)) {
3548            if (saveUpper) {
3549              int numberRows = model2->numberRows();
3550              int numberColumns = model2->numberColumns();
3551              CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3552              CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3553              CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3554              CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3555              delete[] saveLower;
3556              delete[] saveUpper;
3557              saveLower = NULL;
3558              saveUpper = NULL;
3559            }
3560            //int numberRows = model2->numberRows();
3561            //int numberColumns = model2->numberColumns();
3562#ifdef ABC_INHERIT
3563            model2->checkSolution(0);
3564            printf("%d primal infeasibilities summing to %g\n",
3565              model2->numberPrimalInfeasibilities(),
3566              model2->sumPrimalInfeasibilities());
3567            model2->dealWithAbc(1, 1);
3568          }
3569        }
3570#else
3571            // just primal values pass
3572            double saveScale = model2->objectiveScale();
3573            model2->setObjectiveScale(1.0e-3);
3574            model2->primal(2);
3575            model2->setObjectiveScale(saveScale);
3576            // save primal solution and copy back dual
3577            CoinMemcpyN(model2->primalRowSolution(),
3578              numberRows, rowPrimal);
3579            CoinMemcpyN(rowDual,
3580              numberRows, model2->dualRowSolution());
3581            CoinMemcpyN(model2->primalColumnSolution(),
3582              numberColumns, columnPrimal);
3583            CoinMemcpyN(columnDual,
3584              numberColumns, model2->dualColumnSolution());
3585            //model2->primal(1);
3586            // clean up reduced costs and flag variables
3587            {
3588              double *dj = model2->dualColumnSolution();
3589              double *cost = model2->objective();
3590              double *saveCost = new double[numberColumns];
3591              CoinMemcpyN(cost, numberColumns, saveCost);
3592              double *saveLower = new double[numberColumns];
3593              double *lower = model2->columnLower();
3594              CoinMemcpyN(lower, numberColumns, saveLower);
3595              double *saveUpper = new double[numberColumns];
3596              double *upper = model2->columnUpper();
3597              CoinMemcpyN(upper, numberColumns, saveUpper);
3598              int i;
3599              double tolerance = 10.0 * dualTolerance_;
3600              for (i = 0; i < numberColumns; i++) {
3601                if (model2->getStatus(i) == basic) {
3602                  dj[i] = 0.0;
3603                } else if (model2->getStatus(i) == atLowerBound) {
3604                  if (optimizationDirection_ * dj[i] < tolerance) {
3605                    if (optimizationDirection_ * dj[i] < 0.0) {
3606                      //if (dj[i]<-1.0e-3)
3607                      //printf("bad dj at lb %d %g\n",i,dj[i]);
3608                      cost[i] -= dj[i];
3609                      dj[i] = 0.0;
3610                    }
3611                  } else {
3612                    upper[i] = lower[i];
3613                  }
3614                } else if (model2->getStatus(i) == atUpperBound) {
3615                  if (optimizationDirection_ * dj[i] > tolerance) {
3616                    if (optimizationDirection_ * dj[i] > 0.0) {
3617                      //if (dj[i]>1.0e-3)
3618                      //printf("bad dj at ub %d %g\n",i,dj[i]);
3619                      cost[i] -= dj[i];
3620                      dj[i] = 0.0;
3621                    }
3622                  } else {
3623                    lower[i] = upper[i];
3624                  }
3625                }
3626              }
3627              // just dual values pass
3628              //model2->setLogLevel(63);
3629              //model2->setFactorizationFrequency(1);
3630              if (!gap)
3631                model2->dual(2);
3632              CoinMemcpyN(saveCost, numberColumns, cost);
3633              delete[] saveCost;
3634              CoinMemcpyN(saveLower, numberColumns, lower);
3635              delete[] saveLower;
3636              CoinMemcpyN(saveUpper, numberColumns, upper);
3637              delete[] saveUpper;
3638            }
3639          }
3640          // and finish
3641          // move solutions
3642          CoinMemcpyN(rowPrimal,
3643            numberRows, model2->primalRowSolution());
3644          CoinMemcpyN(columnPrimal,
3645            numberColumns, model2->primalColumnSolution());
3646        }
3647        double saveScale = model2->objectiveScale();
3648        model2->setObjectiveScale(1.0e-3);
3649        model2->primal(2);
3650        model2->setObjectiveScale(saveScale);
3651        model2->primal(1);
3652#endif
3653#else
3654        // just primal
3655#ifdef ABC_INHERIT
3656        model2->checkSolution(0);
3657        printf("%d primal infeasibilities summing to %g\n",
3658          model2->numberPrimalInfeasibilities(),
3659          model2->sumPrimalInfeasibilities());
3660        model2->dealWithAbc(1, 1);
3661#else
3662        model2->primal(1);
3663#endif
3664        //model2->primal(1);
3665#endif
3666      } else if (barrierStatus == 4) {
3667        // memory problems
3668        model2->setPerturbation(savePerturbation);
3669        model2->createStatus();
3670        model2->dual();
3671      } else if (maxIts && quadraticObj) {
3672        // make sure no status left
3673        model2->createStatus();
3674        // solve
3675        model2->setPerturbation(100);
3676        model2->reducedGradient(1);
3677      }
3678    }
3679
3680    //model2->setMaximumIterations(saveMaxIts);
3681#ifdef BORROW
3682    model2->setNumberIterations(model2->numberIterations() + saveNumberIterations);
3683    delete[] rowPrimal;
3684    delete[] columnPrimal;
3685    delete[] rowDual;
3686    delete[] columnDual;
3687#endif
3688    if (extraPresolve) {
3689      pinfo2.postsolve(true);
3690      delete model2;
3691      model2 = saveModel2;
3692    }
3693    if (saveUpper) {
3694      if (!forceFixing) {
3695        int numberRows = model2->numberRows();
3696        int numberColumns = model2->numberColumns();
3697        CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3698        CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3699        CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3700        CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3701      }
3702      delete[] saveLower;
3703      delete[] saveUpper;
3704      saveLower = NULL;
3705      saveUpper = NULL;
3706      if (method != ClpSolve::useBarrierNoCross)
3707        model2->primal(1);
3708    }
3709    model2->setPerturbation(savePerturbation);
3710    time2 = CoinCpuTime();
3711    timeCore = time2 - timeX;
3712    handler_->message(CLP_INTERVAL_TIMING, messages_)
3713      << "Crossover" << timeCore << time2 - time1
3714      << CoinMessageEol;
3715    timeX = time2;
3716#else
3717    abort();
3718#endif
3719  } else if (method == ClpSolve::notImplemented) {
3720    printf("done decomposition\n");
3721  } else {
3722    assert(method != ClpSolve::automatic); // later
3723    time2 = 0.0;
3724  }
3725  if (saveMatrix) {
3726    if (model2 == this) {
3727      // delete and replace
3728      delete model2->clpMatrix();
3729      model2->replaceMatrix(saveMatrix);
3730    } else {
3731      delete saveMatrix;
3732    }
3733  }
3734  numberIterations = model2->numberIterations();
3735  finalStatus = model2->status();
3736  int finalSecondaryStatus = model2->secondaryStatus();
3737  if (presolve == ClpSolve::presolveOn) {
3738    int saveLevel = logLevel();
3739    if ((specialOptions_ & 1024) == 0)
3740      setLogLevel(CoinMin(1, saveLevel));
3741    else
3742      setLogLevel(CoinMin(0, saveLevel));
3743    pinfo->postsolve(true);
3744    numberIterations_ = 0;
3745    delete pinfo;
3746    pinfo = NULL;
3747    factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
3748    time2 = CoinCpuTime();
3749    timePresolve += time2 - timeX;
3750    handler_->message(CLP_INTERVAL_TIMING, messages_)
3751      << "Postsolve" << time2 - timeX << time2 - time1
3752      << CoinMessageEol;
3753    timeX = time2;
3754    if (!presolveToFile) {
3755#if 1 //ndef ABC_INHERIT
3756      delete model2;
3757#else
3758      if (model2->abcSimplex())
3759        delete model2->abcSimplex();
3760      else
3761        delete model2;
3762#endif
3763    }
3764    if (interrupt)
3765      currentModel = this;
3766    // checkSolution(); already done by postSolve
3767    setLogLevel(saveLevel);
3768    int oldStatus = problemStatus_;
3769    setProblemStatus(finalStatus);
3770    setSecondaryStatus(finalSecondaryStatus);
3771    /* Code modified so rcode -1 as normal, >0 clean up, 0 say optimal */
3772    int rcode = eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
3773    //#define TREAT_AS_OPTIMAL_TOLERANCE 1.0e-4
3774#ifdef TREAT_AS_OPTIMAL_TOLERANCE
3775    if (rcode == -1 && sumPrimalInfeasibilities_ < TREAT_AS_OPTIMAL_TOLERANCE && sumDualInfeasibilities_ < TREAT_AS_OPTIMAL_TOLERANCE)
3776      rcode = 0;
3777#endif
3778    if (finalStatus != 3 && rcode < 0 && (finalStatus || oldStatus == -1)) {
3779      double sumPrimal = sumPrimalInfeasibilities_;
3780      double sumDual = sumDualInfeasibilities_;
3781      if (sumDual > 1.0e-6 && sumPrimal > 1.0e-6)
3782        moreSpecialOptions_ &= ~2; // be safe and do final solve
3783      // ignore some parts of solution
3784      if (finalStatus == 1) {
3785        // infeasible
3786        sumDual = 0.0;
3787      } else if (finalStatus == 2) {
3788        sumPrimal = 0.0;
3789      }
3790      int savePerturbation = perturbation();
3791      if (savePerturbation == 50)
3792        setPerturbation(51); // small
3793      if (!finalStatus || finalStatus == 2 || (moreSpecialOptions_ & 2) == 0 || fabs(sumDual) + fabs(sumPrimal) < 1.0e-3) {
3794        if (finalStatus == 2) {
3795          if (sumDual > 1.0e-4) {
3796            // unbounded - get feasible first
3797            double save = optimizationDirection_;
3798            optimizationDirection_ = 0.0;
3799            primal(1);
3800            optimizationDirection_ = save;
3801          }
3802          primal(1);
3803        } else if (finalStatus == 1) {
3804          dual();
3805        } else {
3806          if ((moreSpecialOptions_ & 65536) == 0) {
3807            if (numberRows_ < 10000 || true)
3808              setPerturbation(100); // probably better to perturb after n its
3809            else if (savePerturbation < 100)
3810              setPerturbation(51); // probably better to perturb after n its
3811          }
3812#ifndef ABC_INHERIT
3813          // use method thought suitable
3814          int numberSuperBasic = 0;
3815          for (int i = 0; i < numberColumns_; i++) {
3816            if (getColumnStatus(i) == superBasic)
3817              numberSuperBasic++;
3818          }
3819          if (sumDual > 1000.0 * sumPrimal || numberSuperBasic) {
3820            primal(1);
3821          } else if (sumPrimal > 1000.0 * sumDual) {
3822            dual();
3823          } else {
3824            if (method != ClpSolve::useDual)
3825              primal(1);
3826            else
3827              dual();
3828          }
3829#else
3830          dealWithAbc(1, 2, interrupt);
3831#endif
3832        }
3833      } else {
3834        // just set status
3835        problemStatus_ = finalStatus;
3836      }
3837      setPerturbation(savePerturbation);
3838      numberIterations += numberIterations_;
3839      numberIterations_ = numberIterations;
3840      finalStatus = status();
3841      time2 = CoinCpuTime();
3842      handler_->message(CLP_INTERVAL_TIMING, messages_)
3843        << "Cleanup" << time2 - timeX << time2 - time1
3844        << CoinMessageEol;
3845      timeX = time2;
3846    } else if (rcode > 0) { // was >= 0
3847#ifdef ABC_INHERIT
3848      dealWithAbc(1, 2, true);
3849#else
3850      primal(1);
3851#endif
3852    } else {
3853      secondaryStatus_ = finalSecondaryStatus;
3854    }
3855  } else if (model2 != this) {
3856    // not presolved - but different model used (sprint probably)
3857    CoinMemcpyN(model2->primalRowSolution(),
3858      numberRows_, this->primalRowSolution());
3859    CoinMemcpyN(model2->dualRowSolution(),
3860      numberRows_, this->dualRowSolution());
3861    CoinMemcpyN(model2->primalColumnSolution(),
3862      numberColumns_, this->primalColumnSolution());
3863    CoinMemcpyN(model2->dualColumnSolution(),
3864      numberColumns_, this->dualColumnSolution());
3865    CoinMemcpyN(model2->statusArray(),
3866      numberColumns_ + numberRows_, this->statusArray());
3867    objectiveValue_ = model2->objectiveValue_;
3868    numberIterations_ = model2->numberIterations_;
3869    problemStatus_ = model2->problemStatus_;
3870    secondaryStatus_ = model2->secondaryStatus_;
3871    delete model2;
3872  }
3873  if (method != ClpSolve::useBarrierNoCross && method != ClpSolve::useBarrier)
3874    setMaximumIterations(saveMaxIterations);
3875  std::string statusMessage[] = { "Unknown", "Optimal", "PrimalInfeasible", "DualInfeasible", "Stopped",
3876    "Errors", "User stopped" };
3877  assert(finalStatus >= -1 && finalStatus <= 5);
3878  numberIterations_ = numberIterations;
3879  handler_->message(CLP_TIMING, messages_)
3880    << statusMessage[finalStatus + 1] << objectiveValue() << numberIterations << time2 - time1;
3881  handler_->printing(presolve == ClpSolve::presolveOn)
3882    << timePresolve;
3883  handler_->printing(timeIdiot != 0.0)
3884    << timeIdiot;
3885  handler_->message() << CoinMessageEol;
3886  if (interrupt)
3887    signal(SIGINT, saveSignal);
3888  perturbation_ = savePerturbation;
3889  scalingFlag_ = saveScaling;
3890  // If faking objective - put back correct one
3891  if (savedObjective) {
3892    delete objective_;
3893    objective_ = savedObjective;
3894  }
3895  if (options.getSpecialOption(1) == 2 && options.getExtraInfo(1) > 1000000) {
3896    ClpObjective *savedObjective = objective_;
3897    // make up zero objective
3898    double *obj = new double[numberColumns_];
3899    for (int i = 0; i < numberColumns_; i++)
3900      obj[i] = 0.0;
3901    objective_ = new ClpLinearObjective(obj, numberColumns_);
3902    delete[] obj;
3903    primal(1);
3904    delete objective_;
3905    objective_ = savedObjective;
3906    finalStatus = status();
3907  }
3908  eventHandler()->event(ClpEventHandler::presolveEnd);
3909  delete pinfo;
3910  moreSpecialOptions_ = saveMoreOptions;
3911#ifdef CLP_USEFUL_PRINTOUT
3912  debugInt[23] = numberIterations_;
3913#endif
3914  return finalStatus;
3915}
3916// General solve
3917int ClpSimplex::initialSolve()
3918{
3919  // Default so use dual
3920  ClpSolve options;
3921  return initialSolve(options);
3922}
3923// General dual solve
3924int ClpSimplex::initialDualSolve()
3925{
3926  ClpSolve options;
3927  // Use dual
3928  options.setSolveType(ClpSolve::useDual);
3929  return initialSolve(options);
3930}
3931// General primal solve
3932int ClpSimplex::initialPrimalSolve()
3933{
3934  ClpSolve options;
3935  // Use primal
3936  options.setSolveType(ClpSolve::usePrimal);
3937  return initialSolve(options);
3938}
3939// barrier solve, not to be followed by crossover
3940int ClpSimplex::initialBarrierNoCrossSolve()
3941{
3942  ClpSolve options;
3943  // Use primal
3944  options.setSolveType(ClpSolve::useBarrierNoCross);
3945  return initialSolve(options);
3946}
3947
3948// General barrier solve
3949int ClpSimplex::initialBarrierSolve()
3950{
3951  ClpSolve options;
3952  // Use primal
3953  options.setSolveType(ClpSolve::useBarrier);
3954  return initialSolve(options);
3955}
3956
3957// Default constructor
3958ClpSolve::ClpSolve()
3959{
3960  method_ = automatic;
3961  presolveType_ = presolveOn;
3962  numberPasses_ = 5;
3963  int i;
3964  for (i = 0; i < 7; i++)
3965    options_[i] = 0;
3966  // say no +-1 matrix
3967  options_[3] = 1;
3968  for (i = 0; i < 7; i++)
3969    extraInfo_[i] = -1;
3970  independentOptions_[0] = 0;
3971  // But switch off slacks
3972  independentOptions_[1] = 512;
3973  // Substitute up to 3
3974  independentOptions_[2] = 3;
3975}
3976// Constructor when you really know what you are doing
3977ClpSolve::ClpSolve(SolveType method, PresolveType presolveType,
3978  int numberPasses, int options[6],
3979  int extraInfo[6], int independentOptions[3])
3980{
3981  method_ = method;
3982  presolveType_ = presolveType;
3983  numberPasses_ = numberPasses;
3984  int i;
3985  for (i = 0; i < 6; i++)
3986    options_[i] = options[i];
3987  options_[6] = 0;
3988  for (i = 0; i < 6; i++)
3989    extraInfo_[i] = extraInfo[i];
3990  extraInfo_[6] = 0;
3991  for (i = 0; i < 3; i++)
3992    independentOptions_[i] = independentOptions[i];
3993}
3994
3995// Copy constructor.
3996ClpSolve::ClpSolve(const ClpSolve &rhs)
3997{
3998  method_ = rhs.method_;
3999  presolveType_ = rhs.presolveType_;
4000  numberPasses_ = rhs.numberPasses_;
4001  int i;
4002  for (i = 0; i < 7; i++)
4003    options_[i] = rhs.options_[i];
4004  for (i = 0; i < 7; i++)
4005    extraInfo_[i] = rhs.extraInfo_[i];
4006  for (i = 0; i < 3; i++)
4007    independentOptions_[i] = rhs.independentOptions_[i];
4008}
4009// Assignment operator. This copies the data
4010ClpSolve &
4011ClpSolve::operator=(const ClpSolve &rhs)
4012{
4013  if (this != &rhs) {
4014    method_ = rhs.method_;
4015    presolveType_ = rhs.presolveType_;
4016    numberPasses_ = rhs.numberPasses_;
4017    int i;
4018    for (i = 0; i < 7; i++)
4019      options_[i] = rhs.options_[i];
4020    for (i = 0; i < 7; i++)
4021      extraInfo_[i] = rhs.extraInfo_[i];
4022    for (i = 0; i < 3; i++)
4023      independentOptions_[i] = rhs.independentOptions_[i];
4024  }
4025  return *this;
4026}
4027// Destructor
4028ClpSolve::~ClpSolve()
4029{
4030}
4031// See header file for details
4032void ClpSolve::setSpecialOption(int which, int value, int extraInfo)
4033{
4034  options_[which] = value;
4035  extraInfo_[which] = extraInfo;
4036}
4037int ClpSolve::getSpecialOption(int which) const
4038{
4039  return options_[which];
4040}
4041
4042// Solve types
4043void ClpSolve::setSolveType(SolveType method, int /*extraInfo*/)
4044{
4045  method_ = method;
4046}
4047
4048ClpSolve::SolveType
4049ClpSolve::getSolveType()
4050{
4051  return method_;
4052}
4053
4054// Presolve types
4055void ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
4056{
4057  presolveType_ = amount;
4058  numberPasses_ = extraInfo;
4059}
4060ClpSolve::PresolveType
4061ClpSolve::getPresolveType()
4062{
4063  return presolveType_;
4064}
4065// Extra info for idiot (or sprint)
4066int ClpSolve::getExtraInfo(int which) const
4067{
4068  return extraInfo_[which];
4069}
4070int ClpSolve::getPresolvePasses() const
4071{
4072  return numberPasses_;
4073}
4074/* Say to return at once if infeasible,
4075   default is to solve */
4076void ClpSolve::setInfeasibleReturn(bool trueFalse)
4077{
4078  independentOptions_[0] = trueFalse ? 1 : 0;
4079}
4080#include <string>
4081// Generates code for above constructor
4082void ClpSolve::generateCpp(FILE *fp)
4083{
4084  std::string solveType[] = {
4085    "ClpSolve::useDual",
4086    "ClpSolve::usePrimal",
4087    "ClpSolve::usePrimalorSprint",
4088    "ClpSolve::useBarrier",
4089    "ClpSolve::useBarrierNoCross",
4090    "ClpSolve::automatic",
4091    "ClpSolve::notImplemented"
4092  };
4093  std::string presolveType[] = {
4094    "ClpSolve::presolveOn",
4095    "ClpSolve::presolveOff",
4096    "ClpSolve::presolveNumber",
4097    "ClpSolve::presolveNumberCost"
4098  };
4099  fprintf(fp, "3  ClpSolve::SolveType method = %s;\n", solveType[method_].c_str());
4100  fprintf(fp, "3  ClpSolve::PresolveType presolveType = %s;\n",
4101    presolveType[presolveType_].c_str());
4102  fprintf(fp, "3  int numberPasses = %d;\n", numberPasses_);
4103  fprintf(fp, "3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
4104    options_[0], options_[1], options_[2],
4105    options_[3], options_[4], options_[5]);
4106  fprintf(fp, "3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
4107    extraInfo_[0], extraInfo_[1], extraInfo_[2],
4108    extraInfo_[3], extraInfo_[4], extraInfo_[5]);
4109  fprintf(fp, "3  int independentOptions[] = {%d,%d,%d};\n",
4110    independentOptions_[0], independentOptions_[1], independentOptions_[2]);
4111  fprintf(fp, "3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
4112  fprintf(fp, "3                    options,extraInfo,independentOptions);\n");
4113}
4114//#############################################################################
4115#include "ClpNonLinearCost.hpp"
4116
4117ClpSimplexProgress::ClpSimplexProgress()
4118{
4119  int i;
4120  for (i = 0; i < CLP_PROGRESS; i++) {
4121    objective_[i] = COIN_DBL_MAX * 1.0e-50;
4122    infeasibility_[i] = -1.0; // set to an impossible value
4123    realInfeasibility_[i] = COIN_DBL_MAX * 1.0e-50;
4124    numberInfeasibilities_[i] = -1;
4125    iterationNumber_[i] = -1;
4126  }
4127#ifdef CLP_PROGRESS_WEIGHT
4128  for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4129    objectiveWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4130    infeasibilityWeight_[i] = -1.0; // set to an impossible value
4131    realInfeasibilityWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4132    numberInfeasibilitiesWeight_[i] = -1;
4133    iterationNumberWeight_[i] = -1;
4134  }
4135  drop_ = 0.0;
4136  best_ = 0.0;
4137#endif
4138  initialWeight_ = 0.0;
4139  for (i = 0; i < CLP_CYCLE; i++) {
4140    //obj_[i]=COIN_DBL_MAX*1.0e-50;
4141    in_[i] = -1;
4142    out_[i] = -1;
4143    way_[i] = 0;
4144  }
4145  numberTimes_ = 0;
4146  numberBadTimes_ = 0;
4147  numberReallyBadTimes_ = 0;
4148  numberTimesFlagged_ = 0;
4149  model_ = NULL;
4150  oddState_ = 0;
4151}
4152
4153//-----------------------------------------------------------------------------
4154
4155ClpSimplexProgress::~ClpSimplexProgress()
4156{
4157}
4158// Copy constructor.
4159ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
4160{
4161  int i;
4162  for (i = 0; i < CLP_PROGRESS; i++) {
4163    objective_[i] = rhs.objective_[i];
4164    infeasibility_[i] = rhs.infeasibility_[i];
4165    realInfeasibility_[i] = rhs.realInfeasibility_[i];
4166    numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
4167    iterationNumber_[i] = rhs.iterationNumber_[i];
4168  }
4169#ifdef CLP_PROGRESS_WEIGHT
4170  for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4171    objectiveWeight_[i] = rhs.objectiveWeight_[i];
4172    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4173    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4174    numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4175    iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4176  }
4177  drop_ = rhs.drop_;
4178  best_ = rhs.best_;
4179#endif
4180  initialWeight_ = rhs.initialWeight_;
4181  for (i = 0; i < CLP_CYCLE; i++) {
4182    //obj_[i]=rhs.obj_[i];
4183    in_[i] = rhs.in_[i];
4184    out_[i] = rhs.out_[i];
4185    way_[i] = rhs.way_[i];
4186  }
4187  numberTimes_ = rhs.numberTimes_;
4188  numberBadTimes_ = rhs.numberBadTimes_;
4189  numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4190  numberTimesFlagged_ = rhs.numberTimesFlagged_;
4191  model_ = rhs.model_;
4192  oddState_ = rhs.oddState_;
4193}
4194// Copy constructor.from model
4195ClpSimplexProgress::ClpSimplexProgress(ClpSimplex *model)
4196{
4197  model_ = model;
4198  reset();
4199  initialWeight_ = 0.0;
4200}
4201// Fill from model
4202void ClpSimplexProgress::fillFromModel(ClpSimplex *model)
4203{
4204  model_ = model;
4205  reset();
4206  initialWeight_ = 0.0;
4207}
4208// Assignment operator. This copies the data
4209ClpSimplexProgress &
4210ClpSimplexProgress::operator=(const ClpSimplexProgress &rhs)
4211{
4212  if (this != &rhs) {
4213    int i;
4214    for (i = 0; i < CLP_PROGRESS; i++) {
4215      objective_[i] = rhs.objective_[i];
4216      infeasibility_[i] = rhs.infeasibility_[i];
4217      realInfeasibility_[i] = rhs.realInfeasibility_[i];
4218      numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
4219      iterationNumber_[i] = rhs.iterationNumber_[i];
4220    }
4221#ifdef CLP_PROGRESS_WEIGHT
4222    for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4223      objectiveWeight_[i] = rhs.objectiveWeight_[i];
4224      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4225      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4226      numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4227      iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4228    }
4229    drop_ = rhs.drop_;
4230    best_ = rhs.best_;
4231#endif
4232    initialWeight_ = rhs.initialWeight_;
4233    for (i = 0; i < CLP_CYCLE; i++) {
4234      //obj_[i]=rhs.obj_[i];
4235      in_[i] = rhs.in_[i];
4236      out_[i] = rhs.out_[i];
4237      way_[i] = rhs.way_[i];
4238    }
4239    numberTimes_ = rhs.numberTimes_;
4240    numberBadTimes_ = rhs.numberBadTimes_;
4241    numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4242    numberTimesFlagged_ = rhs.numberTimesFlagged_;
4243    model_ = rhs.model_;
4244    oddState_ = rhs.oddState_;
4245  }
4246  return *this;
4247}
4248// Seems to be something odd about exact comparison of doubles on linux
4249static bool equalDouble(double value1, double value2)
4250{
4251
4252  union {
4253    double d;
4254    int i[2];
4255  } v1, v2;
4256  v1.d = value1;
4257  v2.d = value2;
4258  if (sizeof(int) * 2 == sizeof(double))
4259    return (v1.i[0] == v2.i[0] && v1.i[1] == v2.i[1]);
4260  else
4261    return (v1.i[0] == v2.i[0]);
4262}
4263int ClpSimplexProgress::looping()
4264{
4265  if (!model_)
4266    return -1;
4267  double objective;
4268  if (model_->algorithm() < 0) {
4269    objective = model_->rawObjectiveValue();
4270    objective -= model_->bestPossibleImprovement();
4271  } else {
4272    objective = model_->nonLinearCost()->feasibleReportCost();
4273  }
4274  double infeasibility;
4275  double realInfeasibility = 0.0;
4276  int numberInfeasibilities;
4277  int iterationNumber = model_->numberIterations();
4278  //numberTimesFlagged_ = 0;
4279  if (model_->algorithm() < 0) {
4280    // dual
4281    infeasibility = model_->sumPrimalInfeasibilities();
4282    numberInfeasibilities = model_->numberPrimalInfeasibilities();
4283  } else {
4284    //primal
4285    infeasibility = model_->sumDualInfeasibilities();
4286    realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
4287    numberInfeasibilities = model_->numberDualInfeasibilities();
4288  }
4289  int i;
4290  int numberMatched = 0;
4291  int matched = 0;
4292  int nsame = 0;
4293  for (i = 0; i < CLP_PROGRESS; i++) {
4294    bool matchedOnObjective = equalDouble(objective, objective_[i]);
4295    bool matchedOnInfeasibility = equalDouble(infeasibility, infeasibility_[i]);
4296    bool matchedOnInfeasibilities = (numberInfeasibilities == numberInfeasibilities_[i]);
4297
4298    if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
4299      matched |= (1 << i);
4300      // Check not same iteration
4301      if (iterationNumber != iterationNumber_[i]) {
4302        numberMatched++;
4303        // here mainly to get over compiler bug?
4304        if (model_->messageHandler()->logLevel() > 10)
4305          printf("%d %d %d %d %d loop check\n", i, numberMatched,
4306            matchedOnObjective, matchedOnInfeasibility,
4307            matchedOnInfeasibilities);
4308      } else {
4309        // stuck but code should notice
4310        nsame++;
4311      }
4312    }
4313    if (i) {
4314      objective_[i - 1] = objective_[i];
4315      infeasibility_[i - 1] = infeasibility_[i];
4316      realInfeasibility_[i - 1] = realInfeasibility_[i];
4317      numberInfeasibilities_[i - 1] = numberInfeasibilities_[i];
4318      iterationNumber_[i - 1] = iterationNumber_[i];
4319    }
4320  }
4321  objective_[CLP_PROGRESS - 1] = objective;
4322  infeasibility_[CLP_PROGRESS - 1] = infeasibility;
4323  realInfeasibility_[CLP_PROGRESS - 1] = realInfeasibility;
4324  numberInfeasibilities_[CLP_PROGRESS - 1] = numberInfeasibilities;
4325  iterationNumber_[CLP_PROGRESS - 1] = iterationNumber;
4326  if (nsame == CLP_PROGRESS)
4327    numberMatched = CLP_PROGRESS; // really stuck
4328  if (model_->progressFlag())
4329    numberMatched = 0;
4330  numberTimes_++;
4331  if (numberTimes_ < 10)
4332    numberMatched = 0;
4333  // skip if just last time as may be checking something
4334  if (matched == (1 << (CLP_PROGRESS - 1)))
4335    numberMatched = 0;
4336  if (numberMatched && model_->clpMatrix()->type() < 15) {
4337    model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
4338      << numberMatched
4339      << matched
4340      << numberTimes_
4341      << CoinMessageEol;
4342    numberBadTimes_++;
4343    if (numberBadTimes_ < 10) {
4344      // make factorize every iteration
4345      model_->forceFactorization(1);
4346      if (numberBadTimes_ < 2) {
4347        startCheck(); // clear other loop check
4348        if (model_->algorithm() < 0) {
4349          // dual - change tolerance
4350          model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
4351          // if infeasible increase dual bound
4352          if (model_->dualBound() < 1.0e17) {
4353            model_->setDualBound(model_->dualBound() * 1.1);
4354            static_cast< ClpSimplexDual * >(model_)->resetFakeBounds(0);
4355          }
4356        } else {
4357          // primal - change tolerance
4358          if (numberBadTimes_ > 3)
4359            model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
4360          // if infeasible increase infeasibility cost
4361          if (model_->nonLinearCost()->numberInfeasibilities() && model_->infeasibilityCost() < 1.0e17) {
4362            model_->setInfeasibilityCost(model_->infeasibilityCost() * 1.1);
4363          }
4364        }
4365      } else {
4366        // flag
4367        int iSequence;
4368        if (model_->algorithm() < 0) {
4369          // dual
4370          if (model_->dualBound() > 1.0e14)
4371            model_->setDualBound(1.0e14);
4372          iSequence = in_[CLP_CYCLE - 1];
4373        } else {
4374          // primal
4375          //if (model_->infeasibilityCost() > 1.0e14)
4376          //   model_->setInfeasibilityCost(1.0e14);
4377          iSequence = out_[CLP_CYCLE - 1];
4378        }
4379        if (iSequence >= 0) {
4380          char x = model_->isColumn(iSequence) ? 'C' : 'R';
4381          if (model_->messageHandler()->logLevel() >= 63)
4382            model_->messageHandler()->message(CLP_SIMPLEX_FLAG, model_->messages())
4383              << x << model_->sequenceWithin(iSequence)
4384              << CoinMessageEol;
4385          // if Gub then needs to be sequenceIn_
4386          int save = model_->sequenceIn();
4387          model_->setSequenceIn(iSequence);
4388          model_->setFlagged(iSequence);
4389          model_->setSequenceIn(save);
4390          //printf("flagging %d from loop\n",iSequence);
4391          startCheck();
4392        } else {
4393          // Give up
4394          if (model_->messageHandler()->logLevel() >= 63)
4395            printf("***** All flagged?\n");
4396          return 4;
4397        }
4398        // reset
4399        numberBadTimes_ = 2;
4400      }
4401      return -2;
4402    } else {
4403      // look at solution and maybe declare victory
4404      if (infeasibility < 1.0e-4) {
4405        return 0;
4406      } else {
4407        model_->messageHandler()->message(CLP_LOOP, model_->messages())
4408          << CoinMessageEol;
4409#ifndef NDEBUG
4410        printf("debug loop ClpSimplex A\n");
4411        abort();
4412#endif
4413        return 3;
4414      }
4415    }
4416  }
4417  return -1;
4418}
4419// Resets as much as possible
4420void ClpSimplexProgress::reset()
4421{
4422  int i;
4423  for (i = 0; i < CLP_PROGRESS; i++) {
4424    if (model_->algorithm() >= 0)
4425      objective_[i] = COIN_DBL_MAX * 1.0e-50;
4426    else
4427      objective_[i] = -COIN_DBL_MAX * 1.0e-50;
4428    infeasibility_[i] = -1.0; // set to an impossible value
4429    realInfeasibility_[i] = COIN_DBL_MAX * 1.0e-50;
4430    numberInfeasibilities_[i] = -1;
4431    iterationNumber_[i] = -1;
4432  }
4433#ifdef CLP_PROGRESS_WEIGHT
4434  for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4435    objectiveWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4436    infeasibilityWeight_[i] = -1.0; // set to an impossible value
4437    realInfeasibilityWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4438    numberInfeasibilitiesWeight_[i] = -1;
4439    iterationNumberWeight_[i] = -1;
4440  }
4441  drop_ = 0.0;
4442  best_ = 0.0;
4443#endif
4444  for (i = 0; i < CLP_CYCLE; i++) {
4445    //obj_[i]=COIN_DBL_MAX*1.0e-50;
4446    in_[i] = -1;
4447    out_[i] = -1;
4448    way_[i] = 0;
4449  }
4450  numberTimes_ = 0;
4451  numberBadTimes_ = 0;
4452  numberReallyBadTimes_ = 0;
4453  numberTimesFlagged_ = 0;
4454  oddState_ = 0;
4455}
4456// Returns previous objective (if -1) - current if (0)
4457double
4458ClpSimplexProgress::lastObjective(int back) const
4459{
4460  return objective_[CLP_PROGRESS - 1 - back];
4461}
4462// Returns previous infeasibility (if -1) - current if (0)
4463double
4464ClpSimplexProgress::lastInfeasibility(int back) const
4465{
4466  return realInfeasibility_[CLP_PROGRESS - 1 - back];
4467}
4468// Sets real primal infeasibility
4469void ClpSimplexProgress::setInfeasibility(double value)
4470{
4471  for (int i = 1; i < CLP_PROGRESS; i++)
4472    realInfeasibility_[i - 1] = realInfeasibility_[i];
4473  realInfeasibility_[CLP_PROGRESS - 1] = value;
4474}
4475// Returns number of primal infeasibilities (if -1) - current if (0)
4476int ClpSimplexProgress::numberInfeasibilities(int back) const
4477{
4478  return numberInfeasibilities_[CLP_PROGRESS - 1 - back];
4479}
4480// Modify objective e.g. if dual infeasible in dual
4481void ClpSimplexProgress::modifyObjective(double value)
4482{
4483  objective_[CLP_PROGRESS - 1] = value;
4484}
4485// Returns previous iteration number (if -1) - current if (0)
4486int ClpSimplexProgress::lastIterationNumber(int back) const
4487{
4488  return iterationNumber_[CLP_PROGRESS - 1 - back];
4489}
4490// clears iteration numbers (to switch off panic)
4491void ClpSimplexProgress::clearIterationNumbers()
4492{
4493  for (int i = 0; i < CLP_PROGRESS; i++)
4494    iterationNumber_[i] = -1;
4495}
4496// Start check at beginning of whileIterating
4497void ClpSimplexProgress::startCheck()
4498{
4499  int i;
4500  for (i = 0; i < CLP_CYCLE; i++) {
4501    //obj_[i]=COIN_DBL_MAX*1.0e-50;
4502    in_[i] = -1;
4503    out_[i] = -1;
4504    way_[i] = 0;
4505  }
4506}
4507// Returns cycle length in whileIterating
4508int ClpSimplexProgress::cycle(int in, int out, int wayIn, int wayOut)
4509{
4510  int i;
4511#if 0
4512     if (model_->numberIterations() > 206571) {
4513          printf("in %d out %d\n", in, out);
4514          for (i = 0; i < CLP_CYCLE; i++)
4515               printf("cy %d in %d out %d\n", i, in_[i], out_[i]);
4516     }
4517#endif
4518  int matched = 0;
4519  // first see if in matches any out
4520  for (i = 1; i < CLP_CYCLE; i++) {
4521    if (in == out_[i]) {
4522      // even if flip then suspicious
4523      matched = -1;
4524      break;
4525    }
4526  }
4527#if 0
4528     if (!matched || in_[0] < 0) {
4529          // can't be cycle
4530          for (i = 0; i < CLP_CYCLE - 1; i++) {
4531               //obj_[i]=obj_[i+1];
4532               in_[i] = in_[i+1];
4533               out_[i] = out_[i+1];
4534               way_[i] = way_[i+1];
4535          }
4536     } else {
4537          // possible cycle
4538          matched = 0;
4539          for (i = 0; i < CLP_CYCLE - 1; i++) {
4540               int k;
4541               char wayThis = way_[i];
4542               int inThis = in_[i];
4543               int outThis = out_[i];
4544               //double objThis = obj_[i];
4545               for(k = i + 1; k < CLP_CYCLE; k++) {
4546                    if (inThis == in_[k] && outThis == out_[k] && wayThis == way_[k]) {
4547                         int distance = k - i;
4548                         if (k + distance < CLP_CYCLE) {
4549                              // See if repeats
4550                              int j = k + distance;
4551                              if (inThis == in_[j] && outThis == out_[j] && wayThis == way_[j]) {
4552                                   matched = distance;
4553                                   break;
4554                              }
4555                         } else {
4556                              matched = distance;
4557                              break;
4558                         }
4559                    }
4560               }
4561               //obj_[i]=obj_[i+1];
4562               in_[i] = in_[i+1];
4563               out_[i] = out_[i+1];
4564               way_[i] = way_[i+1];
4565          }
4566     }
4567#else
4568  if (matched && in_[0] >= 0) {
4569    // possible cycle - only check [0] against all
4570    matched = 0;
4571    int nMatched = 0;
4572    char way0 = way_[0];
4573    int in0 = in_[0];
4574    int out0 = out_[0];
4575    //double obj0 = obj_[i];
4576    for (int k = 1; k < CLP_CYCLE - 4; k++) {
4577      if (in0 == in_[k] && out0 == out_[k] && way0 == way_[k]) {
4578        nMatched++;
4579        // See if repeats
4580        int end = CLP_CYCLE - k;
4581        int j;
4582        for (j = 1; j < end; j++) {
4583          if (in_[j + k] != in_[j] || out_[j + k] != out_[j] || way_[j + k] != way_[j])
4584            break;
4585        }
4586        if (j == end) {
4587          matched = k;
4588          break;
4589        }
4590      }
4591    }
4592    // If three times then that is too much even if not regular
4593    if (matched <= 0 && nMatched > 1)
4594      matched = 100;
4595  }
4596  for (i = 0; i < CLP_CYCLE - 1; i++) {
4597    //obj_[i]=obj_[i+1];
4598    in_[i] = in_[i + 1];
4599    out_[i] = out_[i + 1];
4600    way_[i] = way_[i + 1];
4601  }
4602#endif
4603  int way = 1 - wayIn + 4 * (1 - wayOut);
4604  //obj_[i]=model_->objectiveValue();
4605  in_[CLP_CYCLE - 1] = in;
4606  out_[CLP_CYCLE - 1] = out;
4607  way_[CLP_CYCLE - 1] = static_cast< char >(way);
4608  return matched;
4609}
4610#include "CoinStructuredModel.hpp"
4611// Solve using structure of model and maybe in parallel
4612int ClpSimplex::solve(CoinStructuredModel *model)
4613{
4614  // analyze structure
4615  int numberRowBlocks = model->numberRowBlocks();
4616  int numberColumnBlocks = model->numberColumnBlocks();
4617  int numberElementBlocks = model->numberElementBlocks();
4618  if (numberElementBlocks == 1) {
4619    loadProblem(*model, false);
4620    return dual();
4621  }
4622  // For now just get top level structure
4623  CoinModelBlockInfo *blockInfo = new CoinModelBlockInfo[numberElementBlocks];
4624  for (int i = 0; i < numberElementBlocks; i++) {
4625    CoinStructuredModel *subModel = dynamic_cast< CoinStructuredModel * >(model->block(i));
4626    CoinModel *thisBlock;
4627    if (subModel) {
4628      thisBlock = subModel->coinModelBlock(blockInfo[i]);
4629      model->setCoinModel(thisBlock, i);
4630    } else {
4631      thisBlock = dynamic_cast< CoinModel * >(model->block(i));
4632      assert(thisBlock);
4633      // just fill in info
4634      CoinModelBlockInfo info = CoinModelBlockInfo();
4635      int whatsSet = thisBlock->whatIsSet();
4636      info.matrix = static_cast< char >(((whatsSet & 1) != 0) ? 1 : 0);
4637      info.rhs = static_cast< char >(((whatsSet & 2) != 0) ? 1 : 0);
4638      info.rowName = static_cast< char >(((whatsSet & 4) != 0) ? 1 : 0);
4639      info.integer = static_cast< char >(((whatsSet & 32) != 0) ? 1 : 0);
4640      info.bounds = static_cast< char >(((whatsSet & 8) != 0) ? 1 : 0);
4641      info.columnName = static_cast< char >(((whatsSet & 16) != 0) ? 1 : 0);
4642      // Which block
4643      int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4644      info.rowBlock = iRowBlock;
4645      int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4646      info.columnBlock = iColumnBlock;
4647      blockInfo[i] = info;
4648    }
4649  }
4650  int *rowCounts = new int[numberRowBlocks];
4651  CoinZeroN(rowCounts, numberRowBlocks);
4652  int *columnCounts = new int[numberColumnBlocks + 1];
4653  CoinZeroN(columnCounts, numberColumnBlocks);
4654  int decomposeType = 0;
4655  for (int i = 0; i < numberElementBlocks; i++) {
4656    int iRowBlock = blockInfo[i].rowBlock;
4657    int iColumnBlock = blockInfo[i].columnBlock;
4658    rowCounts[iRowBlock]++;
4659    columnCounts[iColumnBlock]++;
4660  }
4661  if (numberRowBlocks == numberColumnBlocks || numberRowBlocks == numberColumnBlocks + 1) {
4662    // could be Dantzig-Wolfe
4663    int numberG1 = 0;
4664    for (int i = 0; i < numberRowBlocks; i++) {
4665      if (rowCounts[i] > 1)
4666        numberG1++;
4667    }
4668    bool masterColumns = (numberColumnBlocks == numberRowBlocks);
4669    if ((masterColumns && numberElementBlocks == 2 * numberRowBlocks - 1)
4670      || (!masterColumns && numberElementBlocks == 2 * numberRowBlocks)) {
4671      if (numberG1 < 2)
4672        decomposeType = 1;
4673    }
4674  }
4675  if (!decomposeType && (numberRowBlocks == numberColumnBlocks || numberRowBlocks == numberColumnBlocks - 1)) {
4676    // could be Benders
4677    int numberG1 = 0;
4678    for (int i = 0; i < numberColumnBlocks; i++) {
4679      if (columnCounts[i] > 1)
4680        numberG1++;
4681    }
4682    bool masterRows = (numberColumnBlocks == numberRowBlocks);
4683    if ((masterRows && numberElementBlocks == 2 * numberColumnBlocks - 1)
4684      || (!masterRows && numberElementBlocks == 2 * numberColumnBlocks)) {
4685      if (numberG1 < 2)
4686        decomposeType = 2;
4687    }
4688  }
4689  delete[] rowCounts;
4690  delete[] columnCounts;
4691  delete[] blockInfo;
4692  // decide what to do
4693  ClpSolve options;
4694  options.setIndependentOption(2, 100);
4695  switch (decomposeType) {
4696    // No good
4697  case 0:
4698    loadProblem(*model, false);
4699    return dual();
4700    // DW
4701  case 1:
4702    return solveDW(model, options);
4703    // Benders
4704  case 2:
4705    return solveBenders(model, options);
4706  }
4707  return 0; // to stop compiler warning
4708}
4709/* This loads a model from a CoinStructuredModel object - returns number of errors.
4710   If originalOrder then keep to order stored in blocks,
4711   otherwise first column/rows correspond to first block - etc.
4712   If keepSolution true and size is same as current then
4713   keeps current status and solution
4714*/
4715int ClpSimplex::loadProblem(CoinStructuredModel &coinModel,
4716  bool originalOrder,
4717  bool keepSolution)
4718{
4719  unsigned char *status = NULL;
4720  double *psol = NULL;
4721  double *dsol = NULL;
4722  int numberRows = coinModel.numberRows();
4723  int numberColumns = coinModel.numberColumns();
4724  int numberRowBlocks = coinModel.numberRowBlocks();
4725  int numberColumnBlocks = coinModel.numberColumnBlocks();
4726  int numberElementBlocks = coinModel.numberElementBlocks();
4727  if (status_ && numberRows_ && numberRows_ == numberRows && numberColumns_ == numberColumns && keepSolution) {
4728    status = new unsigned char[numberRows_ + numberColumns_];
4729    CoinMemcpyN(status_, numberRows_ + numberColumns_, status);
4730    psol = new double[numberRows_ + numberColumns_];
4731    CoinMemcpyN(columnActivity_, numberColumns_, psol);
4732    CoinMemcpyN(rowActivity_, numberRows_, psol + numberColumns_);
4733    dsol = new double[numberRows_ + numberColumns_];
4734    CoinMemcpyN(reducedCost_, numberColumns_, dsol);
4735    CoinMemcpyN(dual_, numberRows_, dsol + numberColumns_);
4736  }
4737  int returnCode = 0;
4738  double *rowLower = new double[numberRows];
4739  double *rowUpper = new double[numberRows];
4740  double *columnLower = new double[numberColumns];
4741  double *columnUpper = new double[numberColumns];
4742  double *objective = new double[numberColumns];
4743  int *integerType = new int[numberColumns];
4744  CoinBigIndex numberElements = 0;
4745  // Bases for blocks
4746  int *rowBase = new int[numberRowBlocks];
4747  CoinFillN(rowBase, numberRowBlocks, -1);
4748  // And row to put it
4749  int *whichRow = new int[numberRows + numberRowBlocks];
4750  int *columnBase = new int[numberColumnBlocks];
4751  CoinFillN(columnBase, numberColumnBlocks, -1);
4752  // And column to put it
4753  int *whichColumn = new int[numberColumns + numberColumnBlocks];
4754  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4755    CoinModel *block = coinModel.coinBlock(iBlock);
4756    numberElements += block->numberElements();
4757    //and set up elements etc
4758    double *associated = block->associatedArray();
4759    // If strings then do copies
4760    if (block->stringsExist())
4761      returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
4762        objective, integerType, associated);
4763    const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4764    int iRowBlock = info.rowBlock;
4765    int iColumnBlock = info.columnBlock;
4766    if (rowBase[iRowBlock] < 0) {
4767      rowBase[iRowBlock] = block->numberRows();
4768      // Save block number
4769      whichRow[numberRows + iRowBlock] = iBlock;
4770    } else {
4771      assert(rowBase[iRowBlock] == block->numberRows());
4772    }
4773    if (columnBase[iColumnBlock] < 0) {
4774      columnBase[iColumnBlock] = block->numberColumns();
4775      // Save block number
4776      whichColumn[numberColumns + iColumnBlock] = iBlock;
4777    } else {
4778      assert(columnBase[iColumnBlock] == block->numberColumns());
4779    }
4780  }
4781  // Fill arrays with defaults
4782  CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
4783  CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
4784  CoinFillN(columnLower, numberColumns, 0.0);
4785  CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
4786  CoinFillN(objective, numberColumns, 0.0);
4787  CoinFillN(integerType, numberColumns, 0);
4788  int n = 0;
4789  for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4790    int k = rowBase[iBlock];
4791    rowBase[iBlock] = n;
4792    assert(k >= 0);
4793    // block number
4794    int jBlock = whichRow[numberRows + iBlock];
4795    if (originalOrder) {
4796      memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
4797    } else {
4798      CoinIotaN(whichRow + n, k, n);
4799    }
4800    n += k;
4801  }
4802  assert(n == numberRows);
4803  n = 0;
4804  for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4805    int k = columnBase[iBlock];
4806    columnBase[iBlock] = n;
4807    assert(k >= 0);
4808    if (k) {
4809      // block number
4810      int jBlock = whichColumn[numberColumns + iBlock];
4811      if (originalOrder) {
4812        memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
4813          k * sizeof(int));
4814      } else {
4815        CoinIotaN(whichColumn + n, k, n);
4816      }
4817      n += k;
4818    }
4819  }
4820  assert(n == numberColumns);
4821  bool gotIntegers = false;
4822  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4823    CoinModel *block = coinModel.coinBlock(iBlock);
4824    const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4825    int iRowBlock = info.rowBlock;
4826    int iRowBase = rowBase[iRowBlock];
4827    int iColumnBlock = info.columnBlock;
4828    int iColumnBase = columnBase[iColumnBlock];
4829    if (info.rhs) {
4830      int nRows = block->numberRows();
4831      const double *lower = block->rowLowerArray();
4832      const double *upper = block->rowUpperArray();
4833      for (int i = 0; i < nRows; i++) {
4834        int put = whichRow[i + iRowBase];
4835        rowLower[put] = lower[i];
4836        rowUpper[put] = upper[i];
4837      }
4838    }
4839    if (info.bounds) {
4840      int nColumns = block->numberColumns();
4841      const double *lower = block->columnLowerArray();
4842      const double *upper = block->columnUpperArray();
4843      const double *obj = block->objectiveArray();
4844      for (int i = 0; i < nColumns; i++) {
4845        int put = whichColumn[i + iColumnBase];
4846        columnLower[put] = lower[i];
4847        columnUpper[put] = upper[i];
4848        objective[put] = obj[i];
4849      }
4850    }
4851    if (info.integer) {
4852      gotIntegers = true;
4853      int nColumns = block->numberColumns();
4854      const int *type = block->integerTypeArray();
4855      for (int i = 0; i < nColumns; i++) {
4856        int put = whichColumn[i + iColumnBase];
4857        integerType[put] = type[i];
4858      }
4859    }
4860  }
4861  gutsOfLoadModel(numberRows, numberColumns,
4862    columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
4863  delete[] rowLower;
4864  delete[] rowUpper;
4865  delete[] columnLower;
4866  delete[] columnUpper;
4867  delete[] objective;
4868  // Do integers if wanted
4869  if (gotIntegers) {
4870    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4871      if (integerType[iColumn])
4872        setInteger(iColumn);
4873    }
4874  }
4875  delete[] integerType;
4876  setObjectiveOffset(coinModel.objectiveOffset());
4877  // Space for elements
4878  int *row = new int[numberElements];
4879  int *column = new int[numberElements];
4880  double *element = new double[numberElements];
4881  numberElements = 0;
4882  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4883    CoinModel *block = coinModel.coinBlock(iBlock);
4884    const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4885    int iRowBlock = info.rowBlock;
4886    int iRowBase = rowBase[iRowBlock];
4887    int iColumnBlock = info.columnBlock;
4888    int iColumnBase = columnBase[iColumnBlock];
4889    if (info.rowName) {
4890      int numberItems = block->rowNames()->numberItems();
4891      assert(block->numberRows() >= numberItems);
4892      if (numberItems) {
4893        const char *const *rowNames = block->rowNames()->names();
4894        for (int i = 0; i < numberItems; i++) {
4895          int put = whichRow[i + iRowBase];
4896          std::string name = rowNames[i];
4897          setRowName(put, name);
4898        }
4899      }
4900    }
4901    if (info.columnName) {
4902      int numberItems = block->columnNames()->numberItems();
4903      assert(block->numberColumns() >= numberItems);
4904      if (numberItems) {
4905        const char *const *columnNames = block->columnNames()->names();
4906        for (int i = 0; i < numberItems; i++) {
4907          int put = whichColumn[i + iColumnBase];
4908          std::string name = columnNames[i];
4909          setColumnName(put, name);
4910        }
4911      }
4912    }
4913    if (info.matrix) {
4914      CoinPackedMatrix matrix2;
4915      const CoinPackedMatrix *matrix = block->packedMatrix();
4916      if (!matrix) {
4917        double *associated = block->associatedArray();
4918        block->createPackedMatrix(matrix2, associated);
4919        matrix = &matrix2;
4920      }
4921      // get matrix data pointers
4922      const int *row2 = matrix->getIndices();
4923      const CoinBigIndex *columnStart = matrix->getVectorStarts();
4924      const double *elementByColumn = matrix->getElements();
4925      const int *columnLength = matrix->getVectorLengths();
4926      int n = matrix->getNumCols();
4927      assert(matrix->isColOrdered());
4928      for (int iColumn = 0; iColumn < n; iColumn++) {
4929        CoinBigIndex j;
4930        int jColumn = whichColumn[iColumn + iColumnBase];
4931        for (j = columnStart[iColumn];
4932             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4933          row[numberElements] = whichRow[row2[j] + iRowBase];
4934          column[numberElements] = jColumn;
4935          element[numberElements++] = elementByColumn[j];
4936        }
4937      }
4938    }
4939  }
4940  delete[] whichRow;
4941  delete[] whichColumn;
4942  delete[] rowBase;
4943  delete[] columnBase;
4944  CoinPackedMatrix *matrix = new CoinPackedMatrix(true, row, column, element, numberElements);
4945  matrix_ = new ClpPackedMatrix(matrix);
4946  matrix_->setDimensions(numberRows, numberColumns);
4947  delete[] row;
4948  delete[] column;
4949  delete[] element;
4950  createStatus();
4951  if (status) {
4952    // copy back
4953    CoinMemcpyN(status, numberRows_ + numberColumns_, status_);
4954    CoinMemcpyN(psol, numberColumns_, columnActivity_);
4955    CoinMemcpyN(psol + numberColumns_, numberRows_, rowActivity_);
4956    CoinMemcpyN(dsol, numberColumns_, reducedCost_);
4957    CoinMemcpyN(dsol + numberColumns_, numberRows_, dual_);
4958    delete[] status;
4959    delete[] psol;
4960    delete[] dsol;
4961  }
4962  optimizationDirection_ = coinModel.optimizationDirection();
4963  return returnCode;
4964}
4965/*  If input negative scales objective so maximum <= -value
4966    and returns scale factor used.  If positive unscales and also
4967    redoes dual stuff
4968*/
4969double
4970ClpSimplex::scaleObjective(double value)
4971{
4972  double *obj = objective();
4973  double largest = 0.0;
4974  if (value < 0.0) {
4975    value = -value;
4976    for (int i = 0; i < numberColumns_; i++) {
4977      largest = CoinMax(largest, fabs(obj[i]));
4978    }
4979    if (largest > value) {
4980      double scaleFactor = value / largest;
4981      for (int i = 0; i < numberColumns_; i++) {
4982        obj[i] *= scaleFactor;
4983        reducedCost_[i] *= scaleFactor;
4984      }
4985      for (int i = 0; i < numberRows_; i++) {
4986        dual_[i] *= scaleFactor;
4987      }
4988      largest /= value;
4989    } else {
4990      // no need
4991      largest = 1.0;
4992    }
4993  } else {
4994    // at end
4995    if (value != 1.0) {
4996      for (int i = 0; i < numberColumns_; i++) {
4997        obj[i] *= value;
4998        reducedCost_[i] *= value;
4999      }
5000      for (int i = 0; i < numberRows_; i++) {
5001        dual_[i] *= value;
5002      }
5003      computeObjectiveValue();
5004    }
5005  }
5006  return largest;
5007}
5008#if defined(ABC_INHERIT) || defined(CBC_THREAD) || defined(THREADS_IN_ANALYZE)
5009void *clp_parallelManager(void *stuff)
5010{
5011  CoinPthreadStuff *driver = reinterpret_cast< CoinPthreadStuff * >(stuff);
5012  int whichThread = driver->whichThread();
5013  CoinThreadInfo *threadInfo = driver->threadInfoPointer(whichThread);
5014  threadInfo->status = -1;
5015  int *which = threadInfo->stuff;
5016#ifdef PTHREAD_BARRIER_SERIAL_THREAD
5017  pthread_barrier_wait(driver->barrierPointer());
5018#endif
5019#if 0
5020  int status=-1;
5021  while (status!=100)
5022    status=timedWait(driver,1000,2);
5023  pthread_cond_signal(driver->conditionPointer(1));
5024  pthread_mutex_unlock(driver->mutexPointer(1,whichThread));
5025#endif
5026  // so now mutex_ is locked
5027  int whichLocked = 0;
5028  while (true) {
5029    pthread_mutex_t *mutexPointer = driver->mutexPointer(whichLocked, whichThread);
5030    // wait
5031    //printf("Child waiting for %d - status %d %d %d\n",
5032    //     whichLocked,lockedX[0],lockedX[1],lockedX[2]);
5033#ifdef DETAIL_THREAD
5034    printf("thread %d about to lock mutex %d\n", whichThread, whichLocked);
5035#endif
5036    pthread_mutex_lock(mutexPointer);
5037    whichLocked++;
5038    if (whichLocked == 3)
5039      whichLocked = 0;
5040    int unLock = whichLocked + 1;
5041    if (unLock == 3)
5042      unLock = 0;
5043    //printf("child pointer %x status %d\n",threadInfo,threadInfo->status);
5044    assert(threadInfo->status >= 0);
5045    if (threadInfo->status == 1000)
5046      pthread_exit(NULL);
5047    int type = threadInfo->status;
5048    int &returnCode = which[0];
5049    int iPass = which[1];
5050    ClpSimplex *clpSimplex = reinterpret_cast< ClpSimplex * >(threadInfo->extraInfo);
5051    //CoinIndexedVector * array;
5052    //double dummy;
5053    switch (type) {
5054      // dummy
5055    case 0:
5056      break;
5057    case 1:
5058      if (!clpSimplex->problemStatus() || !iPass)
5059        returnCode = clpSimplex->dual();
5060      else
5061        returnCode = clpSimplex->primal();
5062      break;
5063    case 100:
5064      // initialization
5065      break;
5066    }
5067    threadInfo->status = -1;
5068#ifdef DETAIL_THREAD
5069    printf("thread %d about to unlock mutex %d\n", whichThread, unLock);
5070#endif
5071    pthread_mutex_unlock(driver->mutexPointer(unLock, whichThread));
5072  }
5073}
5074#endif
5075// Solve using Dantzig-Wolfe decomposition and maybe in parallel
5076int ClpSimplex::solveDW(CoinStructuredModel *model, ClpSolve &options)
5077{
5078  double time1 = CoinCpuTime();
5079  int numberColumns = model->numberColumns();
5080  int numberRowBlocks = model->numberRowBlocks();
5081  int numberColumnBlocks = model->numberColumnBlocks();
5082  int numberElementBlocks = model->numberElementBlocks();
5083  // We already have top level structure
5084  CoinModelBlockInfo *blockInfo = new CoinModelBlockInfo[numberElementBlocks];
5085  for (int i = 0; i < numberElementBlocks; i++) {
5086    CoinModel *thisBlock = model->coinBlock(i);
5087    assert(thisBlock);
5088    // just fill in info
5089    CoinModelBlockInfo info = CoinModelBlockInfo();
5090    int whatsSet = thisBlock->whatIsSet();
5091    info.matrix = static_cast< char >(((whatsSet & 1) != 0) ? 1 : 0);
5092    info.rhs = static_cast< char >(((whatsSet & 2) != 0) ? 1 : 0);
5093    info.rowName = static_cast< char >(((whatsSet & 4) != 0) ? 1 : 0);
5094    info.integer = static_cast< char >(((whatsSet & 32) != 0) ? 1 : 0);
5095    info.bounds = static_cast< char >(((whatsSet & 8) != 0) ? 1 : 0);
5096    info.columnName = static_cast< char >(((whatsSet & 16) != 0) ? 1 : 0);
5097    // Which block
5098    int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
5099    info.rowBlock = iRowBlock;
5100    int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
5101    info.columnBlock = iColumnBlock;
5102    blockInfo[i] = info;
5103  }
5104  // make up problems
5105  int numberBlocks = numberRowBlocks - 1;
5106  // Find master rows and columns
5107  int *rowCounts = new int[numberRowBlocks];
5108  CoinZeroN(rowCounts, numberRowBlocks);
5109  int *columnCounts = new int[numberColumnBlocks + 1];
5110  CoinZeroN(columnCounts, numberColumnBlocks);
5111  int iBlock;
5112  for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5113    int iRowBlock = blockInfo[iBlock].rowBlock;
5114    rowCounts[iRowBlock]++;
5115    int iColumnBlock = blockInfo[iBlock].columnBlock;
5116    columnCounts[iColumnBlock]++;
5117  }
5118  int *whichBlock = new int[numberElementBlocks];
5119  int masterRowBlock = -1;
5120  for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5121    if (rowCounts[iBlock] > 1) {
5122      if (masterRowBlock == -1) {
5123        masterRowBlock = iBlock;
5124      } else {
5125        // Can't decode
5126        masterRowBlock = -2;
5127        break;
5128      }
5129    }
5130  }
5131  int masterColumnBlock = -1;
5132  int kBlock = 0;
5133  for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5134    int count = columnCounts[iBlock];
5135    columnCounts[iBlock] = kBlock;
5136    kBlock += count;
5137  }
5138  for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5139    int iColumnBlock = blockInfo[iBlock].columnBlock;
5140    whichBlock[columnCounts[iColumnBlock]] = iBlock;
5141    columnCounts[iColumnBlock]++;
5142  }
5143  for (iBlock = numberColumnBlocks - 1; iBlock >= 0; iBlock--)
5144    columnCounts[iBlock + 1] = columnCounts[iBlock];
5145  columnCounts[0] = 0;
5146  for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5147    int count = columnCounts[iBlock + 1] - columnCounts[iBlock];
5148    if (count == 1) {
5149      int kBlock = whichBlock[columnCounts[iBlock]];
5150      int iRowBlock = blockInfo[kBlock].rowBlock;
5151      if (iRowBlock == masterRowBlock) {
5152        if (masterColumnBlock == -1) {
5153          masterColumnBlock = iBlock;
5154        } else {
5155          // Can't decode
5156          masterColumnBlock = -2;
5157          break;
5158        }
5159      }
5160    }
5161  }
5162  if (masterRowBlock < 0 || masterColumnBlock == -2) {
5163    // What now
5164    abort();
5165  }
5166  delete[] rowCounts;
5167  // create all data
5168  const CoinPackedMatrix **top = new const CoinPackedMatrix *[numberColumnBlocks];
5169  ClpSimplex *sub = new ClpSimplex[numberBlocks];
5170  ClpSimplex master;
5171  // Set offset
5172  master.setObjectiveOffset(model->objectiveOffset());
5173  bool reducePrint = logLevel() == 7;
5174  if (reducePrint) {
5175    // special
5176    setLogLevel(1);
5177    master.setLogLevel(1);
5178  }
5179  kBlock = 0;
5180  int masterBlock = -1;
5181  for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5182    top[kBlock] = NULL;
5183    int start = columnCounts[iBlock];
5184    int end = columnCounts[iBlock + 1];
5185    assert(end - start <= 2);
5186    for (int j = start; j < end; j++) {
5187      int jBlock = whichBlock[j];
5188      int iRowBlock = blockInfo[jBlock].rowBlock;
5189      int iColumnBlock = blockInfo[jBlock].columnBlock;
5190      assert(iColumnBlock == iBlock);
5191      if (iColumnBlock != masterColumnBlock && iRowBlock == masterRowBlock) {
5192        // top matrix
5193        top[kBlock] = model->coinBlock(jBlock)->packedMatrix();
5194      } else {
5195        const CoinPackedMatrix *matrix
5196          = model->coinBlock(jBlock)->packedMatrix();
5197        // Get pointers to arrays
5198        const double *rowLower;
5199        const double *rowUpper;
5200        const double *columnLower;
5201        const double *columnUpper;
5202        const double *objective;
5203        model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
5204          columnLower, columnUpper, objective);
5205        if (iColumnBlock != masterColumnBlock) {
5206          // diagonal block
5207          sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
5208            objective, rowLower, rowUpper);
5209          if (true) {
5210            double *lower = sub[kBlock].columnLower();
5211            double *upper = sub[kBlock].columnUpper();
5212            int n = sub[kBlock].numberColumns();
5213            for (int i = 0; i < n; i++) {
5214              lower[i] = CoinMax(-1.0e8, lower[i]);
5215              upper[i] = CoinMin(1.0e8, upper[i]);
5216            }
5217          }
5218          if (optimizationDirection_ < 0.0) {
5219            double *obj = sub[kBlock].objective();
5220            int n = sub[kBlock].numberColumns();
5221            for (int i = 0; i < n; i++)
5222              obj[i] = -obj[i];
5223          }
5224          if (this->factorizationFrequency() == 200) {
5225            // User did not touch preset
5226            sub[kBlock].defaultFactorizationFrequency();
5227          } else {
5228            // make sure model has correct value
5229            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
5230          }
5231          sub[kBlock].setPerturbation(50);
5232          // Set columnCounts to be diagonal block index for cleanup
5233          columnCounts[kBlock] = jBlock;
5234        } else {
5235          // master
5236          masterBlock = jBlock;
5237          master.loadProblem(*matrix, columnLower, columnUpper,
5238            objective, rowLower, rowUpper);
5239          if (optimizationDirection_ < 0.0) {
5240            double *obj = master.objective();
5241            int n = master.numberColumns();
5242            for (int i = 0; i < n; i++)
5243              obj[i] = -obj[i];
5244          }
5245        }
5246      }
5247    }
5248    if (iBlock != masterColumnBlock)
5249      kBlock++;
5250  }
5251  delete[] whichBlock;
5252  delete[] blockInfo;
5253  // For now master must have been defined (does not have to have columns)
5254  assert(master.numberRows());
5255  assert(masterBlock >= 0);
5256  int numberMasterRows = master.numberRows();
5257  // Overkill in terms of space
5258  int spaceNeeded = CoinMax(numberBlocks * (numberMasterRows + 1),
5259    2 * numberMasterRows);
5260  int *rowAdd = new int[spaceNeeded];
5261  double *elementAdd = new double[spaceNeeded];
5262  spaceNeeded = numberBlocks;
5263  CoinBigIndex *columnAdd = new CoinBigIndex[spaceNeeded + 1];
5264  double *objective = new double[spaceNeeded];
5265  // Add in costed slacks
5266  int firstArtificial = master.numberColumns();
5267  int lastArtificial = firstArtificial;
5268  if (true) {
5269    const double *lower = master.rowLower();
5270    const double *upper = master.rowUpper();
5271    int kCol = 0;
5272    for (int iRow = 0; iRow < numberMasterRows; iRow++) {
5273      if (lower[iRow] > -1.0e10) {
5274        rowAdd[kCol] = iRow;
5275        elementAdd[kCol++] = 1.0;
5276      }
5277      if (upper[iRow] < 1.0e10) {
5278        rowAdd[kCol] = iRow;
5279        elementAdd[kCol++] = -1.0;
5280      }
5281    }
5282    if (kCol > spaceNeeded) {
5283      spaceNeeded = kCol;
5284      delete[] columnAdd;
5285      delete[] objective;
5286      columnAdd = new CoinBigIndex[spaceNeeded + 1];
5287      objective = new double[spaceNeeded];
5288    }
5289    for (int i = 0; i < kCol; i++) {
5290      columnAdd[i] = i;
5291      objective[i] = 1.0e13;
5292    }
5293    columnAdd[kCol] = kCol;
5294    master.addColumns(kCol, NULL, NULL, objective,
5295      columnAdd, rowAdd, elementAdd);
5296    lastArtificial = master.numberColumns();
5297  }
5298  int maxPass = options.independentOption(2);
5299  if (maxPass < 2)
5300    maxPass = 100;
5301  int iPass;
5302  double lastObjective = 1.0e31;
5303  // Create convexity rows for proposals
5304  int numberMasterColumns = master.numberColumns();
5305  master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
5306  if (this->factorizationFrequency() == 200) {
5307    // User did not touch preset
5308    master.defaultFactorizationFrequency();
5309  } else {
5310    // make sure model has correct value
5311    master.setFactorizationFrequency(this->factorizationFrequency());
5312  }
5313  master.setPerturbation(50);
5314  // Arrays to say which block and when created
5315  int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
5316  whichBlock = new int[maximumColumns];
5317  int *when = new int[maximumColumns];
5318  int numberColumnsGenerated = numberBlocks;
5319  // fill in rhs and add in artificials
5320  {
5321    double *rowLower = master.rowLower();
5322    double *rowUpper = master.rowUpper();
5323    int iBlock;
5324    columnAdd[0] = 0;
5325    for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
5326      int iRow = iBlock + numberMasterRows;
5327      ;
5328      rowLower[iRow] = 1.0;
5329      rowUpper[iRow] = 1.0;
5330      rowAdd[iBlock] = iRow;
5331      elementAdd[iBlock] = 1.0;
5332      objective[iBlock] = 1.0e13;
5333      columnAdd[iBlock + 1] = iBlock + 1;
5334      when[iBlock] = -1;
5335      whichBlock[iBlock] = iBlock;
5336    }
5337    master.addColumns(numberBlocks, NULL, NULL, objective,
5338      columnAdd, rowAdd, elementAdd);
5339  }
5340  char generalPrint[200];
5341  // and resize matrix to double check clp will be happy
5342  //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
5343  //                     numberMasterColumns+numberBlocks);
5344  sprintf(generalPrint, "Time to decompose %.2f seconds", CoinCpuTime() - time1);
5345  handler_->message(CLP_GENERAL, messages_)
5346    << generalPrint
5347    << CoinMessageEol;
5348#ifdef ABC_INHERIT
5349  //AbcSimplex abcMaster;
5350  //if (!this->abcState())
5351  //setAbcState(1);
5352  int numberCpu = CoinMin((this->abcState() & 15), 4);
5353  CoinPthreadStuff threadInfo(numberCpu, clp_parallelManager);
5354  master.setAbcState(this->abcState());
5355  //AbcSimplex * tempMaster=master.dealWithAbc(2,10,true);
5356  //abcMaster=*tempMaster;
5357  //delete tempMaster;
5358  //abcMaster.startThreads(numberCpu);
5359  //#define master abcMaster
5360#endif
5361  for (iPass = 0; iPass < maxPass; iPass++) {
5362    sprintf(generalPrint, "Start of pass %d", iPass);
5363    handler_->message(CLP_GENERAL, messages_)
5364      << generalPrint
5365      << CoinMessageEol;
5366    // Solve master - may be infeasible
5367    //master.scaling(0);
5368    if (0) {
5369      master.writeMps("yy.mps");
5370    }
5371    // Correct artificials
5372    double sumArtificials = 0.0;
5373    if (iPass) {
5374      double *upper = master.columnUpper();
5375      double *solution = master.primalColumnSolution();
5376      double *obj = master.objective();
5377      sumArtificials = 0.0;
5378      for (int i = firstArtificial; i < lastArtificial; i++) {
5379        sumArtificials += solution[i];
5380        //assert (solution[i]>-1.0e-2);
5381        if (solution[i] < 1.0e-6) {
5382#if 0
5383                         // Could take out
5384                         obj[i] = 0.0;
5385                         upper[i] = 0.0;
5386#else
5387          obj[i] = 1.0e7;
5388          upper[i] = 1.0e-1;
5389#endif
5390          solution[i] = 0.0;
5391          master.setColumnStatus(i, isFixed);
5392        } else {
5393          upper[i] = solution[i] + 1.0e-5 * (1.0 + solution[i]);
5394        }
5395      }
5396      sprintf(generalPrint, "Sum of artificials before solve is %g", sumArtificials);
5397      handler_->message(CLP_GENERAL, messages_)
5398        << generalPrint
5399        << CoinMessageEol;
5400    }
5401    // scale objective to be reasonable
5402    double scaleFactor = master.scaleObjective(-1.0e9);
5403    {
5404      double *dual = master.dualRowSolution();
5405      int n = master.numberRows();
5406      memset(dual, 0, n * sizeof(double));
5407      double *solution = master.primalColumnSolution();
5408      master.clpMatrix()->times(1.0, solution, dual);
5409      double sum = 0.0;
5410      double *lower = master.rowLower();
5411      double *upper = master.rowUpper();
5412      for (int iRow = 0; iRow < n; iRow++) {
5413        double value = dual[iRow];
5414        if (value > upper[iRow])
5415          sum += value - upper[iRow];
5416        else if (value < lower[iRow])
5417          sum -= value - lower[iRow];
5418      }
5419      printf("** suminf %g\n", sum);
5420      lower = master.columnLower();
5421      upper = master.columnUpper();
5422      n = master.numberColumns();
5423      for (int iColumn = 0; iColumn < n; iColumn++) {
5424        double value = solution[iColumn];
5425        if (value > upper[iColumn] + 1.0e-5)
5426          sum += value - upper[iColumn];
5427        else if (value < lower[iColumn] - 1.0e-5)
5428          sum -= value - lower[iColumn];
5429      }
5430      printf("** suminf %g\n", sum);
5431    }
5432#ifdef ABC_INHERIT
5433    master.dealWithAbc(1, 0, true);
5434#else
5435    master.primal();
5436#endif
5437    //master.primal(1);
5438    // Correct artificials
5439    sumArtificials = 0.0;
5440    {
5441      double *solution = master.primalColumnSolution();
5442      for (int i = firstArtificial; i < lastArtificial; i++) {
5443        sumArtificials += solution[i];
5444      }
5445      printf("** Sum of artificials after solve is %g\n", sumArtificials);
5446    }
5447    master.scaleObjective(scaleFactor);
5448    int problemStatus = master.status(); // do here as can change (delcols)
5449    if (problemStatus == 2 && master.numberColumns()) {
5450#ifdef ABC_INHERIT
5451      master.dealWithAbc(1, 1, true);
5452#else
5453      master.primal(1);
5454#endif
5455      //master.primal(1);
5456      if (problemStatus == 2)