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

Last change on this file since 2241 was 2241, checked in by forrest, 3 years ago

fix overflow

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