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

Last change on this file since 2086 was 2086, checked in by forrest, 6 years ago

cleaner way for thread barrier

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