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

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

fix number of iterations after barrier

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 219.9 KB
Line 
1/* $Id: ClpSolve.cpp 1950 2013-05-01 13:17:25Z 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 ? 1 : 0);
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 ? 1 : 0;
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      int saveOptions=abcModel2->specialOptions();
678      if (startUp==2)
679        abcModel2->setSpecialOptions(8192|saveOptions);
680      abcModel2->ClpSimplex::doAbcPrimal(ifValuesPass);
681      abcModel2->setSpecialOptions(saveOptions);
682    }
683#if ABC_INSTRUMENT
684    double timeCpu=CoinCpuTime()-startTimeCpu;
685    double timeElapsed=CoinGetTimeOfDay()-startTimeElapsed;
686    sprintf(line,"Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
687            this->problemName().c_str(),this->numberRows(),this->numberColumns(),
688            this->getNumElements(),
689            timeCpu,timeElapsed,timeElapsed ? timeCpu/timeElapsed : 1.0,abcModel2->numberIterations());
690    handler_->message(CLP_GENERAL, messages_)
691      << line
692      << CoinMessageEol;
693#if ABC_INSTRUMENT>1
694    {
695      int n;
696      n=0;
697      for (int i=0;i<20;i++) 
698        n+= abcPricing[i];
699      printf("CCSparse pricing done %d times",n);
700      int n2=0;
701      for (int i=0;i<20;i++) 
702        n2+= abcPricingDense[i];
703      if (n2) 
704        printf(" and dense pricing done %d times\n",n2);
705      else
706        printf("\n");
707      n=0;
708      printf ("CCS");
709      for (int i=0;i<19;i++) {
710        if (abcPricing[i]) {
711          if (n==5) {
712            n=0;
713            printf("\nCCS");
714          }
715          n++;
716          printf("(%d els,%d times) ",i,abcPricing[i]);
717        }
718      }
719      if (abcPricing[19]) {
720        if (n==5) {
721          n=0;
722          printf("\nCCS");
723        }
724        n++;
725        printf("(>=19 els,%d times) ",abcPricing[19]);
726      }
727      if (n2) {
728        printf ("CCD");
729        for (int i=0;i<19;i++) {
730          if (abcPricingDense[i]) {
731            if (n==5) {
732              n=0;
733              printf("\nCCD");
734            }
735            n++;
736            int k1=(numberRows_/16)*i;;
737            int k2=CoinMin(numberRows_,k1+(numberRows_/16)-1);
738            printf("(%d-%d els,%d times) ",k1,k2,abcPricingDense[i]);
739          }
740        }
741      }
742      printf("\n");
743    }
744    instrument_print();
745#endif
746#endif
747    abcModel2->moveStatusToClp(this);
748    //ClpModel::stopPermanentArrays();
749    this->setSpecialOptions(this->specialOptions()&~65536);
750#if 0
751    this->writeBasis("a.bas",true);
752    for (int i=0;i<this->numberRows();i++)
753      printf("%d %g\n",i,this->dualRowSolution()[i]);
754    this->dual();
755    this->writeBasis("b.bas",true);
756    for (int i=0;i<this->numberRows();i++)
757      printf("%d %g\n",i,this->dualRowSolution()[i]);
758#endif
759    // switch off initialSolve flag
760    moreSpecialOptions_ &= ~16384;
761    //this->setNumberIterations(abcModel2->numberIterations()+this->numberIterations());
762    delete abcModel2;
763  }
764}
765#endif
766/** General solve algorithm which can do presolve
767    special options (bits)
768    1 - do not perturb
769    2 - do not scale
770    4 - use crash (default allslack in dual, idiot in primal)
771    8 - all slack basis in primal
772    16 - switch off interrupt handling
773    32 - do not try and make plus minus one matrix
774    64 - do not use sprint even if problem looks good
775 */
776int
777ClpSimplex::initialSolve(ClpSolve & options)
778{
779     ClpSolve::SolveType method = options.getSolveType();
780     //ClpSolve::SolveType originalMethod=method;
781     ClpSolve::PresolveType presolve = options.getPresolveType();
782     int saveMaxIterations = maximumIterations();
783     int finalStatus = -1;
784     int numberIterations = 0;
785     double time1 = CoinCpuTime();
786     double timeX = time1;
787     double time2;
788     ClpMatrixBase * saveMatrix = NULL;
789     ClpObjective * savedObjective = NULL;
790     if (!objective_ || !matrix_) {
791          // totally empty
792          handler_->message(CLP_EMPTY_PROBLEM, messages_)
793                    << 0
794                    << 0
795                    << 0
796                    << CoinMessageEol;
797          return -1;
798     } else if (!numberRows_ || !numberColumns_ || !getNumElements()) {
799          presolve = ClpSolve::presolveOff;
800     }
801     if (objective_->type() >= 2 && optimizationDirection_ == 0) {
802          // pretend linear
803          savedObjective = objective_;
804          // make up objective
805          double * obj = new double[numberColumns_];
806          for (int i = 0; i < numberColumns_; i++) {
807               double l = fabs(columnLower_[i]);
808               double u = fabs(columnUpper_[i]);
809               obj[i] = 0.0;
810               if (CoinMin(l, u) < 1.0e20) {
811                    if (l < u)
812                         obj[i] = 1.0 + randomNumberGenerator_.randomDouble() * 1.0e-2;
813                    else
814                         obj[i] = -1.0 - randomNumberGenerator_.randomDouble() * 1.0e-2;
815               }
816          }
817          objective_ = new ClpLinearObjective(obj, numberColumns_);
818          delete [] obj;
819     }
820     ClpSimplex * model2 = this;
821     bool interrupt = (options.getSpecialOption(2) == 0);
822     CoinSighandler_t saveSignal = static_cast<CoinSighandler_t> (0);
823     if (interrupt) {
824          currentModel = model2;
825          // register signal handler
826          saveSignal = signal(SIGINT, signal_handler);
827     }
828     // If no status array - set up basis
829     if (!status_)
830          allSlackBasis();
831     ClpPresolve * pinfo = new ClpPresolve();
832     pinfo->setSubstitution(options.substitution());
833     int presolveOptions = options.presolveActions();
834     bool presolveToFile = (presolveOptions & 0x40000000) != 0;
835     presolveOptions &= ~0x40000000;
836     if ((presolveOptions & 0xffff) != 0)
837          pinfo->setPresolveActions(presolveOptions);
838     // switch off singletons to slacks
839     //pinfo->setDoSingletonColumn(false); // done by bits
840     int printOptions = options.getSpecialOption(5);
841     if ((printOptions & 1) != 0)
842          pinfo->statistics();
843     double timePresolve = 0.0;
844     double timeIdiot = 0.0;
845     double timeCore = 0.0;
846     eventHandler()->event(ClpEventHandler::presolveStart);
847     int savePerturbation = perturbation_;
848     int saveScaling = scalingFlag_;
849#ifndef SLIM_CLP
850#ifndef NO_RTTI
851     if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
852          // network - switch off stuff
853          presolve = ClpSolve::presolveOff;
854     }
855#else
856     if (matrix_->type() == 11) {
857          // network - switch off stuff
858          presolve = ClpSolve::presolveOff;
859     }
860#endif
861#endif
862     if (presolve != ClpSolve::presolveOff) {
863          bool costedSlacks = false;
864#ifdef ABC_INHERIT
865          int numberPasses = 20;
866#else
867          int numberPasses = 5;
868#endif
869          if (presolve == ClpSolve::presolveNumber) {
870               numberPasses = options.getPresolvePasses();
871               presolve = ClpSolve::presolveOn;
872          } else if (presolve == ClpSolve::presolveNumberCost) {
873               numberPasses = options.getPresolvePasses();
874               presolve = ClpSolve::presolveOn;
875               costedSlacks = true;
876               // switch on singletons to slacks
877               pinfo->setDoSingletonColumn(true);
878               // gub stuff for testing
879               //pinfo->setDoGubrow(true);
880          }
881#ifndef CLP_NO_STD
882          if (presolveToFile) {
883               // PreSolve to file - not fully tested
884               printf("Presolving to file - presolve.save\n");
885               pinfo->presolvedModelToFile(*this, "presolve.save", dblParam_[ClpPresolveTolerance],
886                                          false, numberPasses);
887               model2 = this;
888          } else {
889#endif
890               model2 = pinfo->presolvedModel(*this, dblParam_[ClpPresolveTolerance],
891                                             false, numberPasses, true, costedSlacks);
892#ifndef CLP_NO_STD
893          }
894#endif
895          time2 = CoinCpuTime();
896          timePresolve = time2 - timeX;
897          handler_->message(CLP_INTERVAL_TIMING, messages_)
898                    << "Presolve" << timePresolve << time2 - time1
899                    << CoinMessageEol;
900          timeX = time2;
901          if (!model2) {
902               handler_->message(CLP_INFEASIBLE, messages_)
903                         << CoinMessageEol;
904               model2 = this;
905               eventHandler()->event(ClpEventHandler::presolveInfeasible);
906               problemStatus_ = pinfo->presolveStatus(); 
907               if (options.infeasibleReturn() || (moreSpecialOptions_ & 1) != 0) {
908                 delete pinfo;
909                    return -1;
910               }
911               presolve = ClpSolve::presolveOff;
912          } else {
913#if 0 //def ABC_INHERIT
914            {
915              AbcSimplex * abcModel2=new AbcSimplex(*model2);
916              delete model2;
917              model2=abcModel2;
918              pinfo->setPresolvedModel(model2);
919            }
920#else
921            //ClpModel::stopPermanentArrays();
922            //setSpecialOptions(specialOptions()&~65536);
923#endif
924              model2->eventHandler()->setSimplex(model2);
925              int rcode=model2->eventHandler()->event(ClpEventHandler::presolveSize);
926              // see if too big or small
927              if (rcode==2) {
928                  delete model2;
929                 delete pinfo;
930                  return -2;
931              } else if (rcode==3) {
932                  delete model2;
933                 delete pinfo;
934                  return -3;
935              }
936          }
937          model2->setMoreSpecialOptions(model2->moreSpecialOptions()&(~1024));
938          model2->eventHandler()->setSimplex(model2);
939          // We may be better off using original (but if dual leave because of bounds)
940          if (presolve != ClpSolve::presolveOff &&
941                    numberRows_ < 1.01 * model2->numberRows_ && numberColumns_ < 1.01 * model2->numberColumns_
942                    && model2 != this) {
943               if(method != ClpSolve::useDual ||
944                         (numberRows_ == model2->numberRows_ && numberColumns_ == model2->numberColumns_)) {
945                    delete model2;
946                    model2 = this;
947                    presolve = ClpSolve::presolveOff;
948               }
949          }
950     }
951     if (interrupt)
952          currentModel = model2;
953     // For below >0 overrides
954     // 0 means no, -1 means maybe
955     int doIdiot = 0;
956     int doCrash = 0;
957     int doSprint = 0;
958     int doSlp = 0;
959     int primalStartup = 1;
960     model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve);
961     bool tryItSave = false;
962     // switch to primal from automatic if just one cost entry
963     if (method == ClpSolve::automatic && model2->numberColumns() > 5000 &&
964               (specialOptions_ & 1024) != 0) {
965          int numberColumns = model2->numberColumns();
966          int numberRows = model2->numberRows();
967          const double * obj = model2->objective();
968          int nNon = 0;
969          for (int i = 0; i < numberColumns; i++) {
970               if (obj[i])
971                    nNon++;
972          }
973          if (nNon == 1) {
974#ifdef COIN_DEVELOP
975               printf("Forcing primal\n");
976#endif
977               method = ClpSolve::usePrimal;
978               tryItSave = numberRows > 200 && numberColumns > 2000 &&
979                           (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0);
980          }
981     }
982     if (method != ClpSolve::useDual && method != ClpSolve::useBarrier
983               && method != ClpSolve::useBarrierNoCross) {
984          switch (options.getSpecialOption(1)) {
985          case 0:
986               doIdiot = -1;
987               doCrash = -1;
988               doSprint = -1;
989               break;
990          case 1:
991               doIdiot = 0;
992               doCrash = 1;
993               if (options.getExtraInfo(1) > 0)
994                    doCrash = options.getExtraInfo(1);
995               doSprint = 0;
996               break;
997          case 2:
998               doIdiot = 1;
999               if (options.getExtraInfo(1) > 0)
1000                    doIdiot = options.getExtraInfo(1);
1001               doCrash = 0;
1002               doSprint = 0;
1003               break;
1004          case 3:
1005               doIdiot = 0;
1006               doCrash = 0;
1007               doSprint = 1;
1008               break;
1009          case 4:
1010               doIdiot = 0;
1011               doCrash = 0;
1012               doSprint = 0;
1013               break;
1014          case 5:
1015               doIdiot = 0;
1016               doCrash = -1;
1017               doSprint = -1;
1018               break;
1019          case 6:
1020               doIdiot = -1;
1021               doCrash = -1;
1022               doSprint = 0;
1023               break;
1024          case 7:
1025               doIdiot = -1;
1026               doCrash = 0;
1027               doSprint = -1;
1028               break;
1029          case 8:
1030               doIdiot = -1;
1031               doCrash = 0;
1032               doSprint = 0;
1033               break;
1034          case 9:
1035               doIdiot = 0;
1036               doCrash = 0;
1037               doSprint = -1;
1038               break;
1039          case 10:
1040               doIdiot = 0;
1041               doCrash = 0;
1042               doSprint = 0;
1043               if (options.getExtraInfo(1) > 0)
1044                    doSlp = options.getExtraInfo(1);
1045               break;
1046          case 11:
1047               doIdiot = 0;
1048               doCrash = 0;
1049               doSprint = 0;
1050               primalStartup = 0;
1051               break;
1052          default:
1053               abort();
1054          }
1055     } else {
1056          // Dual
1057          switch (options.getSpecialOption(0)) {
1058          case 0:
1059               doIdiot = 0;
1060               doCrash = 0;
1061               doSprint = 0;
1062               break;
1063          case 1:
1064               doIdiot = 0;
1065               doCrash = 1;
1066               if (options.getExtraInfo(0) > 0)
1067                    doCrash = options.getExtraInfo(0);
1068               doSprint = 0;
1069               break;
1070          case 2:
1071               doIdiot = -1;
1072               if (options.getExtraInfo(0) > 0)
1073                    doIdiot = options.getExtraInfo(0);
1074               doCrash = 0;
1075               doSprint = 0;
1076               break;
1077          default:
1078               abort();
1079          }
1080     }
1081#ifndef NO_RTTI
1082     ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objectiveAsObject()));
1083#else
1084     ClpQuadraticObjective * quadraticObj = NULL;
1085     if (objective_->type() == 2)
1086          quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
1087#endif
1088     // If quadratic then primal or barrier or slp
1089     if (quadraticObj) {
1090          doSprint = 0;
1091          doIdiot = 0;
1092          // off
1093          if (method == ClpSolve::useBarrier)
1094               method = ClpSolve::useBarrierNoCross;
1095          else if (method != ClpSolve::useBarrierNoCross)
1096               method = ClpSolve::usePrimal;
1097     }
1098#ifdef COIN_HAS_VOL
1099     // Save number of idiot
1100     int saveDoIdiot = doIdiot;
1101#endif
1102     // Just do this number of passes in Sprint
1103     int maxSprintPass = 100;
1104     // See if worth trying +- one matrix
1105     bool plusMinus = false;
1106     int numberElements = model2->getNumElements();
1107#ifndef SLIM_CLP
1108#ifndef NO_RTTI
1109     if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
1110          // network - switch off stuff
1111          doIdiot = 0;
1112          if (doSprint < 0)
1113               doSprint = 0;
1114     }
1115#else
1116     if (matrix_->type() == 11) {
1117          // network - switch off stuff
1118          doIdiot = 0;
1119          //doSprint=0;
1120     }
1121#endif
1122#endif
1123     int numberColumns = model2->numberColumns();
1124     int numberRows = model2->numberRows();
1125     // If not all slack basis - switch off all except sprint
1126     int numberRowsBasic = 0;
1127     int iRow;
1128     for (iRow = 0; iRow < numberRows; iRow++)
1129          if (model2->getRowStatus(iRow) == basic)
1130               numberRowsBasic++;
1131     if (numberRowsBasic < numberRows) {
1132          doIdiot = 0;
1133          doCrash = 0;
1134          //doSprint=0;
1135     }
1136     if (options.getSpecialOption(3) == 0) {
1137          if(numberElements > 100000)
1138               plusMinus = true;
1139          if(numberElements > 10000 && (doIdiot || doSprint))
1140               plusMinus = true;
1141     } else if ((specialOptions_ & 1024) != 0) {
1142          plusMinus = true;
1143     }
1144#ifndef SLIM_CLP
1145     // Statistics (+1,-1, other) - used to decide on strategy if not +-1
1146     CoinBigIndex statistics[3] = { -1, 0, 0};
1147     if (plusMinus) {
1148          saveMatrix = model2->clpMatrix();
1149#ifndef NO_RTTI
1150          ClpPackedMatrix* clpMatrix =
1151               dynamic_cast< ClpPackedMatrix*>(saveMatrix);
1152#else
1153          ClpPackedMatrix* clpMatrix = NULL;
1154          if (saveMatrix->type() == 1)
1155               clpMatrix =
1156                    static_cast< ClpPackedMatrix*>(saveMatrix);
1157#endif
1158          if (clpMatrix) {
1159               ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
1160               if (newMatrix->getIndices()) {
1161                  // CHECKME This test of specialOptions and the one above
1162                  // don't seem compatible.
1163#ifndef ABC_INHERIT
1164                    if ((specialOptions_ & 1024) == 0) {
1165                         model2->replaceMatrix(newMatrix);
1166                    } else {
1167#endif
1168                         // in integer (or abc) - just use for sprint/idiot
1169                         saveMatrix = NULL;
1170                         delete newMatrix;
1171#ifndef ABC_INHERIT
1172                    }
1173#endif
1174               } else {
1175                    handler_->message(CLP_MATRIX_CHANGE, messages_)
1176                              << "+- 1"
1177                              << CoinMessageEol;
1178                    CoinMemcpyN(newMatrix->startPositive(), 3, statistics);
1179                    saveMatrix = NULL;
1180                    plusMinus = false;
1181                    delete newMatrix;
1182               }
1183          } else {
1184               saveMatrix = NULL;
1185               plusMinus = false;
1186          }
1187     }
1188#endif
1189     if (this->factorizationFrequency() == 200) {
1190          // User did not touch preset
1191          model2->defaultFactorizationFrequency();
1192     } else if (model2 != this) {
1193          // make sure model2 has correct value
1194          model2->setFactorizationFrequency(this->factorizationFrequency());
1195     }
1196     if (method == ClpSolve::automatic) {
1197          if (doSprint == 0 && doIdiot == 0) {
1198               // off
1199               method = ClpSolve::useDual;
1200          } else {
1201               // only do primal if sprint or idiot
1202               if (doSprint > 0) {
1203                    method = ClpSolve::usePrimalorSprint;
1204               } else if (doIdiot > 0) {
1205                    method = ClpSolve::usePrimal;
1206               } else {
1207                    if (numberElements < 500000) {
1208                         // Small problem
1209                         if(numberRows * 10 > numberColumns || numberColumns < 6000
1210                                   || (numberRows * 20 > numberColumns && !plusMinus))
1211                              doSprint = 0; // switch off sprint
1212                    } else {
1213                         // larger problem
1214                         if(numberRows * 8 > numberColumns)
1215                              doSprint = 0; // switch off sprint
1216                    }
1217                    // switch off idiot or sprint if any free variable
1218                    // switch off sprint if very few with costs
1219                    int iColumn;
1220                    const double * columnLower = model2->columnLower();
1221                    const double * columnUpper = model2->columnUpper();
1222                    const double * objective = model2->objective();
1223                    int nObj=0;
1224                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1225                         if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
1226                              doSprint = 0;
1227                              doIdiot = 0;
1228                              break;
1229                         } else if (objective[iColumn]) {
1230                           nObj++;
1231                         }
1232                    }
1233                    if (nObj*10<numberColumns)
1234                      doSprint=0;
1235                    int nPasses = 0;
1236                    // look at rhs
1237                    int iRow;
1238                    double largest = 0.0;
1239                    double smallest = 1.0e30;
1240                    double largestGap = 0.0;
1241                    int numberNotE = 0;
1242                    bool notInteger = false;
1243                    for (iRow = 0; iRow < numberRows; iRow++) {
1244                         double value1 = model2->rowLower_[iRow];
1245                         if (value1 && value1 > -1.0e31) {
1246                              largest = CoinMax(largest, fabs(value1));
1247                              smallest = CoinMin(smallest, fabs(value1));
1248                              if (fabs(value1 - floor(value1 + 0.5)) > 1.0e-8) {
1249                                   notInteger = true;
1250                                   break;
1251                              }
1252                         }
1253                         double value2 = model2->rowUpper_[iRow];
1254                         if (value2 && value2 < 1.0e31) {
1255                              largest = CoinMax(largest, fabs(value2));
1256                              smallest = CoinMin(smallest, fabs(value2));
1257                              if (fabs(value2 - floor(value2 + 0.5)) > 1.0e-8) {
1258                                   notInteger = true;
1259                                   break;
1260                              }
1261                         }
1262                         // CHECKME This next bit can't be right...
1263                         if (value2 > value1) {
1264                              numberNotE++;
1265                              //if (value2 > 1.0e31 || value1 < -1.0e31)
1266                              //   largestGap = COIN_DBL_MAX;
1267                              //else
1268                              //   largestGap = value2 - value1;
1269                         }
1270                    }
1271                    bool tryIt = numberRows > 200 && numberColumns > 2000 &&
1272                                 (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0));
1273                    tryItSave = tryIt;
1274                    if (numberRows < 1000 && numberColumns < 3000)
1275                         tryIt = false;
1276                    if (notInteger)
1277                         tryIt = false;
1278                    if (largest / smallest > 10 || (largest / smallest > 2.0 && largest > 50))
1279                         tryIt = false;
1280                    if (tryIt) {
1281                         if (largest / smallest > 2.0) {
1282                              nPasses = 10 + numberColumns / 100000;
1283                              nPasses = CoinMin(nPasses, 50);
1284                              nPasses = CoinMax(nPasses, 15);
1285                              if (numberRows > 20000 && nPasses > 5) {
1286                                   // Might as well go for it
1287                                   nPasses = CoinMax(nPasses, 71);
1288                              } else if (numberRows > 2000 && nPasses > 5) {
1289                                   nPasses = CoinMax(nPasses, 50);
1290                              } else if (numberElements < 3 * numberColumns) {
1291                                   nPasses = CoinMin(nPasses, 10); // probably not worh it
1292                              }
1293                         } else if (largest / smallest > 1.01 || numberElements <= 3 * numberColumns) {
1294                              nPasses = 10 + numberColumns / 1000;
1295                              nPasses = CoinMin(nPasses, 100);
1296                              nPasses = CoinMax(nPasses, 30);
1297                              if (numberRows > 25000) {
1298                                   // Might as well go for it
1299                                   nPasses = CoinMax(nPasses, 71);
1300                              }
1301                              if (!largestGap)
1302                                   nPasses *= 2;
1303                         } else {
1304                              nPasses = 10 + numberColumns / 1000;
1305                              nPasses = CoinMin(nPasses, 200);
1306                              nPasses = CoinMax(nPasses, 100);
1307                              if (!largestGap)
1308                                   nPasses *= 2;
1309                         }
1310                    }
1311                    //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
1312                    //     numberRows,numberColumns,plusMinus ? 'Y' : 'N',
1313                    //     tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
1314                    //exit(0);
1315                    if (!tryIt || nPasses <= 5)
1316                         doIdiot = 0;
1317                    if (doSprint) {
1318                         method = ClpSolve::usePrimalorSprint;
1319                    } else if (doIdiot) {
1320                         method = ClpSolve::usePrimal;
1321                    } else {
1322                         method = ClpSolve::useDual;
1323                    }
1324               }
1325          }
1326     }
1327     if (method == ClpSolve::usePrimalorSprint) {
1328          if (doSprint < 0) {
1329               if (numberElements < 500000) {
1330                    // Small problem
1331                    if(numberRows * 10 > numberColumns || numberColumns < 6000
1332                              || (numberRows * 20 > numberColumns && !plusMinus))
1333                         method = ClpSolve::usePrimal; // switch off sprint
1334               } else {
1335                    // larger problem
1336                    if(numberRows * 8 > numberColumns)
1337                         method = ClpSolve::usePrimal; // switch off sprint
1338                    // but make lightweight
1339                    if(numberRows * 10 > numberColumns || numberColumns < 6000
1340                              || (numberRows * 20 > numberColumns && !plusMinus))
1341                         maxSprintPass = 10;
1342               }
1343          } else if (doSprint == 0) {
1344               method = ClpSolve::usePrimal; // switch off sprint
1345          }
1346     }
1347     if (method == ClpSolve::useDual) {
1348          double * saveLower = NULL;
1349          double * saveUpper = NULL;
1350          if (presolve == ClpSolve::presolveOn) {
1351               int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0);
1352               if (numberInfeasibilities) {
1353                    handler_->message(CLP_INFEASIBLE, messages_)
1354                              << CoinMessageEol;
1355                    delete model2;
1356                    model2 = this;
1357                    presolve = ClpSolve::presolveOff;
1358               }
1359          } else if (numberRows_ + numberColumns_ > 5000) {
1360               // do anyway
1361               saveLower = new double[numberRows_+numberColumns_];
1362               CoinMemcpyN(model2->columnLower(), numberColumns_, saveLower);
1363               CoinMemcpyN(model2->rowLower(), numberRows_, saveLower + numberColumns_);
1364               saveUpper = new double[numberRows_+numberColumns_];
1365               CoinMemcpyN(model2->columnUpper(), numberColumns_, saveUpper);
1366               CoinMemcpyN(model2->rowUpper(), numberRows_, saveUpper + numberColumns_);
1367               int numberInfeasibilities = model2->tightenPrimalBounds();
1368               if (numberInfeasibilities) {
1369                    handler_->message(CLP_INFEASIBLE, messages_)
1370                              << CoinMessageEol;
1371                    CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
1372                    CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
1373                    delete [] saveLower;
1374                    saveLower = NULL;
1375                    CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
1376                    CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
1377                    delete [] saveUpper;
1378                    saveUpper = NULL;
1379               }
1380          }
1381#ifndef COIN_HAS_VOL
1382          // switch off idiot and volume for now
1383          doIdiot = 0;
1384#endif
1385          // pick up number passes
1386          int nPasses = 0;
1387          int numberNotE = 0;
1388#ifndef SLIM_CLP
1389          if ((doIdiot < 0 && plusMinus) || doIdiot > 0) {
1390               // See if candidate for idiot
1391               nPasses = 0;
1392               Idiot info(*model2);
1393               // Get average number of elements per column
1394               double ratio  = static_cast<double> (numberElements) / static_cast<double> (numberColumns);
1395               // look at rhs
1396               int iRow;
1397               double largest = 0.0;
1398               double smallest = 1.0e30;
1399               for (iRow = 0; iRow < numberRows; iRow++) {
1400                    double value1 = model2->rowLower_[iRow];
1401                    if (value1 && value1 > -1.0e31) {
1402                         largest = CoinMax(largest, fabs(value1));
1403                         smallest = CoinMin(smallest, fabs(value1));
1404                    }
1405                    double value2 = model2->rowUpper_[iRow];
1406                    if (value2 && value2 < 1.0e31) {
1407                         largest = CoinMax(largest, fabs(value2));
1408                         smallest = CoinMin(smallest, fabs(value2));
1409                    }
1410                    if (value2 > value1) {
1411                         numberNotE++;
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/std::pow(10.0,(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          int saveNumberIterations = barrier.numberIterations();
2771          barrier.returnModel(*model2);
2772          double * rowPrimal = new double [numberRows];
2773          double * columnPrimal = new double [numberColumns];
2774          double * rowDual = new double [numberRows];
2775          double * columnDual = new double [numberColumns];
2776          // move solutions other way
2777          CoinMemcpyN(model2->primalRowSolution(),
2778                      numberRows, rowPrimal);
2779          CoinMemcpyN(model2->dualRowSolution(),
2780                      numberRows, rowDual);
2781          CoinMemcpyN(model2->primalColumnSolution(),
2782                      numberColumns, columnPrimal);
2783          CoinMemcpyN(model2->dualColumnSolution(),
2784                      numberColumns, columnDual);
2785#else
2786          double * rowPrimal = barrier.primalRowSolution();
2787          double * columnPrimal = barrier.primalColumnSolution();
2788          double * rowDual = barrier.dualRowSolution();
2789          double * columnDual = barrier.dualColumnSolution();
2790          // move solutions
2791          CoinMemcpyN(rowPrimal,
2792                      numberRows, model2->primalRowSolution());
2793          CoinMemcpyN(rowDual,
2794                      numberRows, model2->dualRowSolution());
2795          CoinMemcpyN(columnPrimal,
2796                      numberColumns, model2->primalColumnSolution());
2797          CoinMemcpyN(columnDual,
2798                      numberColumns, model2->dualColumnSolution());
2799#endif
2800          if (saveModel2) {
2801               // do presolve
2802               model2 = pinfo2.presolvedModel(*model2, dblParam_[ClpPresolveTolerance],
2803                                              false, 5, true);
2804               if (!model2) {
2805                    model2 = saveModel2;
2806                    saveModel2 = NULL;
2807                    int numberRows = model2->numberRows();
2808                    int numberColumns = model2->numberColumns();
2809                    CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
2810                    CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
2811                    delete [] saveLower;
2812                    CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
2813                    CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
2814                    delete [] saveUpper;
2815                    saveLower = NULL;
2816                    saveUpper = NULL;
2817               }
2818          }
2819          if (method == ClpSolve::useBarrier || barrierStatus < 0) {
2820               if (maxIts && barrierStatus < 4 && !quadraticObj) {
2821                    //printf("***** crossover - needs more thought on difficult models\n");
2822#if SAVEIT==1
2823                    model2->ClpSimplex::saveModel("xx.save");
2824#endif
2825                    // make sure no status left
2826                    model2->createStatus();
2827                    // solve
2828                    if (!forceFixing)
2829                         model2->setPerturbation(100);
2830                    if (model2->factorizationFrequency() == 200) {
2831                         // User did not touch preset
2832                         model2->defaultFactorizationFrequency();
2833                    }
2834#if 1 //ndef ABC_INHERIT //#if 1
2835                    // throw some into basis
2836                    if(!forceFixing) {
2837                         int numberRows = model2->numberRows();
2838                         int numberColumns = model2->numberColumns();
2839                         double * dsort = new double[numberColumns];
2840                         int * sort = new int[numberColumns];
2841                         int n = 0;
2842                         const double * columnLower = model2->columnLower();
2843                         const double * columnUpper = model2->columnUpper();
2844                         double * primalSolution = model2->primalColumnSolution();
2845                         const double * dualSolution = model2->dualColumnSolution();
2846                         double tolerance = 10.0 * primalTolerance_;
2847                         int i;
2848                         for ( i = 0; i < numberRows; i++)
2849                              model2->setRowStatus(i, superBasic);
2850                         for ( i = 0; i < numberColumns; i++) {
2851                              double distance = CoinMin(columnUpper[i] - primalSolution[i],
2852                                                        primalSolution[i] - columnLower[i]);
2853                              if (distance > tolerance) {
2854                                   if (fabs(dualSolution[i]) < 1.0e-5)
2855                                        distance *= 100.0;
2856                                   dsort[n] = -distance;
2857                                   sort[n++] = i;
2858                                   model2->setStatus(i, superBasic);
2859                              } else if (distance > primalTolerance_) {
2860                                   model2->setStatus(i, superBasic);
2861                              } else if (primalSolution[i] <= columnLower[i] + primalTolerance_) {
2862                                   model2->setStatus(i, atLowerBound);
2863                                   primalSolution[i] = columnLower[i];
2864                              } else {
2865                                   model2->setStatus(i, atUpperBound);
2866                                   primalSolution[i] = columnUpper[i];
2867                              }
2868                         }
2869                         CoinSort_2(dsort, dsort + n, sort);
2870                         n = CoinMin(numberRows, n);
2871                         for ( i = 0; i < n; i++) {
2872                              int iColumn = sort[i];
2873                              model2->setStatus(iColumn, basic);
2874                         }
2875                         delete [] sort;
2876                         delete [] dsort;
2877                         // model2->allSlackBasis();
2878                         if (gap < 1.0e-3 * static_cast<double> (numberRows + numberColumns)) {
2879                              if (saveUpper) {
2880                                   int numberRows = model2->numberRows();
2881                                   int numberColumns = model2->numberColumns();
2882                                   CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
2883                                   CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
2884                                   CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
2885                                   CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
2886                                   delete [] saveLower;
2887                                   delete [] saveUpper;
2888                                   saveLower = NULL;
2889                                   saveUpper = NULL;
2890                              }
2891                              int numberRows = model2->numberRows();
2892                              int numberColumns = model2->numberColumns();
2893#ifdef ABC_INHERIT
2894                              model2->checkSolution(0);
2895                              printf("%d primal infeasibilities summing to %g\n",
2896                                     model2->numberPrimalInfeasibilities(),
2897                                     model2->sumPrimalInfeasibilities());
2898                              model2->dealWithAbc(1,1);
2899                         }
2900                    }
2901#else
2902                              // just primal values pass
2903                              double saveScale = model2->objectiveScale();
2904                              model2->setObjectiveScale(1.0e-3);
2905                              model2->primal(2);
2906                              model2->setObjectiveScale(saveScale);
2907                              // save primal solution and copy back dual
2908                              CoinMemcpyN(model2->primalRowSolution(),
2909                                          numberRows, rowPrimal);
2910                              CoinMemcpyN(rowDual,
2911                                          numberRows, model2->dualRowSolution());
2912                              CoinMemcpyN(model2->primalColumnSolution(),
2913                                          numberColumns, columnPrimal);
2914                              CoinMemcpyN(columnDual,
2915                                          numberColumns, model2->dualColumnSolution());
2916                              //model2->primal(1);
2917                              // clean up reduced costs and flag variables
2918                              {
2919                                   double * dj = model2->dualColumnSolution();
2920                                   double * cost = model2->objective();
2921                                   double * saveCost = new double[numberColumns];
2922                                   CoinMemcpyN(cost, numberColumns, saveCost);
2923                                   double * saveLower = new double[numberColumns];
2924                                   double * lower = model2->columnLower();
2925                                   CoinMemcpyN(lower, numberColumns, saveLower);
2926                                   double * saveUpper = new double[numberColumns];
2927                                   double * upper = model2->columnUpper();
2928                                   CoinMemcpyN(upper, numberColumns, saveUpper);
2929                                   int i;
2930                                   double tolerance = 10.0 * dualTolerance_;
2931                                   for ( i = 0; i < numberColumns; i++) {
2932                                        if (model2->getStatus(i) == basic) {
2933                                             dj[i] = 0.0;
2934                                        } else if (model2->getStatus(i) == atLowerBound) {
2935                                             if (optimizationDirection_ * dj[i] < tolerance) {
2936                                                  if (optimizationDirection_ * dj[i] < 0.0) {
2937                                                       //if (dj[i]<-1.0e-3)
2938                                                       //printf("bad dj at lb %d %g\n",i,dj[i]);
2939                                                       cost[i] -= dj[i];
2940                                                       dj[i] = 0.0;
2941                                                  }
2942                                             } else {
2943                                                  upper[i] = lower[i];
2944                                             }
2945                                        } else if (model2->getStatus(i) == atUpperBound) {
2946                                             if (optimizationDirection_ * dj[i] > tolerance) {
2947                                                  if (optimizationDirection_ * dj[i] > 0.0) {
2948                                                       //if (dj[i]>1.0e-3)
2949                                                       //printf("bad dj at ub %d %g\n",i,dj[i]);
2950                                                       cost[i] -= dj[i];
2951                                                       dj[i] = 0.0;
2952                                                  }
2953                                             } else {
2954                                                  lower[i] = upper[i];
2955                                             }
2956                                        }
2957                                   }
2958                                   // just dual values pass
2959                                   //model2->setLogLevel(63);
2960                                   //model2->setFactorizationFrequency(1);
2961                                   model2->dual(2);
2962                                   CoinMemcpyN(saveCost, numberColumns, cost);
2963                                   delete [] saveCost;
2964                                   CoinMemcpyN(saveLower, numberColumns, lower);
2965                                   delete [] saveLower;
2966                                   CoinMemcpyN(saveUpper, numberColumns, upper);
2967                                   delete [] saveUpper;
2968                              }
2969                         }
2970                         // and finish
2971                         // move solutions
2972                         CoinMemcpyN(rowPrimal,
2973                                     numberRows, model2->primalRowSolution());
2974                         CoinMemcpyN(columnPrimal,
2975                                     numberColumns, model2->primalColumnSolution());
2976                    }
2977                    double saveScale = model2->objectiveScale();
2978                    model2->setObjectiveScale(1.0e-3);
2979                    model2->primal(2);
2980                    model2->setObjectiveScale(saveScale);
2981                    model2->primal(1);
2982#endif
2983#else
2984                    // just primal
2985#ifdef ABC_INHERIT
2986                    model2->checkSolution(0);
2987                    printf("%d primal infeasibilities summing to %g\n",
2988                           model2->numberPrimalInfeasibilities(),
2989                           model2->sumPrimalInfeasibilities());
2990                    model2->dealWithAbc(1,1);
2991#else
2992                    model2->primal(1);
2993#endif
2994                    //model2->primal(1);
2995#endif
2996               } else if (barrierStatus == 4) {
2997                    // memory problems
2998                    model2->setPerturbation(savePerturbation);
2999                    model2->createStatus();
3000                    model2->dual();
3001               } else if (maxIts && quadraticObj) {
3002                    // make sure no status left
3003                    model2->createStatus();
3004                    // solve
3005                    model2->setPerturbation(100);
3006                    model2->reducedGradient(1);
3007               }
3008          }
3009
3010          //model2->setMaximumIterations(saveMaxIts);
3011#ifdef BORROW
3012          model2->setNumberIterations(model2->numberIterations()+saveNumberIterations);
3013          delete [] rowPrimal;
3014          delete [] columnPrimal;
3015          delete [] rowDual;
3016          delete [] columnDual;
3017#endif
3018          if (extraPresolve) {
3019               pinfo2.postsolve(true);
3020               delete model2;
3021               model2 = saveModel2;
3022          }
3023          if (saveUpper) {
3024               if (!forceFixing) {
3025                    int numberRows = model2->numberRows();
3026                    int numberColumns = model2->numberColumns();
3027                    CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3028                    CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3029                    CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3030                    CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3031               }
3032               delete [] saveLower;
3033               delete [] saveUpper;
3034               saveLower = NULL;
3035               saveUpper = NULL;
3036               if (method != ClpSolve::useBarrierNoCross)
3037                    model2->primal(1);
3038          }
3039          model2->setPerturbation(savePerturbation);
3040          time2 = CoinCpuTime();
3041          timeCore = time2 - timeX;
3042          handler_->message(CLP_INTERVAL_TIMING, messages_)
3043                    << "Crossover" << timeCore << time2 - time1
3044                    << CoinMessageEol;
3045          timeX = time2;
3046#else
3047          abort();
3048#endif
3049     } else {
3050          assert (method != ClpSolve::automatic); // later
3051          time2 = 0.0;
3052     }
3053     if (saveMatrix) {
3054          if (model2 == this) {
3055               // delete and replace
3056               delete model2->clpMatrix();
3057               model2->replaceMatrix(saveMatrix);
3058          } else {
3059               delete saveMatrix;
3060          }
3061     }
3062     numberIterations = model2->numberIterations();
3063     finalStatus = model2->status();
3064     int finalSecondaryStatus = model2->secondaryStatus();
3065     if (presolve == ClpSolve::presolveOn) {
3066          int saveLevel = logLevel();
3067          if ((specialOptions_ & 1024) == 0)
3068               setLogLevel(CoinMin(1, saveLevel));
3069          else
3070               setLogLevel(CoinMin(0, saveLevel));
3071          pinfo->postsolve(true);
3072          numberIterations_ = 0;
3073          delete pinfo;
3074          pinfo = NULL;
3075          factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
3076          time2 = CoinCpuTime();
3077          timePresolve += time2 - timeX;
3078          handler_->message(CLP_INTERVAL_TIMING, messages_)
3079                    << "Postsolve" << time2 - timeX << time2 - time1
3080                    << CoinMessageEol;
3081          timeX = time2;
3082          if (!presolveToFile) {
3083#if 1 //ndef ABC_INHERIT
3084               delete model2;
3085#else
3086               if (model2->abcSimplex())
3087                 delete model2->abcSimplex();
3088               else
3089                 delete model2;
3090#endif
3091          }
3092          if (interrupt)
3093               currentModel = this;
3094          // checkSolution(); already done by postSolve
3095          setLogLevel(saveLevel);
3096          int oldStatus=problemStatus_;
3097          setProblemStatus(finalStatus);
3098          setSecondaryStatus(finalSecondaryStatus);
3099          int rcode=eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
3100          if (finalStatus != 3 && rcode < 0 && (finalStatus || oldStatus == -1)) {
3101               int savePerturbation = perturbation();
3102               if (!finalStatus || (moreSpecialOptions_ & 2) == 0 ||
3103                   fabs(sumDualInfeasibilities_)+
3104                   fabs(sumPrimalInfeasibilities_)<1.0e-3) {
3105                    if (finalStatus == 2) {
3106                         // unbounded - get feasible first
3107                         double save = optimizationDirection_;
3108                         optimizationDirection_ = 0.0;
3109                         primal(1);
3110                         optimizationDirection_ = save;
3111                         primal(1);
3112                    } else if (finalStatus == 1) {
3113                         dual();
3114                    } else {
3115                      if (numberRows_<10000)
3116                        setPerturbation(100); // probably better to perturb after n its
3117                      else if (savePerturbation<100)
3118                        setPerturbation(51); // probably better to perturb after n its
3119#ifndef ABC_INHERIT
3120                        primal(1);
3121#else
3122                        dealWithAbc(1,2,interrupt);
3123#endif
3124                    }
3125               } else {
3126                    // just set status
3127                    problemStatus_ = finalStatus;
3128               }
3129               setPerturbation(savePerturbation);
3130               numberIterations += numberIterations_;
3131               numberIterations_ = numberIterations;
3132               finalStatus = status();
3133               time2 = CoinCpuTime();
3134               handler_->message(CLP_INTERVAL_TIMING, messages_)
3135                         << "Cleanup" << time2 - timeX << time2 - time1
3136                         << CoinMessageEol;
3137               timeX = time2;
3138          } else if (rcode >= 0) {
3139#ifdef ABC_INHERIT
3140            dealWithAbc(1,2,true);
3141#else
3142            primal(1);
3143#endif
3144          } else {
3145               secondaryStatus_ = finalSecondaryStatus;
3146          }
3147     } else if (model2 != this) {
3148          // not presolved - but different model used (sprint probably)
3149          CoinMemcpyN(model2->primalRowSolution(),
3150                      numberRows_, this->primalRowSolution());
3151          CoinMemcpyN(model2->dualRowSolution(),
3152                      numberRows_, this->dualRowSolution());
3153          CoinMemcpyN(model2->primalColumnSolution(),
3154                      numberColumns_, this->primalColumnSolution());
3155          CoinMemcpyN(model2->dualColumnSolution(),
3156                      numberColumns_, this->dualColumnSolution());
3157          CoinMemcpyN(model2->statusArray(),
3158                      numberColumns_ + numberRows_, this->statusArray());
3159          objectiveValue_ = model2->objectiveValue_;
3160          numberIterations_ = model2->numberIterations_;
3161          problemStatus_ = model2->problemStatus_;
3162          secondaryStatus_ = model2->secondaryStatus_;
3163          delete model2;
3164     }
3165     if (method != ClpSolve::useBarrierNoCross &&
3166               method != ClpSolve::useBarrier)
3167          setMaximumIterations(saveMaxIterations);
3168     std::string statusMessage[] = {"Unknown", "Optimal", "PrimalInfeasible", "DualInfeasible", "Stopped",
3169                                    "Errors", "User stopped"
3170                                   };
3171     assert (finalStatus >= -1 && finalStatus <= 5);
3172     numberIterations_ = numberIterations;
3173     handler_->message(CLP_TIMING, messages_)
3174               << statusMessage[finalStatus+1] << objectiveValue() << numberIterations << time2 - time1;
3175     handler_->printing(presolve == ClpSolve::presolveOn)
3176               << timePresolve;
3177     handler_->printing(timeIdiot != 0.0)
3178               << timeIdiot;
3179     handler_->message() << CoinMessageEol;
3180     if (interrupt)
3181          signal(SIGINT, saveSignal);
3182     perturbation_ = savePerturbation;
3183     scalingFlag_ = saveScaling;
3184     // If faking objective - put back correct one
3185     if (savedObjective) {
3186          delete objective_;
3187          objective_ = savedObjective;
3188     }
3189     if (options.getSpecialOption(1) == 2 &&
3190               options.getExtraInfo(1) > 1000000) {
3191          ClpObjective * savedObjective = objective_;
3192          // make up zero objective
3193          double * obj = new double[numberColumns_];
3194          for (int i = 0; i < numberColumns_; i++)
3195               obj[i] = 0.0;
3196          objective_ = new ClpLinearObjective(obj, numberColumns_);
3197          delete [] obj;
3198          primal(1);
3199          delete objective_;
3200          objective_ = savedObjective;
3201          finalStatus = status();
3202     }
3203     eventHandler()->event(ClpEventHandler::presolveEnd);
3204     delete pinfo;
3205     return finalStatus;
3206}
3207// General solve
3208int
3209ClpSimplex::initialSolve()
3210{
3211     // Default so use dual
3212     ClpSolve options;
3213     return initialSolve(options);
3214}
3215// General dual solve
3216int
3217ClpSimplex::initialDualSolve()
3218{
3219     ClpSolve options;
3220     // Use dual
3221     options.setSolveType(ClpSolve::useDual);
3222     return initialSolve(options);
3223}
3224// General primal solve
3225int
3226ClpSimplex::initialPrimalSolve()
3227{
3228     ClpSolve options;
3229     // Use primal
3230     options.setSolveType(ClpSolve::usePrimal);
3231     return initialSolve(options);
3232}
3233// barrier solve, not to be followed by crossover
3234int
3235ClpSimplex::initialBarrierNoCrossSolve()
3236{
3237     ClpSolve options;
3238     // Use primal
3239     options.setSolveType(ClpSolve::useBarrierNoCross);
3240     return initialSolve(options);
3241}
3242
3243// General barrier solve
3244int
3245ClpSimplex::initialBarrierSolve()
3246{
3247     ClpSolve options;
3248     // Use primal
3249     options.setSolveType(ClpSolve::useBarrier);
3250     return initialSolve(options);
3251}
3252
3253// Default constructor
3254ClpSolve::ClpSolve (  )
3255{
3256     method_ = automatic;
3257     presolveType_ = presolveOn;
3258     numberPasses_ = 5;
3259     int i;
3260     for (i = 0; i < 7; i++)
3261          options_[i] = 0;
3262     // say no +-1 matrix
3263     options_[3] = 1;
3264     for (i = 0; i < 7; i++)
3265          extraInfo_[i] = -1;
3266     independentOptions_[0] = 0;
3267     // But switch off slacks
3268     independentOptions_[1] = 512;
3269     // Substitute up to 3
3270     independentOptions_[2] = 3;
3271
3272}
3273// Constructor when you really know what you are doing
3274ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
3275                     int numberPasses, int options[6],
3276                     int extraInfo[6], int independentOptions[3])
3277{
3278     method_ = method;
3279     presolveType_ = presolveType;
3280     numberPasses_ = numberPasses;
3281     int i;
3282     for (i = 0; i < 6; i++)
3283          options_[i] = options[i];
3284     options_[6] = 0;
3285     for (i = 0; i < 6; i++)
3286          extraInfo_[i] = extraInfo[i];
3287     extraInfo_[6] = 0;
3288     for (i = 0; i < 3; i++)
3289          independentOptions_[i] = independentOptions[i];
3290}
3291
3292// Copy constructor.
3293ClpSolve::ClpSolve(const ClpSolve & rhs)
3294{
3295     method_ = rhs.method_;
3296     presolveType_ = rhs.presolveType_;
3297     numberPasses_ = rhs.numberPasses_;
3298     int i;
3299     for ( i = 0; i < 7; i++)
3300          options_[i] = rhs.options_[i];
3301     for ( i = 0; i < 7; i++)
3302          extraInfo_[i] = rhs.extraInfo_[i];
3303     for (i = 0; i < 3; i++)
3304          independentOptions_[i] = rhs.independentOptions_[i];
3305}
3306// Assignment operator. This copies the data
3307ClpSolve &
3308ClpSolve::operator=(const ClpSolve & rhs)
3309{
3310     if (this != &rhs) {
3311          method_ = rhs.method_;
3312          presolveType_ = rhs.presolveType_;
3313          numberPasses_ = rhs.numberPasses_;
3314          int i;
3315          for (i = 0; i < 7; i++)
3316               options_[i] = rhs.options_[i];
3317          for (i = 0; i < 7; i++)
3318               extraInfo_[i] = rhs.extraInfo_[i];
3319          for (i = 0; i < 3; i++)
3320               independentOptions_[i] = rhs.independentOptions_[i];
3321     }
3322     return *this;
3323
3324}
3325// Destructor
3326ClpSolve::~ClpSolve (  )
3327{
3328}
3329// See header file for details
3330void
3331ClpSolve::setSpecialOption(int which, int value, int extraInfo)
3332{
3333     options_[which] = value;
3334     extraInfo_[which] = extraInfo;
3335}
3336int
3337ClpSolve::getSpecialOption(int which) const
3338{
3339     return options_[which];
3340}
3341
3342// Solve types
3343void
3344ClpSolve::setSolveType(SolveType method, int /*extraInfo*/)
3345{
3346     method_ = method;
3347}
3348
3349ClpSolve::SolveType
3350ClpSolve::getSolveType()
3351{
3352     return method_;
3353}
3354
3355// Presolve types
3356void
3357ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
3358{
3359     presolveType_ = amount;
3360     numberPasses_ = extraInfo;
3361}
3362ClpSolve::PresolveType
3363ClpSolve::getPresolveType()
3364{
3365     return presolveType_;
3366}
3367// Extra info for idiot (or sprint)
3368int
3369ClpSolve::getExtraInfo(int which) const
3370{
3371     return extraInfo_[which];
3372}
3373int
3374ClpSolve::getPresolvePasses() const
3375{
3376     return numberPasses_;
3377}
3378/* Say to return at once if infeasible,
3379   default is to solve */
3380void
3381ClpSolve::setInfeasibleReturn(bool trueFalse)
3382{
3383     independentOptions_[0] = trueFalse ? 1 : 0;
3384}
3385#include <string>
3386// Generates code for above constructor
3387void
3388ClpSolve::generateCpp(FILE * fp)
3389{
3390     std::string solveType[] = {
3391          "ClpSolve::useDual",
3392          "ClpSolve::usePrimal",
3393          "ClpSolve::usePrimalorSprint",
3394          "ClpSolve::useBarrier",
3395          "ClpSolve::useBarrierNoCross",
3396          "ClpSolve::automatic",
3397          "ClpSolve::notImplemented"
3398     };
3399     std::string presolveType[] =  {
3400          "ClpSolve::presolveOn",
3401          "ClpSolve::presolveOff",
3402          "ClpSolve::presolveNumber",
3403          "ClpSolve::presolveNumberCost"
3404     };
3405     fprintf(fp, "3  ClpSolve::SolveType method = %s;\n", solveType[method_].c_str());
3406     fprintf(fp, "3  ClpSolve::PresolveType presolveType = %s;\n",
3407             presolveType[presolveType_].c_str());
3408     fprintf(fp, "3  int numberPasses = %d;\n", numberPasses_);
3409     fprintf(fp, "3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
3410             options_[0], options_[1], options_[2],
3411             options_[3], options_[4], options_[5]);
3412     fprintf(fp, "3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
3413             extraInfo_[0], extraInfo_[1], extraInfo_[2],
3414             extraInfo_[3], extraInfo_[4], extraInfo_[5]);
3415     fprintf(fp, "3  int independentOptions[] = {%d,%d,%d};\n",
3416             independentOptions_[0], independentOptions_[1], independentOptions_[2]);
3417     fprintf(fp, "3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
3418     fprintf(fp, "3                    options,extraInfo,independentOptions);\n");
3419}
3420//#############################################################################
3421#include "ClpNonLinearCost.hpp"
3422
3423ClpSimplexProgress::ClpSimplexProgress ()
3424{
3425     int i;
3426     for (i = 0; i < CLP_PROGRESS; i++) {
3427          objective_[i] = COIN_DBL_MAX;
3428          infeasibility_[i] = -1.0; // set to an impossible value
3429          realInfeasibility_[i] = COIN_DBL_MAX;
3430          numberInfeasibilities_[i] = -1;
3431          iterationNumber_[i] = -1;
3432     }
3433#ifdef CLP_PROGRESS_WEIGHT
3434     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3435          objectiveWeight_[i] = COIN_DBL_MAX;
3436          infeasibilityWeight_[i] = -1.0; // set to an impossible value
3437          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3438          numberInfeasibilitiesWeight_[i] = -1;
3439          iterationNumberWeight_[i] = -1;
3440     }
3441     drop_ = 0.0;
3442     best_ = 0.0;
3443#endif
3444     initialWeight_ = 0.0;
3445     for (i = 0; i < CLP_CYCLE; i++) {
3446          //obj_[i]=COIN_DBL_MAX;
3447          in_[i] = -1;
3448          out_[i] = -1;
3449          way_[i] = 0;
3450     }
3451     numberTimes_ = 0;
3452     numberBadTimes_ = 0;
3453     numberReallyBadTimes_ = 0;
3454     numberTimesFlagged_ = 0;
3455     model_ = NULL;
3456     oddState_ = 0;
3457}
3458
3459
3460//-----------------------------------------------------------------------------
3461
3462ClpSimplexProgress::~ClpSimplexProgress ()
3463{
3464}
3465// Copy constructor.
3466ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
3467{
3468     int i;
3469     for (i = 0; i < CLP_PROGRESS; i++) {
3470          objective_[i] = rhs.objective_[i];
3471          infeasibility_[i] = rhs.infeasibility_[i];
3472          realInfeasibility_[i] = rhs.realInfeasibility_[i];
3473          numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
3474          iterationNumber_[i] = rhs.iterationNumber_[i];
3475     }
3476#ifdef CLP_PROGRESS_WEIGHT
3477     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3478          objectiveWeight_[i] = rhs.objectiveWeight_[i];
3479          infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
3480          realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
3481          numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
3482          iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
3483     }
3484     drop_ = rhs.drop_;
3485     best_ = rhs.best_;
3486#endif
3487     initialWeight_ = rhs.initialWeight_;
3488     for (i = 0; i < CLP_CYCLE; i++) {
3489          //obj_[i]=rhs.obj_[i];
3490          in_[i] = rhs.in_[i];
3491          out_[i] = rhs.out_[i];
3492          way_[i] = rhs.way_[i];
3493     }
3494     numberTimes_ = rhs.numberTimes_;
3495     numberBadTimes_ = rhs.numberBadTimes_;
3496     numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
3497     numberTimesFlagged_ = rhs.numberTimesFlagged_;
3498     model_ = rhs.model_;
3499     oddState_ = rhs.oddState_;
3500}
3501// Copy constructor.from model
3502ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model)
3503{
3504     model_ = model;
3505     reset();
3506     initialWeight_ = 0.0;
3507}
3508// Fill from model
3509void
3510ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
3511{
3512     model_ = model;
3513     reset();
3514     initialWeight_ = 0.0;
3515}
3516// Assignment operator. This copies the data
3517ClpSimplexProgress &
3518ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
3519{
3520     if (this != &rhs) {
3521          int i;
3522          for (i = 0; i < CLP_PROGRESS; i++) {
3523               objective_[i] = rhs.objective_[i];
3524               infeasibility_[i] = rhs.infeasibility_[i];
3525               realInfeasibility_[i] = rhs.realInfeasibility_[i];
3526               numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
3527               iterationNumber_[i] = rhs.iterationNumber_[i];
3528          }
3529#ifdef CLP_PROGRESS_WEIGHT
3530          for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3531               objectiveWeight_[i] = rhs.objectiveWeight_[i];
3532               infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
3533               realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
3534               numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
3535               iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
3536          }
3537          drop_ = rhs.drop_;
3538          best_ = rhs.best_;
3539#endif
3540          initialWeight_ = rhs.initialWeight_;
3541          for (i = 0; i < CLP_CYCLE; i++) {
3542               //obj_[i]=rhs.obj_[i];
3543               in_[i] = rhs.in_[i];
3544               out_[i] = rhs.out_[i];
3545               way_[i] = rhs.way_[i];
3546          }
3547          numberTimes_ = rhs.numberTimes_;
3548          numberBadTimes_ = rhs.numberBadTimes_;
3549          numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
3550          numberTimesFlagged_ = rhs.numberTimesFlagged_;
3551          model_ = rhs.model_;
3552          oddState_ = rhs.oddState_;
3553     }
3554     return *this;
3555}
3556// Seems to be something odd about exact comparison of doubles on linux
3557static bool equalDouble(double value1, double value2)
3558{
3559
3560     union {
3561          double d;
3562          int i[2];
3563     } v1, v2;
3564     v1.d = value1;
3565     v2.d = value2;
3566     if (sizeof(int) * 2 == sizeof(double))
3567          return (v1.i[0] == v2.i[0] && v1.i[1] == v2.i[1]);
3568     else
3569          return (v1.i[0] == v2.i[0]);
3570}
3571int
3572ClpSimplexProgress::looping()
3573{
3574     if (!model_)
3575          return -1;
3576     double objective;
3577     if (model_->algorithm() < 0) {
3578       objective = model_->rawObjectiveValue();
3579          objective -= model_->bestPossibleImprovement();
3580     } else {
3581       objective = model_->rawObjectiveValue();
3582     }
3583     double infeasibility;
3584     double realInfeasibility = 0.0;
3585     int numberInfeasibilities;
3586     int iterationNumber = model_->numberIterations();
3587     numberTimesFlagged_ = 0;
3588     if (model_->algorithm() < 0) {
3589          // dual
3590          infeasibility = model_->sumPrimalInfeasibilities();
3591          numberInfeasibilities = model_->numberPrimalInfeasibilities();
3592     } else {
3593          //primal
3594          infeasibility = model_->sumDualInfeasibilities();
3595          realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
3596          numberInfeasibilities = model_->numberDualInfeasibilities();
3597     }
3598     int i;
3599     int numberMatched = 0;
3600     int matched = 0;
3601     int nsame = 0;
3602     for (i = 0; i < CLP_PROGRESS; i++) {
3603          bool matchedOnObjective = equalDouble(objective, objective_[i]);
3604          bool matchedOnInfeasibility = equalDouble(infeasibility, infeasibility_[i]);
3605          bool matchedOnInfeasibilities =
3606               (numberInfeasibilities == numberInfeasibilities_[i]);
3607
3608          if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
3609               matched |= (1 << i);
3610               // Check not same iteration
3611               if (iterationNumber != iterationNumber_[i]) {
3612                    numberMatched++;
3613                    // here mainly to get over compiler bug?
3614                    if (model_->messageHandler()->logLevel() > 10)
3615                         printf("%d %d %d %d %d loop check\n", i, numberMatched,
3616                                matchedOnObjective, matchedOnInfeasibility,
3617                                matchedOnInfeasibilities);
3618               } else {
3619                    // stuck but code should notice
3620                    nsame++;
3621               }
3622          }
3623          if (i) {
3624               objective_[i-1] = objective_[i];
3625               infeasibility_[i-1] = infeasibility_[i];
3626               realInfeasibility_[i-1] = realInfeasibility_[i];
3627               numberInfeasibilities_[i-1] = numberInfeasibilities_[i];
3628               iterationNumber_[i-1] = iterationNumber_[i];
3629          }
3630     }
3631     objective_[CLP_PROGRESS-1] = objective;
3632     infeasibility_[CLP_PROGRESS-1] = infeasibility;
3633     realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
3634     numberInfeasibilities_[CLP_PROGRESS-1] = numberInfeasibilities;
3635     iterationNumber_[CLP_PROGRESS-1] = iterationNumber;
3636     if (nsame == CLP_PROGRESS)
3637          numberMatched = CLP_PROGRESS; // really stuck
3638     if (model_->progressFlag())
3639          numberMatched = 0;
3640     numberTimes_++;
3641     if (numberTimes_ < 10)
3642          numberMatched = 0;
3643     // skip if just last time as may be checking something
3644     if (matched == (1 << (CLP_PROGRESS - 1)))
3645          numberMatched = 0;
3646     if (numberMatched && model_->clpMatrix()->type() < 15) {
3647          model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
3648                    << numberMatched
3649                    << matched
3650                    << numberTimes_
3651                    << CoinMessageEol;
3652          numberBadTimes_++;
3653          if (numberBadTimes_ < 10) {
3654               // make factorize every iteration
3655               model_->forceFactorization(1);
3656               if (numberBadTimes_ < 2) {
3657                    startCheck(); // clear other loop check
3658                    if (model_->algorithm() < 0) {
3659                         // dual - change tolerance
3660                         model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
3661                         // if infeasible increase dual bound
3662                         if (model_->dualBound() < 1.0e17) {
3663                              model_->setDualBound(model_->dualBound() * 1.1);
3664                              static_cast<ClpSimplexDual *>(model_)->resetFakeBounds(0);
3665                         }
3666                    } else {
3667                         // primal - change tolerance
3668                         if (numberBadTimes_ > 3)
3669                              model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
3670                         // if infeasible increase infeasibility cost
3671                         if (model_->nonLinearCost()->numberInfeasibilities() &&
3672                                   model_->infeasibilityCost() < 1.0e17) {
3673                              model_->setInfeasibilityCost(model_->infeasibilityCost() * 1.1);
3674                         }
3675                    }
3676               } else {
3677                    // flag
3678                    int iSequence;
3679                    if (model_->algorithm() < 0) {
3680                         // dual
3681                         if (model_->dualBound() > 1.0e14)
3682                              model_->setDualBound(1.0e14);
3683                         iSequence = in_[CLP_CYCLE-1];
3684                    } else {
3685                         // primal
3686                         if (model_->infeasibilityCost() > 1.0e14)
3687                              model_->setInfeasibilityCost(1.0e14);
3688                         iSequence = out_[CLP_CYCLE-1];
3689                    }
3690                    if (iSequence >= 0) {
3691                         char x = model_->isColumn(iSequence) ? 'C' : 'R';
3692                         if (model_->messageHandler()->logLevel() >= 63)
3693                              model_->messageHandler()->message(CLP_SIMPLEX_FLAG, model_->messages())
3694                                        << x << model_->sequenceWithin(iSequence)
3695                                        << CoinMessageEol;
3696                         // if Gub then needs to be sequenceIn_
3697                         int save = model_->sequenceIn();
3698                         model_->setSequenceIn(iSequence);
3699                         model_->setFlagged(iSequence);
3700                         model_->setSequenceIn(save);
3701                         //printf("flagging %d from loop\n",iSequence);
3702                         startCheck();
3703                    } else {
3704                         // Give up
3705                         if (model_->messageHandler()->logLevel() >= 63)
3706                              printf("***** All flagged?\n");
3707                         return 4;
3708                    }
3709                    // reset
3710                    numberBadTimes_ = 2;
3711               }
3712               return -2;
3713          } else {
3714               // look at solution and maybe declare victory
3715               if (infeasibility < 1.0e-4) {
3716                    return 0;
3717               } else {
3718                    model_->messageHandler()->message(CLP_LOOP, model_->messages())
3719                              << CoinMessageEol;
3720#ifndef NDEBUG
3721                    printf("debug loop ClpSimplex A\n");
3722                    abort();
3723#endif
3724                    return 3;
3725               }
3726          }
3727     }
3728     return -1;
3729}
3730// Resets as much as possible
3731void
3732ClpSimplexProgress::reset()
3733{
3734     int i;
3735     for (i = 0; i < CLP_PROGRESS; i++) {
3736          if (model_->algorithm() >= 0)
3737               objective_[i] = COIN_DBL_MAX;
3738          else
3739               objective_[i] = -COIN_DBL_MAX;
3740          infeasibility_[i] = -1.0; // set to an impossible value
3741          realInfeasibility_[i] = COIN_DBL_MAX;
3742          numberInfeasibilities_[i] = -1;
3743          iterationNumber_[i] = -1;
3744     }
3745#ifdef CLP_PROGRESS_WEIGHT
3746     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3747          objectiveWeight_[i] = COIN_DBL_MAX;
3748          infeasibilityWeight_[i] = -1.0; // set to an impossible value
3749          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3750          numberInfeasibilitiesWeight_[i] = -1;
3751          iterationNumberWeight_[i] = -1;
3752     }
3753     drop_ = 0.0;
3754     best_ = 0.0;
3755#endif
3756     for (i = 0; i < CLP_CYCLE; i++) {
3757          //obj_[i]=COIN_DBL_MAX;
3758          in_[i] = -1;
3759          out_[i] = -1;
3760          way_[i] = 0;
3761     }
3762     numberTimes_ = 0;
3763     numberBadTimes_ = 0;
3764     numberReallyBadTimes_ = 0;
3765     numberTimesFlagged_ = 0;
3766     oddState_ = 0;
3767}
3768// Returns previous objective (if -1) - current if (0)
3769double
3770ClpSimplexProgress::lastObjective(int back) const
3771{
3772     return objective_[CLP_PROGRESS-1-back];
3773}
3774// Returns previous infeasibility (if -1) - current if (0)
3775double
3776ClpSimplexProgress::lastInfeasibility(int back) const
3777{
3778     return realInfeasibility_[CLP_PROGRESS-1-back];
3779}
3780// Sets real primal infeasibility
3781void
3782ClpSimplexProgress::setInfeasibility(double value)
3783{
3784     for (int i = 1; i < CLP_PROGRESS; i++)
3785          realInfeasibility_[i-1] = realInfeasibility_[i];
3786     realInfeasibility_[CLP_PROGRESS-1] = value;
3787}
3788// Modify objective e.g. if dual infeasible in dual
3789void
3790ClpSimplexProgress::modifyObjective(double value)
3791{
3792     objective_[CLP_PROGRESS-1] = value;
3793}
3794// Returns previous iteration number (if -1) - current if (0)
3795int
3796ClpSimplexProgress::lastIterationNumber(int back) const
3797{
3798     return iterationNumber_[CLP_PROGRESS-1-back];
3799}
3800// clears iteration numbers (to switch off panic)
3801void
3802ClpSimplexProgress::clearIterationNumbers()
3803{
3804     for (int i = 0; i < CLP_PROGRESS; i++)
3805          iterationNumber_[i] = -1;
3806}
3807// Start check at beginning of whileIterating
3808void
3809ClpSimplexProgress::startCheck()
3810{
3811     int i;
3812     for (i = 0; i < CLP_CYCLE; i++) {
3813          //obj_[i]=COIN_DBL_MAX;
3814          in_[i] = -1;
3815          out_[i] = -1;
3816          way_[i] = 0;
3817     }
3818}
3819// Returns cycle length in whileIterating
3820int
3821ClpSimplexProgress::cycle(int in, int out, int wayIn, int wayOut)
3822{
3823     int i;
3824#if 0
3825     if (model_->numberIterations() > 206571) {
3826          printf("in %d out %d\n", in, out);
3827          for (i = 0; i < CLP_CYCLE; i++)
3828               printf("cy %d in %d out %d\n", i, in_[i], out_[i]);
3829     }
3830#endif
3831     int matched = 0;
3832     // first see if in matches any out
3833     for (i = 1; i < CLP_CYCLE; i++) {
3834          if (in == out_[i]) {
3835               // even if flip then suspicious
3836               matched = -1;
3837               break;
3838          }
3839     }
3840#if 0
3841     if (!matched || in_[0] < 0) {
3842          // can't be cycle
3843          for (i = 0; i < CLP_CYCLE - 1; i++) {
3844               //obj_[i]=obj_[i+1];
3845               in_[i] = in_[i+1];
3846               out_[i] = out_[i+1];
3847               way_[i] = way_[i+1];
3848          }
3849     } else {
3850          // possible cycle
3851          matched = 0;
3852          for (i = 0; i < CLP_CYCLE - 1; i++) {
3853               int k;
3854               char wayThis = way_[i];
3855               int inThis = in_[i];
3856               int outThis = out_[i];
3857               //double objThis = obj_[i];
3858               for(k = i + 1; k < CLP_CYCLE; k++) {
3859                    if (inThis == in_[k] && outThis == out_[k] && wayThis == way_[k]) {
3860                         int distance = k - i;
3861                         if (k + distance < CLP_CYCLE) {
3862                              // See if repeats
3863                              int j = k + distance;
3864                              if (inThis == in_[j] && outThis == out_[j] && wayThis == way_[j]) {
3865                                   matched = distance;
3866                                   break;
3867                              }
3868                         } else {
3869                              matched = distance;
3870                              break;
3871                         }
3872                    }
3873               }
3874               //obj_[i]=obj_[i+1];
3875               in_[i] = in_[i+1];
3876               out_[i] = out_[i+1];
3877               way_[i] = way_[i+1];
3878          }
3879     }
3880#else
3881     if (matched && in_[0] >= 0) {
3882          // possible cycle - only check [0] against all
3883          matched = 0;
3884          int nMatched = 0;
3885          char way0 = way_[0];
3886          int in0 = in_[0];
3887          int out0 = out_[0];
3888          //double obj0 = obj_[i];
3889          for(int k = 1; k < CLP_CYCLE - 4; k++) {
3890               if (in0 == in_[k] && out0 == out_[k] && way0 == way_[k]) {
3891                    nMatched++;
3892                    // See if repeats
3893                    int end = CLP_CYCLE - k;
3894                    int j;
3895                    for ( j = 1; j < end; j++) {
3896                         if (in_[j+k] != in_[j] || out_[j+k] != out_[j] || way_[j+k] != way_[j])
3897                              break;
3898                    }
3899                    if (j == end) {
3900                         matched = k;
3901                         break;
3902                    }
3903               }
3904          }
3905          // If three times then that is too much even if not regular
3906          if (matched <= 0 && nMatched > 1)
3907               matched = 100;
3908     }
3909     for (i = 0; i < CLP_CYCLE - 1; i++) {
3910          //obj_[i]=obj_[i+1];
3911          in_[i] = in_[i+1];
3912          out_[i] = out_[i+1];
3913          way_[i] = way_[i+1];
3914     }
3915#endif
3916     int way = 1 - wayIn + 4 * (1 - wayOut);
3917     //obj_[i]=model_->objectiveValue();
3918     in_[CLP_CYCLE-1] = in;
3919     out_[CLP_CYCLE-1] = out;
3920     way_[CLP_CYCLE-1] = static_cast<char>(way);
3921     return matched;
3922}
3923#include "CoinStructuredModel.hpp"
3924// Solve using structure of model and maybe in parallel
3925int
3926ClpSimplex::solve(CoinStructuredModel * model)
3927{
3928     // analyze structure
3929     int numberRowBlocks = model->numberRowBlocks();
3930     int numberColumnBlocks = model->numberColumnBlocks();
3931     int numberElementBlocks = model->numberElementBlocks();
3932     if (numberElementBlocks == 1) {
3933          loadProblem(*model, false);
3934          return dual();
3935     }
3936     // For now just get top level structure
3937     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3938     for (int i = 0; i < numberElementBlocks; i++) {
3939          CoinStructuredModel * subModel =
3940               dynamic_cast<CoinStructuredModel *>(model->block(i));
3941          CoinModel * thisBlock;
3942          if (subModel) {
3943               thisBlock = subModel->coinModelBlock(blockInfo[i]);
3944               model->setCoinModel(thisBlock, i);
3945          } else {
3946               thisBlock = dynamic_cast<CoinModel *>(model->block(i));
3947               assert (thisBlock);
3948               // just fill in info
3949               CoinModelBlockInfo info = CoinModelBlockInfo();
3950               int whatsSet = thisBlock->whatIsSet();
3951               info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
3952               info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
3953               info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
3954               info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
3955               info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
3956               info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
3957               // Which block
3958               int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
3959               info.rowBlock = iRowBlock;
3960               int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
3961               info.columnBlock = iColumnBlock;
3962               blockInfo[i] = info;
3963          }
3964     }
3965     int * rowCounts = new int [numberRowBlocks];
3966     CoinZeroN(rowCounts, numberRowBlocks);
3967     int * columnCounts = new int [numberColumnBlocks+1];
3968     CoinZeroN(columnCounts, numberColumnBlocks);
3969     int decomposeType = 0;
3970     for (int i = 0; i < numberElementBlocks; i++) {
3971          int iRowBlock = blockInfo[i].rowBlock;
3972          int iColumnBlock = blockInfo[i].columnBlock;
3973          rowCounts[iRowBlock]++;
3974          columnCounts[iColumnBlock]++;
3975     }
3976     if (numberRowBlocks == numberColumnBlocks ||
3977               numberRowBlocks == numberColumnBlocks + 1) {
3978          // could be Dantzig-Wolfe
3979          int numberG1 = 0;
3980          for (int i = 0; i < numberRowBlocks; i++) {
3981               if (rowCounts[i] > 1)
3982                    numberG1++;
3983          }
3984          bool masterColumns = (numberColumnBlocks == numberRowBlocks);
3985          if ((masterColumns && numberElementBlocks == 2 * numberRowBlocks - 1)
3986                    || (!masterColumns && numberElementBlocks == 2 * numberRowBlocks)) {
3987               if (numberG1 < 2)
3988                    decomposeType = 1;
3989          }
3990     }
3991     if (!decomposeType && (numberRowBlocks == numberColumnBlocks ||
3992                            numberRowBlocks == numberColumnBlocks - 1)) {
3993          // could be Benders
3994          int numberG1 = 0;
3995          for (int i = 0; i < numberColumnBlocks; i++) {
3996               if (columnCounts[i] > 1)
3997                    numberG1++;
3998          }
3999          bool masterRows = (numberColumnBlocks == numberRowBlocks);
4000          if ((masterRows && numberElementBlocks == 2 * numberColumnBlocks - 1)
4001                    || (!masterRows && numberElementBlocks == 2 * numberColumnBlocks)) {
4002               if (numberG1 < 2)
4003                    decomposeType = 2;
4004          }
4005     }
4006     delete [] rowCounts;
4007     delete [] columnCounts;
4008     delete [] blockInfo;
4009     // decide what to do
4010     switch (decomposeType) {
4011          // No good
4012     case 0:
4013          loadProblem(*model, false);
4014          return dual();
4015          // DW
4016     case 1:
4017          return solveDW(model);
4018          // Benders
4019     case 2:
4020          return solveBenders(model);
4021     }
4022     return 0; // to stop compiler warning
4023}
4024/* This loads a model from a CoinStructuredModel object - returns number of errors.
4025   If originalOrder then keep to order stored in blocks,
4026   otherwise first column/rows correspond to first block - etc.
4027   If keepSolution true and size is same as current then
4028   keeps current status and solution
4029*/
4030int
4031ClpSimplex::loadProblem (  CoinStructuredModel & coinModel,
4032                           bool originalOrder,
4033                           bool keepSolution)
4034{
4035     unsigned char * status = NULL;
4036     double * psol = NULL;
4037     double * dsol = NULL;
4038     int numberRows = coinModel.numberRows();
4039     int numberColumns = coinModel.numberColumns();
4040     int numberRowBlocks = coinModel.numberRowBlocks();
4041     int numberColumnBlocks = coinModel.numberColumnBlocks();
4042     int numberElementBlocks = coinModel.numberElementBlocks();
4043     if (status_ && numberRows_ && numberRows_ == numberRows &&
4044               numberColumns_ == numberColumns && keepSolution) {
4045          status = new unsigned char [numberRows_+numberColumns_];
4046          CoinMemcpyN(status_, numberRows_ + numberColumns_, status);
4047          psol = new double [numberRows_+numberColumns_];
4048          CoinMemcpyN(columnActivity_, numberColumns_, psol);
4049          CoinMemcpyN(rowActivity_, numberRows_, psol + numberColumns_);
4050          dsol = new double [numberRows_+numberColumns_];
4051          CoinMemcpyN(reducedCost_, numberColumns_, dsol);
4052          CoinMemcpyN(dual_, numberRows_, dsol + numberColumns_);
4053     }
4054     int returnCode = 0;
4055     double * rowLower = new double [numberRows];
4056     double * rowUpper = new double [numberRows];
4057     double * columnLower = new double [numberColumns];
4058     double * columnUpper = new double [numberColumns];
4059     double * objective = new double [numberColumns];
4060     int * integerType = new int [numberColumns];
4061     CoinBigIndex numberElements = 0;
4062     // Bases for blocks
4063     int * rowBase = new int[numberRowBlocks];
4064     CoinFillN(rowBase, numberRowBlocks, -1);
4065     // And row to put it
4066     int * whichRow = new int [numberRows];
4067     int * columnBase = new int[numberColumnBlocks];
4068     CoinFillN(columnBase, numberColumnBlocks, -1);
4069     // And column to put it
4070     int * whichColumn = new int [numberColumns];
4071     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4072          CoinModel * block = coinModel.coinBlock(iBlock);
4073          numberElements += block->numberElements();
4074          //and set up elements etc
4075          double * associated = block->associatedArray();
4076          // If strings then do copies
4077          if (block->stringsExist())
4078               returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
4079                                                 objective, integerType, associated);
4080          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4081          int iRowBlock = info.rowBlock;
4082          int iColumnBlock = info.columnBlock;
4083          if (rowBase[iRowBlock] < 0) {
4084               rowBase[iRowBlock] = block->numberRows();
4085               // Save block number
4086               whichRow[numberRows-numberRowBlocks+iRowBlock] = iBlock;
4087          } else {
4088               assert(rowBase[iRowBlock] == block->numberRows());
4089          }
4090          if (columnBase[iColumnBlock] < 0) {
4091               columnBase[iColumnBlock] = block->numberColumns();
4092               // Save block number
4093               whichColumn[numberColumns-numberColumnBlocks+iColumnBlock] = iBlock;
4094          } else {
4095               assert(columnBase[iColumnBlock] == block->numberColumns());
4096          }
4097     }
4098     // Fill arrays with defaults
4099     CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
4100     CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
4101     CoinFillN(columnLower, numberColumns, 0.0);
4102     CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
4103     CoinFillN(objective, numberColumns, 0.0);
4104     CoinFillN(integerType, numberColumns, 0);
4105     int n = 0;
4106     for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4107          int k = rowBase[iBlock];
4108          rowBase[iBlock] = n;
4109          assert (k >= 0);
4110          // block number
4111          int jBlock = whichRow[numberRows-numberRowBlocks+iBlock];
4112          if (originalOrder) {
4113               memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
4114          } else {
4115               CoinIotaN(whichRow + n, k, n);
4116          }
4117          n += k;
4118     }
4119     assert (n == numberRows);
4120     n = 0;
4121     for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4122          int k = columnBase[iBlock];
4123          columnBase[iBlock] = n;
4124          assert (k >= 0);
4125          if (k) {
4126               // block number
4127               int jBlock = whichColumn[numberColumns-numberColumnBlocks+iBlock];
4128               if (originalOrder) {
4129                    memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
4130                           k * sizeof(int));
4131               } else {
4132                    CoinIotaN(whichColumn + n, k, n);
4133               }
4134               n += k;
4135          }
4136     }
4137     assert (n == numberColumns);
4138     bool gotIntegers = false;
4139     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4140          CoinModel * block = coinModel.coinBlock(iBlock);
4141          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4142          int iRowBlock = info.rowBlock;
4143          int iRowBase = rowBase[iRowBlock];
4144          int iColumnBlock = info.columnBlock;
4145          int iColumnBase = columnBase[iColumnBlock];
4146          if (info.rhs) {
4147               int nRows = block->numberRows();
4148               const double * lower = block->rowLowerArray();
4149               const double * upper = block->rowUpperArray();
4150               for (int i = 0; i < nRows; i++) {
4151                    int put = whichRow[i+iRowBase];
4152                    rowLower[put] = lower[i];
4153                    rowUpper[put] = upper[i];
4154               }
4155          }
4156          if (info.bounds) {
4157               int nColumns = block->numberColumns();
4158               const double * lower = block->columnLowerArray();
4159               const double * upper = block->columnUpperArray();
4160               const double * obj = block->objectiveArray();
4161               for (int i = 0; i < nColumns; i++) {
4162                    int put = whichColumn[i+iColumnBase];
4163                    columnLower[put] = lower[i];
4164                    columnUpper[put] = upper[i];
4165                    objective[put] = obj[i];
4166               }
4167          }
4168          if (info.integer) {
4169               gotIntegers = true;
4170               int nColumns = block->numberColumns();
4171               const int * type = block->integerTypeArray();
4172               for (int i = 0; i < nColumns; i++) {
4173                    int put = whichColumn[i+iColumnBase];
4174                    integerType[put] = type[i];
4175               }
4176          }
4177     }
4178     gutsOfLoadModel(numberRows, numberColumns,
4179                     columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
4180     delete [] rowLower;
4181     delete [] rowUpper;
4182     delete [] columnLower;
4183     delete [] columnUpper;
4184     delete [] objective;
4185     // Do integers if wanted
4186     if (gotIntegers) {
4187          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4188               if (integerType[iColumn])
4189                    setInteger(iColumn);
4190          }
4191     }
4192     delete [] integerType;
4193     setObjectiveOffset(coinModel.objectiveOffset());
4194     // Space for elements
4195     int * row = new int[numberElements];
4196     int * column = new int[numberElements];
4197     double * element = new double[numberElements];
4198     numberElements = 0;
4199     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4200          CoinModel * block = coinModel.coinBlock(iBlock);
4201          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4202          int iRowBlock = info.rowBlock;
4203          int iRowBase = rowBase[iRowBlock];
4204          int iColumnBlock = info.columnBlock;
4205          int iColumnBase = columnBase[iColumnBlock];
4206          if (info.rowName) {
4207               int numberItems = block->rowNames()->numberItems();
4208               assert( block->numberRows() >= numberItems);
4209               if (numberItems) {
4210                    const char *const * rowNames = block->rowNames()->names();
4211                    for (int i = 0; i < numberItems; i++) {
4212                         int put = whichRow[i+iRowBase];
4213                         std::string name = rowNames[i];
4214                         setRowName(put, name);
4215                    }
4216               }
4217          }
4218          if (info.columnName) {
4219               int numberItems = block->columnNames()->numberItems();
4220               assert( block->numberColumns() >= numberItems);
4221               if (numberItems) {
4222                    const char *const * columnNames = block->columnNames()->names();
4223                    for (int i = 0; i < numberItems; i++) {
4224                         int put = whichColumn[i+iColumnBase];
4225                         std::string name = columnNames[i];
4226                         setColumnName(put, name);
4227                    }
4228               }
4229          }
4230          if (info.matrix) {
4231               CoinPackedMatrix matrix2;
4232               const CoinPackedMatrix * matrix = block->packedMatrix();
4233               if (!matrix) {
4234                    double * associated = block->associatedArray();
4235                    block->createPackedMatrix(matrix2, associated);
4236                    matrix = &matrix2;
4237               }
4238               // get matrix data pointers
4239               const int * row2 = matrix->getIndices();
4240               const CoinBigIndex * columnStart = matrix->getVectorStarts();
4241               const double * elementByColumn = matrix->getElements();
4242               const int * columnLength = matrix->getVectorLengths();
4243               int n = matrix->getNumCols();
4244               assert (matrix->isColOrdered());
4245               for (int iColumn = 0; iColumn < n; iColumn++) {
4246                    CoinBigIndex j;
4247                    int jColumn = whichColumn[iColumn+iColumnBase];
4248                    for (j = columnStart[iColumn];
4249                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4250                         row[numberElements] = whichRow[row2[j] + iRowBase];
4251                         column[numberElements] = jColumn;
4252                         element[numberElements++] = elementByColumn[j];
4253                    }
4254               }
4255          }
4256     }
4257     delete [] whichRow;
4258     delete [] whichColumn;
4259     delete [] rowBase;
4260     delete [] columnBase;
4261     CoinPackedMatrix * matrix =
4262          new CoinPackedMatrix (true, row, column, element, numberElements);
4263     matrix_ = new ClpPackedMatrix(matrix);
4264     matrix_->setDimensions(numberRows, numberColumns);
4265     delete [] row;
4266     delete [] column;
4267     delete [] element;
4268     createStatus();
4269     if (status) {
4270          // copy back
4271          CoinMemcpyN(status, numberRows_ + numberColumns_, status_);
4272          CoinMemcpyN(psol, numberColumns_, columnActivity_);
4273          CoinMemcpyN(psol + numberColumns_, numberRows_, rowActivity_);
4274          CoinMemcpyN(dsol, numberColumns_, reducedCost_);
4275          CoinMemcpyN(dsol + numberColumns_, numberRows_, dual_);
4276          delete [] status;
4277          delete [] psol;
4278          delete [] dsol;
4279     }
4280     optimizationDirection_ = coinModel.optimizationDirection();
4281     return returnCode;
4282}
4283/*  If input negative scales objective so maximum <= -value
4284    and returns scale factor used.  If positive unscales and also
4285    redoes dual stuff
4286*/
4287double
4288ClpSimplex::scaleObjective(double value)
4289{
4290     double * obj = objective();
4291     double largest = 0.0;
4292     if (value < 0.0) {
4293          value = - value;
4294          for (int i = 0; i < numberColumns_; i++) {
4295               largest = CoinMax(largest, fabs(obj[i]));
4296          }
4297          if (largest > value) {
4298               double scaleFactor = value / largest;
4299               for (int i = 0; i < numberColumns_; i++) {
4300                    obj[i] *= scaleFactor;
4301                    reducedCost_[i] *= scaleFactor;
4302               }
4303               for (int i = 0; i < numberRows_; i++) {
4304                    dual_[i] *= scaleFactor;
4305               }
4306               largest /= value;
4307          } else {
4308               // no need
4309               largest = 1.0;
4310          }
4311     } else {
4312          // at end
4313          if (value != 1.0) {
4314               for (int i = 0; i < numberColumns_; i++) {
4315                    obj[i] *= value;
4316                    reducedCost_[i] *= value;
4317               }
4318               for (int i = 0; i < numberRows_; i++) {
4319                    dual_[i] *= value;
4320               }
4321               computeObjectiveValue();
4322          }
4323     }
4324     return largest;
4325}
4326// Solve using Dantzig-Wolfe decomposition and maybe in parallel
4327int
4328ClpSimplex::solveDW(CoinStructuredModel * model)
4329{
4330     double time1 = CoinCpuTime();
4331     int numberColumns = model->numberColumns();
4332     int numberRowBlocks = model->numberRowBlocks();
4333     int numberColumnBlocks = model->numberColumnBlocks();
4334     int numberElementBlocks = model->numberElementBlocks();
4335     // We already have top level structure
4336     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4337     for (int i = 0; i < numberElementBlocks; i++) {
4338          CoinModel * thisBlock = model->coinBlock(i);
4339          assert (thisBlock);
4340          // just fill in info
4341          CoinModelBlockInfo info = CoinModelBlockInfo();
4342          int whatsSet = thisBlock->whatIsSet();
4343          info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
4344          info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
4345          info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
4346          info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
4347          info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
4348          info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
4349          // Which block
4350          int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4351          info.rowBlock = iRowBlock;
4352          int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4353          info.columnBlock = iColumnBlock;
4354          blockInfo[i] = info;
4355     }
4356     // make up problems
4357     int numberBlocks = numberRowBlocks - 1;
4358     // Find master rows and columns
4359     int * rowCounts = new int [numberRowBlocks];
4360     CoinZeroN(rowCounts, numberRowBlocks);
4361     int * columnCounts = new int [numberColumnBlocks+1];
4362     CoinZeroN(columnCounts, numberColumnBlocks);
4363     int iBlock;
4364     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4365          int iRowBlock = blockInfo[iBlock].rowBlock;
4366          rowCounts[iRowBlock]++;
4367          int iColumnBlock = blockInfo[iBlock].columnBlock;
4368          columnCounts[iColumnBlock]++;
4369     }
4370     int * whichBlock = new int [numberElementBlocks];
4371     int masterRowBlock = -1;
4372     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4373          if (rowCounts[iBlock] > 1) {
4374               if (masterRowBlock == -1) {
4375                    masterRowBlock = iBlock;
4376               } else {
4377                    // Can't decode
4378                    masterRowBlock = -2;
4379                    break;
4380               }
4381          }
4382     }
4383     int masterColumnBlock = -1;
4384     int kBlock = 0;
4385     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4386          int count = columnCounts[iBlock];
4387          columnCounts[iBlock] = kBlock;
4388          kBlock += count;
4389     }
4390     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4391          int iColumnBlock = blockInfo[iBlock].columnBlock;
4392          whichBlock[columnCounts[iColumnBlock]] = iBlock;
4393          columnCounts[iColumnBlock]++;
4394     }
4395     for (iBlock = numberColumnBlocks - 1; iBlock >= 0; iBlock--)
4396          columnCounts[iBlock+1] = columnCounts[iBlock];
4397     columnCounts[0] = 0;
4398     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4399          int count = columnCounts[iBlock+1] - columnCounts[iBlock];
4400          if (count == 1) {
4401               int kBlock = whichBlock[columnCounts[iBlock]];
4402               int iRowBlock = blockInfo[kBlock].rowBlock;
4403               if (iRowBlock == masterRowBlock) {
4404                    if (masterColumnBlock == -1) {
4405                         masterColumnBlock = iBlock;
4406                    } else {
4407                         // Can't decode
4408                         masterColumnBlock = -2;
4409                         break;
4410                    }
4411               }
4412          }
4413     }
4414     if (masterRowBlock < 0 || masterColumnBlock == -2) {
4415          // What now
4416          abort();
4417     }
4418     delete [] rowCounts;
4419     // create all data
4420     const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
4421     ClpSimplex * sub = new ClpSimplex [numberBlocks];
4422     ClpSimplex master;
4423     // Set offset
4424     master.setObjectiveOffset(model->objectiveOffset());
4425     kBlock = 0;
4426     int masterBlock = -1;
4427     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4428          top[kBlock] = NULL;
4429          int start = columnCounts[iBlock];
4430          int end = columnCounts[iBlock+1];
4431          assert (end - start <= 2);
4432          for (int j = start; j < end; j++) {
4433               int jBlock = whichBlock[j];
4434               int iRowBlock = blockInfo[jBlock].rowBlock;
4435               int iColumnBlock = blockInfo[jBlock].columnBlock;
4436               assert (iColumnBlock == iBlock);
4437               if (iColumnBlock != masterColumnBlock && iRowBlock == masterRowBlock) {
4438                    // top matrix
4439                    top[kBlock] = model->coinBlock(jBlock)->packedMatrix();
4440               } else {
4441                    const CoinPackedMatrix * matrix
4442                    = model->coinBlock(jBlock)->packedMatrix();
4443                    // Get pointers to arrays
4444                    const double * rowLower;
4445                    const double * rowUpper;
4446                    const double * columnLower;
4447                    const double * columnUpper;
4448                    const double * objective;
4449                    model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
4450                                 columnLower, columnUpper, objective);
4451                    if (iColumnBlock != masterColumnBlock) {
4452                         // diagonal block
4453                         sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
4454                                                 objective, rowLower, rowUpper);
4455                         if (true) {
4456                              double * lower = sub[kBlock].columnLower();
4457                              double * upper = sub[kBlock].columnUpper();
4458                              int n = sub[kBlock].numberColumns();
4459                              for (int i = 0; i < n; i++) {
4460                                   lower[i] = CoinMax(-1.0e8, lower[i]);
4461                                   upper[i] = CoinMin(1.0e8, upper[i]);
4462                              }
4463                         }
4464                         if (optimizationDirection_ < 0.0) {
4465                              double * obj = sub[kBlock].objective();
4466                              int n = sub[kBlock].numberColumns();
4467                              for (int i = 0; i < n; i++)
4468                                   obj[i] = - obj[i];
4469                         }
4470                         if (this->factorizationFrequency() == 200) {
4471                              // User did not touch preset
4472                              sub[kBlock].defaultFactorizationFrequency();
4473                         } else {
4474                              // make sure model has correct value
4475                              sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
4476                         }
4477                         sub[kBlock].setPerturbation(50);
4478                         // Set columnCounts to be diagonal block index for cleanup
4479                         columnCounts[kBlock] = jBlock;
4480                    } else {
4481                         // master
4482                         masterBlock = jBlock;
4483                         master.loadProblem(*matrix, columnLower, columnUpper,
4484                                            objective, rowLower, rowUpper);
4485                         if (optimizationDirection_ < 0.0) {
4486                              double * obj = master.objective();
4487                              int n = master.numberColumns();
4488                              for (int i = 0; i < n; i++)
4489                                   obj[i] = - obj[i];
4490                         }
4491                    }
4492               }
4493          }
4494          if (iBlock != masterColumnBlock)
4495               kBlock++;
4496     }
4497     delete [] whichBlock;
4498     delete [] blockInfo;
4499     // For now master must have been defined (does not have to have columns)
4500     assert (master.numberRows());
4501     assert (masterBlock >= 0);
4502     int numberMasterRows = master.numberRows();
4503     // Overkill in terms of space
4504     int spaceNeeded = CoinMax(numberBlocks * (numberMasterRows + 1),
4505                               2 * numberMasterRows);
4506     int * rowAdd = new int[spaceNeeded];
4507     double * elementAdd = new double[spaceNeeded];
4508     spaceNeeded = numberBlocks;
4509     int * columnAdd = new int[spaceNeeded+1];
4510     double * objective = new double[spaceNeeded];
4511     // Add in costed slacks
4512     int firstArtificial = master.numberColumns();
4513     int lastArtificial = firstArtificial;
4514     if (true) {
4515          const double * lower = master.rowLower();
4516          const double * upper = master.rowUpper();
4517          int kCol = 0;
4518          for (int iRow = 0; iRow < numberMasterRows; iRow++) {
4519               if (lower[iRow] > -1.0e10) {
4520                    rowAdd[kCol] = iRow;
4521                    elementAdd[kCol++] = 1.0;
4522               }
4523               if (upper[iRow] < 1.0e10) {
4524                    rowAdd[kCol] = iRow;
4525                    elementAdd[kCol++] = -1.0;
4526               }
4527          }
4528          if (kCol > spaceNeeded) {
4529               spaceNeeded = kCol;
4530               delete [] columnAdd;
4531               delete [] objective;
4532               columnAdd = new int[spaceNeeded+1];
4533               objective = new double[spaceNeeded];
4534          }
4535          for (int i = 0; i < kCol; i++) {
4536               columnAdd[i] = i;
4537               objective[i] = 1.0e13;
4538          }
4539          columnAdd[kCol] = kCol;
4540          master.addColumns(kCol, NULL, NULL, objective,
4541                            columnAdd, rowAdd, elementAdd);
4542          lastArtificial = master.numberColumns();
4543     }
4544     int maxPass = 500;
4545     int iPass;
4546     double lastObjective = 1.0e31;
4547     // Create convexity rows for proposals
4548     int numberMasterColumns = master.numberColumns();
4549     master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
4550     if (this->factorizationFrequency() == 200) {
4551          // User did not touch preset
4552          master.defaultFactorizationFrequency();
4553     } else {
4554          // make sure model has correct value
4555          master.setFactorizationFrequency(this->factorizationFrequency());
4556     }
4557     master.setPerturbation(50);
4558     // Arrays to say which block and when created
4559     int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
4560     whichBlock = new int[maximumColumns];
4561     int * when = new int[maximumColumns];
4562     int numberColumnsGenerated = numberBlocks;
4563     // fill in rhs and add in artificials
4564     {
4565          double * rowLower = master.rowLower();
4566          double * rowUpper = master.rowUpper();
4567          int iBlock;
4568          columnAdd[0] = 0;
4569          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4570               int iRow = iBlock + numberMasterRows;;
4571               rowLower[iRow] = 1.0;
4572               rowUpper[iRow] = 1.0;
4573               rowAdd[iBlock] = iRow;
4574               elementAdd[iBlock] = 1.0;
4575               objective[iBlock] = 1.0e13;
4576               columnAdd[iBlock+1] = iBlock + 1;
4577               when[iBlock] = -1;
4578               whichBlock[iBlock] = iBlock;
4579          }
4580          master.addColumns(numberBlocks, NULL, NULL, objective,
4581                            columnAdd, rowAdd, elementAdd);
4582     }
4583     // and resize matrix to double check clp will be happy
4584     //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
4585     //                  numberMasterColumns+numberBlocks);
4586     std::cout << "Time to decompose " << CoinCpuTime() - time1 << " seconds" << std::endl;
4587     for (iPass = 0; iPass < maxPass; iPass++) {
4588          printf("Start of pass %d\n", iPass);
4589          // Solve master - may be infeasible
4590          //master.scaling(0);
4591          if (0) {
4592               master.writeMps("yy.mps");
4593          }
4594          // Correct artificials
4595          double sumArtificials = 0.0;
4596          if (iPass) {
4597               double * upper = master.columnUpper();
4598               double * solution = master.primalColumnSolution();
4599               double * obj = master.objective();
4600               sumArtificials = 0.0;
4601               for (int i = firstArtificial; i < lastArtificial; i++) {
4602                    sumArtificials += solution[i];
4603                    //assert (solution[i]>-1.0e-2);
4604                    if (solution[i] < 1.0e-6) {
4605#if 0
4606                         // Could take out
4607                         obj[i] = 0.0;
4608                         upper[i] = 0.0;
4609#else
4610                         obj[i] = 1.0e7;
4611                         upper[i] = 1.0e-1;
4612#endif
4613                         solution[i] = 0.0;
4614                         master.setColumnStatus(i, isFixed);
4615                    } else {
4616                         upper[i] = solution[i] + 1.0e-5 * (1.0 + solution[i]);
4617                    }
4618               }
4619               printf("Sum of artificials before solve is %g\n", sumArtificials);
4620          }
4621          // scale objective to be reasonable
4622          double scaleFactor = master.scaleObjective(-1.0e9);
4623          {
4624               double * dual = master.dualRowSolution();
4625               int n = master.numberRows();
4626               memset(dual, 0, n * sizeof(double));
4627               double * solution = master.primalColumnSolution();
4628               master.clpMatrix()->times(1.0, solution, dual);
4629               double sum = 0.0;
4630               double * lower = master.rowLower();
4631               double * upper = master.rowUpper();
4632               for (int iRow = 0; iRow < n; iRow++) {
4633                    double value = dual[iRow];
4634                    if (value > upper[iRow])
4635                         sum += value - upper[iRow];
4636                    else if (value < lower[iRow])
4637                         sum -= value - lower[iRow];
4638               }
4639               printf("suminf %g\n", sum);
4640               lower = master.columnLower();
4641               upper = master.columnUpper();
4642               n = master.numberColumns();
4643               for (int iColumn = 0; iColumn < n; iColumn++) {
4644                    double value = solution[iColumn];
4645                    if (value > upper[iColumn] + 1.0e-5)
4646                         sum += value - upper[iColumn];
4647                    else if (value < lower[iColumn] - 1.0e-5)
4648                         sum -= value - lower[iColumn];
4649               }
4650               printf("suminf %g\n", sum);
4651          }
4652          master.primal(1);
4653          // Correct artificials
4654          sumArtificials = 0.0;
4655          {
4656               double * solution = master.primalColumnSolution();
4657               for (int i = firstArtificial; i < lastArtificial; i++) {
4658                    sumArtificials += solution[i];
4659               }
4660               printf("Sum of artificials after solve is %g\n", sumArtificials);
4661          }
4662          master.scaleObjective(scaleFactor);
4663          int problemStatus = master.status(); // do here as can change (delcols)
4664          if (master.numberIterations() == 0 && iPass)
4665               break; // finished
4666          if (master.objectiveValue() > lastObjective - 1.0e-7 && iPass > 555)
4667               break; // finished
4668          lastObjective = master.objectiveValue();
4669          // mark basic ones and delete if necessary
4670          int iColumn;
4671          numberColumnsGenerated = master.numberColumns() - numberMasterColumns;
4672          for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
4673               if (master.getStatus(iColumn + numberMasterColumns) == ClpSimplex::basic)
4674                    when[iColumn] = iPass;
4675          }
4676          if (numberColumnsGenerated + numberBlocks > maximumColumns) {
4677               // delete
4678               int numberKeep = 0;
4679               int numberDelete = 0;
4680               int * whichDelete = new int[numberColumnsGenerated];
4681               for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
4682                    if (when[iColumn] > iPass - 7) {
4683                         // keep
4684                         when[numberKeep] = when[iColumn];
4685                         whichBlock[numberKeep++] = whichBlock[iColumn];
4686                    } else {
4687                         // delete
4688                         whichDelete[numberDelete++] = iColumn + numberMasterColumns;
4689                    }
4690               }
4691               numberColumnsGenerated -= numberDelete;
4692               master.deleteColumns(numberDelete, whichDelete);
4693               delete [] whichDelete;
4694          }
4695          const double * dual = NULL;
4696          bool deleteDual = false;
4697          if (problemStatus == 0) {
4698               dual = master.dualRowSolution();
4699          } else if (problemStatus == 1) {
4700               // could do composite objective
4701               dual = master.infeasibilityRay();
4702               deleteDual = true;
4703               printf("The sum of infeasibilities is %g\n",
4704                      master.sumPrimalInfeasibilities());
4705          } else if (!master.numberColumns()) {
4706               assert(!iPass);
4707               dual = master.dualRowSolution();
4708               memset(master.dualRowSolution(),
4709                      0, (numberMasterRows + numberBlocks)*sizeof(double));
4710          } else {
4711               abort();
4712          }
4713          // Scale back on first time
4714          if (!iPass) {
4715               double * dual2 = master.dualRowSolution();
4716               for (int i = 0; i < numberMasterRows + numberBlocks; i++) {
4717                    dual2[i] *= 1.0e-7;
4718               }
4719               dual = master.dualRowSolution();
4720          }
4721          // Create objective for sub problems and solve
4722          columnAdd[0] = 0;
4723          int numberProposals = 0;
4724          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4725               int numberColumns2 = sub[iBlock].numberColumns();
4726               double * saveObj = new double [numberColumns2];
4727               double * objective2 = sub[iBlock].objective();
4728               memcpy(saveObj, objective2, numberColumns2 * sizeof(double));
4729               // new objective
4730               top[iBlock]->transposeTimes(dual, objective2);
4731               int i;
4732               if (problemStatus == 0) {
4733                    for (i = 0; i < numberColumns2; i++)
4734                         objective2[i] = saveObj[i] - objective2[i];
4735               } else {
4736                    for (i = 0; i < numberColumns2; i++)
4737                         objective2[i] = -objective2[i];
4738               }
4739               // scale objective to be reasonable
4740               double scaleFactor =
4741                    sub[iBlock].scaleObjective((sumArtificials > 1.0e-5) ? -1.0e-4 : -1.0e9);
4742               if (iPass) {
4743                    sub[iBlock].primal();
4744               } else {
4745                    sub[iBlock].dual();
4746               }
4747               sub[iBlock].scaleObjective(scaleFactor);
4748               if (!sub[iBlock].isProvenOptimal() &&
4749                         !sub[iBlock].isProvenDualInfeasible()) {
4750                    memset(objective2, 0, numberColumns2 * sizeof(double));
4751                    sub[iBlock].primal();
4752                    if (problemStatus == 0) {
4753                         for (i = 0; i < numberColumns2; i++)
4754                              objective2[i] = saveObj[i] - objective2[i];
4755                    } else {
4756                         for (i = 0; i < numberColumns2; i++)
4757                              objective2[i] = -objective2[i];
4758                    }
4759                    double scaleFactor = sub[iBlock].scaleObjective(-1.0e9);
4760                    sub[iBlock].primal(1);
4761                    sub[iBlock].scaleObjective(scaleFactor);
4762               }
4763               memcpy(objective2, saveObj, numberColumns2 * sizeof(double));
4764               // get proposal
4765               if (sub[iBlock].numberIterations() || !iPass) {
4766                    double objValue = 0.0;
4767                    int start = columnAdd[numberProposals];
4768                    // proposal
4769                    if (sub[iBlock].isProvenOptimal()) {
4770                         const double * solution = sub[iBlock].primalColumnSolution();
4771                         top[iBlock]->times(solution, elementAdd + start);
4772                         for (i = 0; i < numberColumns2; i++)
4773                              objValue += solution[i] * saveObj[i];
4774                         // See if good dj and pack down
4775                         int number = start;
4776                         double dj = objValue;
4777                         if (problemStatus)
4778                              dj = 0.0;
4779                         double smallest = 1.0e100;
4780                         double largest = 0.0;
4781                         for (i = 0; i < numberMasterRows; i++) {
4782                              double value = elementAdd[start+i];
4783                              if (fabs(value) > 1.0e-15) {
4784                                   dj -= dual[i] * value;
4785                                   smallest = CoinMin(smallest, fabs(value));
4786                                   largest = CoinMax(largest, fabs(value));
4787                                   rowAdd[number] = i;
4788                                   elementAdd[number++] = value;
4789                              }
4790                         }
4791                         // and convexity
4792                         dj -= dual[numberMasterRows+iBlock];
4793                         rowAdd[number] = numberMasterRows + iBlock;
4794                         elementAdd[number++] = 1.0;
4795                         // if elements large then scale?
4796                         //if (largest>1.0e8||smallest<1.0e-8)
4797                         printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
4798                                iBlock, smallest, largest, dj);
4799                         if (dj < -1.0e-6 || !iPass) {
4800                              // take
4801                              objective[numberProposals] = objValue;
4802                              columnAdd[++numberProposals] = number;
4803                              when[numberColumnsGenerated] = iPass;
4804                              whichBlock[numberColumnsGenerated++] = iBlock;
4805                         }
4806                    } else if (sub[iBlock].isProvenDualInfeasible()) {
4807                         // use ray
4808                         const double * solution = sub[iBlock].unboundedRay();
4809                         top[iBlock]->times(solution, elementAdd + start);
4810                         for (i = 0; i < numberColumns2; i++)
4811                              objValue += solution[i] * saveObj[i];
4812                         // See if good dj and pack down
4813                         int number = start;
4814                         double dj = objValue;
4815                         double smallest = 1.0e100;
4816                         double largest = 0.0;
4817                         for (i = 0; i < numberMasterRows; i++) {
4818                              double value = elementAdd[start+i];
4819                              if (fabs(value) > 1.0e-15) {
4820                                   dj -= dual[i] * value;
4821                                   smallest = CoinMin(smallest, fabs(value));
4822                                   largest = CoinMax(largest, fabs(value));
4823                                   rowAdd[number] = i;
4824                                   elementAdd[number++] = value;
4825                              }
4826                         }
4827                         // if elements large or small then scale?
4828                         //if (largest>1.0e8||smallest<1.0e-8)
4829                         printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
4830                                iBlock, smallest, largest, dj);
4831                         if (dj < -1.0e-6) {
4832                              // take
4833                              objective[numberProposals] = objValue;
4834                              columnAdd[++numberProposals] = number;
4835                              when[numberColumnsGenerated] = iPass;
4836                              whichBlock[numberColumnsGenerated++] = iBlock;
4837                         }
4838                    } else {
4839                         abort();
4840                    }
4841               }
4842               delete [] saveObj;
4843          }
4844          if (deleteDual)
4845               delete [] dual;
4846          if (numberProposals)
4847               master.addColumns(numberProposals, NULL, NULL, objective,
4848                                 columnAdd, rowAdd, elementAdd);
4849     }
4850     std::cout << "Time at end of D-W " << CoinCpuTime() - time1 << " seconds" << std::endl;
4851     //master.scaling(0);
4852     //master.primal(1);
4853     loadProblem(*model);
4854     // now put back a good solution
4855     double * lower = new double[numberMasterRows+numberBlocks];
4856     double * upper = new double[numberMasterRows+numberBlocks];
4857     numberColumnsGenerated  += numberMasterColumns;
4858     double * sol = new double[numberColumnsGenerated];
4859     const double * solution = master.primalColumnSolution();
4860     const double * masterLower = master.rowLower();
4861     const double * masterUpper = master.rowUpper();
4862     double * fullSolution = primalColumnSolution();
4863     const double * fullLower = columnLower();
4864     const double * fullUpper = columnUpper();
4865     const double * rowSolution = master.primalRowSolution();
4866     double * fullRowSolution = primalRowSolution();
4867     const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4868     int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4869     const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4870     int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4871     for (int iRow = 0; iRow < numberRows2; iRow++) {
4872          int kRow = rowBack[iRow];
4873          setRowStatus(kRow, master.getRowStatus(iRow));
4874          fullRowSolution[kRow] = rowSolution[iRow];
4875     }
4876     for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4877          int kColumn = columnBack[iColumn];
4878          setStatus(kColumn, master.getStatus(iColumn));
4879          fullSolution[kColumn] = solution[iColumn];
4880     }
4881     for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4882          // move basis
4883          int kBlock = columnCounts[iBlock];
4884          const int * rowBack = model->coinBlock(kBlock)->originalRows();
4885          int numberRows2 = model->coinBlock(kBlock)->numberRows();
4886          const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4887          int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4888          for (int iRow = 0; iRow < numberRows2; iRow++) {
4889               int kRow = rowBack[iRow];
4890               setRowStatus(kRow, sub[iBlock].getRowStatus(iRow));
4891          }
4892          for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4893               int kColumn = columnBack[iColumn];
4894               setStatus(kColumn, sub[iBlock].getStatus(iColumn));
4895          }
4896          // convert top bit to by rows
4897          CoinPackedMatrix topMatrix = *top[iBlock];
4898          topMatrix.reverseOrdering();
4899          // zero solution
4900          memset(sol, 0, numberColumnsGenerated * sizeof(double));
4901
4902          for (int i = numberMasterColumns; i < numberColumnsGenerated; i++) {
4903               if (whichBlock[i-numberMasterColumns] == iBlock)
4904                    sol[i] = solution[i];
4905          }
4906          memset(lower, 0, (numberMasterRows + numberBlocks)*sizeof(double));
4907          master.clpMatrix()->times(1.0, sol, lower);
4908          for (int iRow = 0; iRow < numberMasterRows; iRow++) {
4909               double value = lower[iRow];
4910               if (masterUpper[iRow] < 1.0e20)
4911                    upper[iRow] = value;
4912               else
4913                    upper[iRow] = COIN_DBL_MAX;
4914               if (masterLower[iRow] > -1.0e20)
4915                    lower[iRow] = value;
4916               else
4917                    lower[iRow] = -COIN_DBL_MAX;
4918          }
4919          sub[iBlock].addRows(numberMasterRows, lower, upper,
4920                              topMatrix.getVectorStarts(),
4921                              topMatrix.getVectorLengths(),
4922                              topMatrix.getIndices(),
4923                              topMatrix.getElements());
4924          sub[iBlock].primal(1);
4925          const double * subSolution = sub[iBlock].primalColumnSolution();
4926          const double * subRowSolution = sub[iBlock].primalRowSolution();
4927          // move solution
4928          for (int iRow = 0; iRow < numberRows2; iRow++) {
4929               int kRow = rowBack[iRow];
4930               fullRowSolution[kRow] = subRowSolution[iRow];
4931          }
4932          for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4933               int kColumn = columnBack[iColumn];
4934               fullSolution[kColumn] = subSolution[iColumn];
4935          }
4936     }
4937     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4938          if (fullSolution[iColumn] < fullUpper[iColumn] - 1.0e-8 &&
4939                    fullSolution[iColumn] > fullLower[iColumn] + 1.0e-8) {
4940               if (getStatus(iColumn) != ClpSimplex::basic) {
4941                    if (columnLower_[iColumn] > -1.0e30 ||
4942                              columnUpper_[iColumn] < 1.0e30)
4943                         setStatus(iColumn, ClpSimplex::superBasic);
4944                    else
4945                         setStatus(iColumn, ClpSimplex::isFree);
4946               }
4947          } else if (fullSolution[iColumn] >= fullUpper[iColumn] - 1.0e-8) {
4948               // may help to make rest non basic
4949               if (getStatus(iColumn) != ClpSimplex::basic)
4950                    setStatus(iColumn, ClpSimplex::atUpperBound);
4951          } else if (fullSolution[iColumn] <= fullLower[iColumn] + 1.0e-8) {
4952               // may help to make rest non basic
4953               if (getStatus(iColumn) != ClpSimplex::basic)
4954                    setStatus(iColumn, ClpSimplex::atLowerBound);
4955          }
4956     }
4957     //int numberRows=model->numberRows();
4958     //for (int iRow=0;iRow<numberRows;iRow++)
4959     //setRowStatus(iRow,ClpSimplex::superBasic);
4960     std::cout << "Time before cleanup of full model " << CoinCpuTime() - time1 << " seconds" << std::endl;
4961     primal(1);
4962     std::cout << "Total time " << CoinCpuTime() - time1 << " seconds" << std::endl;
4963     delete [] columnCounts;
4964     delete [] sol;
4965     delete [] lower;
4966     delete [] upper;
4967     delete [] whichBlock;
4968     delete [] when;
4969     delete [] columnAdd;
4970     delete [] rowAdd;
4971     delete [] elementAdd;
4972     delete [] objective;
4973     delete [] top;
4974     delete [] sub;
4975     return 0;
4976}
4977// Solve using Benders decomposition and maybe in parallel
4978int
4979ClpSimplex::solveBenders(CoinStructuredModel * model)
4980{
4981     double time1 = CoinCpuTime();
4982     //int numberColumns = model->numberColumns();
4983     int numberRowBlocks = model->numberRowBlocks();
4984     int numberColumnBlocks = model->numberColumnBlocks();
4985     int numberElementBlocks = model->numberElementBlocks();
4986     // We already have top level structure
4987     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4988     for (int i = 0; i < numberElementBlocks; i++) {
4989          CoinModel * thisBlock = model->coinBlock(i);
4990          assert (thisBlock);
4991          // just fill in info
4992          CoinModelBlockInfo info = CoinModelBlockInfo();
4993          int whatsSet = thisBlock->whatIsSet();
4994          info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
4995          info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
4996          info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
4997          info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
4998          info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
4999          info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
5000          // Which block
5001          int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
5002          info.rowBlock = iRowBlock;
5003          int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
5004          info.columnBlock = iColumnBlock;
5005          blockInfo[i] = info;
5006     }
5007     // make up problems
5008     int numberBlocks = numberColumnBlocks - 1;
5009     // Find master columns and rows
5010     int * columnCounts = new int [numberColumnBlocks];
5011     CoinZeroN(columnCounts, numberColumnBlocks);
5012     int * rowCounts = new int [numberRowBlocks+1];
5013     CoinZeroN(rowCounts, numberRowBlocks);
5014     int iBlock;
5015     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5016          int iColumnBlock = blockInfo[iBlock].columnBlock;
5017          columnCounts[iColumnBlock]++;
5018          int iRowBlock = blockInfo[iBlock].rowBlock;
5019          rowCounts[iRowBlock]++;
5020     }
5021     int * whichBlock = new int [numberElementBlocks];
5022     int masterColumnBlock = -1;
5023     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5024          if (columnCounts[iBlock] > 1) {
5025               if (masterColumnBlock == -1) {
5026                    masterColumnBlock = iBlock;
5027               } else {
5028                    // Can't decode
5029                    masterColumnBlock = -2;
5030                    break;
5031               }
5032          }
5033     }
5034     int masterRowBlock = -1;
5035     int kBlock = 0;
5036     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5037          int count = rowCounts[iBlock];
5038          rowCounts[iBlock] = kBlock;
5039          kBlock += count;
5040     }
5041     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5042          int iRowBlock = blockInfo[iBlock].rowBlock;
5043          whichBlock[rowCounts[iRowBlock]] = iBlock;
5044          rowCounts[iRowBlock]++;
5045     }
5046     for (iBlock = numberRowBlocks - 1; iBlock >= 0; iBlock--)
5047          rowCounts[iBlock+1] = rowCounts[iBlock];
5048     rowCounts[0] = 0;
5049     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5050          int count = rowCounts[iBlock+1] - rowCounts[iBlock];
5051          if (count == 1) {
5052               int kBlock = whichBlock[rowCounts[iBlock]];
5053               int iColumnBlock = blockInfo[kBlock].columnBlock;
5054               if (iColumnBlock == masterColumnBlock) {
5055                    if (masterRowBlock == -1) {
5056                         masterRowBlock = iBlock;
5057                    } else {
5058                         // Can't decode
5059                         masterRowBlock = -2;
5060                         break;
5061                    }
5062               }
5063          }
5064     }
5065     if (masterColumnBlock < 0 || masterRowBlock == -2) {
5066          // What now
5067          abort();
5068     }
5069     delete [] columnCounts;
5070     // create all data
5071     const CoinPackedMatrix ** first = new const CoinPackedMatrix * [numberRowBlocks];
5072     ClpSimplex * sub = new ClpSimplex [numberBlocks];
5073     ClpSimplex master;
5074     // Set offset
5075     master.setObjectiveOffset(model->objectiveOffset());
5076     kBlock = 0;
5077     int masterBlock = -1;
5078     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5079          first[kBlock] = NULL;
5080          int start = rowCounts[iBlock];
5081          int end = rowCounts[iBlock+1];
5082          assert (end - start <= 2);
5083          for (int j = start; j < end; j++) {
5084               int jBlock = whichBlock[j];
5085               int iColumnBlock = blockInfo[jBlock].columnBlock;
5086               int iRowBlock = blockInfo[jBlock].rowBlock;
5087               assert (iRowBlock == iBlock);
5088               if (iRowBlock != masterRowBlock && iColumnBlock == masterColumnBlock) {
5089                    // first matrix
5090                    first[kBlock] = model->coinBlock(jBlock)->packedMatrix();
5091               } else {
5092                    const CoinPackedMatrix * matrix
5093                    = model->coinBlock(jBlock)->packedMatrix();
5094                    // Get pointers to arrays
5095                    const double * columnLower;
5096                    const double * columnUpper;
5097                    const double * rowLower;
5098                    const double * rowUpper;
5099                    const double * objective;
5100                    model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
5101                                 columnLower, columnUpper, objective);
5102                    if (iRowBlock != masterRowBlock) {
5103                         // diagonal block
5104                         sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
5105                                                 objective, rowLower, rowUpper);
5106                         if (optimizationDirection_ < 0.0) {
5107                              double * obj = sub[kBlock].objective();
5108                              int n = sub[kBlock].numberColumns();
5109                              for (int i = 0; i < n; i++)
5110                                   obj[i] = - obj[i];
5111                         }
5112                         if (this->factorizationFrequency() == 200) {
5113                              // User did not touch preset
5114                              sub[kBlock].defaultFactorizationFrequency();
5115                         } else {
5116                              // make sure model has correct value
5117                              sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
5118                         }
5119                         sub[kBlock].setPerturbation(50);
5120                         // Set rowCounts to be diagonal block index for cleanup
5121                         rowCounts[kBlock] = jBlock;
5122                    } else {
5123                         // master
5124                         masterBlock = jBlock;
5125                         master.loadProblem(*matrix, columnLower, columnUpper,
5126                                            objective, rowLower, rowUpper);
5127