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

Last change on this file since 1878 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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