source: stable/1.16/Clp/src/ClpSolve.cpp @ 2183

Last change on this file since 2183 was 2183, checked in by forrest, 5 years ago

out flaky perturbation

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 294.1 KB
Line 
1/* $Id: ClpSolve.cpp 2183 2015-10-28 10:18:34Z forrest $ */
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#include "ClpCholeskyMumps.hpp"
73#ifdef COIN_HAS_VOL
74#include "VolVolume.hpp"
75#include "CoinHelperFunctions.hpp"
76#include "CoinPackedMatrix.hpp"
77#include "CoinMpsIO.hpp"
78//#############################################################################
79
80class lpHook : public VOL_user_hooks {
81private:
82     lpHook(const lpHook&);
83     lpHook& operator=(const lpHook&);
84private:
85     /// Pointer to dense vector of structural variable upper bounds
86     double  *colupper_;
87     /// Pointer to dense vector of structural variable lower bounds
88     double  *collower_;
89     /// Pointer to dense vector of objective coefficients
90     double  *objcoeffs_;
91     /// Pointer to dense vector of right hand sides
92     double  *rhs_;
93     /// Pointer to dense vector of senses
94     char    *sense_;
95
96     /// The problem matrix in a row ordered form
97     CoinPackedMatrix rowMatrix_;
98     /// The problem matrix in a column ordered form
99     CoinPackedMatrix colMatrix_;
100
101public:
102     lpHook(double* clb, double* cub, double* obj,
103            double* rhs, char* sense, const CoinPackedMatrix& mat);
104     virtual ~lpHook();
105
106public:
107     // for all hooks: return value of -1 means that volume should quit
108     /** compute reduced costs
109         @param u (IN) the dual variables
110         @param rc (OUT) the reduced cost with respect to the dual values
111     */
112     virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc);
113
114     /** Solve the subproblem for the subgradient step.
115         @param dual (IN) the dual variables
116         @param rc (IN) the reduced cost with respect to the dual values
117         @param lcost (OUT) the lagrangean cost with respect to the dual values
118         @param x (OUT) the primal result of solving the subproblem
119         @param v (OUT) b-Ax for the relaxed constraints
120         @param pcost (OUT) the primal objective value of <code>x</code>
121     */
122     virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
123                                  double& lcost, VOL_dvector& x, VOL_dvector& v,
124                                  double& pcost);
125     /** Starting from the primal vector x, run a heuristic to produce
126         an integer solution
127         @param x (IN) the primal vector
128         @param heur_val (OUT) the value of the integer solution (return
129         <code>DBL_MAX</code> here if no feas sol was found
130     */
131     virtual int heuristics(const VOL_problem& p,
132                            const VOL_dvector& x, double& heur_val) {
133          return 0;
134     }
135};
136
137//#############################################################################
138
139lpHook::lpHook(double* clb, double* cub, double* obj,
140               double* rhs, char* sense,
141               const CoinPackedMatrix& mat)
142{
143     colupper_ = cub;
144     collower_ = clb;
145     objcoeffs_ = obj;
146     rhs_ = rhs;
147     sense_ = sense;
148     assert (mat.isColOrdered());
149     colMatrix_.copyOf(mat);
150     rowMatrix_.reverseOrderedCopyOf(mat);
151}
152
153//-----------------------------------------------------------------------------
154
155lpHook::~lpHook()
156{
157}
158
159//#############################################################################
160
161int
162lpHook::compute_rc(const VOL_dvector& u, VOL_dvector& rc)
163{
164     rowMatrix_.transposeTimes(u.v, rc.v);
165     const int psize = rowMatrix_.getNumCols();
166
167     for (int i = 0; i < psize; ++i)
168          rc[i] = objcoeffs_[i] - rc[i];
169     return 0;
170}
171
172//-----------------------------------------------------------------------------
173
174int
175lpHook::solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
176                         double& lcost, VOL_dvector& x, VOL_dvector& v,
177                         double& pcost)
178{
179     int i;
180     const int psize = x.size();
181     const int dsize = v.size();
182
183     // compute the lagrangean solution corresponding to the reduced costs
184     for (i = 0; i < psize; ++i)
185          x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
186
187     // compute the lagrangean value (rhs*dual + primal*rc)
188     lcost = 0;
189     for (i = 0; i < dsize; ++i)
190          lcost += rhs_[i] * dual[i];
191     for (i = 0; i < psize; ++i)
192          lcost += x[i] * rc[i];
193
194     // compute the rhs - lhs
195     colMatrix_.times(x.v, v.v);
196     for (i = 0; i < dsize; ++i)
197          v[i] = rhs_[i] - v[i];
198
199     // compute the lagrangean primal objective
200     pcost = 0;
201     for (i = 0; i < psize; ++i)
202          pcost += x[i] * objcoeffs_[i];
203
204     return 0;
205}
206
207//#############################################################################
208/** A quick inlined function to convert from lb/ub style constraint
209    definition to sense/rhs/range style */
210inline void
211convertBoundToSense(const double lower, const double upper,
212                    char& sense, double& right,
213                    double& range)
214{
215     range = 0.0;
216     if (lower > -1.0e20) {
217          if (upper < 1.0e20) {
218               right = upper;
219               if (upper == lower) {
220                    sense = 'E';
221               } else {
222                    sense = 'R';
223                    range = upper - lower;
224               }
225          } else {
226               sense = 'G';
227               right = lower;
228          }
229     } else {
230          if (upper < 1.0e20) {
231               sense = 'L';
232               right = upper;
233          } else {
234               sense = 'N';
235               right = 0.0;
236          }
237     }
238}
239
240static int
241solveWithVolume(ClpSimplex * model, int numberPasses, int doIdiot)
242{
243     VOL_problem volprob;
244     volprob.parm.gap_rel_precision = 0.00001;
245     volprob.parm.maxsgriters = 3000;
246     if(numberPasses > 3000) {
247          volprob.parm.maxsgriters = numberPasses;
248          volprob.parm.primal_abs_precision = 0.0;
249          volprob.parm.minimum_rel_ascent = 0.00001;
250     } else if (doIdiot > 0) {
251          volprob.parm.maxsgriters = doIdiot;
252     }
253     if (model->logLevel() < 2)
254          volprob.parm.printflag = 0;
255     else
256          volprob.parm.printflag = 3;
257     const CoinPackedMatrix* mat = model->matrix();
258     int psize = model->numberColumns();
259     int dsize = model->numberRows();
260     char * sense = new char[dsize];
261     double * rhs = new double [dsize];
262
263     // Set the lb/ub on the duals
264     volprob.dsize = dsize;
265     volprob.psize = psize;
266     volprob.dual_lb.allocate(dsize);
267     volprob.dual_ub.allocate(dsize);
268     int i;
269     const double * rowLower = model->rowLower();
270     const double * rowUpper = model->rowUpper();
271     for (i = 0; i < dsize; ++i) {
272          double range;
273          convertBoundToSense(rowLower[i], rowUpper[i],
274                              sense[i], rhs[i], range);
275          switch (sense[i]) {
276          case 'E':
277               volprob.dual_lb[i] = -1.0e31;
278               volprob.dual_ub[i] = 1.0e31;
279               break;
280          case 'L':
281               volprob.dual_lb[i] = -1.0e31;
282               volprob.dual_ub[i] = 0.0;
283               break;
284          case 'G':
285               volprob.dual_lb[i] = 0.0;
286               volprob.dual_ub[i] = 1.0e31;
287               break;
288          default:
289               printf("Volume Algorithm can't work if there is a non ELG row\n");
290               return 1;
291          }
292     }
293     // Check bounds
294     double * saveLower = model->columnLower();
295     double * saveUpper = model->columnUpper();
296     bool good = true;
297     for (i = 0; i < psize; i++) {
298          if (saveLower[i] < -1.0e20 || saveUpper[i] > 1.0e20) {
299               good = false;
300               break;
301          }
302     }
303     if (!good) {
304          saveLower = CoinCopyOfArray(model->columnLower(), psize);
305          saveUpper = CoinCopyOfArray(model->columnUpper(), psize);
306          for (i = 0; i < psize; i++) {
307               if (saveLower[i] < -1.0e20)
308                    saveLower[i] = -1.0e20;
309               if(saveUpper[i] > 1.0e20)
310                    saveUpper[i] = 1.0e20;
311          }
312     }
313     lpHook myHook(saveLower, saveUpper,
314                   model->objective(),
315                   rhs, sense, *mat);
316
317     volprob.solve(myHook, false /* no warmstart */);
318
319     if (saveLower != model->columnLower()) {
320          delete [] saveLower;
321          delete [] saveUpper;
322     }
323     //------------- extract the solution ---------------------------
324
325     //printf("Best lagrangean value: %f\n", volprob.value);
326
327     double avg = 0;
328     for (i = 0; i < dsize; ++i) {
329          switch (sense[i]) {
330          case 'E':
331               avg += CoinAbs(volprob.viol[i]);
332               break;
333          case 'L':
334               if (volprob.viol[i] < 0)
335                    avg +=  (-volprob.viol[i]);
336               break;
337          case 'G':
338               if (volprob.viol[i] > 0)
339                    avg +=  volprob.viol[i];
340               break;
341          }
342     }
343
344     //printf("Average primal constraint violation: %f\n", avg/dsize);
345
346     // volprob.dsol contains the dual solution (dual feasible)
347     // volprob.psol contains the primal solution
348     //              (NOT necessarily primal feasible)
349     CoinMemcpyN(volprob.dsol.v, dsize, model->dualRowSolution());
350     CoinMemcpyN(volprob.psol.v, psize, model->primalColumnSolution());
351     return 0;
352}
353#endif
354static ClpInterior * currentModel2 = NULL;
355#endif
356//#############################################################################
357// Allow for interrupts
358// But is this threadsafe ? (so switched off by option)
359
360#include "CoinSignal.hpp"
361static ClpSimplex * currentModel = NULL;
362#ifdef ABC_INHERIT
363static AbcSimplex * currentAbcModel = NULL;
364#endif
365
366extern "C" {
367     static void
368#if defined(_MSC_VER)
369     __cdecl
370#endif // _MSC_VER
371     signal_handler(int /*whichSignal*/)
372     {
373          if (currentModel != NULL)
374               currentModel->setMaximumIterations(0); // stop at next iterations
375#ifdef ABC_INHERIT
376          if (currentAbcModel != NULL)
377               currentAbcModel->setMaximumIterations(0); // stop at next iterations
378#endif
379#ifndef SLIM_CLP
380          if (currentModel2 != NULL)
381               currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
382#endif
383          return;
384     }
385}
386#if ABC_INSTRUMENT>1
387int abcPricing[20];
388int abcPricingDense[20];
389static int trueNumberRows;
390static int numberTypes;
391#define MAX_TYPES 25
392#define MAX_COUNT 20
393#define MAX_FRACTION 101
394static char * types[MAX_TYPES];
395static double counts[MAX_TYPES][MAX_COUNT];
396static double countsFraction[MAX_TYPES][MAX_FRACTION];
397static double * currentCounts;
398static double * currentCountsFraction;
399static int currentType;
400static double workMultiplier[MAX_TYPES];
401static double work[MAX_TYPES];
402static double currentWork;
403static double otherWork[MAX_TYPES];
404static int timesCalled[MAX_TYPES];
405static int timesStarted[MAX_TYPES];
406static int fractionDivider;
407void instrument_initialize(int numberRows)
408{
409  trueNumberRows=numberRows;
410  numberTypes=0;
411  memset(counts,0,sizeof(counts));
412  currentCounts=NULL;
413  memset(countsFraction,0,sizeof(countsFraction));
414  currentCountsFraction=NULL;
415  memset(workMultiplier,0,sizeof(workMultiplier));
416  memset(work,0,sizeof(work));
417  memset(otherWork,0,sizeof(otherWork));
418  memset(timesCalled,0,sizeof(timesCalled));
419  memset(timesStarted,0,sizeof(timesStarted));
420  currentType=-1;
421  fractionDivider=(numberRows+MAX_FRACTION-2)/(MAX_FRACTION-1);
422}
423void instrument_start(const char * type,int numberRowsEtc)
424{
425  if (currentType>=0)
426    instrument_end();
427  currentType=-1;
428  currentWork=0.0;
429  for (int i=0;i<numberTypes;i++) {
430    if (!strcmp(types[i],type)) {
431      currentType=i;
432      break;
433    }
434  }
435  if (currentType==-1) {
436    assert (numberTypes<MAX_TYPES);
437    currentType=numberTypes;
438    types[numberTypes++]=strdup(type);
439  }
440  currentCounts = &counts[currentType][0];
441  currentCountsFraction = &countsFraction[currentType][0];
442  timesStarted[currentType]++;
443  assert (trueNumberRows);
444  workMultiplier[currentType]+=static_cast<double>(numberRowsEtc)/static_cast<double>(trueNumberRows);
445}
446void instrument_add(int count)
447{
448  assert (currentType>=0);
449  currentWork+=count;
450  timesCalled[currentType]++;
451  if (count<MAX_COUNT-1)
452    currentCounts[count]++;
453  else
454    currentCounts[MAX_COUNT-1]++;
455  assert(count/fractionDivider>=0&&count/fractionDivider<MAX_FRACTION);
456  currentCountsFraction[count/fractionDivider]++;
457}
458void instrument_do(const char * type,double count)
459{
460  int iType=-1;
461  for (int i=0;i<numberTypes;i++) {
462    if (!strcmp(types[i],type)) {
463      iType=i;
464      break;
465    }
466  }
467  if (iType==-1) {
468    assert (numberTypes<MAX_TYPES);
469    iType=numberTypes;
470    types[numberTypes++]=strdup(type);
471  }
472  timesStarted[iType]++;
473  otherWork[iType]+=count;
474}
475void instrument_end()
476{
477  work[currentType]+=currentWork;
478  currentType=-1;
479}
480void instrument_end_and_adjust(double factor)
481{
482  work[currentType]+=currentWork*factor;
483  currentType=-1;
484}
485void instrument_print()
486{
487  for (int iType=0;iType<numberTypes;iType++) {
488    currentCounts = &counts[iType][0];
489    currentCountsFraction = &countsFraction[iType][0];
490    if (!otherWork[iType]) {
491      printf("%s started %d times, used %d times, work %g (average length %.1f) multiplier %g\n",
492             types[iType],timesStarted[iType],timesCalled[iType],
493             work[iType],work[iType]/(timesCalled[iType]+1.0e-100),workMultiplier[iType]/(timesStarted[iType]+1.0e-100));
494      int n=0;
495      for (int i=0;i<MAX_COUNT-1;i++) {
496        if (currentCounts[i]) {
497          if (n==5) {
498            n=0;
499            printf("\n");
500          }
501          n++;
502          printf("(%d els,%.0f times) ",i,currentCounts[i]);
503        }
504      }
505      if (currentCounts[MAX_COUNT-1]) {
506        if (n==5) {
507          n=0;
508          printf("\n");
509        }
510        n++;
511        printf("(>=%d els,%.0f times) ",MAX_COUNT-1,currentCounts[MAX_COUNT-1]);
512      }
513      printf("\n");
514      int largestFraction;
515      int nBig=0;
516      for (largestFraction=MAX_FRACTION-1;largestFraction>=10;largestFraction--) {
517        double count = currentCountsFraction[largestFraction];
518        if (count&&largestFraction>10)
519          nBig++;
520        if (nBig>4)
521          break;
522      }
523      int chunk=(largestFraction+5)/10;
524      int lo=0;
525      for (int iChunk=0;iChunk<largestFraction;iChunk+=chunk) {
526        int hi=CoinMin(lo+chunk*fractionDivider,trueNumberRows);
527        double sum=0.0;
528        for (int i=iChunk;i<CoinMin(iChunk+chunk,MAX_FRACTION);i++)
529          sum += currentCountsFraction[i];
530        if (sum)
531          printf("(%d-%d %.0f) ",lo,hi,sum);
532        lo=hi;
533      }
534      for (int i=lo/fractionDivider;i<MAX_FRACTION;i++) {
535        if (currentCountsFraction[i])
536          printf("(%d %.0f) ",i*fractionDivider,currentCountsFraction[i]);
537      }
538      printf("\n");
539    } else {
540      printf("%s started %d times, used %d times, work %g multiplier %g other work %g\n",
541             types[iType],timesStarted[iType],timesCalled[iType],
542             work[iType],workMultiplier[iType],otherWork[iType]);
543    }
544    free(types[iType]);
545  }
546}
547#endif
548#if ABC_PARALLEL==2
549#ifndef FAKE_CILK
550int number_cilk_workers=0;
551#include <cilk/cilk_api.h>
552#endif
553#endif
554#ifdef ABC_INHERIT
555AbcSimplex *
556ClpSimplex::dealWithAbc(int solveType, int startUp,
557                        bool interrupt)
558{
559  bool keepAbc = startUp>=10;
560  startUp=startUp%10;
561  AbcSimplex * abcModel2 = NULL;
562  if (!this->abcState()||!numberRows_||!numberColumns_) {
563    //this->readBasis("aaa.bas");
564    if (!solveType) {
565      this->dual(0);
566    } else if (solveType==1) {
567      int ifValuesPass=startUp ? 1 : 0;
568      if (startUp==3)
569        ifValuesPass=2;
570      this->primal(ifValuesPass);
571    }
572    //this->writeBasis("a.bas",true);
573  } else {
574    abcModel2=new AbcSimplex(*this);
575    if (interrupt)
576      currentAbcModel = abcModel2;
577    //if (abcSimplex_) {
578    // move factorization stuff
579    abcModel2->factorization()->synchronize(this->factorization(),abcModel2);
580    //}
581    //abcModel2->startPermanentArrays();
582    int crashState=abcModel2->abcState()&(256+512+1024);
583    abcModel2->setAbcState(CLP_ABC_WANTED|crashState|(abcModel2->abcState()&15));
584    int ifValuesPass=startUp ? 1 : 0;
585    if (startUp==3)
586      ifValuesPass=2;
587    // temp
588    if (fabs(this->primalTolerance()-1.001e-6)<0.999e-9) {
589      int type=1;
590      double diff=this->primalTolerance()-1.001e-6;
591      printf("Diff %g\n",diff);
592      if (fabs(diff-1.0e-10)<1.0e-13)
593        type=2;
594      else if (fabs(diff-1.0e-11)<1.0e-13)
595        type=3;
596#if 0
597      ClpSimplex * thisModel = static_cast<ClpSimplexOther *> (this)->dualOfModel(1.0,1.0);
598      if (thisModel) {
599        printf("Dual of model has %d rows and %d columns\n",
600               thisModel->numberRows(), thisModel->numberColumns());
601        thisModel->setOptimizationDirection(1.0);
602        Idiot info(*thisModel);
603        info.setStrategy(512 | info.getStrategy());
604        // Allow for scaling
605        info.setStrategy(32 | info.getStrategy());
606        info.setStartingWeight(1.0e3);
607        info.setReduceIterations(6);
608        info.crash(50, this->messageHandler(), this->messagesPointer(),false);
609        // make sure later that primal solution in correct place
610        // and has correct sign
611        abcModel2->setupDualValuesPass(thisModel->primalColumnSolution(),
612                                       thisModel->dualRowSolution(),type);
613        //thisModel->dual();
614        delete thisModel;
615      }
616#else
617      if (!solveType) {
618        this->dual(0);
619        abcModel2->setupDualValuesPass(this->dualRowSolution(),
620                                       this->primalColumnSolution(),type);
621      } else if (solveType==1) {
622        ifValuesPass=1;
623        abcModel2->setStateOfProblem(abcModel2->stateOfProblem() | VALUES_PASS);
624        Idiot info(*abcModel2);
625        info.setStrategy(512 | info.getStrategy());
626        // Allow for scaling
627        info.setStrategy(32 | info.getStrategy());
628        info.setStartingWeight(1.0e3);
629        info.setReduceIterations(6);
630        info.crash(200, abcModel2->messageHandler(), abcModel2->messagesPointer(),false);
631        //memcpy(abcModel2->primalColumnSolution(),this->primalColumnSolution(),
632        //     this->numberColumns()*sizeof(double));
633      }
634#endif
635    }
636    int numberCpu=this->abcState()&15;
637    if (numberCpu==9) {
638      numberCpu=1;
639#if ABC_PARALLEL==2
640#ifndef FAKE_CILK
641      if (number_cilk_workers>1)
642      numberCpu=CoinMin(2*number_cilk_workers,8);
643#endif
644#endif
645    } else if (numberCpu==10) {
646      // maximum
647      numberCpu=4;
648    } else if (numberCpu==10) {
649      // decide
650      if (abcModel2->getNumElements()<5000)
651        numberCpu=1;
652#if ABC_PARALLEL==2
653#ifndef FAKE_CILK
654      else if (number_cilk_workers>1)
655        numberCpu=CoinMin(2*number_cilk_workers,8);
656#endif
657#endif
658      else
659        numberCpu=1;
660    } else {
661#if ABC_PARALLEL==2
662#ifndef FAKE_CILK
663      char temp[3];
664      sprintf(temp,"%d",numberCpu);
665      __cilkrts_set_param("nworkers",temp);
666      printf("setting cilk workers to %d\n",numberCpu);
667      number_cilk_workers=numberCpu;
668
669#endif
670#endif
671    }
672    char line[200];
673#if ABC_PARALLEL
674#if ABC_PARALLEL==2
675#ifndef FAKE_CILK
676    if (!number_cilk_workers) {
677      number_cilk_workers=__cilkrts_get_nworkers();
678      sprintf(line,"%d cilk workers",number_cilk_workers);
679      handler_->message(CLP_GENERAL, messages_)
680        << line
681        << CoinMessageEol;
682    }
683#endif
684#endif
685    abcModel2->setParallelMode(numberCpu-1);
686#endif
687    //if (abcState()==3||abcState()==4) {
688    //abcModel2->setMoreSpecialOptions((131072*2)|abcModel2->moreSpecialOptions());
689    //}
690    //if (processTune>0&&processTune<8)
691    //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|131072*processTune);
692#if ABC_INSTRUMENT
693    double startTimeCpu=CoinCpuTime();
694    double startTimeElapsed=CoinGetTimeOfDay();
695#if ABC_INSTRUMENT>1
696    memset(abcPricing,0,sizeof(abcPricing));
697    memset(abcPricingDense,0,sizeof(abcPricing));
698    instrument_initialize(abcModel2->numberRows());
699#endif
700#endif
701    if (!solveType) {
702      abcModel2->ClpSimplex::doAbcDual();
703    } else if (solveType==1) {
704      int saveOptions=abcModel2->specialOptions();
705      //if (startUp==2)
706      //abcModel2->setSpecialOptions(8192|saveOptions);
707      abcModel2->ClpSimplex::doAbcPrimal(ifValuesPass);
708      abcModel2->setSpecialOptions(saveOptions);
709    }
710#if ABC_INSTRUMENT
711    if (solveType<2) {
712      double timeCpu=CoinCpuTime()-startTimeCpu;
713      double timeElapsed=CoinGetTimeOfDay()-startTimeElapsed;
714      sprintf(line,"Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
715              this->problemName().c_str(),this->numberRows(),this->numberColumns(),
716              this->getNumElements(),
717              timeCpu,timeElapsed,timeElapsed ? timeCpu/timeElapsed : 1.0,abcModel2->numberIterations());
718      handler_->message(CLP_GENERAL, messages_)
719        << line
720        << CoinMessageEol;
721#if ABC_INSTRUMENT>1
722      {
723        int n;
724        n=0;
725        for (int i=0;i<20;i++) 
726          n+= abcPricing[i];
727        printf("CCSparse pricing done %d times",n);
728        int n2=0;
729        for (int i=0;i<20;i++) 
730          n2+= abcPricingDense[i];
731        if (n2) 
732          printf(" and dense pricing done %d times\n",n2);
733        else
734          printf("\n");
735        n=0;
736        printf ("CCS");
737        for (int i=0;i<19;i++) {
738          if (abcPricing[i]) {
739            if (n==5) {
740              n=0;
741              printf("\nCCS");
742            }
743            n++;
744            printf("(%d els,%d times) ",i,abcPricing[i]);
745          }
746        }
747        if (abcPricing[19]) {
748          if (n==5) {
749            n=0;
750            printf("\nCCS");
751          }
752          n++;
753          printf("(>=19 els,%d times) ",abcPricing[19]);
754        }
755        if (n2) {
756          printf ("CCD");
757          for (int i=0;i<19;i++) {
758            if (abcPricingDense[i]) {
759              if (n==5) {
760                n=0;
761                printf("\nCCD");
762              }
763              n++;
764              int k1=(numberRows_/16)*i;;
765              int k2=CoinMin(numberRows_,k1+(numberRows_/16)-1);
766              printf("(%d-%d els,%d times) ",k1,k2,abcPricingDense[i]);
767            }
768          }
769        }
770        printf("\n");
771      }
772      instrument_print();
773#endif
774    }
775#endif
776    abcModel2->moveStatusToClp(this);
777#if 1
778    if (!problemStatus_) {
779      double offset;
780      CoinMemcpyN(this->objectiveAsObject()->gradient(this,
781                                                      this->primalColumnSolution(), offset, true),
782                  numberColumns_, this->dualColumnSolution());
783      this->clpMatrix()->transposeTimes(-1.0,
784                                        this->dualRowSolution(),
785                                        this->dualColumnSolution());
786      memset(this->primalRowSolution(), 0, numberRows_ * sizeof(double));
787      this->clpMatrix()->times(1.0, 
788                               this->primalColumnSolution(),
789                               this->primalRowSolution());
790      this->checkSolutionInternal();
791      if (sumDualInfeasibilities_>100.0*dualTolerance_) {
792#if CBC_USEFUL_PRINTING>0
793        printf("internal check on duals failed %d %g\n",
794               numberDualInfeasibilities_,sumDualInfeasibilities_);
795#endif
796      } else {
797        sumDualInfeasibilities_=0.0;
798        numberDualInfeasibilities_=0;
799      }
800      if (sumPrimalInfeasibilities_>100.0*primalTolerance_) {
801#if CBC_USEFUL_PRINTING>0
802        printf("internal check on primals failed %d %g\n",
803               numberPrimalInfeasibilities_,sumPrimalInfeasibilities_);
804#endif
805      } else {
806        sumPrimalInfeasibilities_=0.0;
807        numberPrimalInfeasibilities_=0;
808      }
809      problemStatus_=0;
810    }
811#endif
812    //ClpModel::stopPermanentArrays();
813    this->setSpecialOptions(this->specialOptions()&~65536);
814#if 0
815    this->writeBasis("a.bas",true);
816    for (int i=0;i<this->numberRows();i++)
817      printf("%d %g\n",i,this->dualRowSolution()[i]);
818    this->dual();
819    this->writeBasis("b.bas",true);
820    for (int i=0;i<this->numberRows();i++)
821      printf("%d %g\n",i,this->dualRowSolution()[i]);
822#endif
823    // switch off initialSolve flag
824    moreSpecialOptions_ &= ~16384;
825    //this->setNumberIterations(abcModel2->numberIterations()+this->numberIterations());
826    if (!keepAbc) {
827      delete abcModel2;
828      abcModel2=NULL;
829    }
830  }
831  return abcModel2;
832}
833#endif
834/** General solve algorithm which can do presolve
835    special options (bits)
836    1 - do not perturb
837    2 - do not scale
838    4 - use crash (default allslack in dual, idiot in primal)
839    8 - all slack basis in primal
840    16 - switch off interrupt handling
841    32 - do not try and make plus minus one matrix
842    64 - do not use sprint even if problem looks good
843 */
844int
845ClpSimplex::initialSolve(ClpSolve & options)
846{
847     ClpSolve::SolveType method = options.getSolveType();
848     //ClpSolve::SolveType originalMethod=method;
849     ClpSolve::PresolveType presolve = options.getPresolveType();
850     int saveMaxIterations = maximumIterations();
851     int finalStatus = -1;
852     int numberIterations = 0;
853     double time1 = CoinCpuTime();
854     double timeX = time1;
855     double time2=0.0;
856     ClpMatrixBase * saveMatrix = NULL;
857     ClpObjective * savedObjective = NULL;
858     int idiotOptions=0;
859     if(options.getSpecialOption(6))
860       idiotOptions=options.getExtraInfo(6)*32768;
861#ifdef CLP_USEFUL_PRINTOUT
862     debugInt[0]=numberRows();
863     debugInt[1]=numberColumns();
864     debugInt[2]=matrix()->getNumElements();
865#endif
866     if (!objective_ || !matrix_) {
867          // totally empty
868          handler_->message(CLP_EMPTY_PROBLEM, messages_)
869                    << 0
870                    << 0
871                    << 0
872                    << CoinMessageEol;
873          return -1;
874     } else if (!numberRows_ || !numberColumns_ || !getNumElements()) {
875          presolve = ClpSolve::presolveOff;
876     }
877     if (objective_->type() >= 2 && optimizationDirection_ == 0) {
878          // pretend linear
879          savedObjective = objective_;
880          // make up objective
881          double * obj = new double[numberColumns_];
882          for (int i = 0; i < numberColumns_; i++) {
883               double l = fabs(columnLower_[i]);
884               double u = fabs(columnUpper_[i]);
885               obj[i] = 0.0;
886               if (CoinMin(l, u) < 1.0e20) {
887                    if (l < u)
888                         obj[i] = 1.0 + randomNumberGenerator_.randomDouble() * 1.0e-2;
889                    else
890                         obj[i] = -1.0 - randomNumberGenerator_.randomDouble() * 1.0e-2;
891               }
892          }
893          objective_ = new ClpLinearObjective(obj, numberColumns_);
894          delete [] obj;
895     }
896     ClpSimplex * model2 = this;
897     bool interrupt = (options.getSpecialOption(2) == 0);
898     CoinSighandler_t saveSignal = static_cast<CoinSighandler_t> (0);
899     if (interrupt) {
900          currentModel = model2;
901          // register signal handler
902          saveSignal = signal(SIGINT, signal_handler);
903     }
904     // If no status array - set up basis
905     if (!status_)
906          allSlackBasis();
907     ClpPresolve * pinfo = new ClpPresolve();
908     pinfo->setSubstitution(options.substitution());
909     int presolveOptions = options.presolveActions();
910     bool presolveToFile = (presolveOptions & 0x40000000) != 0;
911     presolveOptions &= ~0x40000000;
912     if ((presolveOptions & 0xffff) != 0)
913          pinfo->setPresolveActions(presolveOptions);
914     // switch off singletons to slacks
915     //pinfo->setDoSingletonColumn(false); // done by bits
916     int printOptions = options.getSpecialOption(5);
917     if ((printOptions & 1) != 0)
918          pinfo->statistics();
919     double timePresolve = 0.0;
920     double timeIdiot = 0.0;
921     double timeCore = 0.0;
922     eventHandler()->event(ClpEventHandler::presolveStart);
923     int savePerturbation = perturbation_;
924     int saveScaling = scalingFlag_;
925#ifndef SLIM_CLP
926#ifndef NO_RTTI
927     if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
928          // network - switch off stuff
929          presolve = ClpSolve::presolveOff;
930     }
931#else
932     if (matrix_->type() == 11) {
933          // network - switch off stuff
934          presolve = ClpSolve::presolveOff;
935     }
936#endif
937#endif
938#ifndef CLPSOLVE_ACTIONS
939#define CLPSOLVE_ACTIONS 2
940#endif
941#if CLPSOLVE_ACTIONS
942     bool wasAutomatic = (method==ClpSolve::automatic);
943#endif
944     if (presolve != ClpSolve::presolveOff) {
945          bool costedSlacks = false;
946#if CLP_INHERIT_MODE>1
947          int numberPasses = 20;
948#else
949          int numberPasses = 5;
950#endif
951          if (presolve == ClpSolve::presolveNumber) {
952               numberPasses = options.getPresolvePasses();
953               presolve = ClpSolve::presolveOn;
954          } else if (presolve == ClpSolve::presolveNumberCost) {
955               numberPasses = options.getPresolvePasses();
956               presolve = ClpSolve::presolveOn;
957               costedSlacks = true;
958               // switch on singletons to slacks
959               pinfo->setDoSingletonColumn(true);
960               // gub stuff for testing
961               //pinfo->setDoGubrow(true);
962          }
963#ifndef CLP_NO_STD
964          if (presolveToFile) {
965               // PreSolve to file - not fully tested
966               printf("Presolving to file - presolve.save\n");
967               pinfo->presolvedModelToFile(*this, "presolve.save", dblParam_[ClpPresolveTolerance],
968                                          false, numberPasses);
969               model2 = this;
970          } else {
971#endif
972               model2 = pinfo->presolvedModel(*this, dblParam_[ClpPresolveTolerance],
973                                             false, numberPasses, true, costedSlacks);
974#ifndef CLP_NO_STD
975          }
976#endif
977          time2 = CoinCpuTime();
978          timePresolve = time2 - timeX;
979          handler_->message(CLP_INTERVAL_TIMING, messages_)
980                    << "Presolve" << timePresolve << time2 - time1
981                    << CoinMessageEol;
982          timeX = time2;
983          if (!model2) {
984               handler_->message(CLP_INFEASIBLE, messages_)
985                         << CoinMessageEol;
986               model2 = this;
987               eventHandler()->event(ClpEventHandler::presolveInfeasible);
988               problemStatus_ = pinfo->presolveStatus(); 
989               if (options.infeasibleReturn() || (moreSpecialOptions_ & 1) != 0) {
990                 delete pinfo;
991                    return -1;
992               }
993               presolve = ClpSolve::presolveOff;
994          } else {
995#if 0 //def ABC_INHERIT
996            {
997              AbcSimplex * abcModel2=new AbcSimplex(*model2);
998              delete model2;
999              model2=abcModel2;
1000              pinfo->setPresolvedModel(model2);
1001            }
1002#else
1003            //ClpModel::stopPermanentArrays();
1004            //setSpecialOptions(specialOptions()&~65536);
1005            // try setting tolerances up
1006#if CLPSOLVE_ACTIONS
1007            bool changeTolerances=wasAutomatic;
1008#if CLPSOLVE_ACTIONS>1
1009            changeTolerances=true;
1010#endif
1011            if (changeTolerances && model2 != this) {
1012#define CLP_NEW_TOLERANCE 1.0e-6
1013              if (model2->primalTolerance()==1.0e-7&&model2->dualTolerance()==1.0e-7) {
1014                model2->setPrimalTolerance(CLP_NEW_TOLERANCE);
1015                model2->setDualTolerance(CLP_NEW_TOLERANCE);
1016              }
1017            }
1018#endif
1019#endif
1020              model2->eventHandler()->setSimplex(model2);
1021              int rcode=model2->eventHandler()->event(ClpEventHandler::presolveSize);
1022              // see if too big or small
1023              if (rcode==2) {
1024                  delete model2;
1025                 delete pinfo;
1026                  return -2;
1027              } else if (rcode==3) {
1028                  delete model2;
1029                 delete pinfo;
1030                  return -3;
1031              }
1032          }
1033          model2->setMoreSpecialOptions(model2->moreSpecialOptions()&(~1024));
1034          model2->eventHandler()->setSimplex(model2);
1035          // We may be better off using original (but if dual leave because of bounds)
1036          if (presolve != ClpSolve::presolveOff &&
1037                    numberRows_ < 1.01 * model2->numberRows_ && numberColumns_ < 1.01 * model2->numberColumns_
1038                    && model2 != this) {
1039               if(method != ClpSolve::useDual ||
1040                         (numberRows_ == model2->numberRows_ && numberColumns_ == model2->numberColumns_)) {
1041                    delete model2;
1042                    model2 = this;
1043                    presolve = ClpSolve::presolveOff;
1044               }
1045          }
1046     }
1047#ifdef CLP_USEFUL_PRINTOUT
1048     debugInt[3]=model2->numberRows();
1049     debugInt[4]=model2->numberColumns();
1050     debugInt[5]=model2->matrix()->getNumElements();
1051     // analyze
1052     {
1053       double time1 =CoinCpuTime();
1054       int numberColumns = model2->numberColumns();
1055       const double * columnLower = model2->columnLower();
1056       const double * columnUpper = model2->columnUpper();
1057       int numberRows = model2->numberRows();
1058       const double * rowLower = model2->rowLower();
1059       const double * rowUpper = model2->rowUpper();
1060       const double * objective = model2->objective();
1061       CoinPackedMatrix * matrix = model2->matrix();
1062       CoinBigIndex numberElements = matrix->getNumElements();
1063       const int * columnLength = matrix->getVectorLengths();
1064       //const CoinBigIndex * columnStart = matrix->getVectorStarts();
1065       const double * elementByColumn = matrix->getElements();
1066       const int * row = matrix->getIndices();
1067       int * rowCount = new int [numberRows];
1068       memset(rowCount,0,numberRows*sizeof(int));
1069       int n=CoinMax (2*numberRows,numberElements);
1070       n= CoinMax(2*numberColumns,n);
1071       double * check = new double[n];
1072       memcpy(check,elementByColumn,numberElements*sizeof(double));
1073       for (int i=0;i<numberElements;i++) {
1074         check[i]=fabs(check[i]);
1075         rowCount[row[i]]++;
1076       }
1077       int largestIndex=0;
1078       for (int i=0;i<numberColumns;i++) {
1079         largestIndex=CoinMax(largestIndex,columnLength[i]);
1080       }
1081       debugInt[12]=largestIndex;
1082       largestIndex=0;
1083       for (int i=0;i<numberRows;i++) {
1084         largestIndex=CoinMax(largestIndex,rowCount[i]);
1085       }
1086       n=numberElements;
1087       delete [] rowCount;
1088       debugInt[11]=largestIndex;
1089       std::sort(check,check+n);
1090       debugDouble[4]=check[0];
1091       debugDouble[5]=check[n-1];
1092       int nAtOne=0;
1093       int nDifferent=0;
1094       double last=-COIN_DBL_MAX;
1095       for (int i=0;i<n;i++) {
1096         if (fabs(last-check[i])>1.0e-12) {
1097           nDifferent++;
1098           last=check[i];
1099         }
1100         if (check[i]==1.0)
1101           nAtOne++;
1102       }
1103       debugInt[10]=nDifferent;
1104       debugInt[15]=nAtOne;
1105       int nInf=0;
1106       int nZero=0;
1107       n=0;
1108       for (int i=0;i<numberRows;i++) {
1109         double value=fabs(rowLower[i]);
1110         if (!value)
1111           nZero++;
1112         else if (value!=COIN_DBL_MAX)
1113           check[n++]=value;
1114         else
1115           nInf++;
1116       }
1117       for (int i=0;i<numberRows;i++) {
1118         double value=fabs(rowUpper[i]);
1119         if (!value)
1120           nZero++;
1121         else if (value!=COIN_DBL_MAX)
1122           check[n++]=value;
1123         else
1124           nInf++;
1125       }
1126       debugInt[16]=nInf;
1127       debugInt[20]=nZero;
1128       if (n) {
1129         std::sort(check,check+n);
1130         debugDouble[0]=check[0];
1131         debugDouble[1]=check[n-1];
1132       } else {
1133         debugDouble[0]=0.0;
1134         debugDouble[1]=0.0;
1135       }
1136       nAtOne=0;
1137       nDifferent=0;
1138       last=-COIN_DBL_MAX;
1139       for (int i=0;i<n;i++) {
1140         if (fabs(last-check[i])>1.0e-12) {
1141           nDifferent++;
1142           last=check[i];
1143         }
1144         if (check[i]==1.0)
1145           nAtOne++;
1146       }
1147       debugInt[8]=nDifferent;
1148       debugInt[13]=nAtOne;
1149       nZero=0;
1150       n=0;
1151       for (int i=0;i<numberColumns;i++) {
1152         double value=fabs(objective[i]);
1153         if (value)
1154           check[n++]=value;
1155         else
1156           nZero++;
1157       }
1158       debugInt[21]=nZero;
1159       if (n) {
1160         std::sort(check,check+n);
1161         debugDouble[2]=check[0];
1162         debugDouble[3]=check[n-1];
1163       } else {
1164         debugDouble[2]=0.0;
1165         debugDouble[3]=0.0;
1166       }
1167       nAtOne=0;
1168       nDifferent=0;
1169       last=-COIN_DBL_MAX;
1170       for (int i=0;i<n;i++) {
1171         if (fabs(last-check[i])>1.0e-12) {
1172           nDifferent++;
1173           last=check[i];
1174         }
1175         if (check[i]==1.0)
1176           nAtOne++;
1177       }
1178       debugInt[9]=nDifferent;
1179       debugInt[14]=nAtOne;
1180       nInf=0;
1181       nZero=0;
1182       n=0;
1183       for (int i=0;i<numberColumns;i++) {
1184         double value=fabs(columnLower[i]);
1185         if (!value)
1186           nZero++;
1187         else if (value!=COIN_DBL_MAX)
1188           check[n++]=value;
1189         else
1190           nInf++;
1191       }
1192       for (int i=0;i<numberColumns;i++) {
1193         double value=fabs(columnUpper[i]);
1194         if (!value)
1195           nZero++;
1196         else if (value!=COIN_DBL_MAX)
1197           check[n++]=value;
1198         else
1199           nInf++;
1200       }
1201       debugInt[17]=nInf;
1202       double smallestColBound;
1203       double largestColBound;
1204       if (n) {
1205         std::sort(check,check+n);
1206         smallestColBound=check[0];
1207         largestColBound=check[n-1];
1208       } else {
1209         smallestColBound=0.0;
1210         largestColBound=0.0;
1211       }
1212       nAtOne=0;
1213       nDifferent=0;
1214       last=-COIN_DBL_MAX;
1215       for (int i=0;i<n;i++) {
1216         if (fabs(last-check[i])>1.0e-12) {
1217           nDifferent++;
1218           last=check[i];
1219         }
1220         if (check[i]==1.0)
1221           nAtOne++;
1222       }
1223       //debugInt[8]=nDifferent;
1224       //debugInt[13]=nAtOne;
1225       printf("BENCHMARK_STATS rhs %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
1226              debugInt[8],debugDouble[0],debugDouble[1],debugInt[13],debugInt[16],debugInt[20]);
1227       printf("BENCHMARK_STATS col %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
1228              nDifferent,smallestColBound,largestColBound,nAtOne,nInf,nZero);
1229       printf("BENCHMARK_STATS els %d different - %g -> %g (%d at one) - longest r,c %d,%d\n",
1230              debugInt[10],debugDouble[4],debugDouble[5],debugInt[15],
1231              debugInt[11],debugInt[12]);
1232       printf("BENCHMARK_STATS obj %d different - %g -> %g (%d at one, %d zero) - time %g\n",
1233              debugInt[9],debugDouble[2],debugDouble[3],debugInt[14],debugInt[21],
1234              CoinCpuTime()-time1);
1235       delete [] check;
1236     }
1237#endif
1238     if (interrupt)
1239          currentModel = model2;
1240     int saveMoreOptions = moreSpecialOptions_;
1241     // For below >0 overrides
1242     // 0 means no, -1 means maybe
1243     int doIdiot = 0;
1244     int doCrash = 0;
1245     int doSprint = 0;
1246     int doSlp = 0;
1247     int primalStartup = 1;
1248     model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve);
1249     int tryItSave = 0;
1250#if CLPSOLVE_ACTIONS
1251     if (method == ClpSolve::automatic )
1252       model2->moreSpecialOptions_ |= 8192; // stop switch over
1253#endif
1254     // switch to primal from automatic if just one cost entry
1255     if (method == ClpSolve::automatic && model2->numberColumns() > 5000
1256#ifndef CLPSOLVE_ACTIONS
1257         && (specialOptions_ & 1024) != 0
1258#endif
1259         ) {
1260          // look at original model for objective
1261          int numberColumns = model2->numberColumns();
1262          int numberRows = model2->numberRows();
1263          int numberColumnsOrig = this->numberColumns();
1264          const double * obj = this->objective();
1265          int nNon = 0;
1266          bool allOnes=true;
1267          for (int i = 0; i < numberColumnsOrig; i++) {
1268            if (obj[i]) {
1269              nNon++;
1270              if (fabs(obj[i])!=1.0)
1271                allOnes=false;
1272            }
1273          }
1274          if (nNon <= 1 || allOnes || 
1275              (options.getExtraInfo(1) > 0 && options.getSpecialOption(1)==2)) {
1276#ifdef COIN_DEVELOP
1277               printf("Forcing primal\n");
1278#endif
1279               method = ClpSolve::usePrimal;
1280#ifndef CLPSOLVE_ACTIONS
1281               tryItSave = (numberRows > 200 && numberColumns > 2000 &&
1282                            (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0)) ? 3 : 0;
1283#else
1284               if (numberRows > 200 && numberColumns > 2000 &&
1285                   numberColumns > 2 * numberRows) {
1286                 tryItSave= 3;
1287               } else {
1288                 // If rhs also rubbish then maybe
1289                 int numberRowsOrig = this->numberRows();
1290                 const double * rowLower = this->rowLower();
1291                 const double * rowUpper = this->rowUpper();
1292                 double last=COIN_DBL_MAX;
1293                 int nDifferent=0;
1294                 for (int i=0;i<numberRowsOrig;i++) {
1295                   double value=fabs(rowLower[i]);
1296                   if (value&&value!=COIN_DBL_MAX) {
1297                     if (value!=last) {
1298                       nDifferent++;
1299                       last=value;
1300                     }
1301                   }
1302                   value=fabs(rowUpper[i]);
1303                   if (value&&value!=COIN_DBL_MAX) {
1304                     if (value!=last) {
1305                       nDifferent++;
1306                       last=value;
1307                     }
1308                   }
1309                 }
1310                 if (nDifferent<2)
1311                   tryItSave=1;
1312               }
1313#endif
1314          }
1315     }
1316     if (method != ClpSolve::useDual && method != ClpSolve::useBarrier
1317         && method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe
1318               && method != ClpSolve::useBarrierNoCross) {
1319          switch (options.getSpecialOption(1)) {
1320          case 0:
1321               doIdiot = -1;
1322               doCrash = -1;
1323               doSprint = -1;
1324               break;
1325          case 1:
1326               doIdiot = 0;
1327               doCrash = 1;
1328               if (options.getExtraInfo(1) > 0)
1329                    doCrash = options.getExtraInfo(1);
1330               doSprint = 0;
1331               break;
1332          case 2:
1333               doIdiot = 1;
1334               if (options.getExtraInfo(1) > 0)
1335                    doIdiot = options.getExtraInfo(1);
1336               doCrash = 0;
1337               doSprint = 0;
1338               break;
1339          case 3:
1340               doIdiot = 0;
1341               doCrash = 0;
1342               doSprint = 1;
1343               break;
1344          case 4:
1345               doIdiot = 0;
1346               doCrash = 0;
1347               doSprint = 0;
1348               break;
1349          case 5:
1350               doIdiot = 0;
1351               doCrash = -1;
1352               doSprint = -1;
1353               break;
1354          case 6:
1355               doIdiot = -1;
1356               doCrash = -1;
1357               doSprint = 0;
1358               break;
1359          case 7:
1360               doIdiot = -1;
1361               doCrash = 0;
1362               doSprint = -1;
1363               break;
1364          case 8:
1365               doIdiot = -1;
1366               doCrash = 0;
1367               doSprint = 0;
1368               break;
1369          case 9:
1370               doIdiot = 0;
1371               doCrash = 0;
1372               doSprint = -1;
1373               break;
1374          case 10:
1375               doIdiot = 0;
1376               doCrash = 0;
1377               doSprint = 0;
1378               if (options.getExtraInfo(1) )
1379                    doSlp = options.getExtraInfo(1);
1380               break;
1381          case 11:
1382               doIdiot = 0;
1383               doCrash = 0;
1384               doSprint = 0;
1385               primalStartup = 0;
1386               break;
1387          default:
1388               abort();
1389          }
1390     } else if (method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe) {
1391          // Dual
1392          switch (options.getSpecialOption(0)) {
1393          case 0:
1394               doIdiot = 0;
1395               doCrash = 0;
1396               doSprint = 0;
1397               break;
1398          case 1:
1399               doIdiot = 0;
1400               doCrash = 1;
1401               if (options.getExtraInfo(0) > 0)
1402                    doCrash = options.getExtraInfo(0);
1403               doSprint = 0;
1404               break;
1405          case 2:
1406               doIdiot = -1;
1407               if (options.getExtraInfo(0) > 0)
1408                    doIdiot = options.getExtraInfo(0);
1409               doCrash = 0;
1410               doSprint = 0;
1411               break;
1412          default:
1413               abort();
1414          }
1415     } else {
1416       // decomposition
1417     }
1418#ifndef NO_RTTI
1419     ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objectiveAsObject()));
1420#else
1421     ClpQuadraticObjective * quadraticObj = NULL;
1422     if (objective_->type() == 2)
1423          quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1424#endif
1425     // If quadratic then primal or barrier or slp
1426     if (quadraticObj) {
1427          doSprint = 0;
1428          //doIdiot = 0;
1429          // off
1430          if (method == ClpSolve::useBarrier)
1431               method = ClpSolve::useBarrierNoCross;
1432          else if (method != ClpSolve::useBarrierNoCross)
1433               method = ClpSolve::usePrimal;
1434     }
1435#ifdef COIN_HAS_VOL
1436     // Save number of idiot
1437     int saveDoIdiot = doIdiot;
1438#endif
1439     // Just do this number of passes in Sprint
1440     int maxSprintPass = 100;
1441     // See if worth trying +- one matrix
1442     bool plusMinus = false;
1443     int numberElements = model2->getNumElements();
1444#ifndef SLIM_CLP
1445#ifndef NO_RTTI
1446     if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
1447          // network - switch off stuff
1448          doIdiot = 0;
1449          if (doSprint < 0)
1450               doSprint = 0;
1451     }
1452#else
1453     if (matrix_->type() == 11) {
1454          // network - switch off stuff
1455          doIdiot = 0;
1456          //doSprint=0;
1457     }
1458#endif
1459#endif
1460     int numberColumns = model2->numberColumns();
1461     int numberRows = model2->numberRows();
1462     // If not all slack basis - switch off all except sprint
1463     int numberRowsBasic = 0;
1464     int iRow;
1465     for (iRow = 0; iRow < numberRows; iRow++)
1466          if (model2->getRowStatus(iRow) == basic)
1467               numberRowsBasic++;
1468     if (numberRowsBasic < numberRows && objective_->type()<2) {
1469          doIdiot = 0;
1470          doCrash = 0;
1471          //doSprint=0;
1472     }
1473     if (options.getSpecialOption(3) == 0) {
1474          if(numberElements > 100000)
1475               plusMinus = true;
1476          if(numberElements > 10000 && (doIdiot || doSprint))
1477               plusMinus = true;
1478     } else if ((specialOptions_ & 1024) != 0) {
1479          plusMinus = true;
1480     }
1481#ifndef SLIM_CLP
1482     // Statistics (+1,-1, other) - used to decide on strategy if not +-1
1483     CoinBigIndex statistics[3] = { -1, 0, 0};
1484     if (plusMinus) {
1485          saveMatrix = model2->clpMatrix();
1486#ifndef NO_RTTI
1487          ClpPackedMatrix* clpMatrix =
1488               dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1489#else
1490          ClpPackedMatrix* clpMatrix = NULL;
1491          if (saveMatrix->type() == 1)
1492               clpMatrix =
1493                    static_cast< ClpPackedMatrix*>(saveMatrix);
1494#endif
1495          if (clpMatrix) {
1496               ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
1497               if (newMatrix->getIndices()) {
1498                  // CHECKME This test of specialOptions and the one above
1499                  // don't seem compatible.
1500#ifndef ABC_INHERIT
1501                    if ((specialOptions_ & 1024) == 0) {
1502                         model2->replaceMatrix(newMatrix);
1503                    } else {
1504#endif
1505                         // in integer (or abc) - just use for sprint/idiot
1506                         saveMatrix = NULL;
1507                         delete newMatrix;
1508#ifndef ABC_INHERIT
1509                    }
1510#endif
1511               } else {
1512                    handler_->message(CLP_MATRIX_CHANGE, messages_)
1513                              << "+- 1"
1514                              << CoinMessageEol;
1515                    CoinMemcpyN(newMatrix->startPositive(), 3, statistics);
1516                    saveMatrix = NULL;
1517                    plusMinus = false;
1518                    delete newMatrix;
1519               }
1520          } else {
1521               saveMatrix = NULL;
1522               plusMinus = false;
1523          }
1524     }
1525#endif
1526     if (this->factorizationFrequency() == 200) {
1527          // User did not touch preset
1528          model2->defaultFactorizationFrequency();
1529     } else if (model2 != this) {
1530          // make sure model2 has correct value
1531          model2->setFactorizationFrequency(this->factorizationFrequency());
1532     }
1533     if (method == ClpSolve::automatic) {
1534          if (doSprint == 0 && doIdiot == 0) {
1535               // off
1536               method = ClpSolve::useDual;
1537          } else {
1538               // only do primal if sprint or idiot
1539               if (doSprint > 0) {
1540                    method = ClpSolve::usePrimalorSprint;
1541               } else if (doIdiot > 0) {
1542                    method = ClpSolve::usePrimal;
1543               } else {
1544                    if (numberElements < 500000) {
1545                         // Small problem
1546                         if(numberRows * 10 > numberColumns || numberColumns < 6000
1547                                   || (numberRows * 20 > numberColumns && !plusMinus))
1548                              doSprint = 0; // switch off sprint
1549                    } else {
1550                         // larger problem
1551                         if(numberRows * 8 > numberColumns)
1552                              doSprint = 0; // switch off sprint
1553                    }
1554                    // switch off idiot or sprint if any free variable
1555                    // switch off sprint if very few with costs
1556                    // or great variation in cost
1557                    int iColumn;
1558                    const double * columnLower = model2->columnLower();
1559                    const double * columnUpper = model2->columnUpper();
1560                    const double * objective = model2->objective();
1561                    int nObj=0;
1562                    int nFree=0;
1563                    double smallestObj=COIN_DBL_MAX;
1564                    double largestObj=0.0;
1565                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1566                         if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
1567                           nFree++;
1568                         } else if (objective[iColumn]) {
1569                           nObj++;
1570                           smallestObj=CoinMin(smallestObj,objective[iColumn]);
1571                           largestObj=CoinMax(largestObj,objective[iColumn]);
1572                         }
1573                    }
1574                    if (nObj*10<numberColumns ||
1575                        smallestObj*10.0<largestObj)
1576                      doSprint=0;
1577                    if (nFree)
1578                      doIdiot=0;
1579                    if (nFree*10>numberRows)
1580                      doSprint=0;
1581                    int nPasses = 0;
1582                    // look at rhs
1583                    int iRow;
1584                    double largest = 0.0;
1585                    double smallest = 1.0e30;
1586                    double largestGap = 0.0;
1587                    int numberNotE = 0;
1588                    bool notInteger = false;
1589                    for (iRow = 0; iRow < numberRows; iRow++) {
1590                         double value1 = model2->rowLower_[iRow];
1591                         if (value1 && value1 > -1.0e31) {
1592                              largest = CoinMax(largest, fabs(value1));
1593                              smallest = CoinMin(smallest, fabs(value1));
1594                              if (fabs(value1 - floor(value1 + 0.5)) > 1.0e-8) {
1595                                   notInteger = true;
1596                                   break;
1597                              }
1598                         }
1599                         double value2 = model2->rowUpper_[iRow];
1600                         if (value2 && value2 < 1.0e31) {
1601                              largest = CoinMax(largest, fabs(value2));
1602                              smallest = CoinMin(smallest, fabs(value2));
1603                              if (fabs(value2 - floor(value2 + 0.5)) > 1.0e-8) {
1604                                   notInteger = true;
1605                                   break;
1606                              }
1607                         }
1608                         // CHECKME This next bit can't be right...
1609                         if (value2 > value1) {
1610                              numberNotE++;
1611                              //if (value2 > 1.0e31 || value1 < -1.0e31)
1612                              //   largestGap = COIN_DBL_MAX;
1613                              //else
1614                              //   largestGap = value2 - value1;
1615                         }
1616                    }
1617                    int tryIt = (numberRows > 200 && numberColumns > 2000 &&
1618                                 (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0))) ? 3 : 0;
1619                    tryItSave = tryIt;
1620                    if (numberRows < 1000 && numberColumns < 3000)
1621                         tryIt = 0;
1622                    if (notInteger)
1623                         tryIt = 0;
1624                    if (largest / smallest > 10 || (largest / smallest > 2.0 && largest > 50))
1625                         tryIt = 0;
1626                    if (tryIt) {
1627                         if (largest / smallest > 2.0) {
1628                              nPasses = 10 + numberColumns / 100000;
1629                              nPasses = CoinMin(nPasses, 50);
1630                              nPasses = CoinMax(nPasses, 15);
1631                              if (numberRows > 20000 && nPasses > 5) {
1632                                   // Might as well go for it
1633                                   nPasses = CoinMax(nPasses, 71);
1634                              } else if (numberRows > 2000 && nPasses > 5) {
1635                                   nPasses = CoinMax(nPasses, 50);
1636                              } else if (numberElements < 3 * numberColumns) {
1637                                   nPasses = CoinMin(nPasses, 10); // probably not worh it
1638                              }
1639                         } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) {
1640                              nPasses = 10 + numberColumns / 1000;
1641                              nPasses = CoinMin(nPasses, 100);
1642                              nPasses = CoinMax(nPasses, 30);
1643                              if (numberRows > 25000) {
1644                                   // Might as well go for it
1645                                   nPasses = CoinMax(nPasses, 71);
1646                              }
1647                              if (!largestGap)
1648                                   nPasses *= 2;
1649                         } else {
1650                              nPasses = 10 + numberColumns / 1000;
1651                              nPasses = CoinMax(nPasses, 100);
1652                              if (!largestGap)
1653                                   nPasses *= 2;
1654                              nPasses = CoinMin(nPasses, 200);
1655                         }
1656                    }
1657                    //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
1658                    //     numberRows,numberColumns,plusMinus ? 'Y' : 'N',
1659                    //     tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
1660                    //exit(0);
1661                    if (!tryIt || nPasses <= 5)
1662                         doIdiot = 0;
1663#if CLPSOLVE_ACTIONS
1664                    if (doIdiot&&doSprint<0&&wasAutomatic &&
1665                        20*model2->numberRows()>model2->numberColumns())
1666                      doSprint=0; // switch off sprint
1667#endif
1668                    if (doSprint) {
1669                         method = ClpSolve::usePrimalorSprint;
1670                    } else if (doIdiot) {
1671                         method = ClpSolve::usePrimal;
1672                    } else {
1673                         method = ClpSolve::useDual;
1674                    }
1675               }
1676          }
1677     }
1678     if (method == ClpSolve::tryBenders) {
1679       // Now build model
1680       int lengthNames=model2->lengthNames();
1681       model2->setLengthNames(0);
1682       CoinModel * build = model2->createCoinModel();
1683       model2->setLengthNames(lengthNames);
1684       CoinStructuredModel benders;
1685       build->convertMatrix();
1686       int numberBlocks = options.independentOption(0);
1687       benders.setMessageHandler(handler_);
1688       numberBlocks=benders.decompose(*build,2,numberBlocks,NULL);
1689       delete build;
1690       //exit(0);
1691       if (numberBlocks) {
1692         options.setIndependentOption(1,1); // don't do final clean up
1693         model2->solveBenders(&benders,options);
1694         //move solution
1695         method=ClpSolve::notImplemented;
1696         time2 = CoinCpuTime();
1697         timeCore = time2 - timeX;
1698         handler_->message(CLP_INTERVAL_TIMING, messages_)
1699           << "Crossover" << timeCore << time2 - time1
1700           << CoinMessageEol;
1701         timeX = time2;
1702       } else {
1703         printf("No structure\n");
1704         method=ClpSolve::useDual;
1705       }
1706     } else if (method == ClpSolve::tryDantzigWolfe) {
1707         abort();
1708     }
1709     if (method == ClpSolve::usePrimalorSprint) {
1710          if (doSprint < 0) {
1711               if (numberElements < 500000) {
1712                    // Small problem
1713                    if(numberRows * 10 > numberColumns || numberColumns < 6000
1714                              || (numberRows * 20 > numberColumns && !plusMinus))
1715                         method = ClpSolve::usePrimal; // switch off sprint
1716               } else {
1717                    // larger problem
1718                    if(numberRows * 8 > numberColumns)
1719                         method = ClpSolve::usePrimal; // switch off sprint
1720                    // but make lightweight
1721                    if(numberRows * 10 > numberColumns || numberColumns < 6000
1722                              || (numberRows * 20 > numberColumns && !plusMinus))
1723                         maxSprintPass = 10;
1724               }
1725          } else if (doSprint == 0) {
1726               method = ClpSolve::usePrimal; // switch off sprint
1727          }
1728     }
1729     if (method == ClpSolve::useDual) {
1730#ifdef CLP_USEFUL_PRINTOUT
1731       debugInt[6]=1;
1732#endif
1733          double * saveLower = NULL;
1734          double * saveUpper = NULL;
1735          if (presolve == ClpSolve::presolveOn) {
1736               int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0);
1737               if (numberInfeasibilities) {
1738                    handler_->message(CLP_INFEASIBLE, messages_)
1739                              << CoinMessageEol;
1740                    delete model2;
1741                    model2 = this;
1742                    presolve = ClpSolve::presolveOff;
1743               }
1744          } else if (numberRows_ + numberColumns_ > 5000) {
1745               // do anyway
1746               saveLower = new double[numberRows_+numberColumns_];
1747               CoinMemcpyN(model2->columnLower(), numberColumns_, saveLower);
1748               CoinMemcpyN(model2->rowLower(), numberRows_, saveLower + numberColumns_);
1749               saveUpper = new double[numberRows_+numberColumns_];
1750               CoinMemcpyN(model2->columnUpper(), numberColumns_, saveUpper);
1751               CoinMemcpyN(model2->rowUpper(), numberRows_, saveUpper + numberColumns_);
1752               int numberInfeasibilities = model2->tightenPrimalBounds();
1753               if (numberInfeasibilities) {
1754                    handler_->message(CLP_INFEASIBLE, messages_)
1755                              << CoinMessageEol;
1756                    CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
1757                    CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
1758                    delete [] saveLower;
1759                    saveLower = NULL;
1760                    CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
1761                    CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
1762                    delete [] saveUpper;
1763                    saveUpper = NULL;
1764               }
1765          }
1766#ifndef COIN_HAS_VOL
1767          // switch off idiot and volume for now
1768          doIdiot = 0;
1769#endif
1770          // pick up number passes
1771          int nPasses = 0;
1772          int numberNotE = 0;
1773#ifndef SLIM_CLP
1774          if ((doIdiot < 0 && plusMinus) || doIdiot > 0) {
1775               // See if candidate for idiot
1776               nPasses = 0;
1777               Idiot info(*model2);
1778               info.setStrategy(idiotOptions | info.getStrategy());
1779               // Get average number of elements per column
1780               double ratio  = static_cast<double> (numberElements) / static_cast<double> (numberColumns);
1781               // look at rhs
1782               int iRow;
1783               double largest = 0.0;
1784               double smallest = 1.0e30;
1785               for (iRow = 0; iRow < numberRows; iRow++) {
1786                    double value1 = model2->rowLower_[iRow];
1787                    if (value1 && value1 > -1.0e31) {
1788                         largest = CoinMax(largest, fabs(value1));
1789                         smallest = CoinMin(smallest, fabs(value1));
1790                    }
1791                    double value2 = model2->rowUpper_[iRow];
1792                    if (value2 && value2 < 1.0e31) {
1793                         largest = CoinMax(largest, fabs(value2));
1794                         smallest = CoinMin(smallest, fabs(value2));
1795                    }
1796                    if (value2 > value1) {
1797                         numberNotE++;
1798                    }
1799               }
1800               if (doIdiot < 0) {
1801                    if (numberRows > 200 && numberColumns > 5000 && ratio >= 3.0 &&
1802                              largest / smallest < 1.1 && !numberNotE) {
1803                         nPasses = 71;
1804                    }
1805               }
1806               if (doIdiot > 0) {
1807                    nPasses = CoinMax(nPasses, doIdiot);
1808                    if (nPasses > 70) {
1809                         info.setStartingWeight(1.0e3);
1810                         info.setDropEnoughFeasibility(0.01);
1811                    }
1812               }
1813               if (nPasses > 20) {
1814#ifdef COIN_HAS_VOL
1815                    int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
1816                    if (!returnCode) {
1817                         time2 = CoinCpuTime();
1818                         timeIdiot = time2 - timeX;
1819                         handler_->message(CLP_INTERVAL_TIMING, messages_)
1820                                   << "Idiot Crash" << timeIdiot << time2 - time1
1821                                   << CoinMessageEol;
1822                         timeX = time2;
1823                    } else {
1824                         nPasses = 0;
1825                    }
1826#else
1827                    nPasses = 0;
1828#endif
1829               } else {
1830                    nPasses = 0;
1831               }
1832          }
1833#endif
1834          if (doCrash) {
1835#ifdef ABC_INHERIT
1836            if (!model2->abcState()) {
1837#endif
1838               switch(doCrash) {
1839                    // standard
1840               case 1:
1841                    model2->crash(1000, 1);
1842                    break;
1843                    // As in paper by Solow and Halim (approx)
1844               case 2:
1845               case 3:
1846                    model2->crash(model2->dualBound(), 0);
1847                    break;
1848                    // Just put free in basis
1849               case 4:
1850                    model2->crash(0.0, 3);
1851                    break;
1852               }
1853#ifdef ABC_INHERIT
1854            } else if (doCrash>=0) {
1855               model2->setAbcState(model2->abcState()|256*doCrash);
1856            }
1857#endif
1858          }
1859          if (!nPasses) {
1860               int saveOptions = model2->specialOptions();
1861               if (model2->numberRows() > 100)
1862                    model2->setSpecialOptions(saveOptions | 64); // go as far as possible
1863               //int numberRows = model2->numberRows();
1864               //int numberColumns = model2->numberColumns();
1865               if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1866                    // See if original wanted vector
1867                    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1868                    ClpMatrixBase * matrix = model2->clpMatrix();
1869                    if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
1870                         ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1871                         clpMatrix->makeSpecialColumnCopy();
1872                         //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
1873                         model2->dual(0);
1874                         clpMatrix->releaseSpecialColumnCopy();
1875                    } else {
1876#ifndef ABC_INHERIT
1877                      model2->dual(0);
1878#else
1879                      model2->dealWithAbc(0,0,interrupt);
1880#endif
1881                    }
1882               } else {
1883                    model2->dual(0);
1884               }
1885          } else if (!numberNotE && 0) {
1886               // E so we can do in another way
1887               double * pi = model2->dualRowSolution();
1888               int i;
1889               int numberColumns = model2->numberColumns();
1890               int numberRows = model2->numberRows();
1891               double * saveObj = new double[numberColumns];
1892               CoinMemcpyN(model2->objective(), numberColumns, saveObj);
1893               CoinMemcpyN(model2->objective(),
1894                           numberColumns, model2->dualColumnSolution());
1895               model2->clpMatrix()->transposeTimes(-1.0, pi, model2->dualColumnSolution());
1896               CoinMemcpyN(model2->dualColumnSolution(),
1897                           numberColumns, model2->objective());
1898               const double * rowsol = model2->primalRowSolution();
1899               double offset = 0.0;
1900               for (i = 0; i < numberRows; i++) {
1901                    offset += pi[i] * rowsol[i];
1902               }
1903               double value2;
1904               model2->getDblParam(ClpObjOffset, value2);
1905               //printf("Offset %g %g\n",offset,value2);
1906               model2->setDblParam(ClpObjOffset, value2 - offset);
1907               model2->setPerturbation(51);
1908               //model2->setRowObjective(pi);
1909               // zero out pi
1910               //memset(pi,0,numberRows*sizeof(double));
1911               // Could put some in basis - only partially tested
1912               model2->allSlackBasis();
1913               //model2->factorization()->maximumPivots(200);
1914               //model2->setLogLevel(63);
1915               // solve
1916               model2->dual(0);
1917               model2->setDblParam(ClpObjOffset, value2);
1918               CoinMemcpyN(saveObj, numberColumns, model2->objective());
1919               // zero out pi
1920               //memset(pi,0,numberRows*sizeof(double));
1921               //model2->setRowObjective(pi);
1922               delete [] saveObj;
1923               //model2->dual(0);
1924               model2->setPerturbation(50);
1925               model2->primal();
1926          } else {
1927               // solve
1928               model2->setPerturbation(100);
1929               model2->dual(2);
1930               model2->setPerturbation(50);
1931               model2->dual(0);
1932          }
1933          if (saveLower) {
1934               CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
1935               CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
1936               delete [] saveLower;
1937               saveLower = NULL;
1938               CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
1939               CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
1940               delete [] saveUpper;
1941               saveUpper = NULL;
1942          }
1943          time2 = CoinCpuTime();
1944          timeCore = time2 - timeX;
1945          handler_->message(CLP_INTERVAL_TIMING, messages_)
1946                    << "Dual" << timeCore << time2 - time1
1947                    << CoinMessageEol;
1948          timeX = time2;
1949     } else if (method == ClpSolve::usePrimal) {
1950#ifdef CLP_USEFUL_PRINTOUT
1951       debugInt[6]=2;
1952#endif
1953#ifndef SLIM_CLP
1954          if (doIdiot) {
1955               int nPasses = 0;
1956               Idiot info(*model2);
1957               info.setStrategy(idiotOptions | info.getStrategy());
1958               // Get average number of elements per column
1959               double ratio  = static_cast<double> (numberElements) / static_cast<double> (numberColumns);
1960               // look at rhs
1961               int iRow;
1962               double largest = 0.0;
1963               double smallest = 1.0e30;
1964               double largestGap = 0.0;
1965               int numberNotE = 0;
1966               for (iRow = 0; iRow < numberRows; iRow++) {
1967                    double value1 = model2->rowLower_[iRow];
1968                    if (value1 && value1 > -1.0e31) {
1969                         largest = CoinMax(largest, fabs(value1));
1970                         smallest = CoinMin(smallest, fabs(value1));
1971                    }
1972                    double value2 = model2->rowUpper_[iRow];
1973                    if (value2 && value2 < 1.0e31) {
1974                         largest = CoinMax(largest, fabs(value2));
1975                         smallest = CoinMin(smallest, fabs(value2));
1976                    }
1977                    if (value2 > value1) {
1978                         numberNotE++;
1979                         if (value2 > 1.0e31 || value1 < -1.0e31)
1980                              largestGap = COIN_DBL_MAX;
1981                         else
1982                              largestGap = value2 - value1;
1983                    }
1984               }
1985               bool increaseSprint = plusMinus;
1986               if ((specialOptions_ & 1024) != 0)
1987                    increaseSprint = false;
1988               if (!plusMinus) {
1989                    // If 90% +- 1 then go for sprint
1990                    if (statistics[0] >= 0 && 10 * statistics[2] < statistics[0] + statistics[1])
1991                         increaseSprint = true;
1992               }
1993               int tryIt = tryItSave;
1994               if (numberRows < 1000 && numberColumns < 3000)
1995                    tryIt = 0;
1996               if (tryIt) {
1997                    if (increaseSprint) {
1998                         info.setStartingWeight(1.0e3);
1999                         info.setReduceIterations(6);
2000                         // also be more lenient on infeasibilities
2001                         info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
2002                         info.setDropEnoughWeighted(-2.0);
2003                         if (largest / smallest > 2.0) {
2004                              nPasses = 10 + numberColumns / 100000;
2005                              nPasses = CoinMin(nPasses, 50);
2006                              nPasses = CoinMax(nPasses, 15);
2007                              if (numberRows > 20000 && nPasses > 5) {
2008                                   // Might as well go for it
2009                                   nPasses = CoinMax(nPasses, 71);
2010                              } else if (numberRows > 2000 && nPasses > 5) {
2011                                   nPasses = CoinMax(nPasses, 50);
2012                              } else if (numberElements < 3 * numberColumns) {
2013                                   nPasses = CoinMin(nPasses, 10); // probably not worh it
2014                                   if (doIdiot < 0)
2015                                        info.setLightweight(1); // say lightweight idiot
2016                              } else {
2017                                   if (doIdiot < 0)
2018                                        info.setLightweight(1); // say lightweight idiot
2019                              }
2020                         } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) {
2021                              nPasses = 10 + numberColumns / 1000;
2022                              nPasses = CoinMin(nPasses, 100);
2023                              nPasses = CoinMax(nPasses, 30);
2024                              if (numberRows > 25000) {
2025                                   // Might as well go for it
2026                                   nPasses = CoinMax(nPasses, 71);
2027                              }
2028                              if (!largestGap)
2029                                   nPasses *= 2;
2030                         } else {
2031                              nPasses = 10 + numberColumns / 1000;
2032                              nPasses = CoinMin(nPasses, 200);
2033                              nPasses = CoinMax(nPasses, 100);
2034                              info.setStartingWeight(1.0e-1);
2035                              info.setReduceIterations(6);
2036                              if (!largestGap)
2037                                   nPasses *= 2;
2038                              //info.setFeasibilityTolerance(1.0e-7);
2039                         }
2040                         // If few passes - don't bother
2041                         if (nPasses <= 5 && !plusMinus)
2042                              nPasses = 0;
2043                    } else {
2044                         if (doIdiot < 0)
2045                              info.setLightweight(1); // say lightweight idiot
2046                         if (largest / smallest > 1.01 || numberNotE || statistics[2] > statistics[0] + statistics[1]) {
2047                              if (numberRows > 25000 || numberColumns > 5 * numberRows) {
2048                                   nPasses = 50;
2049                              } else if (numberColumns > 4 * numberRows) {
2050                                   nPasses = 20;
2051                              } else {
2052                                   nPasses = 5;
2053                              }
2054                         } else {
2055                              if (numberRows > 25000 || numberColumns > 5 * numberRows) {
2056                                   nPasses = 50;
2057                                   info.setLightweight(0); // say not lightweight idiot
2058                              } else if (numberColumns > 4 * numberRows) {
2059                                   nPasses = 20;
2060                              } else {
2061                                   nPasses = 15;
2062                              }
2063                         }
2064                         if (ratio < 3.0) {
2065                              nPasses = static_cast<int> (ratio * static_cast<double> (nPasses) / 4.0); // probably not worth it
2066                         } else {
2067                              nPasses = CoinMax(nPasses, 5);
2068                         }
2069                         if (numberRows > 25000 && nPasses > 5) {
2070                              // Might as well go for it
2071                              nPasses = CoinMax(nPasses, 71);
2072                         } else if (increaseSprint) {
2073                              nPasses *= 2;
2074                              nPasses = CoinMin(nPasses, 71);
2075                         } else if (nPasses == 5 && ratio > 5.0) {
2076                              nPasses = static_cast<int> (static_cast<double> (nPasses) * (ratio / 5.0)); // increase if lots of elements per column
2077                         }
2078                         if (nPasses <= 5 && !plusMinus)
2079                              nPasses = 0;
2080                         //info.setStartingWeight(1.0e-1);
2081                    }
2082                    if (tryIt==1) {
2083                      idiotOptions |= 262144;
2084                      info.setStrategy(idiotOptions | info.getStrategy());
2085                      //model2->setSpecialOptions(model2->specialOptions()
2086                      //                        |8388608);
2087                    }
2088               }
2089               if (doIdiot > 0) {
2090                    // pick up number passes
2091                    nPasses = options.getExtraInfo(1) % 1000000;
2092#ifdef COIN_HAS_VOL
2093                    int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
2094                    nPasses=0;
2095                    if (!returnCode) {
2096                      time2 = CoinCpuTime();
2097                      timeIdiot = time2 - timeX;
2098                      handler_->message(CLP_INTERVAL_TIMING, messages_)
2099                        << "Idiot Crash" << timeIdiot << time2 - time1
2100                        << CoinMessageEol;
2101                      timeX = time2;
2102                    }
2103#endif
2104#ifdef CLP_USEFUL_PRINTOUT
2105                    debugInt[6]=4;
2106                    debugInt[7]=nPasses;
2107#endif
2108                    if (nPasses > 70) {
2109                         info.setStartingWeight(1.0e3);
2110                         info.setReduceIterations(6);
2111                         //if (nPasses > 200)
2112                         //info.setFeasibilityTolerance(1.0e-9);
2113                         //if (nPasses > 1900)
2114                         //info.setWeightFactor(0.93);
2115                         if (nPasses > 900) {
2116                           double reductions=nPasses/6.0;
2117                           if (nPasses<5000) {
2118                             reductions /= 12.0;
2119                           } else {
2120                             reductions /= 13.0;
2121                             info.setStartingWeight(1.0e4);
2122                           }
2123                           double ratio=1.0/std::pow(10.0,(1.0/reductions));
2124                           printf("%d passes reduction factor %g\n",nPasses,ratio);
2125                           info.setWeightFactor(ratio);
2126                         } else if (nPasses > 500) {
2127                           info.setWeightFactor(0.7);
2128                         } else if (nPasses > 200) {
2129                           info.setWeightFactor(0.5);
2130                         }
2131                         if (maximumIterations()<nPasses) {
2132                           printf("Presuming maximumIterations is just for Idiot\n");
2133                           nPasses=maximumIterations();
2134                           setMaximumIterations(COIN_INT_MAX);
2135                           model2->setMaximumIterations(COIN_INT_MAX);
2136                         }
2137                         if (nPasses >= 10000&&nPasses<100000) {
2138                              int k = nPasses % 100;
2139                              nPasses /= 200;
2140                              info.setReduceIterations(3);
2141                              if (k)
2142                                   info.setStartingWeight(1.0e2);
2143                         }
2144                         // also be more lenient on infeasibilities
2145                         info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
2146                         info.setDropEnoughWeighted(-2.0);
2147                    } else if (nPasses >= 50) {
2148                         info.setStartingWeight(1.0e3);
2149                         //info.setReduceIterations(6);
2150                    }
2151                    // For experimenting
2152                    if (nPasses < 70 && (nPasses % 10) > 0 && (nPasses % 10) < 4) {
2153                         info.setStartingWeight(1.0e3);
2154                         info.setLightweight(nPasses % 10); // special testing
2155#ifdef COIN_DEVELOP
2156                         printf("warning - odd lightweight %d\n", nPasses % 10);
2157                         //info.setReduceIterations(6);
2158#endif
2159                    }
2160               }
2161               if (options.getExtraInfo(1) > 1000000)
2162                    nPasses += 1000000;
2163               if (nPasses) {
2164                    doCrash = 0;
2165#if 0
2166                    double * solution = model2->primalColumnSolution();
2167                    int iColumn;
2168                    double * saveLower = new double[numberColumns];
2169                    CoinMemcpyN(model2->columnLower(), numberColumns, saveLower);
2170                    double * saveUpper = new double[numberColumns];
2171                    CoinMemcpyN(model2->columnUpper(), numberColumns, saveUpper);
2172                    printf("doing tighten before idiot\n");
2173                    model2->tightenPrimalBounds();
2174                    // Move solution
2175                    double * columnLower = model2->columnLower();
2176                    double * columnUpper = model2->columnUpper();
2177                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2178                         if (columnLower[iColumn] > 0.0)
2179                              solution[iColumn] = columnLower[iColumn];
2180                         else if (columnUpper[iColumn] < 0.0)
2181                              solution[iColumn] = columnUpper[iColumn];
2182                         else
2183                              solution[iColumn] = 0.0;
2184                    }
2185                    CoinMemcpyN(saveLower, numberColumns, columnLower);
2186                    CoinMemcpyN(saveUpper, numberColumns, columnUpper);
2187                    delete [] saveLower;
2188                    delete [] saveUpper;
2189#else
2190                    // Allow for crossover
2191                    //#define LACI_TRY
2192#ifndef LACI_TRY
2193                    //if (doIdiot>0)
2194#if 0 //def ABC_INHERIT
2195                    if (!model2->abcState())
2196#endif
2197                    info.setStrategy(512 | info.getStrategy());
2198#endif
2199                    // Allow for scaling
2200                    info.setStrategy(32 | info.getStrategy());
2201                    int saveScalingFlag=model2->scalingFlag();
2202                    bool linearObjective = objective_->type()<2;
2203                    if (!linearObjective)
2204                      model2->scaling(0);
2205#define CLP_DUAL_IDIOT
2206#ifdef CLP_DUAL_IDIOT
2207                    bool doubleIdiot=false;
2208                    if (nPasses==99&&linearObjective)
2209                      doubleIdiot=true;
2210                    if (doubleIdiot) {
2211                      ClpSimplex * dualModel2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(1.0,1.0);
2212                      if (dualModel2) {
2213                        printf("Dual of model has %d rows and %d columns\n",
2214                               dualModel2->numberRows(), dualModel2->numberColumns());
2215                        dualModel2->setOptimizationDirection(1.0);
2216                        Idiot infoDual(info);
2217                        infoDual.setModel(dualModel2);
2218#if 0
2219                        info.setStrategy(512 | info.getStrategy());
2220                        // Allow for scaling
2221                        info.setStrategy(32 | info.getStrategy());
2222                        info.setStartingWeight(1.0e3);
2223                        info.setReduceIterations(6);
2224#endif
2225                        info.crash(nPasses, model2->messageHandler(),
2226                                   model2->messagesPointer(),false);
2227                        infoDual.crash(nPasses, model2->messageHandler(),
2228                                   model2->messagesPointer(),false);
2229                        // two copies of solutions
2230                        ClpSimplex temp(*model2);
2231#if 0
2232                        static_cast<ClpSimplexOther *> (&temp)->restoreFromDual(dualModel2);
2233#else
2234                        // move duals and just copy primal
2235                        memcpy(temp.dualRowSolution(),dualModel2->primalColumnSolution(),
2236                               numberRows*sizeof(double));
2237                        memcpy(temp.primalColumnSolution(),model2->primalColumnSolution(),
2238                               numberColumns*sizeof(double));
2239#endif
2240                        delete dualModel2;
2241                        int numberRows=model2->numberRows();
2242                        int numberColumns=model2->numberColumns();
2243                        ClpSimplex * tempModel[2];
2244                        tempModel[0]=model2;
2245                        tempModel[1]=&temp;
2246                        const double * primalColumn[2];
2247                        const double * dualRow[2];
2248                        double * dualColumn[2];
2249                        double * primalRow[2];
2250                        for (int i=0;i<2;i++) {
2251                          primalColumn[i]=tempModel[i]->primalColumnSolution();
2252                          dualRow[i]=tempModel[i]->dualRowSolution();
2253                          dualColumn[i]=tempModel[i]->dualColumnSolution();
2254                          primalRow[i]=tempModel[i]->primalRowSolution();
2255                          memcpy(dualColumn[i],model2->objective(),
2256                                 numberColumns*sizeof(double));
2257                          memset(primalRow[i], 0, numberRows * sizeof(double));
2258                          tempModel[i]->clpMatrix()->transposeTimes(-1.0,
2259                                                              dualRow[i],
2260                                                              dualColumn[i]);
2261                          tempModel[i]->clpMatrix()->times(1.0,
2262                                                     primalColumn[i],
2263                                                     primalRow[i]);
2264                          tempModel[i]->checkSolutionInternal();
2265                          printf("model %d - dual inf %g primal inf %g\n",
2266                                 i,tempModel[i]->sumDualInfeasibilities(),
2267                                 tempModel[i]->sumPrimalInfeasibilities());
2268                        }
2269                        printf("What now\n");
2270                      } else {
2271                        doubleIdiot=false;
2272                      }
2273                    }
2274                    if (!doubleIdiot)
2275                      info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),
2276                                (objective_->type() <2));
2277#else
2278                    info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),(objective_->type() <2));
2279#endif
2280                      model2->scaling(saveScalingFlag);
2281#endif
2282                    time2 = CoinCpuTime();
2283                    timeIdiot = time2 - timeX;
2284                    handler_->message(CLP_INTERVAL_TIMING, messages_)
2285                              << "Idiot Crash" << timeIdiot << time2 - time1
2286                              << CoinMessageEol;
2287                    timeX = time2;
2288                    if (nPasses>100000&&nPasses<100500) {
2289                      // make sure no status left
2290                      model2->createStatus();
2291                      // solve
2292                      if (model2->factorizationFrequency() == 200) {
2293                        // User did not touch preset
2294                        model2->defaultFactorizationFrequency();
2295                      }
2296                      //int numberRows = model2->numberRows();
2297                      int numberColumns = model2->numberColumns();
2298                      // save duals
2299                      //double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows);
2300                      // for moment this only works on nug etc (i.e. all ==)
2301                      // needs row objective
2302                      double * saveObj = CoinCopyOfArray(model2->objective(),numberColumns);
2303                      double * pi = model2->dualRowSolution();
2304                      model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective());
2305                      // just primal values pass
2306                      double saveScale = model2->objectiveScale();
2307                      model2->setObjectiveScale(1.0e-3);
2308                      model2->primal(2);
2309                      model2->writeMps("xx.mps");
2310                      double * solution = model2->primalColumnSolution();
2311                      double * upper = model2->columnUpper();
2312                      for (int i=0;i<numberColumns;i++) {
2313                        if (solution[i]<100.0)
2314                          upper[i]=1000.0;
2315                      }
2316                      model2->setProblemStatus(-1);
2317                      model2->setObjectiveScale(saveScale);
2318#ifdef ABC_INHERIT
2319                      AbcSimplex * abcModel2=new AbcSimplex(*model2);
2320                      if (interrupt)
2321                        currentAbcModel = abcModel2;
2322                      if (abcSimplex_) {
2323                        // move factorization stuff
2324                        abcModel2->factorization()->synchronize(model2->factorization(),abcModel2);
2325                      }
2326                      abcModel2->startPermanentArrays();
2327                      abcModel2->setAbcState(CLP_ABC_WANTED);
2328#if ABC_PARALLEL
2329                      int parallelMode=1;
2330                      printf("Parallel mode %d\n",parallelMode);
2331                      abcModel2->setParallelMode(parallelMode);
2332#endif
2333                      //if (processTune>0&&processTune<8)
2334                      //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|65536*processTune);
2335                      abcModel2->doAbcDual();
2336                      abcModel2->moveStatusToClp(model2);
2337                      //ClpModel::stopPermanentArrays();
2338                      model2->setSpecialOptions(model2->specialOptions()&~65536);
2339                      //model2->dual();
2340                      //model2->setNumberIterations(abcModel2->numberIterations()+model2->numberIterations());
2341                      delete abcModel2;
2342#endif
2343                      memcpy(model2->objective(),saveObj,numberColumns*sizeof(double));
2344                      //delete [] saveDuals;
2345                      delete [] saveObj;
2346                      model2->dual(2);
2347                    } // end dubious idiot
2348               }
2349          }
2350#endif
2351          // ?
2352          if (doCrash) {
2353               switch(doCrash) {
2354                    // standard
2355               case 1:
2356                    model2->crash(1000, 1);
2357                    break;
2358                    // As in paper by Solow and Halim (approx)
2359               case 2:
2360                    model2->crash(model2->dualBound(), 0);
2361                    break;
2362                    // My take on it
2363               case 3:
2364                    model2->crash(model2->dualBound(), -1);
2365                    break;
2366                    // Just put free in basis
2367               case 4:
2368                    model2->crash(0.0, 3);
2369                    break;
2370               }
2371          }
2372#ifndef SLIM_CLP
2373          if (doSlp && objective_->type() == 2) {
2374               model2->nonlinearSLP(doSlp, 1.0e-5);
2375          }
2376#endif
2377#ifndef LACI_TRY
2378          if (options.getSpecialOption(1) != 2 ||
2379                    options.getExtraInfo(1) < 1000000) {
2380               if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
2381                    // See if original wanted vector
2382                    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
2383                    ClpMatrixBase * matrix = model2->clpMatrix();
2384                    if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
2385                         ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
2386                         clpMatrix->makeSpecialColumnCopy();
2387                         //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
2388                         model2->primal(primalStartup);
2389                         clpMatrix->releaseSpecialColumnCopy();
2390                    } else {
2391#ifndef ABC_INHERIT
2392                        model2->primal(primalStartup);
2393#else
2394                        model2->dealWithAbc(1,primalStartup,interrupt);
2395#endif
2396                    }
2397               } else {
2398#ifndef ABC_INHERIT
2399                    model2->primal(primalStartup);
2400#else
2401                    model2->dealWithAbc(1,primalStartup,interrupt);
2402#endif
2403               }
2404          }
2405#endif
2406          time2 = CoinCpuTime();
2407          timeCore = time2 - timeX;
2408          handler_->message(CLP_INTERVAL_TIMING, messages_)
2409                    << "Primal" << timeCore << time2 - time1
2410                    << CoinMessageEol;
2411          timeX = time2;
2412     } else if (method == ClpSolve::usePrimalorSprint) {
2413          // Sprint
2414          /*
2415            This driver implements what I called Sprint when I introduced the idea
2416            many years ago.  Cplex calls it "sifting" which I think is just as silly.
2417            When I thought of this trivial idea
2418            it reminded me of an LP code of the 60's called sprint which after
2419            every factorization took a subset of the matrix into memory (all
2420            64K words!) and then iterated very fast on that subset.  On the
2421            problems of those days it did not work very well, but it worked very
2422            well on aircrew scheduling problems where there were very large numbers
2423            of columns all with the same flavor.
2424          */
2425
2426          /* The idea works best if you can get feasible easily.  To make it
2427             more general we can add in costed slacks */
2428
2429          int originalNumberColumns = model2->numberColumns();
2430          int numberRows = model2->numberRows();
2431          ClpSimplex * originalModel2 = model2;
2432
2433          // We will need arrays to choose variables.  These are too big but ..
2434          double * weight = new double [numberRows+originalNumberColumns];
2435          int * sort = new int [numberRows+originalNumberColumns];
2436          int numberSort = 0;
2437          // We are going to add slacks to get feasible.
2438          // initial list will just be artificials
2439          int iColumn;
2440          const double * columnLower = model2->columnLower();
2441          const double * columnUpper = model2->columnUpper();
2442          double * columnSolution = model2->primalColumnSolution();
2443
2444          // See if we have costed slacks
2445          int * negSlack = new int[numberRows];
2446          int * posSlack = new int[numberRows];
2447          int iRow;
2448          for (iRow = 0; iRow < numberRows; iRow++) {
2449               negSlack[iRow] = -1;
2450               posSlack[iRow] = -1;
2451          }
2452          const double * element = model2->matrix()->getElements();
2453          const int * row = model2->matrix()->getIndices();
2454          const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
2455          const int * columnLength = model2->matrix()->getVectorLengths();
2456          //bool allSlack = (numberRowsBasic==numberRows);
2457          for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2458               if (!columnSolution[iColumn] || fabs(columnSolution[iColumn]) > 1.0e20) {
2459                    double value = 0.0;
2460                    if (columnLower[iColumn] > 0.0)
2461                         value = columnLower[iColumn];
2462                    else if (columnUpper[iColumn] < 0.0)
2463                         value = columnUpper[iColumn];
2464                    columnSolution[iColumn] = value;
2465               }
2466               if (columnLength[iColumn] == 1) {
2467                    int jRow = row[columnStart[iColumn]];
2468                    if (!columnLower[iColumn]) {
2469                         if (element[columnStart[iColumn]] > 0.0 && posSlack[jRow] < 0)
2470                              posSlack[jRow] = iColumn;
2471                         else if (element[columnStart[iColumn]] < 0.0 && negSlack[jRow] < 0)
2472                              negSlack[jRow] = iColumn;
2473                    } else if (!columnUpper[iColumn]) {
2474                         if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0)
2475                              posSlack[jRow] = iColumn;
2476                         else if (element[columnStart[iColumn]] > 0.0 && negSlack[jRow] < 0)
2477                              negSlack[jRow] = iColumn;
2478                    }
2479               }
2480          }
2481          // now see what that does to row solution
2482          double * rowSolution = model2->primalRowSolution();
2483          CoinZeroN (rowSolution, numberRows);
2484          model2->clpMatrix()->times(1.0, columnSolution, rowSolution);
2485          // See if we can adjust using costed slacks
2486          double penalty = CoinMax(1.0e5, CoinMin(infeasibilityCost_ * 0.01, 1.0e10)) * optimizationDirection_;
2487          const double * lower = model2->rowLower();
2488          const double * upper = model2->rowUpper();
2489          for (iRow = 0; iRow < numberRows; iRow++) {
2490               if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
2491                    int jColumn = posSlack[iRow];
2492                    if (jColumn >= 0) {
2493                         if (columnSolution[jColumn])
2494                              continue;
2495                         double difference = lower[iRow] - rowSolution[iRow];
2496                         double elementValue = element[columnStart[jColumn]];
2497                         if (elementValue > 0.0) {
2498                              double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
2499                              columnSolution[jColumn] = movement;
2500                              rowSolution[iRow] += movement * elementValue;
2501                         } else {
2502                              double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
2503                              columnSolution[jColumn] = movement;
2504                              rowSolution[iRow] += movement * elementValue;
2505                         }
2506                    }
2507               } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
2508                    int jColumn = negSlack[iRow];
2509                    if (jColumn >= 0) {
2510                         if (columnSolution[jColumn])
2511                              continue;
2512                         double difference = upper[iRow] - rowSolution[iRow];
2513                         double elementValue = element[columnStart[jColumn]];
2514                         if (elementValue < 0.0) {
2515                              double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
2516                              columnSolution[jColumn] = movement;
2517                              rowSolution[iRow] += movement * elementValue;
2518                         } else {
2519                              double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
2520                              columnSolution[jColumn] = movement;
2521                              rowSolution[iRow] += movement * elementValue;
2522                         }
2523                    }
2524               }
2525          }
2526          delete [] negSlack;
2527          delete [] posSlack;
2528          int nRow = numberRows;
2529          bool network = false;
2530          if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
2531               network = true;
2532               nRow *= 2;
2533          }
2534          int * addStarts = new int [nRow+1];
2535          int * addRow = new int[nRow];
2536          double * addElement = new double[nRow];
2537          addStarts[0] = 0;
2538          int numberArtificials = 0;
2539          int numberAdd = 0;
2540          double * addCost = new double [numberRows];
2541          for (iRow = 0; iRow < numberRows; iRow++) {
2542               if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
2543                    addRow[numberAdd] = iRow;
2544                    addElement[numberAdd++] = 1.0;
2545                    if (network) {
2546                         addRow[numberAdd] = numberRows;
2547                         addElement[numberAdd++] = -1.0;
2548                    }
2549                    addCost[numberArtificials] = penalty;
2550                    numberArtificials++;
2551                    addStarts[numberArtificials] = numberAdd;
2552               } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
2553                    addRow[numberAdd] = iRow;
2554                    addElement[numberAdd++] = -1.0;
2555                    if (network) {
2556                         addRow[numberAdd] = numberRows;
2557                         addElement[numberAdd++] = 1.0;
2558                    }
2559                    addCost[numberArtificials] = penalty;
2560                    numberArtificials++;
2561                    addStarts[numberArtificials] = numberAdd;
2562               }
2563          }
2564          if (numberArtificials) {
2565               // need copy so as not to disturb original
2566               model2 = new ClpSimplex(*model2);
2567               if (network) {
2568                    // network - add a null row
2569                    model2->addRow(0, NULL, NULL, -COIN_DBL_MAX, COIN_DBL_MAX);
2570                    numberRows++;
2571               }
2572               model2->addColumns(numberArtificials, NULL, NULL, addCost,
2573                                  addStarts, addRow, addElement);
2574          }
2575          delete [] addStarts;
2576          delete [] addRow;
2577          delete [] addElement;
2578          delete [] addCost;
2579          // look at rhs to see if to perturb
2580          double largest = 0.0;
2581          double smallest = 1.0e30;
2582          for (iRow = 0; iRow < numberRows; iRow++) {
2583               double value;
2584               value = fabs(model2->rowLower_[iRow]);
2585               if (value && value < 1.0e30) {
2586                    largest = CoinMax(largest, value);
2587                    smallest = CoinMin(smallest, value);
2588               }
2589               value = fabs(model2->rowUpper_[iRow]);
2590               if (value && value < 1.0e30) {
2591                    largest = CoinMax(largest, value);
2592                    smallest = CoinMin(smallest, value);
2593               }
2594          }
2595          double * saveLower = NULL;
2596          double * saveUpper = NULL;
2597          if (largest < -2.01 * smallest) { // switch off for now
2598               // perturb - so switch off standard
2599               model2->setPerturbation(100);
2600               saveLower = new double[numberRows];
2601               CoinMemcpyN(model2->rowLower_, numberRows, saveLower);
2602               saveUpper = new double[numberRows];
2603               CoinMemcpyN(model2->rowUpper_, numberRows, saveUpper);
2604               double * lower = model2->rowLower();
2605               double * upper = model2->rowUpper();
2606               for (iRow = 0; iRow < numberRows; iRow++) {
2607                    double lowerValue = lower[iRow], upperValue = upper[iRow];
2608                    double value = randomNumberGenerator_.randomDouble();
2609                    if (upperValue > lowerValue + primalTolerance_) {
2610                         if (lowerValue > -1.0e20 && lowerValue)
2611                              lowerValue -= value * 1.0e-6 * fabs(lowerValue);
2612                         if (upperValue < 1.0e20 && upperValue)
2613                              upperValue += value * 1.0e-6 * fabs(upperValue);
2614                    } else if (upperValue > 0.0) {
2615                         upperValue -= value * 1.0e-6 * fabs(lowerValue);
2616                         lowerValue -= value * 1.0e-6 * fabs(lowerValue);
2617                    } else if (upperValue < 0.0) {
2618                         upperValue += value * 1.0e-6 * fabs(lowerValue);
2619                         lowerValue += value * 1.0e-6 * fabs(lowerValue);
2620                    } else {
2621                    }
2622                    lower[iRow] = lowerValue;
2623                    upper[iRow] = upperValue;
2624               }
2625          }
2626          int i;
2627          // Just do this number of passes in Sprint
2628          if (doSprint > 0)
2629               maxSprintPass = options.getExtraInfo(1);
2630          // but if big use to get ratio
2631          double ratio = 3;
2632#ifdef CLP_USEFUL_PRINTOUT
2633          debugInt[6]=3;
2634          debugInt[7]=maxSprintPass;
2635#endif
2636          if (maxSprintPass > 1000) {
2637               ratio = static_cast<double> (maxSprintPass) * 0.0001;
2638               ratio = CoinMax(ratio, 1.1);
2639               maxSprintPass = maxSprintPass % 1000;
2640#ifdef COIN_DEVELOP
2641               printf("%d passes wanted with ratio of %g\n", maxSprintPass, ratio);
2642#endif
2643          }
2644          // Just take this number of columns in small problem
2645          int smallNumberColumns = static_cast<int> (CoinMin(ratio * numberRows, static_cast<double> (numberColumns)));
2646          smallNumberColumns = CoinMax(smallNumberColumns, 3000);
2647          smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
2648          //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
2649          //smallNumberColumns = CoinMax(smallNumberColumns,3000);
2650          //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
2651          // redo as may have changed
2652          columnLower = model2->columnLower();
2653          columnUpper = model2->columnUpper();
2654          columnSolution = model2->primalColumnSolution();
2655          // Set up initial list
2656          numberSort = 0;
2657          if (numberArtificials) {
2658               numberSort = numberArtificials;
2659               for (i = 0; i < numberSort; i++)
2660                    sort[i] = i + originalNumberColumns;
2661          }
2662          // put in free
2663          // maybe a solution there already
2664          for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2665               if (model2->getColumnStatus(iColumn) == basic||
2666                   (columnLower[iColumn]<-1.0e30&&columnUpper[iColumn]>1.0e30))
2667                    sort[numberSort++] = iColumn;
2668          }
2669          for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2670               if (model2->getColumnStatus(iColumn) != basic) {
2671                    if (columnSolution[iColumn] > columnLower[iColumn] &&
2672                              columnSolution[iColumn] < columnUpper[iColumn] &&
2673                              columnSolution[iColumn])
2674                         sort[numberSort++] = iColumn;
2675               }
2676          }
2677          numberSort = CoinMin(numberSort, smallNumberColumns);
2678
2679          int numberColumns = model2->numberColumns();
2680          double * fullSolution = model2->primalColumnSolution();
2681
2682
2683          int iPass;
2684          double lastObjective[] = {1.0e31,1.0e31};
2685          // It will be safe to allow dense
2686          model2->setInitialDenseFactorization(true);
2687
2688          // We will be using all rows
2689          int * whichRows = new int [numberRows];
2690          for (iRow = 0; iRow < numberRows; iRow++)
2691               whichRows[iRow] = iRow;
2692          double originalOffset;
2693          model2->getDblParam(ClpObjOffset, originalOffset);
2694          int totalIterations = 0;
2695          double lastSumArtificials = COIN_DBL_MAX;
2696          int originalMaxSprintPass = maxSprintPass;
2697          maxSprintPass = 20; // so we do that many if infeasible
2698          for (iPass = 0; iPass < maxSprintPass; iPass++) {
2699               //printf("Bug until submodel new version\n");
2700               //CoinSort_2(sort,sort+numberSort,weight);
2701               // Create small problem
2702               ClpSimplex small(model2, numberRows, whichRows, numberSort, sort);
2703               small.setPerturbation(model2->perturbation());
2704               small.setInfeasibilityCost(model2->infeasibilityCost());
2705               if (model2->factorizationFrequency() == 200) {
2706                    // User did not touch preset
2707                    small.defaultFactorizationFrequency();
2708               }
2709               // now see what variables left out do to row solution
2710               double * rowSolution = model2->primalRowSolution();
2711               double * sumFixed = new double[numberRows];
2712               CoinZeroN (sumFixed, numberRows);
2713               int iRow, iColumn;
2714               // zero out ones in small problem
2715               for (iColumn = 0; iColumn < numberSort; iColumn++) {
2716                    int kColumn = sort[iColumn];
2717                    fullSolution[kColumn] = 0.0;
2718               }
2719               // Get objective offset
2720               const double * objective = model2->objective();
2721               double offset = 0.0;
2722               for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
2723                    offset += fullSolution[iColumn] * objective[iColumn];
2724#if 0
2725               // Set artificials to zero if first time close to zero
2726               for (iColumn = originalNumberColumns; iColumn < numberColumns; iColumn++) {
2727                 if (fullSolution[iColumn]<primalTolerance_&&objective[iColumn]==penalty) {
2728                   model2->objective()[iColumn]=2.0*penalty;
2729                   fullSolution[iColumn]=0.0;
2730                 }
2731               }
2732#endif
2733               small.setDblParam(ClpObjOffset, originalOffset - offset);
2734               model2->clpMatrix()->times(1.0, fullSolution, sumFixed);
2735
2736               double * lower = small.rowLower();
2737               double * upper = small.rowUpper();
2738               for (iRow = 0; iRow < numberRows; iRow++) {
2739                    if (lower[iRow] > -1.0e50)
2740                         lower[iRow] -= sumFixed[iRow];
2741                    if (upper[iRow] < 1.0e50)
2742                         upper[iRow] -= sumFixed[iRow];
2743                    rowSolution[iRow] -= sumFixed[iRow];
2744               }
2745               delete [] sumFixed;
2746               // Solve
2747               if (interrupt)
2748                    currentModel = &small;
2749               small.defaultFactorizationFrequency();
2750               if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
2751                    // See if original wanted vector
2752                    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
2753                    ClpMatrixBase * matrix = small.clpMatrix();
2754                    if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
2755                         ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
2756                         clpMatrix->makeSpecialColumnCopy();
2757                         small.primal(1);
2758                         clpMatrix->releaseSpecialColumnCopy();
2759                    } else {
2760#if 1
2761#ifdef ABC_INHERIT
2762                      //small.writeMps("try.mps");
2763                      if (iPass||!numberArtificials) 
2764                         small.dealWithAbc(1,1);
2765                       else 
2766                         small.dealWithAbc(0,0);
2767#else
2768                      if (iPass||!numberArtificials)
2769                         small.primal(1);
2770                      else
2771                         small.dual(0);
2772#endif
2773                      if (small.problemStatus())
2774                        small.dual(0);
2775#else
2776                         int numberColumns = small.numberColumns();
2777                         int numberRows = small.numberRows();
2778                         // Use dual region
2779                         double * rhs = small.dualRowSolution();
2780                         int * whichRow = new int[3*numberRows];
2781                         int * whichColumn = new int[2*numberColumns];
2782                         int nBound;
2783                         ClpSimplex * small2 = ((ClpSimplexOther *) (&small))->crunch(rhs, whichRow, whichColumn,
2784                                               nBound, false, false);
2785                         if (small2) {
2786#ifdef ABC_INHERIT
2787                              small2->dealWithAbc(1,1);
2788#else
2789                              small.primal(1);
2790#endif
2791                              if (small2->problemStatus() == 0) {
2792                                   small.setProblemStatus(0);
2793                                   ((ClpSimplexOther *) (&small))->afterCrunch(*small2, whichRow, whichColumn, nBound);
2794                              } else {
2795#ifdef ABC_INHERIT
2796                                   small2->dealWithAbc(1,1);
2797#else
2798                                   small.primal(1);
2799#endif
2800                                   if (small2->problemStatus())
2801                                        small.primal(1);
2802                              }
2803                              delete small2;
2804                         } else {
2805                              small.primal(1);
2806                         }
2807                         delete [] whichRow;
2808                         delete [] whichColumn;
2809#endif
2810                    }
2811               } else {
2812                    small.primal(1);
2813               }
2814               totalIterations += small.numberIterations();
2815               // move solution back
2816               const double * solution = small.primalColumnSolution();
2817               for (iColumn = 0; iColumn < numberSort; iColumn++) {
2818                    int kColumn = sort[iColumn];
2819                    model2->setColumnStatus(kColumn, small.getColumnStatus(iColumn));
2820                    fullSolution[kColumn] = solution[iColumn];
2821               }
2822               for (iRow = 0; iRow < numberRows; iRow++)
2823                    model2->setRowStatus(iRow, small.getRowStatus(iRow));
2824               CoinMemcpyN(small.primalRowSolution(),
2825                           numberRows, model2->primalRowSolution());
2826               double sumArtificials = 0.0;
2827               for (i = 0; i < numberArtificials; i++)
2828                    sumArtificials += fullSolution[i + originalNumberColumns];
2829               if (sumArtificials && iPass > 5 && sumArtificials >= lastSumArtificials) {
2830                    // increase costs
2831                    double * cost = model2->objective() + originalNumberColumns;
2832                    double newCost = CoinMin(1.0e10, cost[0] * 1.5);
2833                    for (i = 0; i < numberArtificials; i++)
2834                         cost[i] = newCost;
2835               }
2836               lastSumArtificials = sumArtificials;
2837               // get reduced cost for large problem
2838               double * djs = model2->dualColumnSolution();
2839               CoinMemcpyN(model2->objective(), numberColumns, djs);
2840               model2->clpMatrix()->transposeTimes(-1.0, small.dualRowSolution(), djs);
2841               int numberNegative = 0;
2842               double sumNegative = 0.0;
2843               // now massage weight so all basic in plus good djs
2844               // first count and do basic
2845               numberSort = 0;
2846               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2847                    double dj = djs[iColumn] * optimizationDirection_;
2848                    double value = fullSolution[iColumn];
2849                    if (model2->getColumnStatus(iColumn) == ClpSimplex::basic) {
2850                         sort[numberSort++] = iColumn;
2851                    } else if (dj < -dualTolerance_ && value < columnUpper[iColumn]) {
2852                         numberNegative++;
2853                         sumNegative -= dj;
2854                    } else if (dj > dualTolerance_ && value > columnLower[iColumn]) {
2855                         numberNegative++;
2856                         sumNegative += dj;
2857                    }
2858               }
2859               handler_->message(CLP_SPRINT, messages_)
2860                         << iPass + 1 << small.numberIterations() << small.objectiveValue() << sumNegative
2861                         << numberNegative
2862                         << CoinMessageEol;
2863               if (sumArtificials < 1.0e-8 && originalMaxSprintPass >= 0) {
2864                    maxSprintPass = iPass + originalMaxSprintPass;
2865                    originalMaxSprintPass = -1;
2866               }
2867               if (iPass > 20)
2868                    sumArtificials = 0.0;
2869               if ((small.objectiveValue()*optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 5 && sumArtificials < 1.0e-8) ||
2870                         (!small.numberIterations() && iPass) ||
2871                         iPass == maxSprintPass - 1 || small.status() == 3) {
2872
2873                    break; // finished
2874               } else {
2875                    lastObjective[1] = lastObjective[0];
2876                    lastObjective[0] = small.objectiveValue() * optimizationDirection_;
2877                    double tolerance;
2878                    double averageNegDj = sumNegative / static_cast<double> (numberNegative + 1);
2879                    if (numberNegative + numberSort > smallNumberColumns)
2880                         tolerance = -dualTolerance_;
2881                    else
2882                         tolerance = 10.0 * averageNegDj;
2883                    int saveN = numberSort;
2884                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2885                         double dj = djs[iColumn] * optimizationDirection_;
2886                         double value = fullSolution[iColumn];
2887                         if (model2->getColumnStatus(iColumn) != ClpSimplex::basic) {
2888                              if (dj < -dualTolerance_ && value < columnUpper[iColumn])
2889                                   dj = dj;
2890                              else if (dj > dualTolerance_ && value > columnLower[iColumn])
2891                                   dj = -dj;
2892                              else if (columnUpper[iColumn] > columnLower[iColumn])
2893                                   dj = fabs(dj);
2894                              else
2895                                   dj = 1.0e50;
2896                              if (dj < tolerance) {
2897                                   weight[numberSort] = dj;
2898                                   sort[numberSort++] = iColumn;
2899                              }
2900                         }
2901                    }
2902                    // sort
2903                    CoinSort_2(weight + saveN, weight + numberSort, sort + saveN);
2904                    numberSort = CoinMin(smallNumberColumns, numberSort);
2905               }
2906          }
2907          if (interrupt)
2908               currentModel = model2;
2909          for (i = 0; i < numberArtificials; i++)
2910               sort[i] = i + originalNumberColumns;
2911          model2->deleteColumns(numberArtificials, sort);
2912          if (network) {
2913               int iRow = numberRows - 1;
2914               model2->deleteRows(1, &iRow);
2915          }
2916          delete [] weight;
2917          delete [] sort;
2918          delete [] whichRows;
2919          if (saveLower) {
2920               // unperturb and clean
2921               for (iRow = 0; iRow < numberRows; iRow++) {
2922                    model2->rowLower_[iRow] = saveLower[iRow];
2923                    model2->rowUpper_[iRow] = saveUpper[iRow];
2924               }
2925               delete [] saveLower;
2926               delete [] saveUpper;
2927          }
2928#ifdef ABC_INHERIT
2929          model2->dealWithAbc(1,1);
2930#else
2931          model2->primal(1);
2932#endif
2933          model2->setPerturbation(savePerturbation);
2934          if (model2 != originalModel2) {
2935               originalModel2->moveInfo(*model2);
2936               delete model2;
2937               model2 = originalModel2;
2938          }
2939          time2 = CoinCpuTime();
2940          timeCore = time2 - timeX;
2941          handler_->message(CLP_INTERVAL_TIMING, messages_)
2942                    << "Sprint" << timeCore << time2 - time1
2943                    << CoinMessageEol;
2944          timeX = time2;
2945          model2->setNumberIterations(model2->numberIterations() + totalIterations);
2946     } else if (method == ClpSolve::useBarrier || method == ClpSolve::useBarrierNoCross) {
2947#ifndef SLIM_CLP
2948          //printf("***** experimental pretty crude barrier\n");
2949          //#define SAVEIT 2
2950#ifndef SAVEIT
2951#define BORROW
2952#endif
2953#ifdef BORROW
2954          ClpInterior barrier;
2955          barrier.borrowModel(*model2);
2956#else
2957          ClpInterior barrier(*model2);
2958#endif
2959          if (interrupt)
2960               currentModel2 = &barrier;
2961          if (barrier.numberRows()+barrier.numberColumns()>10000)
2962            barrier.setMaximumBarrierIterations(1000);
2963          int barrierOptions = options.getSpecialOption(4);
2964          int aggressiveGamma = 0;
2965          bool presolveInCrossover = false;
2966          bool scale = false;
2967          bool doKKT = false;
2968          bool forceFixing = false;
2969          int speed = 0;
2970          if (barrierOptions & 16) {
2971               barrierOptions &= ~16;
2972               doKKT = true;
2973          }
2974          if (barrierOptions&(32 + 64 + 128)) {
2975               aggressiveGamma = (barrierOptions & (32 + 64 + 128)) >> 5;
2976               barrierOptions &= ~(32 + 64 + 128);
2977          }
2978          if (barrierOptions & 256) {
2979               barrierOptions &= ~256;
2980               presolveInCrossover = true;
2981          }
2982          if (barrierOptions & 512) {
2983               barrierOptions &= ~512;
2984               forceFixing = true;
2985          }
2986          if (barrierOptions & 1024) {
2987               barrierOptions &= ~1024;
2988               barrier.setProjectionTolerance(1.0e-9);
2989          }
2990          if (barrierOptions&(2048 | 4096)) {
2991               speed = (barrierOptions & (2048 | 4096)) >> 11;
2992               barrierOptions &= ~(2048 | 4096);
2993          }
2994          if (barrierOptions & 8) {
2995               barrierOptions &= ~8;
2996               scale = true;
2997          }
2998          // If quadratic force KKT
2999          if (quadraticObj) {
3000               doKKT = true;
3001          }
3002          switch (barrierOptions) {
3003          case 0:
3004          default:
3005               if (!doKKT) {
3006                    ClpCholeskyBase * cholesky = new ClpCholeskyBase(options.getExtraInfo(1));
3007                    cholesky->setIntegerParameter(0, speed);
3008                    barrier.setCholesky(cholesky);
3009               } else {
3010                    ClpCholeskyBase * cholesky = new ClpCholeskyBase();
3011                    cholesky->setKKT(true);
3012                    barrier.setCholesky(cholesky);
3013               }
3014               break;
3015          case 1:
3016               if (!doKKT) {
3017                    ClpCholeskyDense * cholesky = new ClpCholeskyDense();
3018                    barrier.setCholesky(cholesky);
3019               } else {
3020                    ClpCholeskyDense * cholesky = new ClpCholeskyDense();
3021                    cholesky->setKKT(true);
3022                    barrier.setCholesky(cholesky);
3023               }
3024               break;
3025#ifdef COIN_HAS_WSMP
3026          case 2: {
3027               if (!doKKT) {
3028                    ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100, model2->numberRows() / 10));
3029                    barrier.setCholesky(cholesky);
3030               } else {
3031                    ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3032                    barrier.setCholesky(cholesky);
3033               }
3034          }
3035          break;
3036          case 3:
3037               if (!doKKT) {
3038                    ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
3039                    barrier.setCholesky(cholesky);
3040               } else {
3041                    ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3042                    barrier.setCholesky(cholesky);
3043               }
3044               break;
3045#endif
3046#ifdef UFL_BARRIER
3047          case 4:
3048               if (!doKKT) {
3049                    ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
3050                    barrier.setCholesky(cholesky);
3051               } else {
3052                    ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
3053                    cholesky->setKKT(true);
3054                    barrier.setCholesky(cholesky);
3055               }
3056               break;
3057#endif
3058#ifdef TAUCS_BARRIER
3059          case 5: {
3060               ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
3061               barrier.setCholesky(cholesky);
3062               assert (!doKKT);
3063          }
3064          break;
3065#endif
3066#ifdef COIN_HAS_MUMPS
3067          case 6: {
3068               if (!doKKT) {
3069                    ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
3070                    barrier.setCholesky(cholesky);
3071               } else {
3072                    ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
3073                    cholesky->setKKT(true);
3074                    barrier.setCholesky(cholesky);
3075               }
3076          }
3077          break;
3078#endif
3079          }
3080          int numberRows = model2->numberRows();
3081          int numberColumns = model2->numberColumns();
3082          int saveMaxIts = model2->maximumIterations();
3083          if (saveMaxIts < 1000) {
3084               barrier.setMaximumBarrierIterations(saveMaxIts);
3085               model2->setMaximumIterations(10000000);
3086          }
3087#ifndef SAVEIT
3088          //barrier.setDiagonalPerturbation(1.0e-25);
3089          if (aggressiveGamma) {
3090               switch (aggressiveGamma) {
3091               case 1:
3092                    barrier.setGamma(1.0e-5);
3093                    barrier.setDelta(1.0e-5);
3094                    break;
3095               case 2:
3096                    barrier.setGamma(1.0e-7);
3097                    break;
3098               case 3:
3099                    barrier.setDelta(1.0e-5);
3100                    break;
3101               case 4:
3102                    barrier.setGamma(1.0e-3);
3103                    barrier.setDelta(1.0e-3);
3104                    break;
3105               case 5:
3106                    barrier.setGamma(1.0e-3);
3107                    break;
3108               case 6:
3109                    barrier.setDelta(1.0e-3);
3110                    break;
3111               }
3112          }
3113          if (scale)
3114               barrier.scaling(1);
3115          else
3116               barrier.scaling(0);
3117          barrier.primalDual();
3118#elif SAVEIT==1
3119          barrier.primalDual();
3120#else
3121          model2->restoreModel("xx.save");
3122          // move solutions
3123          CoinMemcpyN(model2->primalRowSolution(),
3124                      numberRows, barrier.primalRowSolution());
3125          CoinMemcpyN(model2->dualRowSolution(),
3126                      numberRows, barrier.dualRowSolution());
3127          CoinMemcpyN(model2->primalColumnSolution(),
3128                      numberColumns, barrier.primalColumnSolution());
3129          CoinMemcpyN(model2->dualColumnSolution(),
3130                      numberColumns, barrier.dualColumnSolution());
3131#endif
3132          time2 = CoinCpuTime();
3133          timeCore = time2 - timeX;
3134          handler_->message(CLP_INTERVAL_TIMING, messages_)
3135                    << "Barrier" << timeCore << time2 - time1
3136                    << CoinMessageEol;
3137          timeX = time2;
3138          int maxIts = barrier.maximumBarrierIterations();
3139          int barrierStatus = barrier.status();
3140          double gap = barrier.complementarityGap();
3141          // get which variables are fixed
3142          double * saveLower = NULL;
3143          double * saveUpper = NULL;
3144          ClpPresolve pinfo2;
3145          ClpSimplex * saveModel2 = NULL;
3146          bool extraPresolve = false;
3147          int numberFixed = barrier.numberFixed();
3148          if (numberFixed) {
3149               int numberRows = barrier.numberRows();
3150               int numberColumns = barrier.numberColumns();
3151               int numberTotal = numberRows + numberColumns;
3152               saveLower = new double [numberTotal];
3153               saveUpper = new double [numberTotal];
3154               CoinMemcpyN(barrier.columnLower(), numberColumns, saveLower);
3155               CoinMemcpyN(barrier.rowLower(), numberRows, saveLower + numberColumns);
3156               CoinMemcpyN(barrier.columnUpper(), numberColumns, saveUpper);
3157               CoinMemcpyN(barrier.rowUpper(), numberRows, saveUpper + numberColumns);
3158          }
3159          if (((numberFixed * 20 > barrier.numberRows() && numberFixed > 5000) || forceFixing) &&
3160                    presolveInCrossover) {
3161               // may as well do presolve
3162               if (!forceFixing) {
3163                    barrier.fixFixed();
3164               } else {
3165                    // Fix
3166                    int n = barrier.numberColumns();
3167                    double * lower = barrier.columnLower();
3168                    double * upper = barrier.columnUpper();
3169                    double * solution = barrier.primalColumnSolution();
3170                    int nFix = 0;
3171                    for (int i = 0; i < n; i++) {
3172                         if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3173                              double value = solution[i];
3174                              if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3175                                   solution[i] = lower[i];
3176                                   upper[i] = lower[i];
3177                                   nFix++;
3178                              } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3179                                   solution[i] = upper[i];
3180                                   lower[i] = upper[i];
3181                                   nFix++;
3182                              }
3183                         }
3184                    }
3185#ifdef CLP_INVESTIGATE
3186                    printf("%d columns fixed\n", nFix);
3187#endif
3188                    int nr = barrier.numberRows();
3189                    lower = barrier.rowLower();
3190                    upper = barrier.rowUpper();
3191                    solution = barrier.primalRowSolution();
3192                    nFix = 0;
3193                    for (int i = 0; i < nr; i++) {
3194                         if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3195                              double value = solution[i];
3196                              if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3197                                   solution[i] = lower[i];
3198                                   upper[i] = lower[i];
3199                                   nFix++;
3200                              } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3201                                   solution[i] = upper[i];
3202                                   lower[i] = upper[i];
3203                                   nFix++;
3204                              }
3205                         }
3206                    }
3207#ifdef CLP_INVESTIGATE
3208                    printf("%d row slacks fixed\n", nFix);
3209#endif
3210               }
3211               saveModel2 = model2;
3212               extraPresolve = true;
3213          } else if (numberFixed) {
3214               // Set fixed to bounds (may have restored earlier solution)
3215               if (!forceFixing) {
3216                    barrier.fixFixed(false);
3217               } else {
3218                    // Fix
3219                    int n = barrier.numberColumns();
3220                    double * lower = barrier.columnLower();
3221                    double * upper = barrier.columnUpper();
3222                    double * solution = barrier.primalColumnSolution();
3223                    int nFix = 0;
3224                    for (int i = 0; i < n; i++) {
3225                         if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3226                              double value = solution[i];
3227                              if (value < lower[i] + 1.0e-8 && value - lower[i] < upper[i] - value) {
3228                                   solution[i] = lower[i];
3229                                   upper[i] = lower[i];
3230                                   nFix++;
3231                              } else if (value > upper[i] - 1.0e-8 && value - lower[i] > upper[i] - value) {
3232                                   solution[i] = upper[i];
3233                                   lower[i] = upper[i];
3234                                   nFix++;
3235                              } else {
3236                                   //printf("fixcol %d %g <= %g <= %g\n",
3237                                   //     i,lower[i],solution[i],upper[i]);
3238                              }
3239                         }
3240                    }
3241#ifdef CLP_INVESTIGATE
3242                    printf("%d columns fixed\n", nFix);
3243#endif
3244                    int nr = barrier.numberRows();
3245                    lower = barrier.rowLower();
3246                    upper = barrier.rowUpper();
3247                    solution = barrier.primalRowSolution();
3248                    nFix = 0;
3249                    for (int i = 0; i < nr; i++) {
3250                         if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3251                              double value = solution[i];
3252                              if (value < lower[i] + 1.0e-5 && value - lower[i] < upper[i] - value) {
3253                                   solution[i] = lower[i];
3254                                   upper[i] = lower[i];
3255                                   nFix++;
3256                              } else if (value > upper[i] - 1.0e-5 && value - lower[i] > upper[i] - value) {
3257                                   solution[i] = upper[i];
3258                                   lower[i] = upper[i];
3259                                   nFix++;
3260                              } else {
3261                                   //printf("fixrow %d %g <= %g <= %g\n",
3262                                   //     i,lower[i],solution[i],upper[i]);
3263                              }
3264                         }
3265                    }
3266#ifdef CLP_INVESTIGATE
3267                    printf("%d row slacks fixed\n", nFix);
3268#endif
3269               }
3270          }
3271#ifdef BORROW
3272          int saveNumberIterations = barrier.numberIterations();
3273          barrier.returnModel(*model2);
3274          double * rowPrimal = new double [numberRows];
3275          double * columnPrimal = new double [numberColumns];
3276          double * rowDual = new double [numberRows];
3277          double * columnDual = new double [numberColumns];
3278          // move solutions other way
3279          CoinMemcpyN(model2->primalRowSolution(),
3280                      numberRows, rowPrimal);
3281          CoinMemcpyN(model2->dualRowSolution(),
3282                      numberRows, rowDual);
3283          CoinMemcpyN(model2->primalColumnSolution(),
3284                      numberColumns, columnPrimal);
3285          CoinMemcpyN(model2->dualColumnSolution(),
3286                      numberColumns, columnDual);
3287#else
3288          double * rowPrimal = barrier.primalRowSolution();
3289          double * columnPrimal = barrier.primalColumnSolution();
3290          double * rowDual = barrier.dualRowSolution();
3291          double * columnDual = barrier.dualColumnSolution();
3292          // move solutions
3293          CoinMemcpyN(rowPrimal,
3294                      numberRows, model2->primalRowSolution());
3295          CoinMemcpyN(rowDual,
3296                      numberRows, model2->dualRowSolution());
3297          CoinMemcpyN(columnPrimal,
3298                      numberColumns, model2->primalColumnSolution());
3299          CoinMemcpyN(columnDual,
3300                      numberColumns, model2->dualColumnSolution());
3301#endif
3302          if (saveModel2) {
3303               // do presolve
3304               model2 = pinfo2.presolvedModel(*model2, dblParam_[ClpPresolveTolerance],
3305                                              false, 5, true);
3306               if (!model2) {
3307                    model2 = saveModel2;
3308                    saveModel2 = NULL;
3309                    int numberRows = model2->numberRows();
3310                    int numberColumns = model2->numberColumns();
3311                    CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3312                    CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3313                    delete [] saveLower;
3314                    CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3315                    CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3316                    delete [] saveUpper;
3317                    saveLower = NULL;
3318                    saveUpper = NULL;
3319               }
3320          }
3321          if (method == ClpSolve::useBarrier || barrierStatus < 0) {
3322               if (maxIts && barrierStatus < 4 && !quadraticObj) {
3323                    //printf("***** crossover - needs more thought on difficult models\n");
3324#if SAVEIT==1
3325                    model2->ClpSimplex::saveModel("xx.save");
3326#endif
3327                    // make sure no status left
3328                    model2->createStatus();
3329                    // solve
3330                    if (!forceFixing)
3331                         model2->setPerturbation(100);
3332                    if (model2->factorizationFrequency() == 200) {
3333                         // User did not touch preset
3334                         model2->defaultFactorizationFrequency();
3335                    }
3336#if 1 //ndef ABC_INHERIT //#if 1
3337                    // throw some into basis
3338                    if(!forceFixing) {
3339                         int numberRows = model2->numberRows();
3340                         int numberColumns = model2->numberColumns();
3341                         double * dsort = new double[numberColumns];
3342                         int * sort = new int[numberColumns];
3343                         int n = 0;
3344                         const double * columnLower = model2->columnLower();
3345                         const double * columnUpper = model2->columnUpper();
3346                         double * primalSolution = model2->primalColumnSolution();
3347                         const double * dualSolution = model2->dualColumnSolution();
3348                         double tolerance = 10.0 * primalTolerance_;
3349                         int i;
3350                         for ( i = 0; i < numberRows; i++)
3351                              model2->setRowStatus(i, superBasic);
3352                         for ( i = 0; i < numberColumns; i++) {
3353                              double distance = CoinMin(columnUpper[i] - primalSolution[i],
3354                                                        primalSolution[i] - columnLower[i]);
3355                              if (distance > tolerance) {
3356                                   if (fabs(dualSolution[i]) < 1.0e-5)
3357                                        distance *= 100.0;
3358                                   dsort[n] = -distance;
3359                                   sort[n++] = i;
3360                                   model2->setStatus(i, superBasic);
3361                              } else if (distance > primalTolerance_) {
3362                                   model2->setStatus(i, superBasic);
3363                              } else if (primalSolution[i] <= columnLower[i] + primalTolerance_) {
3364                                   model2->setStatus(i, atLowerBound);
3365                                   primalSolution[i] = columnLower[i];
3366                              } else {
3367                                   model2->setStatus(i, atUpperBound);
3368                                   primalSolution[i] = columnUpper[i];
3369                              }
3370                         }
3371                         CoinSort_2(dsort, dsort + n, sort);
3372                         n = CoinMin(numberRows, n);
3373                         for ( i = 0; i < n; i++) {
3374                              int iColumn = sort[i];
3375                              model2->setStatus(iColumn, basic);
3376                         }
3377                         delete [] sort;
3378                         delete [] dsort;
3379                         // model2->allSlackBasis();
3380                         if (gap < 1.0e-3 * static_cast<double> (numberRows + numberColumns)) {
3381                              if (saveUpper) {
3382                                   int numberRows = model2->numberRows();
3383                                   int numberColumns = model2->numberColumns();
3384                                   CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3385                                   CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3386                                   CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3387                                   CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3388                                   delete [] saveLower;
3389                                   delete [] saveUpper;
3390                                   saveLower = NULL;
3391                                   saveUpper = NULL;
3392                              }
3393                              //int numberRows = model2->numberRows();
3394                              //int numberColumns = model2->numberColumns();
3395#ifdef ABC_INHERIT
3396                              model2->checkSolution(0);
3397                              printf("%d primal infeasibilities summing to %g\n",
3398                                     model2->numberPrimalInfeasibilities(),
3399                                     model2->sumPrimalInfeasibilities());
3400                              model2->dealWithAbc(1,1);
3401                         }
3402                    }
3403#else
3404                              // just primal values pass
3405                              double saveScale = model2->objectiveScale();
3406                              model2->setObjectiveScale(1.0e-3);
3407                              model2->primal(2);
3408                              model2->setObjectiveScale(saveScale);
3409                              // save primal solution and copy back dual
3410                              CoinMemcpyN(model2->primalRowSolution(),
3411                                          numberRows, rowPrimal);
3412                              CoinMemcpyN(rowDual,
3413                                          numberRows, model2->dualRowSolution());
3414                              CoinMemcpyN(model2->primalColumnSolution(),
3415                                          numberColumns, columnPrimal);
3416                              CoinMemcpyN(columnDual,
3417                                          numberColumns, model2->dualColumnSolution());
3418                              //model2->primal(1);
3419                              // clean up reduced costs and flag variables
3420                              {
3421                                   double * dj = model2->dualColumnSolution();
3422                                   double * cost = model2->objective();
3423                                   double * saveCost = new double[numberColumns];
3424                                   CoinMemcpyN(cost, numberColumns, saveCost);
3425                                   double * saveLower = new double[numberColumns];
3426                                   double * lower = model2->columnLower();
3427                                   CoinMemcpyN(lower, numberColumns, saveLower);
3428                                   double * saveUpper = new double[numberColumns];
3429                                   double * upper = model2->columnUpper();
3430                                   CoinMemcpyN(upper, numberColumns, saveUpper);
3431                                   int i;
3432                                   double tolerance = 10.0 * dualTolerance_;
3433                                   for ( i = 0; i < numberColumns; i++) {
3434                                        if (model2->getStatus(i) == basic) {
3435                                             dj[i] = 0.0;
3436                                        } else if (model2->getStatus(i) == atLowerBound) {
3437                                             if (optimizationDirection_ * dj[i] < tolerance) {
3438                                                  if (optimizationDirection_ * dj[i] < 0.0) {
3439                                                       //if (dj[i]<-1.0e-3)
3440                                                       //printf("bad dj at lb %d %g\n",i,dj[i]);
3441                                                       cost[i] -= dj[i];
3442                                                       dj[i] = 0.0;
3443                                                  }
3444                                             } else {
3445                                                  upper[i] = lower[i];
3446                                             }
3447                                        } else if (model2->getStatus(i) == atUpperBound) {
3448                                             if (optimizationDirection_ * dj[i] > tolerance) {
3449                                                  if (optimizationDirection_ * dj[i] > 0.0) {
3450                                                       //if (dj[i]>1.0e-3)
3451                                                       //printf("bad dj at ub %d %g\n",i,dj[i]);
3452                                                       cost[i] -= dj[i];
3453                                                       dj[i] = 0.0;
3454                                                  }
3455                                             } else {
3456                                                  lower[i] = upper[i];
3457                                             }
3458                                        }
3459                                   }
3460                                   // just dual values pass
3461                                   //model2->setLogLevel(63);
3462                                   //model2->setFactorizationFrequency(1);
3463                                   model2->dual(2);
3464                                   CoinMemcpyN(saveCost, numberColumns, cost);
3465                                   delete [] saveCost;
3466                                   CoinMemcpyN(saveLower, numberColumns, lower);
3467                                   delete [] saveLower;
3468                                   CoinMemcpyN(saveUpper, numberColumns, upper);
3469                                   delete [] saveUpper;
3470                              }
3471                         }
3472                         // and finish
3473                         // move solutions
3474                         CoinMemcpyN(rowPrimal,
3475                                     numberRows, model2->primalRowSolution());
3476                         CoinMemcpyN(columnPrimal,
3477                                     numberColumns, model2->primalColumnSolution());
3478                    }
3479                    double saveScale = model2->objectiveScale();
3480                    model2->setObjectiveScale(1.0e-3);
3481                    model2->primal(2);
3482                    model2->setObjectiveScale(saveScale);
3483                    model2->primal(1);
3484#endif
3485#else
3486                    // just primal
3487#ifdef ABC_INHERIT
3488                    model2->checkSolution(0);
3489                    printf("%d primal infeasibilities summing to %g\n",
3490                           model2->numberPrimalInfeasibilities(),
3491                           model2->sumPrimalInfeasibilities());
3492                    model2->dealWithAbc(1,1);
3493#else
3494                    model2->primal(1);
3495#endif
3496                    //model2->primal(1);
3497#endif
3498               } else if (barrierStatus == 4) {
3499                    // memory problems
3500                    model2->setPerturbation(savePerturbation);
3501                    model2->createStatus();
3502                    model2->dual();
3503               } else if (maxIts && quadraticObj) {
3504                    // make sure no status left
3505                    model2->createStatus();
3506                    // solve
3507                    model2->setPerturbation(100);
3508                    model2->reducedGradient(1);
3509               }
3510          }
3511
3512          //model2->setMaximumIterations(saveMaxIts);
3513#ifdef BORROW
3514          model2->setNumberIterations(model2->numberIterations()+saveNumberIterations);
3515          delete [] rowPrimal;
3516          delete [] columnPrimal;
3517          delete [] rowDual;
3518          delete [] columnDual;
3519#endif
3520          if (extraPresolve) {
3521               pinfo2.postsolve(true);
3522               delete model2;
3523               model2 = saveModel2;
3524          }
3525          if (saveUpper) {
3526               if (!forceFixing) {
3527                    int numberRows = model2->numberRows();
3528                    int numberColumns = model2->numberColumns();
3529                    CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3530                    CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3531                    CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3532                    CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3533               }
3534               delete [] saveLower;
3535               delete [] saveUpper;
3536               saveLower = NULL;
3537               saveUpper = NULL;
3538               if (method != ClpSolve::useBarrierNoCross)
3539                    model2->primal(1);
3540          }
3541          model2->setPerturbation(savePerturbation);
3542          time2 = CoinCpuTime();
3543          timeCore = time2 - timeX;
3544          handler_->message(CLP_INTERVAL_TIMING, messages_)
3545                    << "Crossover" << timeCore << time2 - time1
3546                    << CoinMessageEol;
3547          timeX = time2;
3548#else
3549          abort();
3550#endif
3551     } else if (method == ClpSolve::notImplemented) {
3552         printf("done decomposition\n");
3553     } else {
3554          assert (method != ClpSolve::automatic); // later
3555          time2 = 0.0;
3556     }
3557     if (saveMatrix) {
3558          if (model2 == this) {
3559               // delete and replace
3560               delete model2->clpMatrix();
3561               model2->replaceMatrix(saveMatrix);
3562          } else {
3563               delete saveMatrix;
3564          }
3565     }
3566     numberIterations = model2->numberIterations();
3567     finalStatus = model2->status();
3568     int finalSecondaryStatus = model2->secondaryStatus();
3569     if (presolve == ClpSolve::presolveOn) {
3570          int saveLevel = logLevel();
3571          if ((specialOptions_ & 1024) == 0)
3572               setLogLevel(CoinMin(1, saveLevel));
3573          else
3574               setLogLevel(CoinMin(0, saveLevel));
3575          pinfo->postsolve(true);
3576          numberIterations_ = 0;
3577          delete pinfo;
3578          pinfo = NULL;
3579          factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
3580          time2 = CoinCpuTime();
3581          timePresolve += time2 - timeX;
3582          handler_->message(CLP_INTERVAL_TIMING, messages_)
3583                    << "Postsolve" << time2 - timeX << time2 - time1
3584                    << CoinMessageEol;
3585          timeX = time2;
3586          if (!presolveToFile) {
3587#if 1 //ndef ABC_INHERIT
3588               delete model2;
3589#else
3590               if (model2->abcSimplex())
3591                 delete model2->abcSimplex();
3592               else
3593                 delete model2;
3594#endif
3595          }
3596          if (interrupt)
3597               currentModel = this;
3598          // checkSolution(); already done by postSolve
3599          setLogLevel(saveLevel);
3600          int oldStatus=problemStatus_;
3601          setProblemStatus(finalStatus);
3602          setSecondaryStatus(finalSecondaryStatus);
3603          int rcode=eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
3604          if (finalStatus != 3 && rcode < 0 && (finalStatus || oldStatus == -1)) {
3605               double sumPrimal=sumPrimalInfeasibilities_;
3606               double sumDual=sumDualInfeasibilities_;
3607               // ignore some parts of solution
3608               if (finalStatus == 1) {
3609                 // infeasible
3610                 sumDual=0.0;
3611               } else if (finalStatus == 2) {
3612                 sumPrimal=0.0;
3613               }
3614               int savePerturbation = perturbation();
3615               if (savePerturbation==50)
3616                 setPerturbation(51); // small
3617               if (!finalStatus || finalStatus == 2 ||(moreSpecialOptions_ & 2) == 0 ||
3618                   fabs(sumDual)+
3619                   fabs(sumPrimal)<1.0e-3) {
3620                    if (finalStatus == 2) {
3621                      if (sumDual>1.0e-4) {
3622                         // unbounded - get feasible first
3623                         double save = optimizationDirection_;
3624                         optimizationDirection_ = 0.0;
3625                         primal(1);
3626                         optimizationDirection_ = save;
3627                      }
3628                      primal(1);
3629                    } else if (finalStatus == 1) {
3630                         dual();
3631                    } else {
3632                        if ((moreSpecialOptions_&65536)==0) {
3633                          if (numberRows_<10000||true) 
3634                            setPerturbation(100); // probably better to perturb after n its
3635                          else if (savePerturbation<100)
3636                            setPerturbation(51); // probably better to perturb after n its
3637                        }
3638#ifndef ABC_INHERIT
3639                        primal(1);
3640#else
3641                        dealWithAbc(1,2,interrupt);
3642#endif
3643                    }
3644               } else {
3645                    // just set status
3646                    problemStatus_ = finalStatus;
3647               }
3648               setPerturbation(savePerturbation);
3649               numberIterations += numberIterations_;
3650               numberIterations_ = numberIterations;
3651               finalStatus = status();
3652               time2 = CoinCpuTime();
3653               handler_->message(CLP_INTERVAL_TIMING, messages_)
3654                         << "Cleanup" << time2 - timeX << time2 - time1
3655                         << CoinMessageEol;
3656               timeX = time2;
3657          } else if (rcode >= 0) {
3658#ifdef ABC_INHERIT
3659            dealWithAbc(1,2,true);
3660#else
3661            primal(1);
3662#endif
3663          } else {
3664               secondaryStatus_ = finalSecondaryStatus;
3665          }
3666     } else if (model2 != this) {
3667          // not presolved - but different model used (sprint probably)
3668          CoinMemcpyN(model2->primalRowSolution(),
3669                      numberRows_, this->primalRowSolution());
3670          CoinMemcpyN(model2->dualRowSolution(),
3671                      numberRows_, this->dualRowSolution());
3672          CoinMemcpyN(model2->primalColumnSolution(),
3673                      numberColumns_, this->primalColumnSolution());
3674          CoinMemcpyN(model2->dualColumnSolution(),
3675                      numberColumns_, this->dualColumnSolution());
3676          CoinMemcpyN(model2->statusArray(),
3677                      numberColumns_ + numberRows_, this->statusArray());
3678          objectiveValue_ = model2->objectiveValue_;
3679          numberIterations_ = model2->numberIterations_;
3680          problemStatus_ = model2->problemStatus_;
3681          secondaryStatus_ = model2->secondaryStatus_;
3682          delete model2;
3683     }
3684     if (method != ClpSolve::useBarrierNoCross &&
3685               method != ClpSolve::useBarrier)
3686          setMaximumIterations(saveMaxIterations);
3687     std::string statusMessage[] = {"Unknown", "Optimal", "PrimalInfeasible", "DualInfeasible", "Stopped",
3688                                    "Errors", "User stopped"
3689                                   };
3690     assert (finalStatus >= -1 && finalStatus <= 5);
3691     numberIterations_ = numberIterations;
3692     handler_->message(CLP_TIMING, messages_)
3693               << statusMessage[finalStatus+1] << objectiveValue() << numberIterations << time2 - time1;
3694     handler_->printing(presolve == ClpSolve::presolveOn)
3695               << timePresolve;
3696     handler_->printing(timeIdiot != 0.0)
3697               << timeIdiot;
3698     handler_->message() << CoinMessageEol;
3699     if (interrupt)
3700          signal(SIGINT, saveSignal);
3701     perturbation_ = savePerturbation;
3702     scalingFlag_ = saveScaling;
3703     // If faking objective - put back correct one
3704     if (savedObjective) {
3705          delete objective_;
3706          objective_ = savedObjective;
3707     }
3708     if (options.getSpecialOption(1) == 2 &&
3709               options.getExtraInfo(1) > 1000000) {
3710          ClpObjective * savedObjective = objective_;
3711          // make up zero objective
3712          double * obj = new double[numberColumns_];
3713          for (int i = 0; i < numberColumns_; i++)
3714               obj[i] = 0.0;
3715          objective_ = new ClpLinearObjective(obj, numberColumns_);
3716          delete [] obj;
3717          primal(1);
3718          delete objective_;
3719          objective_ = savedObjective;
3720          finalStatus = status();
3721     }
3722     eventHandler()->event(ClpEventHandler::presolveEnd);
3723     delete pinfo;
3724     moreSpecialOptions_= saveMoreOptions;
3725#ifdef CLP_USEFUL_PRINTOUT
3726     debugInt[23]=numberIterations_;
3727#endif
3728     return finalStatus;
3729}
3730// General solve
3731int
3732ClpSimplex::initialSolve()
3733{
3734     // Default so use dual
3735     ClpSolve options;
3736     return initialSolve(options);
3737}
3738// General dual solve
3739int
3740ClpSimplex::initialDualSolve()
3741{
3742     ClpSolve options;
3743     // Use dual
3744     options.setSolveType(ClpSolve::useDual);
3745     return initialSolve(options);
3746}
3747// General primal solve
3748int
3749ClpSimplex::initialPrimalSolve()
3750{
3751     ClpSolve options;
3752     // Use primal
3753     options.setSolveType(ClpSolve::usePrimal);
3754     return initialSolve(options);
3755}
3756// barrier solve, not to be followed by crossover
3757int
3758ClpSimplex::initialBarrierNoCrossSolve()
3759{
3760     ClpSolve options;
3761     // Use primal
3762     options.setSolveType(ClpSolve::useBarrierNoCross);
3763     return initialSolve(options);
3764}
3765
3766// General barrier solve
3767int
3768ClpSimplex::initialBarrierSolve()
3769{
3770     ClpSolve options;
3771     // Use primal
3772     options.setSolveType(ClpSolve::useBarrier);
3773     return initialSolve(options);
3774}
3775
3776// Default constructor
3777ClpSolve::ClpSolve (  )
3778{
3779     method_ = automatic;
3780     presolveType_ = presolveOn;
3781     numberPasses_ = 5;
3782     int i;
3783     for (i = 0; i < 7; i++)
3784          options_[i] = 0;
3785     // say no +-1 matrix
3786     options_[3] = 1;
3787     for (i = 0; i < 7; i++)
3788          extraInfo_[i] = -1;
3789     independentOptions_[0] = 0;
3790     // But switch off slacks
3791     independentOptions_[1] = 512;
3792     // Substitute up to 3
3793     independentOptions_[2] = 3;
3794
3795}
3796// Constructor when you really know what you are doing
3797ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
3798                     int numberPasses, int options[6],
3799                     int extraInfo[6], int independentOptions[3])
3800{
3801     method_ = method;
3802     presolveType_ = presolveType;
3803     numberPasses_ = numberPasses;
3804     int i;
3805     for (i = 0; i < 6; i++)
3806          options_[i] = options[i];
3807     options_[6] = 0;
3808     for (i = 0; i < 6; i++)
3809          extraInfo_[i] = extraInfo[i];
3810     extraInfo_[6] = 0;
3811     for (i = 0; i < 3; i++)
3812          independentOptions_[i] = independentOptions[i];
3813}
3814
3815// Copy constructor.
3816ClpSolve::ClpSolve(const ClpSolve & rhs)
3817{
3818     method_ = rhs.method_;
3819     presolveType_ = rhs.presolveType_;
3820     numberPasses_ = rhs.numberPasses_;
3821     int i;
3822     for ( i = 0; i < 7; i++)
3823          options_[i] = rhs.options_[i];
3824     for ( i = 0; i < 7; i++)
3825          extraInfo_[i] = rhs.extraInfo_[i];
3826     for (i = 0; i < 3; i++)
3827          independentOptions_[i] = rhs.independentOptions_[i];
3828}
3829// Assignment operator. This copies the data
3830ClpSolve &
3831ClpSolve::operator=(const ClpSolve & rhs)
3832{
3833     if (this != &rhs) {
3834          method_ = rhs.method_;
3835          presolveType_ = rhs.presolveType_;
3836          numberPasses_ = rhs.numberPasses_;
3837          int i;
3838          for (i = 0; i < 7; i++)
3839               options_[i] = rhs.options_[i];
3840          for (i = 0; i < 7; i++)
3841               extraInfo_[i] = rhs.extraInfo_[i];
3842          for (i = 0; i < 3; i++)
3843               independentOptions_[i] = rhs.independentOptions_[i];
3844     }
3845     return *this;
3846
3847}
3848// Destructor
3849ClpSolve::~ClpSolve (  )
3850{
3851}
3852// See header file for details
3853void
3854ClpSolve::setSpecialOption(int which, int value, int extraInfo)
3855{
3856     options_[which] = value;
3857     extraInfo_[which] = extraInfo;
3858}
3859int
3860ClpSolve::getSpecialOption(int which) const
3861{
3862     return options_[which];
3863}
3864
3865// Solve types
3866void
3867ClpSolve::setSolveType(SolveType method, int /*extraInfo*/)
3868{
3869     method_ = method;
3870}
3871
3872ClpSolve::SolveType
3873ClpSolve::getSolveType()
3874{
3875     return method_;
3876}
3877
3878// Presolve types
3879void
3880ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
3881{
3882     presolveType_ = amount;
3883     numberPasses_ = extraInfo;
3884}
3885ClpSolve::PresolveType
3886ClpSolve::getPresolveType()
3887{
3888     return presolveType_;
3889}
3890// Extra info for idiot (or sprint)
3891int
3892ClpSolve::getExtraInfo(int which) const
3893{
3894     return extraInfo_[which];
3895}
3896int
3897ClpSolve::getPresolvePasses() const
3898{
3899     return numberPasses_;
3900}
3901/* Say to return at once if infeasible,
3902   default is to solve */
3903void
3904ClpSolve::setInfeasibleReturn(bool trueFalse)
3905{
3906     independentOptions_[0] = trueFalse ? 1 : 0;
3907}
3908#include <string>
3909// Generates code for above constructor
3910void
3911ClpSolve::generateCpp(FILE * fp)
3912{
3913     std::string solveType[] = {
3914          "ClpSolve::useDual",
3915          "ClpSolve::usePrimal",
3916          "ClpSolve::usePrimalorSprint",
3917          "ClpSolve::useBarrier",
3918          "ClpSolve::useBarrierNoCross",
3919          "ClpSolve::automatic",
3920          "ClpSolve::notImplemented"
3921     };
3922     std::string presolveType[] =  {
3923          "ClpSolve::presolveOn",
3924          "ClpSolve::presolveOff",
3925          "ClpSolve::presolveNumber",
3926          "ClpSolve::presolveNumberCost"
3927     };
3928     fprintf(fp, "3  ClpSolve::SolveType method = %s;\n", solveType[method_].c_str());
3929     fprintf(fp, "3  ClpSolve::PresolveType presolveType = %s;\n",
3930             presolveType[presolveType_].c_str());
3931     fprintf(fp, "3  int numberPasses = %d;\n", numberPasses_);
3932     fprintf(fp, "3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
3933             options_[0], options_[1], options_[2],
3934             options_[3], options_[4], options_[5]);
3935     fprintf(fp, "3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
3936             extraInfo_[0], extraInfo_[1], extraInfo_[2],
3937             extraInfo_[3], extraInfo_[4], extraInfo_[5]);
3938     fprintf(fp, "3  int independentOptions[] = {%d,%d,%d};\n",
3939             independentOptions_[0], independentOptions_[1], independentOptions_[2]);
3940     fprintf(fp, "3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
3941     fprintf(fp, "3                    options,extraInfo,independentOptions);\n");
3942}
3943//#############################################################################
3944#include "ClpNonLinearCost.hpp"
3945
3946ClpSimplexProgress::ClpSimplexProgress ()
3947{
3948     int i;
3949     for (i = 0; i < CLP_PROGRESS; i++) {
3950          objective_[i] = COIN_DBL_MAX;
3951          infeasibility_[i] = -1.0; // set to an impossible value
3952          realInfeasibility_[i] = COIN_DBL_MAX;
3953          numberInfeasibilities_[i] = -1;
3954          iterationNumber_[i] = -1;
3955     }
3956#ifdef CLP_PROGRESS_WEIGHT
3957     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3958          objectiveWeight_[i] = COIN_DBL_MAX;
3959          infeasibilityWeight_[i] = -1.0; // set to an impossible value
3960          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3961          numberInfeasibilitiesWeight_[i] = -1;
3962          iterationNumberWeight_[i] = -1;
3963     }
3964     drop_ = 0.0;
3965     best_ = 0.0;
3966#endif
3967     initialWeight_ = 0.0;
3968     for (i = 0; i < CLP_CYCLE; i++) {
3969          //obj_[i]=COIN_DBL_MAX;
3970          in_[i] = -1;
3971          out_[i] = -1;
3972          way_[i] = 0;
3973     }
3974     numberTimes_ = 0;
3975     numberBadTimes_ = 0;
3976     numberReallyBadTimes_ = 0;
3977     numberTimesFlagged_ = 0;
3978     model_ = NULL;
3979     oddState_ = 0;
3980}
3981
3982
3983//-----------------------------------------------------------------------------
3984
3985ClpSimplexProgress::~ClpSimplexProgress ()
3986{
3987}
3988// Copy constructor.
3989ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
3990{
3991     int i;
3992     for (i = 0; i < CLP_PROGRESS; i++) {
3993          objective_[i] = rhs.objective_[i];
3994          infeasibility_[i] = rhs.infeasibility_[i];
3995          realInfeasibility_[i] = rhs.realInfeasibility_[i];
3996          numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
3997          iterationNumber_[i] = rhs.iterationNumber_[i];
3998     }
3999#ifdef CLP_PROGRESS_WEIGHT
4000     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4001          objectiveWeight_[i] = rhs.objectiveWeight_[i];
4002          infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4003          realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4004          numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4005          iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4006     }
4007     drop_ = rhs.drop_;
4008     best_ = rhs.best_;
4009#endif
4010     initialWeight_ = rhs.initialWeight_;
4011     for (i = 0; i < CLP_CYCLE; i++) {
4012          //obj_[i]=rhs.obj_[i];
4013          in_[i] = rhs.in_[i];
4014          out_[i] = rhs.out_[i];
4015          way_[i] = rhs.way_[i];
4016     }
4017     numberTimes_ = rhs.numberTimes_;
4018     numberBadTimes_ = rhs.numberBadTimes_;
4019     numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4020     numberTimesFlagged_ = rhs.numberTimesFlagged_;
4021     model_ = rhs.model_;
4022     oddState_ = rhs.oddState_;
4023}
4024// Copy constructor.from model
4025ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model)
4026{
4027     model_ = model;
4028     reset();
4029     initialWeight_ = 0.0;
4030}
4031// Fill from model
4032void
4033ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
4034{
4035     model_ = model;
4036     reset();
4037     initialWeight_ = 0.0;
4038}
4039// Assignment operator. This copies the data
4040ClpSimplexProgress &
4041ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
4042{
4043     if (this != &rhs) {
4044          int i;
4045          for (i = 0; i < CLP_PROGRESS; i++) {
4046               objective_[i] = rhs.objective_[i];
4047               infeasibility_[i] = rhs.infeasibility_[i];
4048               realInfeasibility_[i] = rhs.realInfeasibility_[i];
4049               numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
4050               iterationNumber_[i] = rhs.iterationNumber_[i];
4051          }
4052#ifdef CLP_PROGRESS_WEIGHT
4053          for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4054               objectiveWeight_[i] = rhs.objectiveWeight_[i];
4055               infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4056               realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4057               numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4058               iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4059          }
4060          drop_ = rhs.drop_;
4061          best_ = rhs.best_;
4062#endif
4063          initialWeight_ = rhs.initialWeight_;
4064          for (i = 0; i < CLP_CYCLE; i++) {
4065               //obj_[i]=rhs.obj_[i];
4066               in_[i] = rhs.in_[i];
4067               out_[i] = rhs.out_[i];
4068               way_[i] = rhs.way_[i];
4069          }
4070          numberTimes_ = rhs.numberTimes_;
4071          numberBadTimes_ = rhs.numberBadTimes_;
4072          numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4073          numberTimesFlagged_ = rhs.numberTimesFlagged_;
4074          model_ = rhs.model_;
4075          oddState_ = rhs.oddState_;
4076     }
4077     return *this;
4078}
4079// Seems to be something odd about exact comparison of doubles on linux
4080static bool equalDouble(double value1, double value2)
4081{
4082
4083     union {
4084          double d;
4085          int i[2];
4086     } v1, v2;
4087     v1.d = value1;
4088     v2.d = value2;
4089     if (sizeof(int) * 2 == sizeof(double))
4090          return (v1.i[0] == v2.i[0] && v1.i[1] == v2.i[1]);
4091     else
4092          return (v1.i[0] == v2.i[0]);
4093}
4094int
4095ClpSimplexProgress::looping()
4096{
4097     if (!model_)
4098          return -1;
4099     double objective;
4100     if (model_->algorithm() < 0) {
4101       objective = model_->rawObjectiveValue();
4102       objective -= model_->bestPossibleImprovement();
4103     } else {
4104       objective = model_->nonLinearCost()->feasibleReportCost();
4105     }
4106     double infeasibility;
4107     double realInfeasibility = 0.0;
4108     int numberInfeasibilities;
4109     int iterationNumber = model_->numberIterations();
4110     //numberTimesFlagged_ = 0;
4111     if (model_->algorithm() < 0) {
4112          // dual
4113          infeasibility = model_->sumPrimalInfeasibilities();
4114          numberInfeasibilities = model_->numberPrimalInfeasibilities();
4115     } else {
4116          //primal
4117          infeasibility = model_->sumDualInfeasibilities();
4118          realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
4119          numberInfeasibilities = model_->numberDualInfeasibilities();
4120     }
4121     int i;
4122     int numberMatched = 0;
4123     int matched = 0;
4124     int nsame = 0;
4125     for (i = 0; i < CLP_PROGRESS; i++) {
4126          bool matchedOnObjective = equalDouble(objective, objective_[i]);
4127          bool matchedOnInfeasibility = equalDouble(infeasibility, infeasibility_[i]);
4128          bool matchedOnInfeasibilities =
4129               (numberInfeasibilities == numberInfeasibilities_[i]);
4130
4131          if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
4132               matched |= (1 << i);
4133               // Check not same iteration
4134               if (iterationNumber != iterationNumber_[i]) {
4135                    numberMatched++;
4136                    // here mainly to get over compiler bug?
4137                    if (model_->messageHandler()->logLevel() > 10)
4138                         printf("%d %d %d %d %d loop check\n", i, numberMatched,
4139                                matchedOnObjective, matchedOnInfeasibility,
4140                                matchedOnInfeasibilities);
4141               } else {
4142                    // stuck but code should notice
4143                    nsame++;
4144               }
4145          }
4146          if (i) {
4147               objective_[i-1] = objective_[i];
4148               infeasibility_[i-1] = infeasibility_[i];
4149               realInfeasibility_[i-1] = realInfeasibility_[i];
4150               numberInfeasibilities_[i-1] = numberInfeasibilities_[i];
4151               iterationNumber_[i-1] = iterationNumber_[i];
4152          }
4153     }
4154     objective_[CLP_PROGRESS-1] = objective;
4155     infeasibility_[CLP_PROGRESS-1] = infeasibility;
4156     realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
4157     numberInfeasibilities_[CLP_PROGRESS-1] = numberInfeasibilities;
4158     iterationNumber_[CLP_PROGRESS-1] = iterationNumber;
4159     if (nsame == CLP_PROGRESS)
4160          numberMatched = CLP_PROGRESS; // really stuck
4161     if (model_->progressFlag())
4162          numberMatched = 0;
4163     numberTimes_++;
4164     if (numberTimes_ < 10)
4165          numberMatched = 0;
4166     // skip if just last time as may be checking something
4167     if (matched == (1 << (CLP_PROGRESS - 1)))
4168          numberMatched = 0;
4169     if (numberMatched && model_->clpMatrix()->type() < 15) {
4170          model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
4171                    << numberMatched
4172                    << matched
4173                    << numberTimes_
4174                    << CoinMessageEol;
4175          numberBadTimes_++;
4176          if (numberBadTimes_ < 10) {
4177               // make factorize every iteration
4178               model_->forceFactorization(1);
4179               if (numberBadTimes_ < 2) {
4180                    startCheck(); // clear other loop check
4181                    if (model_->algorithm() < 0) {
4182                         // dual - change tolerance
4183                         model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
4184                         // if infeasible increase dual bound
4185                         if (model_->dualBound() < 1.0e17) {
4186                              model_->setDualBound(model_->dualBound() * 1.1);
4187                              static_cast<ClpSimplexDual *>(model_)->resetFakeBounds(0);
4188                         }
4189                    } else {
4190                         // primal - change tolerance
4191                         if (numberBadTimes_ > 3)
4192                              model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
4193                         // if infeasible increase infeasibility cost
4194                         if (model_->nonLinearCost()->numberInfeasibilities() &&
4195                                   model_->infeasibilityCost() < 1.0e17) {
4196                              model_->setInfeasibilityCost(model_->infeasibilityCost() * 1.1);
4197                         }
4198                    }
4199               } else {
4200                    // flag
4201                    int iSequence;
4202                    if (model_->algorithm() < 0) {
4203                         // dual
4204                         if (model_->dualBound() > 1.0e14)
4205                              model_->setDualBound(1.0e14);
4206                         iSequence = in_[CLP_CYCLE-1];
4207                    } else {
4208                         // primal
4209                         if (model_->infeasibilityCost() > 1.0e14)
4210                              model_->setInfeasibilityCost(1.0e14);
4211                         iSequence = out_[CLP_CYCLE-1];
4212                    }
4213                    if (iSequence >= 0) {
4214                         char x = model_->isColumn(iSequence) ? 'C' : 'R';
4215                         if (model_->messageHandler()->logLevel() >= 63)
4216                              model_->messageHandler()->message(CLP_SIMPLEX_FLAG, model_->messages())
4217                                        << x << model_->sequenceWithin(iSequence)
4218                                        << CoinMessageEol;
4219                         // if Gub then needs to be sequenceIn_
4220                         int save = model_->sequenceIn();
4221                         model_->setSequenceIn(iSequence);
4222                         model_->setFlagged(iSequence);
4223                         model_->setSequenceIn(save);
4224                         //printf("flagging %d from loop\n",iSequence);
4225                         startCheck();
4226                    } else {
4227                         // Give up
4228                         if (model_->messageHandler()->logLevel() >= 63)
4229                              printf("***** All flagged?\n");
4230                         return 4;
4231                    }
4232                    // reset
4233                    numberBadTimes_ = 2;
4234               }
4235               return -2;
4236          } else {
4237               // look at solution and maybe declare victory
4238               if (infeasibility < 1.0e-4) {
4239                    return 0;
4240               } else {
4241                    model_->messageHandler()->message(CLP_LOOP, model_->messages())
4242                              << CoinMessageEol;
4243#ifndef NDEBUG
4244                    printf("debug loop ClpSimplex A\n");
4245                    abort();
4246#endif
4247                    return 3;
4248               }
4249          }
4250     }
4251     return -1;
4252}
4253// Resets as much as possible
4254void
4255ClpSimplexProgress::reset()
4256{
4257     int i;
4258     for (i = 0; i < CLP_PROGRESS; i++) {
4259          if (model_->algorithm() >= 0)
4260               objective_[i] = COIN_DBL_MAX;
4261          else
4262               objective_[i] = -COIN_DBL_MAX;
4263          infeasibility_[i] = -1.0; // set to an impossible value
4264          realInfeasibility_[i] = COIN_DBL_MAX;
4265          numberInfeasibilities_[i] = -1;
4266          iterationNumber_[i] = -1;
4267     }
4268#ifdef CLP_PROGRESS_WEIGHT
4269     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4270          objectiveWeight_[i] = COIN_DBL_MAX;
4271          infeasibilityWeight_[i] = -1.0; // set to an impossible value
4272          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
4273          numberInfeasibilitiesWeight_[i] = -1;
4274          iterationNumberWeight_[i] = -1;
4275     }
4276     drop_ = 0.0;
4277     best_ = 0.0;
4278#endif
4279     for (i = 0; i < CLP_CYCLE; i++) {
4280          //obj_[i]=COIN_DBL_MAX;
4281          in_[i] = -1;
4282          out_[i] = -1;
4283          way_[i] = 0;
4284     }
4285     numberTimes_ = 0;
4286     numberBadTimes_ = 0;
4287     numberReallyBadTimes_ = 0;
4288     numberTimesFlagged_ = 0;
4289     oddState_ = 0;
4290}
4291// Returns previous objective (if -1) - current if (0)
4292double
4293ClpSimplexProgress::lastObjective(int back) const
4294{
4295     return objective_[CLP_PROGRESS-1-back];
4296}
4297// Returns previous infeasibility (if -1) - current if (0)
4298double
4299ClpSimplexProgress::lastInfeasibility(int back) const
4300{
4301     return realInfeasibility_[CLP_PROGRESS-1-back];
4302}
4303// Sets real primal infeasibility
4304void
4305ClpSimplexProgress::setInfeasibility(double value)
4306{
4307     for (int i = 1; i < CLP_PROGRESS; i++)
4308          realInfeasibility_[i-1] = realInfeasibility_[i];
4309     realInfeasibility_[CLP_PROGRESS-1] = value;
4310}
4311// Returns number of primal infeasibilities (if -1) - current if (0)
4312int 
4313ClpSimplexProgress::numberInfeasibilities(int back) const
4314{
4315    return numberInfeasibilities_[CLP_PROGRESS-1-back];
4316}
4317// Modify objective e.g. if dual infeasible in dual
4318void
4319ClpSimplexProgress::modifyObjective(double value)
4320{
4321     objective_[CLP_PROGRESS-1] = value;
4322}
4323// Returns previous iteration number (if -1) - current if (0)
4324int
4325ClpSimplexProgress::lastIterationNumber(int back) const
4326{
4327     return iterationNumber_[CLP_PROGRESS-1-back];
4328}
4329// clears iteration numbers (to switch off panic)
4330void
4331ClpSimplexProgress::clearIterationNumbers()
4332{
4333     for (int i = 0; i < CLP_PROGRESS; i++)
4334          iterationNumber_[i] = -1;
4335}
4336// Start check at beginning of whileIterating
4337void
4338ClpSimplexProgress::startCheck()
4339{
4340     int i;
4341     for (i = 0; i < CLP_CYCLE; i++) {
4342          //obj_[i]=COIN_DBL_MAX;
4343          in_[i] = -1;
4344          out_[i] = -1;
4345          way_[i] = 0;
4346     }
4347}
4348// Returns cycle length in whileIterating
4349int
4350ClpSimplexProgress::cycle(int in, int out, int wayIn, int wayOut)
4351{
4352     int i;
4353#if 0
4354     if (model_->numberIterations() > 206571) {
4355          printf("in %d out %d\n", in, out);
4356          for (i = 0; i < CLP_CYCLE; i++)
4357               printf("cy %d in %d out %d\n", i, in_[i], out_[i]);
4358     }
4359#endif
4360     int matched = 0;
4361     // first see if in matches any out
4362     for (i = 1; i < CLP_CYCLE; i++) {
4363          if (in == out_[i]) {
4364               // even if flip then suspicious
4365               matched = -1;
4366               break;
4367          }
4368     }
4369#if 0
4370     if (!matched || in_[0] < 0) {
4371          // can't be cycle
4372          for (i = 0; i < CLP_CYCLE - 1; i++) {
4373               //obj_[i]=obj_[i+1];
4374               in_[i] = in_[i+1];
4375               out_[i] = out_[i+1];
4376               way_[i] = way_[i+1];
4377          }
4378     } else {
4379          // possible cycle
4380          matched = 0;
4381          for (i = 0; i < CLP_CYCLE - 1; i++) {
4382               int k;
4383               char wayThis = way_[i];
4384               int inThis = in_[i];
4385               int outThis = out_[i];
4386               //double objThis = obj_[i];
4387               for(k = i + 1; k < CLP_CYCLE; k++) {
4388                    if (inThis == in_[k] && outThis == out_[k] && wayThis == way_[k]) {
4389                         int distance = k - i;
4390                         if (k + distance < CLP_CYCLE) {
4391                              // See if repeats
4392                              int j = k + distance;
4393                              if (inThis == in_[j] && outThis == out_[j] && wayThis == way_[j]) {
4394                                   matched = distance;
4395                                   break;
4396                              }
4397                         } else {
4398                              matched = distance;
4399                              break;
4400                         }
4401                    }
4402               }
4403               //obj_[i]=obj_[i+1];
4404               in_[i] = in_[i+1];
4405               out_[i] = out_[i+1];
4406               way_[i] = way_[i+1];
4407          }
4408     }
4409#else
4410     if (matched && in_[0] >= 0) {
4411          // possible cycle - only check [0] against all
4412          matched = 0;
4413          int nMatched = 0;
4414          char way0 = way_[0];
4415          int in0 = in_[0];
4416          int out0 = out_[0];
4417          //double obj0 = obj_[i];
4418          for(int k = 1; k < CLP_CYCLE - 4; k++) {
4419               if (in0 == in_[k] && out0 == out_[k] && way0 == way_[k]) {
4420                    nMatched++;
4421                    // See if repeats
4422                    int end = CLP_CYCLE - k;
4423                    int j;
4424                    for ( j = 1; j < end; j++) {
4425                         if (in_[j+k] != in_[j] || out_[j+k] != out_[j] || way_[j+k] != way_[j])
4426                              break;
4427                    }
4428                    if (j == end) {
4429                         matched = k;
4430                         break;
4431                    }
4432               }
4433          }
4434          // If three times then that is too much even if not regular
4435          if (matched <= 0 && nMatched > 1)
4436               matched = 100;
4437     }
4438     for (i = 0; i < CLP_CYCLE - 1; i++) {
4439          //obj_[i]=obj_[i+1];
4440          in_[i] = in_[i+1];
4441          out_[i] = out_[i+1];
4442          way_[i] = way_[i+1];
4443     }
4444#endif
4445     int way = 1 - wayIn + 4 * (1 - wayOut);
4446     //obj_[i]=model_->objectiveValue();
4447     in_[CLP_CYCLE-1] = in;
4448     out_[CLP_CYCLE-1] = out;
4449     way_[CLP_CYCLE-1] = static_cast<char>(way);
4450     return matched;
4451}
4452#include "CoinStructuredModel.hpp"
4453// Solve using structure of model and maybe in parallel
4454int
4455ClpSimplex::solve(CoinStructuredModel * model)
4456{
4457     // analyze structure
4458     int numberRowBlocks = model->numberRowBlocks();
4459     int numberColumnBlocks = model->numberColumnBlocks();
4460     int numberElementBlocks = model->numberElementBlocks();
4461     if (numberElementBlocks == 1) {
4462          loadProblem(*model, false);
4463          return dual();
4464     }
4465     // For now just get top level structure
4466     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4467     for (int i = 0; i < numberElementBlocks; i++) {
4468          CoinStructuredModel * subModel =
4469               dynamic_cast<CoinStructuredModel *>(model->block(i));
4470          CoinModel * thisBlock;
4471          if (subModel) {
4472               thisBlock = subModel->coinModelBlock(blockInfo[i]);
4473               model->setCoinModel(thisBlock, i);
4474          } else {
4475               thisBlock = dynamic_cast<CoinModel *>(model->block(i));
4476               assert (thisBlock);
4477               // just fill in info
4478               CoinModelBlockInfo info = CoinModelBlockInfo();
4479               int whatsSet = thisBlock->whatIsSet();
4480               info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
4481               info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
4482               info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
4483               info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
4484               info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
4485               info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
4486               // Which block
4487               int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4488               info.rowBlock = iRowBlock;
4489               int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4490               info.columnBlock = iColumnBlock;
4491               blockInfo[i] = info;
4492          }
4493     }
4494     int * rowCounts = new int [numberRowBlocks];
4495     CoinZeroN(rowCounts, numberRowBlocks);
4496     int * columnCounts = new int [numberColumnBlocks+1];
4497     CoinZeroN(columnCounts, numberColumnBlocks);
4498     int decomposeType = 0;
4499     for (int i = 0; i < numberElementBlocks; i++) {
4500          int iRowBlock = blockInfo[i].rowBlock;
4501          int iColumnBlock = blockInfo[i].columnBlock;
4502          rowCounts[iRowBlock]++;
4503          columnCounts[iColumnBlock]++;
4504     }
4505     if (numberRowBlocks == numberColumnBlocks ||
4506               numberRowBlocks == numberColumnBlocks + 1) {
4507          // could be Dantzig-Wolfe
4508          int numberG1 = 0;
4509          for (int i = 0; i < numberRowBlocks; i++) {
4510               if (rowCounts[i] > 1)
4511                    numberG1++;
4512          }
4513          bool masterColumns = (numberColumnBlocks == numberRowBlocks);
4514          if ((masterColumns && numberElementBlocks == 2 * numberRowBlocks - 1)
4515                    || (!masterColumns && numberElementBlocks == 2 * numberRowBlocks)) {
4516               if (numberG1 < 2)
4517                    decomposeType = 1;
4518          }
4519     }
4520     if (!decomposeType && (numberRowBlocks == numberColumnBlocks ||
4521                            numberRowBlocks == numberColumnBlocks - 1)) {
4522          // could be Benders
4523          int numberG1 = 0;
4524          for (int i = 0; i < numberColumnBlocks; i++) {
4525               if (columnCounts[i] > 1)
4526                    numberG1++;
4527          }
4528          bool masterRows = (numberColumnBlocks == numberRowBlocks);
4529          if ((masterRows && numberElementBlocks == 2 * numberColumnBlocks - 1)
4530                    || (!masterRows && numberElementBlocks == 2 * numberColumnBlocks)) {
4531               if (numberG1 < 2)
4532                    decomposeType = 2;
4533          }
4534     }
4535     delete [] rowCounts;
4536     delete [] columnCounts;
4537     delete [] blockInfo;
4538     // decide what to do
4539     ClpSolve options;
4540     options.setIndependentOption(2,100);
4541     switch (decomposeType) {
4542          // No good
4543     case 0:
4544          loadProblem(*model, false);
4545          return dual();
4546          // DW
4547     case 1:
4548       return solveDW(model,options);
4549          // Benders
4550     case 2:
4551       return solveBenders(model,options);
4552     }
4553     return 0; // to stop compiler warning
4554}
4555/* This loads a model from a CoinStructuredModel object - returns number of errors.
4556   If originalOrder then keep to order stored in blocks,
4557   otherwise first column/rows correspond to first block - etc.
4558   If keepSolution true and size is same as current then
4559   keeps current status and solution
4560*/
4561int
4562ClpSimplex::loadProblem (  CoinStructuredModel & coinModel,
4563                           bool originalOrder,
4564                           bool keepSolution)
4565{
4566     unsigned char * status = NULL;
4567     double * psol = NULL;
4568     double * dsol = NULL;
4569     int numberRows = coinModel.numberRows();
4570     int numberColumns = coinModel.numberColumns();
4571     int numberRowBlocks = coinModel.numberRowBlocks();
4572     int numberColumnBlocks = coinModel.numberColumnBlocks();
4573     int numberElementBlocks = coinModel.numberElementBlocks();
4574     if (status_ && numberRows_ && numberRows_ == numberRows &&
4575               numberColumns_ == numberColumns && keepSolution) {
4576          status = new unsigned char [numberRows_+numberColumns_];
4577          CoinMemcpyN(status_, numberRows_ + numberColumns_, status);
4578          psol = new double [numberRows_+numberColumns_];
4579          CoinMemcpyN(columnActivity_, numberColumns_, psol);
4580          CoinMemcpyN(rowActivity_, numberRows_, psol + numberColumns_);
4581          dsol = new double [numberRows_+numberColumns_];
4582          CoinMemcpyN(reducedCost_, numberColumns_, dsol);
4583          CoinMemcpyN(dual_, numberRows_, dsol + numberColumns_);
4584     }
4585     int returnCode = 0;
4586     double * rowLower = new double [numberRows];
4587     double * rowUpper = new double [numberRows];
4588     double * columnLower = new double [numberColumns];
4589     double * columnUpper = new double [numberColumns];
4590     double * objective = new double [numberColumns];
4591     int * integerType = new int [numberColumns];
4592     CoinBigIndex numberElements = 0;
4593     // Bases for blocks
4594     int * rowBase = new int[numberRowBlocks];
4595     CoinFillN(rowBase, numberRowBlocks, -1);
4596     // And row to put it
4597     int * whichRow = new int [numberRows+numberRowBlocks];
4598     int * columnBase = new int[numberColumnBlocks];
4599     CoinFillN(columnBase, numberColumnBlocks, -1);
4600     // And column to put it
4601     int * whichColumn = new int [numberColumns+numberColumnBlocks];
4602     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4603          CoinModel * block = coinModel.coinBlock(iBlock);
4604          numberElements += block->numberElements();
4605          //and set up elements etc
4606          double * associated = block->associatedArray();
4607          // If strings then do copies
4608          if (block->stringsExist())
4609               returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
4610                                                 objective, integerType, associated);
4611          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4612          int iRowBlock = info.rowBlock;
4613          int iColumnBlock = info.columnBlock;
4614          if (rowBase[iRowBlock] < 0) {
4615               rowBase[iRowBlock] = block->numberRows();
4616               // Save block number
4617               whichRow[numberRows+iRowBlock] = iBlock;
4618          } else {
4619               assert(rowBase[iRowBlock] == block->numberRows());
4620          }
4621          if (columnBase[iColumnBlock] < 0) {
4622               columnBase[iColumnBlock] = block->numberColumns();
4623               // Save block number
4624               whichColumn[numberColumns+iColumnBlock] = iBlock;
4625          } else {
4626               assert(columnBase[iColumnBlock] == block->numberColumns());
4627          }
4628     }
4629     // Fill arrays with defaults
4630     CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
4631     CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
4632     CoinFillN(columnLower, numberColumns, 0.0);
4633     CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
4634     CoinFillN(objective, numberColumns, 0.0);
4635     CoinFillN(integerType, numberColumns, 0);
4636     int n = 0;
4637     for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4638          int k = rowBase[iBlock];
4639          rowBase[iBlock] = n;
4640          assert (k >= 0);
4641          // block number
4642          int jBlock = whichRow[numberRows+iBlock];
4643          if (originalOrder) {
4644               memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
4645          } else {
4646               CoinIotaN(whichRow + n, k, n);
4647          }
4648          n += k;
4649     }
4650     assert (n == numberRows);
4651     n = 0;
4652     for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4653          int k = columnBase[iBlock];
4654          columnBase[iBlock] = n;
4655          assert (k >= 0);
4656          if (k) {
4657               // block number
4658               int jBlock = whichColumn[numberColumns+iBlock];
4659               if (originalOrder) {
4660                    memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
4661                           k * sizeof(int));
4662               } else {
4663                    CoinIotaN(whichColumn + n, k, n);
4664               }
4665               n += k;
4666          }
4667     }
4668     assert (n == numberColumns);
4669     bool gotIntegers = false;
4670     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4671          CoinModel * block = coinModel.coinBlock(iBlock);
4672          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4673          int iRowBlock = info.rowBlock;
4674          int iRowBase = rowBase[iRowBlock];
4675          int iColumnBlock = info.columnBlock;
4676          int iColumnBase = columnBase[iColumnBlock];
4677          if (info.rhs) {
4678               int nRows = block->numberRows();
4679               const double * lower = block->rowLowerArray();
4680               const double * upper = block->rowUpperArray();
4681               for (int i = 0; i < nRows; i++) {
4682                    int put = whichRow[i+iRowBase];
4683                    rowLower[put] = lower[i];
4684                    rowUpper[put] = upper[i];
4685               }
4686          }
4687          if (info.bounds) {
4688               int nColumns = block->numberColumns();
4689               const double * lower = block->columnLowerArray();
4690               const double * upper = block->columnUpperArray();
4691               const double * obj = block->objectiveArray();
4692               for (int i = 0; i < nColumns; i++) {
4693                    int put = whichColumn[i+iColumnBase];
4694                    columnLower[put] = lower[i];
4695                    columnUpper[put] = upper[i];
4696                    objective[put] = obj[i];
4697               }
4698          }
4699          if (info.integer) {
4700               gotIntegers = true;
4701               int nColumns = block->numberColumns();
4702               const int * type = block->integerTypeArray();
4703               for (int i = 0; i < nColumns; i++) {
4704                    int put = whichColumn[i+iColumnBase];
4705                    integerType[put] = type[i];
4706               }
4707          }
4708     }
4709     gutsOfLoadModel(numberRows, numberColumns,
4710                     columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
4711     delete [] rowLower;
4712     delete [] rowUpper;
4713     delete [] columnLower;
4714     delete [] columnUpper;
4715     delete [] objective;
4716     // Do integers if wanted
4717     if (gotIntegers) {
4718          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4719               if (integerType[iColumn])
4720                    setInteger(iColumn);
4721          }
4722     }
4723     delete [] integerType;
4724     setObjectiveOffset(coinModel.objectiveOffset());
4725     // Space for elements
4726     int * row = new int[numberElements];
4727     int * column = new int[numberElements];
4728     double * element = new double[numberElements];
4729     numberElements = 0;
4730     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4731          CoinModel * block = coinModel.coinBlock(iBlock);
4732          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4733          int iRowBlock = info.rowBlock;
4734          int iRowBase = rowBase[iRowBlock];
4735          int iColumnBlock = info.columnBlock;
4736          int iColumnBase = columnBase[iColumnBlock];
4737          if (info.rowName) {
4738               int numberItems = block->rowNames()->numberItems();
4739               assert( block->numberRows() >= numberItems);
4740               if (numberItems) {
4741                    const char *const * rowNames = block->rowNames()->names();
4742                    for (int i = 0; i < numberItems; i++) {
4743                         int put = whichRow[i+iRowBase];
4744                         std::string name = rowNames[i];
4745                         setRowName(put, name);
4746                    }
4747               }
4748          }
4749          if (info.columnName) {
4750               int numberItems = block->columnNames()->numberItems();
4751               assert( block->numberColumns() >= numberItems);
4752               if (numberItems) {
4753                    const char *const * columnNames = block->columnNames()->names();
4754                    for (int i = 0; i < numberItems; i++) {
4755                         int put = whichColumn[i+iColumnBase];
4756                         std::string name = columnNames[i];
4757                         setColumnName(put, name);
4758                    }
4759               }
4760          }
4761          if (info.matrix) {
4762               CoinPackedMatrix matrix2;
4763               const CoinPackedMatrix * matrix = block->packedMatrix();
4764               if (!matrix) {
4765                    double * associated = block->associatedArray();
4766                    block->createPackedMatrix(matrix2, associated);
4767                    matrix = &matrix2;
4768               }
4769               // get matrix data pointers
4770               const int * row2 = matrix->getIndices();
4771               const CoinBigIndex * columnStart = matrix->getVectorStarts();
4772               const double * elementByColumn = matrix->getElements();
4773               const int * columnLength = matrix->getVectorLengths();
4774               int n = matrix->getNumCols();
4775               assert (matrix->isColOrdered());
4776               for (int iColumn = 0; iColumn < n; iColumn++) {
4777                    CoinBigIndex j;
4778                    int jColumn = whichColumn[iColumn+iColumnBase];
4779                    for (j = columnStart[iColumn];
4780                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4781                         row[numberElements] = whichRow[row2[j] + iRowBase];
4782                         column[numberElements] = jColumn;
4783                         element[numberElements++] = elementByColumn[j];
4784                    }
4785               }
4786          }
4787     }
4788     delete [] whichRow;
4789     delete [] whichColumn;
4790     delete [] rowBase;
4791     delete [] columnBase;
4792     CoinPackedMatrix * matrix =
4793          new CoinPackedMatrix (true, row, column, element, numberElements);
4794     matrix_ = new ClpPackedMatrix(matrix);
4795     matrix_->setDimensions(numberRows, numberColumns);
4796     delete [] row;
4797     delete [] column;
4798     delete [] element;
4799     createStatus();
4800     if (status) {
4801          // copy back
4802          CoinMemcpyN(status, numberRows_ + numberColumns_, status_);
4803          CoinMemcpyN(psol, numberColumns_, columnActivity_);
4804          CoinMemcpyN(psol + numberColumns_, numberRows_, rowActivity_);
4805          CoinMemcpyN(dsol, numberColumns_, reducedCost_);
4806          CoinMemcpyN(dsol + numberColumns_, numberRows_, dual_);
4807          delete [] status;
4808          delete [] psol;
4809          delete [] dsol;
4810     }
4811     optimizationDirection_ = coinModel.optimizationDirection();
4812     return returnCode;
4813}
4814/*  If input negative scales objective so maximum <= -value
4815    and returns scale factor used.  If positive unscales and also
4816    redoes dual stuff
4817*/
4818double
4819ClpSimplex::scaleObjective(double value)
4820{
4821     double * obj = objective();
4822     double largest = 0.0;
4823     if (value < 0.0) {
4824          value = - value;
4825          for (int i = 0; i < numberColumns_; i++) {
4826               largest = CoinMax(largest, fabs(obj[i]));
4827          }
4828          if (largest > value) {
4829               double scaleFactor = value / largest;
4830               for (int i = 0; i < numberColumns_; i++) {
4831                    obj[i] *= scaleFactor;
4832                    reducedCost_[i] *= scaleFactor;
4833               }
4834               for (int i = 0; i < numberRows_; i++) {
4835                    dual_[i] *= scaleFactor;
4836               }
4837               largest /= value;
4838          } else {
4839               // no need
4840               largest = 1.0;
4841          }
4842     } else {
4843          // at end
4844          if (value != 1.0) {
4845               for (int i = 0; i < numberColumns_; i++) {
4846                    obj[i] *= value;
4847                    reducedCost_[i] *= value;
4848               }
4849               for (int i = 0; i < numberRows_; i++) {
4850                    dual_[i] *= value;
4851               }
4852               computeObjectiveValue();
4853          }
4854     }
4855     return largest;
4856}
4857// Solve using Dantzig-Wolfe decomposition and maybe in parallel
4858int
4859ClpSimplex::solveDW(CoinStructuredModel * model,ClpSolve & options)
4860{
4861     double time1 = CoinCpuTime();
4862     int numberColumns = model->numberColumns();
4863     int numberRowBlocks = model->numberRowBlocks();
4864     int numberColumnBlocks = model->numberColumnBlocks();
4865     int numberElementBlocks = model->numberElementBlocks();
4866     // We already have top level structure
4867     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4868     for (int i = 0; i < numberElementBlocks; i++) {
4869          CoinModel * thisBlock = model->coinBlock(i);
4870          assert (thisBlock);
4871          // just fill in info
4872          CoinModelBlockInfo info = CoinModelBlockInfo();
4873          int whatsSet = thisBlock->whatIsSet();
4874          info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
4875          info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
4876          info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
4877          info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
4878          info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
4879          info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
4880          // Which block
4881          int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4882          info.rowBlock = iRowBlock;
4883          int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4884          info.columnBlock = iColumnBlock;
4885          blockInfo[i] = info;
4886     }
4887     // make up problems
4888     int numberBlocks = numberRowBlocks - 1;
4889     // Find master rows and columns
4890     int * rowCounts = new int [numberRowBlocks];
4891     CoinZeroN(rowCounts, numberRowBlocks);
4892     int * columnCounts = new int [numberColumnBlocks+1];
4893     CoinZeroN(columnCounts, numberColumnBlocks);
4894     int iBlock;
4895     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4896          int iRowBlock = blockInfo[iBlock].rowBlock;
4897          rowCounts[iRowBlock]++;
4898          int iColumnBlock = blockInfo[iBlock].columnBlock;
4899          columnCounts[iColumnBlock]++;
4900     }
4901     int * whichBlock = new int [numberElementBlocks];
4902     int masterRowBlock = -1;
4903     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4904          if (rowCounts[iBlock] > 1) {
4905               if (masterRowBlock == -1) {
4906                    masterRowBlock = iBlock;
4907               } else {
4908                    // Can't decode
4909                    masterRowBlock = -2;
4910                    break;
4911               }
4912          }
4913     }
4914     int masterColumnBlock = -1;
4915     int kBlock = 0;
4916     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4917          int count = columnCounts[iBlock];
4918          columnCounts[iBlock] = kBlock;
4919          kBlock += count;
4920     }
4921     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4922          int iColumnBlock = blockInfo[iBlock].columnBlock;
4923          whichBlock[columnCounts[iColumnBlock]] = iBlock;
4924          columnCounts[iColumnBlock]++;
4925     }
4926     for (iBlock = numberColumnBlocks - 1; iBlock >= 0; iBlock--)
4927          columnCounts[iBlock+1] = columnCounts[iBlock];
4928     columnCounts[0] = 0;
4929     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4930          int count = columnCounts[iBlock+1] - columnCounts[iBlock];
4931          if (count == 1) {
4932               int kBlock = whichBlock[columnCounts[iBlock]];
4933               int iRowBlock = blockInfo[kBlock].rowBlock;
4934               if (iRowBlock == masterRowBlock) {
4935                    if (masterColumnBlock == -1) {
4936                         masterColumnBlock = iBlock;
4937                    } else {
4938                         // Can't decode
4939                         masterColumnBlock = -2;
4940                         break;
4941                    }
4942               }
4943          }
4944     }
4945     if (masterRowBlock < 0 || masterColumnBlock == -2) {
4946          // What now
4947          abort();
4948     }
4949     delete [] rowCounts;
4950     // create all data
4951     const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
4952     ClpSimplex * sub = new ClpSimplex [numberBlocks];
4953     ClpSimplex master;
4954     // Set offset
4955     master.setObjectiveOffset(model->objectiveOffset());
4956     bool reducePrint = logLevel()==7;
4957     if (reducePrint) {
4958       // special
4959       setLogLevel(1);
4960       master.setLogLevel(1);
4961     }
4962     kBlock = 0;
4963     int masterBlock = -1;
4964     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4965          top[kBlock] = NULL;
4966          int start = columnCounts[iBlock];
4967          int end = columnCounts[iBlock+1];
4968          assert (end - start <= 2);
4969          for (int j = start; j < end; j++) {
4970               int jBlock = whichBlock[j];
4971               int iRowBlock = blockInfo[jBlock].rowBlock;
4972               int iColumnBlock = blockInfo[jBlock].columnBlock;
4973               assert (iColumnBlock == iBlock);
4974               if (iColumnBlock != masterColumnBlock && iRowBlock == masterRowBlock) {
4975                    // top matrix
4976                    top[kBlock] = model->coinBlock(jBlock)->packedMatrix();
4977               } else {
4978                    const CoinPackedMatrix * matrix
4979                    = model->coinBlock(jBlock)->packedMatrix();
4980                    // Get pointers to arrays
4981                    const double * rowLower;
4982                    const double * rowUpper;
4983                    const double * columnLower;
4984                    const double * columnUpper;
4985                    const double * objective;
4986                    model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
4987                                 columnLower, columnUpper, objective);
4988                    if (iColumnBlock != masterColumnBlock) {
4989                         // diagonal block
4990                         sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
4991                                                 objective, rowLower, rowUpper);
4992                         if (true) {
4993                              double * lower = sub[kBlock].columnLower();
4994                              double * upper = sub[kBlock].columnUpper();
4995                              int n = sub[kBlock].numberColumns();
4996                              for (int i = 0; i < n; i++) {
4997                                   lower[i] = CoinMax(-1.0e8, lower[i]);
4998                                   upper[i] = CoinMin(1.0e8, upper[i]);
4999                              }
5000                         }
5001                         if (optimizationDirection_ < 0.0) {
5002                              double * obj = sub[kBlock].objective();
5003                              int n = sub[kBlock].numberColumns();
5004                              for (int i = 0; i < n; i++)
5005                                   obj[i] = - obj[i];
5006                         }
5007                         if (this->factorizationFrequency() == 200) {
5008                              // User did not touch preset
5009                              sub[kBlock].defaultFactorizationFrequency();
5010                         } else {
5011                              // make sure model has correct value
5012                              sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
5013                         }
5014                         sub[kBlock].setPerturbation(50);
5015                         // Set columnCounts to be diagonal block index for cleanup
5016                         columnCounts[kBlock] = jBlock;
5017                    } else {
5018                         // master
5019                         masterBlock = jBlock;
5020                         master.loadProblem(*matrix, columnLower, columnUpper,
5021                                            objective, rowLower, rowUpper);
5022                         if (optimizationDirection_ < 0.0) {
5023                              double * obj = master.objective();
5024                              int n = master.numberColumns();
5025                              for (int i = 0; i < n; i++)
5026                                   obj[i] = - obj[i];
5027                         }
5028                    }
5029               }
5030          }
5031          if (iBlock != masterColumnBlock)
5032               kBlock++;
5033     }
5034     delete [] whichBlock;
5035     delete [] blockInfo;
5036     // For now master must have been defined (does not have to have columns)
5037     assert (master.numberRows());
5038     assert (masterBlock >= 0);
5039     int numberMasterRows = master.numberRows();
5040     // Overkill in terms of space
5041     int spaceNeeded = CoinMax(numberBlocks * (numberMasterRows + 1),
5042                               2 * numberMasterRows);
5043     int * rowAdd = new int[spaceNeeded];
5044     double * elementAdd = new double[spaceNeeded];
5045     spaceNeeded = numberBlocks;
5046     int * columnAdd = new int[spaceNeeded+1];
5047     double * objective = new double[spaceNeeded];
5048     // Add in costed slacks
5049     int firstArtificial = master.numberColumns();
5050     int lastArtificial = firstArtificial;
5051     if (true) {
5052          const double * lower = master.rowLower();
5053          const double * upper = master.rowUpper();
5054          int kCol = 0;
5055          for (int iRow = 0; iRow < numberMasterRows; iRow++) {
5056               if (lower[iRow] > -1.0e10) {
5057                    rowAdd[kCol] = iRow;
5058                    elementAdd[kCol++] = 1.0;
5059               }
5060               if (upper[iRow] < 1.0e10) {
5061                    rowAdd[kCol] = iRow;
5062                    elementAdd[kCol++] = -1.0;
5063               }
5064          }
5065          if (kCol > spaceNeeded) {
5066               spaceNeeded = kCol;
5067               delete [] columnAdd;
5068               delete [] objective;
5069               columnAdd = new int[spaceNeeded+1];
5070               objective = new double[spaceNeeded];
5071          }
5072          for (int i = 0; i < kCol; i++) {
5073               columnAdd[i] = i;
5074               objective[i] = 1.0e13;
5075          }
5076          columnAdd[kCol] = kCol;
5077          master.addColumns(kCol, NULL, NULL, objective,
5078                            columnAdd, rowAdd, elementAdd);
5079          lastArtificial = master.numberColumns();
5080     }
5081     int maxPass = options.independentOption(2);
5082     if (maxPass<2)
5083       maxPass=100;
5084     int iPass;
5085     double lastObjective = 1.0e31;
5086     // Create convexity rows for proposals
5087     int numberMasterColumns = master.numberColumns();
5088     master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
5089     if (this->factorizationFrequency() == 200) {
5090          // User did not touch preset
5091          master.defaultFactorizationFrequency();
5092     } else {
5093          // make sure model has correct value
5094          master.setFactorizationFrequency(this->factorizationFrequency());
5095     }
5096     master.setPerturbation(50);
5097     // Arrays to say which block and when created
5098     int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
5099     whichBlock = new int[maximumColumns];
5100     int * when = new int[maximumColumns];
5101     int numberColumnsGenerated = numberBlocks;
5102     // fill in rhs and add in artificials
5103     {
5104          double * rowLower = master.rowLower();
5105          double * rowUpper = master.rowUpper();
5106          int iBlock;
5107          columnAdd[0] = 0;
5108          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
5109               int iRow = iBlock + numberMasterRows;;
5110               rowLower[iRow] = 1.0;
5111               rowUpper[iRow] = 1.0;
5112               rowAdd[iBlock] = iRow;
5113               elementAdd[iBlock] = 1.0;
5114               objective[iBlock] = 1.0e13;
5115               columnAdd[iBlock+1] = iBlock + 1;
5116               when[iBlock] = -1;
5117               whichBlock[iBlock] = iBlock;
5118          }
5119          master.addColumns(numberBlocks, NULL, NULL, objective,
5120                            columnAdd, rowAdd, elementAdd);
5121     }
5122     char generalPrint[200];
5123     // and resize matrix to double check clp will be happy
5124     //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
5125     //                  numberMasterColumns+numberBlocks);
5126     sprintf(generalPrint,"Time to decompose %.2f seconds",CoinCpuTime() - time1);
5127     handler_->message(CLP_GENERAL, messages_)
5128       << generalPrint
5129       << CoinMessageEol;
5130     for (iPass = 0; iPass < maxPass; iPass++) {
5131          sprintf(generalPrint,"Start of pass %d", iPass);
5132          handler_->message(CLP_GENERAL, messages_)
5133            << generalPrint
5134            << CoinMessageEol;
5135          // Solve master - may be infeasible
5136          //master.scaling(0);
5137          if (0) {
5138               master.writeMps("yy.mps");
5139          }
5140          // Correct artificials
5141          double sumArtificials = 0.0;
5142          if (iPass) {
5143               double * upper = master.columnUpper();
5144               double * solution = master.primalColumnSolution();
5145               double * obj = master.objective();
5146               sumArtificials = 0.0;
5147               for (int i = firstArtificial; i < lastArtificial; i++) {
5148                    sumArtificials += solution[i];
5149                    //assert (solution[i]>-1.0e-2);
5150                    if (solution[i] < 1.0e-6) {
5151#if 0
5152                         // Could take out
5153                         obj[i] = 0.0;
5154                         upper[i] = 0.0;
5155#else
5156                         obj[i] = 1.0e7;
5157                         upper[i] = 1.0e-1;
5158#endif
5159                         solution[i] = 0.0;
5160                         master.setColumnStatus(i, isFixed);
5161                    } else {
5162                         upper[i] = solution[i] + 1.0e-5 * (1.0 + solution[i]);
5163                    }
5164               }
5165               sprintf(generalPrint,"Sum of artificials before solve is %g", sumArtificials);
5166               handler_->message(CLP_GENERAL, messages_)
5167                 << generalPrint
5168                 << CoinMessageEol;
5169          }
5170          // scale objective to be reasonable
5171          double scaleFactor = master.scaleObjective(-1.0e9);
5172          {
5173               double * dual = master.dualRowSolution();
5174               int n = master.numberRows();
5175               memset(dual, 0, n * sizeof(double));
5176               double * solution = master.primalColumnSolution();
5177               master.clpMatrix()->times(1.0, solution, dual);
5178               double sum = 0.0;
5179               double * lower = master.rowLower();
5180               double * upper = master.rowUpper();
5181               for (int iRow = 0; iRow < n; iRow++) {
5182                    double value = dual[iRow];
5183                    if (value > upper[iRow])
5184                         sum += value - upper[iRow];
5185                    else if (value < lower[iRow])
5186                         sum -= value - lower[iRow];
5187               }
5188               printf("** suminf %g\n", sum);
5189               lower = master.columnLower();
5190               upper = master.columnUpper();
5191               n = master.numberColumns();
5192               for (int iColumn = 0; iColumn < n; iColumn++) {
5193                    double value = solution[iColumn];
5194                    if (value > upper[iColumn] + 1.0e-5)
5195                         sum += value - upper[iColumn];
5196                    else if (value < lower[iColumn] - 1.0e-5)
5197                         sum -= value - lower[iColumn];
5198               }
5199               printf("** suminf %g\n", sum);
5200          }
5201          master.primal(1);
5202          // Correct artificials
5203          sumArtificials = 0.0;
5204          {
5205               double * solution = master.primalColumnSolution();
5206               for (int i = firstArtificial; i < lastArtificial; i++) {
5207                    sumArtificials += solution[i];
5208               }
5209               printf("** Sum of artificials after solve is %g\n", sumArtificials);
5210          }
5211          master.scaleObjective(scaleFactor);
5212          int problemStatus = master.status(); // do here as can change (delcols)
5213          if (problemStatus == 2 && master.numberColumns()) {
5214            master.primal(1);
5215            if (problemStatus==2) {
5216              int numberColumns = master.numberColumns();
5217              double * lower = master.columnLower();
5218              double * upper = master.columnUpper();
5219              for (int i=0;i<numberColumns;i++) {
5220<