source: stable/1.15/Clp/src/ClpSolve.cpp @ 2018

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

changes for parallel and idiot

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 220.5 KB
Line 
1/* $Id: ClpSolve.cpp 2006 2013-12-12 15:40:41Z stefan $ */
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()||!numberRows_||!numberColumns_) {
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#ifdef COIN_HAS_VOL
1697                    int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
1698                    nPasses=0;
1699                    if (!returnCode) {
1700                      time2 = CoinCpuTime();
1701                      timeIdiot = time2 - timeX;
1702                      handler_->message(CLP_INTERVAL_TIMING, messages_)
1703                        << "Idiot Crash" << timeIdiot << time2 - time1
1704                        << CoinMessageEol;
1705                      timeX = time2;
1706                    }
1707#endif
1708                    if (nPasses > 70) {
1709                         info.setStartingWeight(1.0e3);
1710                         info.setReduceIterations(6);
1711                         //if (nPasses > 200)
1712                         //info.setFeasibilityTolerance(1.0e-9);
1713                         //if (nPasses > 1900)
1714                         //info.setWeightFactor(0.93);
1715                         if (nPasses > 900) {
1716                           double reductions=nPasses/6.0;
1717                           if (nPasses<5000) {
1718                             reductions /= 12.0;
1719                           } else {
1720                             reductions /= 13.0;
1721                             info.setStartingWeight(1.0e4);
1722                           }
1723                           double ratio=1.0/std::pow(10.0,(1.0/reductions));
1724                           printf("%d passes reduction factor %g\n",nPasses,ratio);
1725                           info.setWeightFactor(ratio);
1726                         } else if (nPasses > 500) {
1727                           info.setWeightFactor(0.7);
1728                         } else if (nPasses > 200) {
1729                           info.setWeightFactor(0.5);
1730                         }
1731                         if (maximumIterations()<nPasses) {
1732                           printf("Presuming maximumIterations is just for Idiot\n");
1733                           nPasses=maximumIterations();
1734                           setMaximumIterations(COIN_INT_MAX);
1735                           model2->setMaximumIterations(COIN_INT_MAX);
1736                         }
1737                         if (nPasses >= 10000&&nPasses<100000) {
1738                              int k = nPasses % 100;
1739                              nPasses /= 200;
1740                              info.setReduceIterations(3);
1741                              if (k)
1742                                   info.setStartingWeight(1.0e2);
1743                         }
1744                         // also be more lenient on infeasibilities
1745                         info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
1746                         info.setDropEnoughWeighted(-2.0);
1747                    } else if (nPasses >= 50) {
1748                         info.setStartingWeight(1.0e3);
1749                         //info.setReduceIterations(6);
1750                    }
1751                    // For experimenting
1752                    if (nPasses < 70 && (nPasses % 10) > 0 && (nPasses % 10) < 4) {
1753                         info.setStartingWeight(1.0e3);
1754                         info.setLightweight(nPasses % 10); // special testing
1755#ifdef COIN_DEVELOP
1756                         printf("warning - odd lightweight %d\n", nPasses % 10);
1757                         //info.setReduceIterations(6);
1758#endif
1759                    }
1760               }
1761               if (options.getExtraInfo(1) > 1000000)
1762                    nPasses += 1000000;
1763               if (nPasses) {
1764                    doCrash = 0;
1765#if 0
1766                    double * solution = model2->primalColumnSolution();
1767                    int iColumn;
1768                    double * saveLower = new double[numberColumns];
1769                    CoinMemcpyN(model2->columnLower(), numberColumns, saveLower);
1770                    double * saveUpper = new double[numberColumns];
1771                    CoinMemcpyN(model2->columnUpper(), numberColumns, saveUpper);
1772                    printf("doing tighten before idiot\n");
1773                    model2->tightenPrimalBounds();
1774                    // Move solution
1775                    double * columnLower = model2->columnLower();
1776                    double * columnUpper = model2->columnUpper();
1777                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1778                         if (columnLower[iColumn] > 0.0)
1779                              solution[iColumn] = columnLower[iColumn];
1780                         else if (columnUpper[iColumn] < 0.0)
1781                              solution[iColumn] = columnUpper[iColumn];
1782                         else
1783                              solution[iColumn] = 0.0;
1784                    }
1785                    CoinMemcpyN(saveLower, numberColumns, columnLower);
1786                    CoinMemcpyN(saveUpper, numberColumns, columnUpper);
1787                    delete [] saveLower;
1788                    delete [] saveUpper;
1789#else
1790                    // Allow for crossover
1791                    //#define LACI_TRY
1792#ifndef LACI_TRY
1793                    //if (doIdiot>0)
1794#ifdef ABC_INHERIT
1795                    if (!model2->abcState())
1796#endif
1797                    info.setStrategy(512 | info.getStrategy());
1798#endif
1799                    // Allow for scaling
1800                    info.setStrategy(32 | info.getStrategy());
1801                    info.crash(nPasses, model2->messageHandler(), model2->messagesPointer());
1802#endif
1803                    time2 = CoinCpuTime();
1804                    timeIdiot = time2 - timeX;
1805                    handler_->message(CLP_INTERVAL_TIMING, messages_)
1806                              << "Idiot Crash" << timeIdiot << time2 - time1
1807                              << CoinMessageEol;
1808                    timeX = time2;
1809                    if (nPasses>100000&&nPasses<100500) {
1810                      // make sure no status left
1811                      model2->createStatus();
1812                      // solve
1813                      if (model2->factorizationFrequency() == 200) {
1814                        // User did not touch preset
1815                        model2->defaultFactorizationFrequency();
1816                      }
1817                      //int numberRows = model2->numberRows();
1818                      int numberColumns = model2->numberColumns();
1819                      // save duals
1820                      //double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows);
1821                      // for moment this only works on nug etc (i.e. all ==)
1822                      // needs row objective
1823                      double * saveObj = CoinCopyOfArray(model2->objective(),numberColumns);
1824                      double * pi = model2->dualRowSolution();
1825                      model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective());
1826                      // just primal values pass
1827                      double saveScale = model2->objectiveScale();
1828                      model2->setObjectiveScale(1.0e-3);
1829                      model2->primal(2);
1830                      model2->writeMps("xx.mps");
1831                      double * solution = model2->primalColumnSolution();
1832                      double * upper = model2->columnUpper();
1833                      for (int i=0;i<numberColumns;i++) {
1834                        if (solution[i]<100.0)
1835                          upper[i]=1000.0;
1836                      }
1837                      model2->setProblemStatus(-1);
1838                      model2->setObjectiveScale(saveScale);
1839#ifdef ABC_INHERIT
1840                      AbcSimplex * abcModel2=new AbcSimplex(*model2);
1841                      if (interrupt)
1842                        currentAbcModel = abcModel2;
1843                      if (abcSimplex_) {
1844                        // move factorization stuff
1845                        abcModel2->factorization()->synchronize(model2->factorization(),abcModel2);
1846                      }
1847                      abcModel2->startPermanentArrays();
1848                      abcModel2->setAbcState(CLP_ABC_WANTED);
1849#if ABC_PARALLEL
1850                      int parallelMode=1;
1851                      printf("Parallel mode %d\n",parallelMode);
1852                      abcModel2->setParallelMode(parallelMode);
1853#endif
1854                      //if (processTune>0&&processTune<8)
1855                      //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|65536*processTune);
1856                      abcModel2->doAbcDual();
1857                      abcModel2->moveStatusToClp(model2);
1858                      //ClpModel::stopPermanentArrays();
1859                      model2->setSpecialOptions(model2->specialOptions()&~65536);
1860                      //model2->dual();
1861                      //model2->setNumberIterations(abcModel2->numberIterations()+model2->numberIterations());
1862                      delete abcModel2;
1863#endif
1864                      memcpy(model2->objective(),saveObj,numberColumns*sizeof(double));
1865                      //delete [] saveDuals;
1866                      delete [] saveObj;
1867                      model2->dual(2);
1868                    } // end dubious idiot
1869               }
1870          }
1871#endif
1872          // ?
1873          if (doCrash) {
1874               switch(doCrash) {
1875                    // standard
1876               case 1:
1877                    model2->crash(1000, 1);
1878                    break;
1879                    // As in paper by Solow and Halim (approx)
1880               case 2:
1881                    model2->crash(model2->dualBound(), 0);
1882                    break;
1883                    // My take on it
1884               case 3:
1885                    model2->crash(model2->dualBound(), -1);
1886                    break;
1887                    // Just put free in basis
1888               case 4:
1889                    model2->crash(0.0, 3);
1890                    break;
1891               }
1892          }
1893#ifndef SLIM_CLP
1894          if (doSlp > 0 && objective_->type() == 2) {
1895               model2->nonlinearSLP(doSlp, 1.0e-5);
1896          }
1897#endif
1898#ifndef LACI_TRY
1899          if (options.getSpecialOption(1) != 2 ||
1900                    options.getExtraInfo(1) < 1000000) {
1901               if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1902                    // See if original wanted vector
1903                    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1904                    ClpMatrixBase * matrix = model2->clpMatrix();
1905                    if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
1906                         ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1907                         clpMatrix->makeSpecialColumnCopy();
1908                         //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
1909                         model2->primal(primalStartup);
1910                         clpMatrix->releaseSpecialColumnCopy();
1911                    } else {
1912#ifndef ABC_INHERIT
1913                        model2->primal(primalStartup);
1914#else
1915                        model2->dealWithAbc(1,primalStartup,interrupt);
1916#endif
1917                    }
1918               } else {
1919#ifndef ABC_INHERIT
1920                    model2->primal(primalStartup);
1921#else
1922                    model2->dealWithAbc(1,primalStartup,interrupt);
1923#endif
1924               }
1925          }
1926#endif
1927          time2 = CoinCpuTime();
1928          timeCore = time2 - timeX;
1929          handler_->message(CLP_INTERVAL_TIMING, messages_)
1930                    << "Primal" << timeCore << time2 - time1
1931                    << CoinMessageEol;
1932          timeX = time2;
1933     } else if (method == ClpSolve::usePrimalorSprint) {
1934          // Sprint
1935          /*
1936            This driver implements what I called Sprint when I introduced the idea
1937            many years ago.  Cplex calls it "sifting" which I think is just as silly.
1938            When I thought of this trivial idea
1939            it reminded me of an LP code of the 60's called sprint which after
1940            every factorization took a subset of the matrix into memory (all
1941            64K words!) and then iterated very fast on that subset.  On the
1942            problems of those days it did not work very well, but it worked very
1943            well on aircrew scheduling problems where there were very large numbers
1944            of columns all with the same flavor.
1945          */
1946
1947          /* The idea works best if you can get feasible easily.  To make it
1948             more general we can add in costed slacks */
1949
1950          int originalNumberColumns = model2->numberColumns();
1951          int numberRows = model2->numberRows();
1952          ClpSimplex * originalModel2 = model2;
1953
1954          // We will need arrays to choose variables.  These are too big but ..
1955          double * weight = new double [numberRows+originalNumberColumns];
1956          int * sort = new int [numberRows+originalNumberColumns];
1957          int numberSort = 0;
1958          // We are going to add slacks to get feasible.
1959          // initial list will just be artificials
1960          int iColumn;
1961          const double * columnLower = model2->columnLower();
1962          const double * columnUpper = model2->columnUpper();
1963          double * columnSolution = model2->primalColumnSolution();
1964
1965          // See if we have costed slacks
1966          int * negSlack = new int[numberRows];
1967          int * posSlack = new int[numberRows];
1968          int iRow;
1969          for (iRow = 0; iRow < numberRows; iRow++) {
1970               negSlack[iRow] = -1;
1971               posSlack[iRow] = -1;
1972          }
1973          const double * element = model2->matrix()->getElements();
1974          const int * row = model2->matrix()->getIndices();
1975          const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
1976          const int * columnLength = model2->matrix()->getVectorLengths();
1977          //bool allSlack = (numberRowsBasic==numberRows);
1978          for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
1979               if (!columnSolution[iColumn] || fabs(columnSolution[iColumn]) > 1.0e20) {
1980                    double value = 0.0;
1981                    if (columnLower[iColumn] > 0.0)
1982                         value = columnLower[iColumn];
1983                    else if (columnUpper[iColumn] < 0.0)
1984                         value = columnUpper[iColumn];
1985                    columnSolution[iColumn] = value;
1986               }
1987               if (columnLength[iColumn] == 1) {
1988                    int jRow = row[columnStart[iColumn]];
1989                    if (!columnLower[iColumn]) {
1990                         if (element[columnStart[iColumn]] > 0.0 && posSlack[jRow] < 0)
1991                              posSlack[jRow] = iColumn;
1992                         else if (element[columnStart[iColumn]] < 0.0 && negSlack[jRow] < 0)
1993                              negSlack[jRow] = iColumn;
1994                    } else if (!columnUpper[iColumn]) {
1995                         if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0)
1996                              posSlack[jRow] = iColumn;
1997                         else if (element[columnStart[iColumn]] > 0.0 && negSlack[jRow] < 0)
1998                              negSlack[jRow] = iColumn;
1999                    }
2000               }
2001          }
2002          // now see what that does to row solution
2003          double * rowSolution = model2->primalRowSolution();
2004          CoinZeroN (rowSolution, numberRows);
2005          model2->clpMatrix()->times(1.0, columnSolution, rowSolution);
2006          // See if we can adjust using costed slacks
2007          double penalty = CoinMax(1.0e5, CoinMin(infeasibilityCost_ * 0.01, 1.0e10)) * optimizationDirection_;
2008          const double * lower = model2->rowLower();
2009          const double * upper = model2->rowUpper();
2010          for (iRow = 0; iRow < numberRows; iRow++) {
2011               if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
2012                    int jColumn = posSlack[iRow];
2013                    if (jColumn >= 0) {
2014                         if (columnSolution[jColumn])
2015                              continue;
2016                         double difference = lower[iRow] - rowSolution[iRow];
2017                         double elementValue = element[columnStart[jColumn]];
2018                         if (elementValue > 0.0) {
2019                              double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
2020                              columnSolution[jColumn] = movement;
2021                              rowSolution[iRow] += movement * elementValue;
2022                         } else {
2023                              double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
2024                              columnSolution[jColumn] = movement;
2025                              rowSolution[iRow] += movement * elementValue;
2026                         }
2027                    }
2028               } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
2029                    int jColumn = negSlack[iRow];
2030                    if (jColumn >= 0) {
2031                         if (columnSolution[jColumn])
2032                              continue;
2033                         double difference = upper[iRow] - rowSolution[iRow];
2034                         double elementValue = element[columnStart[jColumn]];
2035                         if (elementValue < 0.0) {
2036                              double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
2037                              columnSolution[jColumn] = movement;
2038                              rowSolution[iRow] += movement * elementValue;
2039                         } else {
2040                              double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
2041                              columnSolution[jColumn] = movement;
2042                              rowSolution[iRow] += movement * elementValue;
2043                         }
2044                    }
2045               }
2046          }
2047          delete [] negSlack;
2048          delete [] posSlack;
2049          int nRow = numberRows;
2050          bool network = false;
2051          if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
2052               network = true;
2053               nRow *= 2;
2054          }
2055          int * addStarts = new int [nRow+1];
2056          int * addRow = new int[nRow];
2057          double * addElement = new double[nRow];
2058          addStarts[0] = 0;
2059          int numberArtificials = 0;
2060          int numberAdd = 0;
2061          double * addCost = new double [numberRows];
2062          for (iRow = 0; iRow < numberRows; iRow++) {
2063               if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
2064                    addRow[numberAdd] = iRow;
2065                    addElement[numberAdd++] = 1.0;
2066                    if (network) {
2067                         addRow[numberAdd] = numberRows;
2068                         addElement[numberAdd++] = -1.0;
2069                    }
2070                    addCost[numberArtificials] = penalty;
2071                    numberArtificials++;
2072                    addStarts[numberArtificials] = numberAdd;
2073               } else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
2074                    addRow[numberAdd] = iRow;
2075                    addElement[numberAdd++] = -1.0;
2076                    if (network) {
2077                         addRow[numberAdd] = numberRows;
2078                         addElement[numberAdd++] = 1.0;
2079                    }
2080                    addCost[numberArtificials] = penalty;
2081                    numberArtificials++;
2082                    addStarts[numberArtificials] = numberAdd;
2083               }
2084          }
2085          if (numberArtificials) {
2086               // need copy so as not to disturb original
2087               model2 = new ClpSimplex(*model2);
2088               if (network) {
2089                    // network - add a null row
2090                    model2->addRow(0, NULL, NULL, -COIN_DBL_MAX, COIN_DBL_MAX);
2091                    numberRows++;
2092               }
2093               model2->addColumns(numberArtificials, NULL, NULL, addCost,
2094                                  addStarts, addRow, addElement);
2095          }
2096          delete [] addStarts;
2097          delete [] addRow;
2098          delete [] addElement;
2099          delete [] addCost;
2100          // look at rhs to see if to perturb
2101          double largest = 0.0;
2102          double smallest = 1.0e30;
2103          for (iRow = 0; iRow < numberRows; iRow++) {
2104               double value;
2105               value = fabs(model2->rowLower_[iRow]);
2106               if (value && value < 1.0e30) {
2107                    largest = CoinMax(largest, value);
2108                    smallest = CoinMin(smallest, value);
2109               }
2110               value = fabs(model2->rowUpper_[iRow]);
2111               if (value && value < 1.0e30) {
2112                    largest = CoinMax(largest, value);
2113                    smallest = CoinMin(smallest, value);
2114               }
2115          }
2116          double * saveLower = NULL;
2117          double * saveUpper = NULL;
2118          if (largest < 2.01 * smallest) {
2119               // perturb - so switch off standard
2120               model2->setPerturbation(100);
2121               saveLower = new double[numberRows];
2122               CoinMemcpyN(model2->rowLower_, numberRows, saveLower);
2123               saveUpper = new double[numberRows];
2124               CoinMemcpyN(model2->rowUpper_, numberRows, saveUpper);
2125               double * lower = model2->rowLower();
2126               double * upper = model2->rowUpper();
2127               for (iRow = 0; iRow < numberRows; iRow++) {
2128                    double lowerValue = lower[iRow], upperValue = upper[iRow];
2129                    double value = randomNumberGenerator_.randomDouble();
2130                    if (upperValue > lowerValue + primalTolerance_) {
2131                         if (lowerValue > -1.0e20 && lowerValue)
2132                              lowerValue -= value * 1.0e-4 * fabs(lowerValue);
2133                         if (upperValue < 1.0e20 && upperValue)
2134                              upperValue += value * 1.0e-4 * fabs(upperValue);
2135                    } else if (upperValue > 0.0) {
2136                         upperValue -= value * 1.0e-4 * fabs(lowerValue);
2137                         lowerValue -= value * 1.0e-4 * fabs(lowerValue);
2138                    } else if (upperValue < 0.0) {
2139                         upperValue += value * 1.0e-4 * fabs(lowerValue);
2140                         lowerValue += value * 1.0e-4 * fabs(lowerValue);
2141                    } else {
2142                    }
2143                    lower[iRow] = lowerValue;
2144                    upper[iRow] = upperValue;
2145               }
2146          }
2147          int i;
2148          // Just do this number of passes in Sprint
2149          if (doSprint > 0)
2150               maxSprintPass = options.getExtraInfo(1);
2151          // but if big use to get ratio
2152          double ratio = 3;
2153          if (maxSprintPass > 1000) {
2154               ratio = static_cast<double> (maxSprintPass) * 0.0001;
2155               ratio = CoinMax(ratio, 1.1);
2156               maxSprintPass = maxSprintPass % 1000;
2157#ifdef COIN_DEVELOP
2158               printf("%d passes wanted with ratio of %g\n", maxSprintPass, ratio);
2159#endif
2160          }
2161          // Just take this number of columns in small problem
2162          int smallNumberColumns = static_cast<int> (CoinMin(ratio * numberRows, static_cast<double> (numberColumns)));
2163          smallNumberColumns = CoinMax(smallNumberColumns, 3000);
2164          smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
2165          //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
2166          //smallNumberColumns = CoinMax(smallNumberColumns,3000);
2167          //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
2168          // redo as may have changed
2169          columnLower = model2->columnLower();
2170          columnUpper = model2->columnUpper();
2171          columnSolution = model2->primalColumnSolution();
2172          // Set up initial list
2173          numberSort = 0;
2174          if (numberArtificials) {
2175               numberSort = numberArtificials;
2176               for (i = 0; i < numberSort; i++)
2177                    sort[i] = i + originalNumberColumns;
2178          }
2179          // maybe a solution there already
2180          for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2181               if (model2->getColumnStatus(iColumn) == basic)
2182                    sort[numberSort++] = iColumn;
2183          }
2184          for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2185               if (model2->getColumnStatus(iColumn) != basic) {
2186                    if (columnSolution[iColumn] > columnLower[iColumn] &&
2187                              columnSolution[iColumn] < columnUpper[iColumn] &&
2188                              columnSolution[iColumn])
2189                         sort[numberSort++] = iColumn;
2190               }
2191          }
2192          numberSort = CoinMin(numberSort, smallNumberColumns);
2193
2194          int numberColumns = model2->numberColumns();
2195          double * fullSolution = model2->primalColumnSolution();
2196
2197
2198          int iPass;
2199          double lastObjective[] = {1.0e31,1.0e31};
2200          // It will be safe to allow dense
2201          model2->setInitialDenseFactorization(true);
2202
2203          // We will be using all rows
2204          int * whichRows = new int [numberRows];
2205          for (iRow = 0; iRow < numberRows; iRow++)
2206               whichRows[iRow] = iRow;
2207          double originalOffset;
2208          model2->getDblParam(ClpObjOffset, originalOffset);
2209          int totalIterations = 0;
2210          double lastSumArtificials = COIN_DBL_MAX;
2211          int originalMaxSprintPass = maxSprintPass;
2212          maxSprintPass = 20; // so we do that many if infeasible
2213          for (iPass = 0; iPass < maxSprintPass; iPass++) {
2214               //printf("Bug until submodel new version\n");
2215               //CoinSort_2(sort,sort+numberSort,weight);
2216               // Create small problem
2217               ClpSimplex small(model2, numberRows, whichRows, numberSort, sort);
2218               small.setPerturbation(model2->perturbation());
2219               small.setInfeasibilityCost(model2->infeasibilityCost());
2220               if (model2->factorizationFrequency() == 200) {
2221                    // User did not touch preset
2222                    small.defaultFactorizationFrequency();
2223               }
2224               // now see what variables left out do to row solution
2225               double * rowSolution = model2->primalRowSolution();
2226               double * sumFixed = new double[numberRows];
2227               CoinZeroN (sumFixed, numberRows);
2228               int iRow, iColumn;
2229               // zero out ones in small problem
2230               for (iColumn = 0; iColumn < numberSort; iColumn++) {
2231                    int kColumn = sort[iColumn];
2232                    fullSolution[kColumn] = 0.0;
2233               }
2234               // Get objective offset
2235               const double * objective = model2->objective();
2236               double offset = 0.0;
2237               for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
2238                    offset += fullSolution[iColumn] * objective[iColumn];
2239#if 0
2240               // Set artificials to zero if first time close to zero
2241               for (iColumn = originalNumberColumns; iColumn < numberColumns; iColumn++) {
2242                 if (fullSolution[iColumn]<primalTolerance_&&objective[iColumn]==penalty) {
2243                   model2->objective()[iColumn]=2.0*penalty;
2244                   fullSolution[iColumn]=0.0;
2245                 }
2246               }
2247#endif
2248               small.setDblParam(ClpObjOffset, originalOffset - offset);
2249               model2->clpMatrix()->times(1.0, fullSolution, sumFixed);
2250
2251               double * lower = small.rowLower();
2252               double * upper = small.rowUpper();
2253               for (iRow = 0; iRow < numberRows; iRow++) {
2254                    if (lower[iRow] > -1.0e50)
2255                         lower[iRow] -= sumFixed[iRow];
2256                    if (upper[iRow] < 1.0e50)
2257                         upper[iRow] -= sumFixed[iRow];
2258                    rowSolution[iRow] -= sumFixed[iRow];
2259               }
2260               delete [] sumFixed;
2261               // Solve
2262               if (interrupt)
2263                    currentModel = &small;
2264               small.defaultFactorizationFrequency();
2265               if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
2266                    // See if original wanted vector
2267                    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
2268                    ClpMatrixBase * matrix = small.clpMatrix();
2269                    if (dynamic_cast< ClpPackedMatrix*>(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
2270                         ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
2271                         clpMatrix->makeSpecialColumnCopy();
2272                         small.primal(1);
2273                         clpMatrix->releaseSpecialColumnCopy();
2274                    } else {
2275#if 1
2276#ifdef ABC_INHERIT
2277                      //small.writeMps("try.mps");
2278                      if (iPass||!numberArtificials) 
2279                         small.dealWithAbc(1,1);
2280                       else 
2281                         small.dealWithAbc(0,0);
2282#else
2283                      if (iPass||!numberArtificials) 
2284                         small.primal(1);
2285                      else
2286                         small.dual(0);
2287#endif
2288                      if (small.problemStatus())
2289                        small.dual(0);
2290#else
2291                         int numberColumns = small.numberColumns();
2292                         int numberRows = small.numberRows();
2293                         // Use dual region
2294                         double * rhs = small.dualRowSolution();
2295                         int * whichRow = new int[3*numberRows];
2296                         int * whichColumn = new int[2*numberColumns];
2297                         int nBound;
2298                         ClpSimplex * small2 = ((ClpSimplexOther *) (&small))->crunch(rhs, whichRow, whichColumn,
2299                                               nBound, false, false);
2300                         if (small2) {
2301#ifdef ABC_INHERIT
2302                              small2->dealWithAbc(1,1);
2303#else
2304                              small.primal(1);
2305#endif
2306                              if (small2->problemStatus() == 0) {
2307                                   small.setProblemStatus(0);
2308                                   ((ClpSimplexOther *) (&small))->afterCrunch(*small2, whichRow, whichColumn, nBound);
2309                              } else {
2310#ifdef ABC_INHERIT
2311                                   small2->dealWithAbc(1,1);
2312#else
2313                                   small.primal(1);
2314#endif
2315                                   if (small2->problemStatus())
2316                                        small.primal(1);
2317                              }
2318                              delete small2;
2319                         } else {
2320                              small.primal(1);
2321                         }
2322                         delete [] whichRow;
2323                         delete [] whichColumn;
2324#endif
2325                    }
2326               } else {
2327                    small.primal(1);
2328               }
2329               totalIterations += small.numberIterations();
2330               // move solution back
2331               const double * solution = small.primalColumnSolution();
2332               for (iColumn = 0; iColumn < numberSort; iColumn++) {
2333                    int kColumn = sort[iColumn];
2334                    model2->setColumnStatus(kColumn, small.getColumnStatus(iColumn));
2335                    fullSolution[kColumn] = solution[iColumn];
2336               }
2337               for (iRow = 0; iRow < numberRows; iRow++)
2338                    model2->setRowStatus(iRow, small.getRowStatus(iRow));
2339               CoinMemcpyN(small.primalRowSolution(),
2340                           numberRows, model2->primalRowSolution());
2341               double sumArtificials = 0.0;
2342               for (i = 0; i < numberArtificials; i++)
2343                    sumArtificials += fullSolution[i + originalNumberColumns];
2344               if (sumArtificials && iPass > 5 && sumArtificials >= lastSumArtificials) {
2345                    // increase costs
2346                    double * cost = model2->objective() + originalNumberColumns;
2347                    double newCost = CoinMin(1.0e10, cost[0] * 1.5);
2348                    for (i = 0; i < numberArtificials; i++)
2349                         cost[i] = newCost;
2350               }
2351               lastSumArtificials = sumArtificials;
2352               // get reduced cost for large problem
2353               double * djs = model2->dualColumnSolution();
2354               CoinMemcpyN(model2->objective(), numberColumns, djs);
2355               model2->clpMatrix()->transposeTimes(-1.0, small.dualRowSolution(), djs);
2356               int numberNegative = 0;
2357               double sumNegative = 0.0;
2358               // now massage weight so all basic in plus good djs
2359               // first count and do basic
2360               numberSort = 0;
2361               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2362                    double dj = djs[iColumn] * optimizationDirection_;
2363                    double value = fullSolution[iColumn];
2364                    if (model2->getColumnStatus(iColumn) == ClpSimplex::basic) {
2365                         sort[numberSort++] = iColumn;
2366                    } else if (dj < -dualTolerance_ && value < columnUpper[iColumn]) {
2367                         numberNegative++;
2368                         sumNegative -= dj;
2369                    } else if (dj > dualTolerance_ && value > columnLower[iColumn]) {
2370                         numberNegative++;
2371                         sumNegative += dj;
2372                    }
2373               }
2374               handler_->message(CLP_SPRINT, messages_)
2375                         << iPass + 1 << small.numberIterations() << small.objectiveValue() << sumNegative
2376                         << numberNegative
2377                         << CoinMessageEol;
2378               if (sumArtificials < 1.0e-8 && originalMaxSprintPass >= 0) {
2379                    maxSprintPass = iPass + originalMaxSprintPass;
2380                    originalMaxSprintPass = -1;
2381               }
2382               if (iPass > 20)
2383                    sumArtificials = 0.0;
2384               if ((small.objectiveValue()*optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 5 && sumArtificials < 1.0e-8) ||
2385                         (!small.numberIterations() && iPass) ||
2386                         iPass == maxSprintPass - 1 || small.status() == 3) {
2387
2388                    break; // finished
2389               } else {
2390                    lastObjective[1] = lastObjective[0];
2391                    lastObjective[0] = small.objectiveValue() * optimizationDirection_;
2392                    double tolerance;
2393                    double averageNegDj = sumNegative / static_cast<double> (numberNegative + 1);
2394                    if (numberNegative + numberSort > smallNumberColumns)
2395                         tolerance = -dualTolerance_;
2396                    else
2397                         tolerance = 10.0 * averageNegDj;
2398                    int saveN = numberSort;
2399                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2400                         double dj = djs[iColumn] * optimizationDirection_;
2401                         double value = fullSolution[iColumn];
2402                         if (model2->getColumnStatus(iColumn) != ClpSimplex::basic) {
2403                              if (dj < -dualTolerance_ && value < columnUpper[iColumn])
2404                                   dj = dj;
2405                              else if (dj > dualTolerance_ && value > columnLower[iColumn])
2406                                   dj = -dj;
2407                              else if (columnUpper[iColumn] > columnLower[iColumn])
2408                                   dj = fabs(dj);
2409                              else
2410                                   dj = 1.0e50;
2411                              if (dj < tolerance) {
2412                                   weight[numberSort] = dj;
2413                                   sort[numberSort++] = iColumn;
2414                              }
2415                         }
2416                    }
2417                    // sort
2418                    CoinSort_2(weight + saveN, weight + numberSort, sort + saveN);
2419                    numberSort = CoinMin(smallNumberColumns, numberSort);
2420               }
2421          }
2422          if (interrupt)
2423               currentModel = model2;
2424          for (i = 0; i < numberArtificials; i++)
2425               sort[i] = i + originalNumberColumns;
2426          model2->deleteColumns(numberArtificials, sort);
2427          if (network) {
2428               int iRow = numberRows - 1;
2429               model2->deleteRows(1, &iRow);
2430          }
2431          delete [] weight;
2432          delete [] sort;
2433          delete [] whichRows;
2434          if (saveLower) {
2435               // unperturb and clean
2436               for (iRow = 0; iRow < numberRows; iRow++) {
2437                    double diffLower = saveLower[iRow] - model2->rowLower_[iRow];
2438                    double diffUpper = saveUpper[iRow] - model2->rowUpper_[iRow];
2439                    model2->rowLower_[iRow] = saveLower[iRow];
2440                    model2->rowUpper_[iRow] = saveUpper[iRow];
2441                    if (diffLower)
2442                         assert (!diffUpper || fabs(diffLower - diffUpper) < 1.0e-5);
2443                    else
2444                         diffLower = diffUpper;
2445                    model2->rowActivity_[iRow] += diffLower;
2446               }
2447               delete [] saveLower;
2448               delete [] saveUpper;
2449          }
2450#ifdef ABC_INHERIT
2451          model2->dealWithAbc(1,1);
2452#else
2453          model2->primal(1);
2454#endif
2455          model2->setPerturbation(savePerturbation);
2456          if (model2 != originalModel2) {
2457               originalModel2->moveInfo(*model2);
2458               delete model2;
2459               model2 = originalModel2;
2460          }
2461          time2 = CoinCpuTime();
2462          timeCore = time2 - timeX;
2463          handler_->message(CLP_INTERVAL_TIMING, messages_)
2464                    << "Sprint" << timeCore << time2 - time1
2465                    << CoinMessageEol;
2466          timeX = time2;
2467          model2->setNumberIterations(model2->numberIterations() + totalIterations);
2468     } else if (method == ClpSolve::useBarrier || method == ClpSolve::useBarrierNoCross) {
2469#ifndef SLIM_CLP
2470          //printf("***** experimental pretty crude barrier\n");
2471          //#define SAVEIT 2
2472#ifndef SAVEIT
2473#define BORROW
2474#endif
2475#ifdef BORROW
2476          ClpInterior barrier;
2477          barrier.borrowModel(*model2);
2478#else
2479          ClpInterior barrier(*model2);
2480#endif
2481          if (interrupt)
2482               currentModel2 = &barrier;
2483          if (barrier.numberRows()+barrier.numberColumns()>10000)
2484            barrier.setMaximumBarrierIterations(1000);
2485          int barrierOptions = options.getSpecialOption(4);
2486          int aggressiveGamma = 0;
2487          bool presolveInCrossover = false;
2488          bool scale = false;
2489          bool doKKT = false;
2490          bool forceFixing = false;
2491          int speed = 0;
2492          if (barrierOptions & 16) {
2493               barrierOptions &= ~16;
2494               doKKT = true;
2495          }
2496          if (barrierOptions&(32 + 64 + 128)) {
2497               aggressiveGamma = (barrierOptions & (32 + 64 + 128)) >> 5;
2498               barrierOptions &= ~(32 + 64 + 128);
2499          }
2500          if (barrierOptions & 256) {
2501               barrierOptions &= ~256;
2502               presolveInCrossover = true;
2503          }
2504          if (barrierOptions & 512) {
2505               barrierOptions &= ~512;
2506               forceFixing = true;
2507          }
2508          if (barrierOptions & 1024) {
2509               barrierOptions &= ~1024;
2510               barrier.setProjectionTolerance(1.0e-9);
2511          }
2512          if (barrierOptions&(2048 | 4096)) {
2513               speed = (barrierOptions & (2048 | 4096)) >> 11;
2514               barrierOptions &= ~(2048 | 4096);
2515          }
2516          if (barrierOptions & 8) {
2517               barrierOptions &= ~8;
2518               scale = true;
2519          }
2520          // If quadratic force KKT
2521          if (quadraticObj) {
2522               doKKT = true;
2523          }
2524          switch (barrierOptions) {
2525          case 0:
2526          default:
2527               if (!doKKT) {
2528                    ClpCholeskyBase * cholesky = new ClpCholeskyBase(options.getExtraInfo(1));
2529                    cholesky->setIntegerParameter(0, speed);
2530                    barrier.setCholesky(cholesky);
2531               } else {
2532                    ClpCholeskyBase * cholesky = new ClpCholeskyBase();
2533                    cholesky->setKKT(true);
2534                    barrier.setCholesky(cholesky);
2535               }
2536               break;
2537          case 1:
2538               if (!doKKT) {
2539                    ClpCholeskyDense * cholesky = new ClpCholeskyDense();
2540                    barrier.setCholesky(cholesky);
2541               } else {
2542                    ClpCholeskyDense * cholesky = new ClpCholeskyDense();
2543                    cholesky->setKKT(true);
2544                    barrier.setCholesky(cholesky);
2545               }
2546               break;
2547#ifdef COIN_HAS_WSMP
2548          case 2: {
2549               ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100, model2->numberRows() / 10));
2550               barrier.setCholesky(cholesky);
2551               assert (!doKKT);
2552          }
2553          break;
2554          case 3:
2555               if (!doKKT) {
2556                    ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
2557                    barrier.setCholesky(cholesky);
2558               } else {
2559                    ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
2560                    barrier.setCholesky(cholesky);
2561               }
2562               break;
2563#endif
2564#ifdef UFL_BARRIER
2565          case 4:
2566               if (!doKKT) {
2567                    ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
2568                    barrier.setCholesky(cholesky);
2569               } else {
2570                    ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
2571                    cholesky->setKKT(true);
2572                    barrier.setCholesky(cholesky);
2573               }
2574               break;
2575#endif
2576#ifdef TAUCS_BARRIER
2577          case 5: {
2578               ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
2579               barrier.setCholesky(cholesky);
2580               assert (!doKKT);
2581          }
2582          break;
2583#endif
2584#ifdef COIN_HAS_MUMPS
2585          case 6: {
2586               ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
2587               barrier.setCholesky(cholesky);
2588               assert (!doKKT);
2589          }
2590          break;
2591#endif
2592          }
2593          int numberRows = model2->numberRows();
2594          int numberColumns = model2->numberColumns();
2595          int saveMaxIts = model2->maximumIterations();
2596          if (saveMaxIts < 1000) {
2597               barrier.setMaximumBarrierIterations(saveMaxIts);
2598               model2->setMaximumIterations(10000000);
2599          }
2600#ifndef SAVEIT
2601          //barrier.setDiagonalPerturbation(1.0e-25);
2602          if (aggressiveGamma) {
2603               switch (aggressiveGamma) {
2604               case 1:
2605                    barrier.setGamma(1.0e-5);
2606                    barrier.setDelta(1.0e-5);
2607                    break;
2608               case 2:
2609                    barrier.setGamma(1.0e-7);
2610                    break;
2611               case 3:
2612                    barrier.setDelta(1.0e-5);
2613                    break;
2614               case 4:
2615                    barrier.setGamma(1.0e-3);
2616                    barrier.setDelta(1.0e-3);
2617                    break;
2618               case 5:
2619                    barrier.setGamma(1.0e-3);
2620                    break;
2621               case 6:
2622                    barrier.setDelta(1.0e-3);
2623                    break;
2624               }
2625          }
2626          if (scale)
2627               barrier.scaling(1);
2628          else
2629               barrier.scaling(0);
2630          barrier.primalDual();
2631#elif SAVEIT==1
2632          barrier.primalDual();
2633#else
2634          model2->restoreModel("xx.save");
2635          // move solutions
2636          CoinMemcpyN(model2->primalRowSolution(),
2637                      numberRows, barrier.primalRowSolution());
2638          CoinMemcpyN(model2->dualRowSolution(),
2639                      numberRows, barrier.dualRowSolution());
2640          CoinMemcpyN(model2->primalColumnSolution(),
2641                      numberColumns, barrier.primalColumnSolution());
2642          CoinMemcpyN(model2->dualColumnSolution(),
2643                      numberColumns, barrier.dualColumnSolution());
2644#endif
2645          time2 = CoinCpuTime();
2646          timeCore = time2 - timeX;
2647          handler_->message(CLP_INTERVAL_TIMING, messages_)
2648                    << "Barrier" << timeCore << time2 - time1
2649                    << CoinMessageEol;
2650          timeX = time2;
2651          int maxIts = barrier.maximumBarrierIterations();
2652          int barrierStatus = barrier.status();
2653          double gap = barrier.complementarityGap();
2654          // get which variables are fixed
2655          double * saveLower = NULL;
2656          double * saveUpper = NULL;
2657          ClpPresolve pinfo2;
2658          ClpSimplex * saveModel2 = NULL;
2659          bool extraPresolve = false;
2660          int numberFixed = barrier.numberFixed();
2661          if (numberFixed) {
2662               int numberRows = barrier.numberRows();
2663               int numberColumns = barrier.numberColumns();
2664               int numberTotal = numberRows + numberColumns;
2665               saveLower = new double [numberTotal];
2666               saveUpper = new double [numberTotal];
2667               CoinMemcpyN(barrier.columnLower(), numberColumns, saveLower);
2668               CoinMemcpyN(barrier.rowLower(), numberRows, saveLower + numberColumns);
2669               CoinMemcpyN(barrier.columnUpper(), numberColumns, saveUpper);
2670               CoinMemcpyN(barrier.rowUpper(), numberRows, saveUpper + numberColumns);
2671          }
2672          if (((numberFixed * 20 > barrier.numberRows() && numberFixed > 5000) || forceFixing) &&
2673                    presolveInCrossover) {
2674               // may as well do presolve
2675               if (!forceFixing) {
2676                    barrier.fixFixed();
2677               } else {
2678                    // Fix
2679                    int n = barrier.numberColumns();
2680                    double * lower = barrier.columnLower();
2681                    double * upper = barrier.columnUpper();
2682                    double * solution = barrier.primalColumnSolution();
2683                    int nFix = 0;
2684                    for (int i = 0; i < n; i++) {
2685                         if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
2686                              double value = solution[i];
2687                              if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
2688                                   solution[i] = lower[i];
2689                                   upper[i] = lower[i];
2690                                   nFix++;
2691                              } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
2692                                   solution[i] = upper[i];
2693                                   lower[i] = upper[i];
2694                                   nFix++;
2695                              }
2696                         }
2697                    }
2698#ifdef CLP_INVESTIGATE
2699                    printf("%d columns fixed\n", nFix);
2700#endif
2701                    int nr = barrier.numberRows();
2702                    lower = barrier.rowLower();
2703                    upper = barrier.rowUpper();
2704                    solution = barrier.primalRowSolution();
2705                    nFix = 0;
2706                    for (int i = 0; i < nr; i++) {
2707                         if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
2708                              double value = solution[i];
2709                              if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
2710                                   solution[i] = lower[i];
2711                                   upper[i] = lower[i];
2712                                   nFix++;
2713                              } else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
2714                                   solution[i] = upper[i];
2715                                   lower[i] = upper[i];
2716                                   nFix++;
2717                              }
2718                         }
2719                    }
2720#ifdef CLP_INVESTIGATE
2721                    printf("%d row slacks fixed\n", nFix);
2722#endif
2723               }
2724               saveModel2 = model2;
2725               extraPresolve = true;
2726          } else if (numberFixed) {
2727               // Set fixed to bounds (may have restored earlier solution)
2728               if (!forceFixing) {
2729                    barrier.fixFixed(false);
2730               } else {
2731                    // Fix
2732                    int n = barrier.numberColumns();
2733                    double * lower = barrier.columnLower();
2734                    double * upper = barrier.columnUpper();
2735                    double * solution = barrier.primalColumnSolution();
2736                    int nFix = 0;
2737                    for (int i = 0; i < n; i++) {
2738                         if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
2739                              double value = solution[i];
2740                              if (value < lower[i] + 1.0e-8 && value - lower[i] < upper[i] - value) {
2741                                   solution[i] = lower[i];
2742                                   upper[i] = lower[i];
2743                                   nFix++;
2744                              } else if (value > upper[i] - 1.0e-8 && value - lower[i] > upper[i] - value) {
2745                                   solution[i] = upper[i];
2746                                   lower[i] = upper[i];
2747                                   nFix++;
2748                              } else {
2749                                   //printf("fixcol %d %g <= %g <= %g\n",
2750                                   //     i,lower[i],solution[i],upper[i]);
2751                              }
2752                         }
2753                    }
2754#ifdef CLP_INVESTIGATE
2755                    printf("%d columns fixed\n", nFix);
2756#endif
2757                    int nr = barrier.numberRows();
2758                    lower = barrier.rowLower();
2759                    upper = barrier.rowUpper();
2760                    solution = barrier.primalRowSolution();
2761                    nFix = 0;
2762                    for (int i = 0; i < nr; i++) {
2763                         if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
2764                              double value = solution[i];
2765                              if (value < lower[i] + 1.0e-5 && value - lower[i] < upper[i] - value) {
2766                                   solution[i] = lower[i];
2767                                   upper[i] = lower[i];
2768                                   nFix++;
2769                              } else if (value > upper[i] - 1.0e-5 && value - lower[i] > upper[i] - value) {
2770                                   solution[i] = upper[i];
2771                                   lower[i] = upper[i];
2772                                   nFix++;
2773                              } else {
2774                                   //printf("fixrow %d %g <= %g <= %g\n",
2775                                   //     i,lower[i],solution[i],upper[i]);
2776                              }
2777                         }
2778                    }
2779#ifdef CLP_INVESTIGATE
2780                    printf("%d row slacks fixed\n", nFix);
2781#endif
2782               }
2783          }
2784#ifdef BORROW
2785          int saveNumberIterations = barrier.numberIterations();
2786          barrier.returnModel(*model2);
2787          double * rowPrimal = new double [numberRows];
2788          double * columnPrimal = new double [numberColumns];
2789          double * rowDual = new double [numberRows];
2790          double * columnDual = new double [numberColumns];
2791          // move solutions other way
2792          CoinMemcpyN(model2->primalRowSolution(),
2793                      numberRows, rowPrimal);
2794          CoinMemcpyN(model2->dualRowSolution(),
2795                      numberRows, rowDual);
2796          CoinMemcpyN(model2->primalColumnSolution(),
2797                      numberColumns, columnPrimal);
2798          CoinMemcpyN(model2->dualColumnSolution(),
2799                      numberColumns, columnDual);
2800#else
2801          double * rowPrimal = barrier.primalRowSolution();
2802          double * columnPrimal = barrier.primalColumnSolution();
2803          double * rowDual = barrier.dualRowSolution();
2804          double * columnDual = barrier.dualColumnSolution();
2805          // move solutions
2806          CoinMemcpyN(rowPrimal,
2807                      numberRows, model2->primalRowSolution());
2808          CoinMemcpyN(rowDual,
2809                      numberRows, model2->dualRowSolution());
2810          CoinMemcpyN(columnPrimal,
2811                      numberColumns, model2->primalColumnSolution());
2812          CoinMemcpyN(columnDual,
2813                      numberColumns, model2->dualColumnSolution());
2814#endif
2815          if (saveModel2) {
2816               // do presolve
2817               model2 = pinfo2.presolvedModel(*model2, dblParam_[ClpPresolveTolerance],
2818                                              false, 5, true);
2819               if (!model2) {
2820                    model2 = saveModel2;
2821                    saveModel2 = NULL;
2822                    int numberRows = model2->numberRows();
2823                    int numberColumns = model2->numberColumns();
2824                    CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
2825                    CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
2826                    delete [] saveLower;
2827                    CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
2828                    CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
2829                    delete [] saveUpper;
2830                    saveLower = NULL;
2831                    saveUpper = NULL;
2832               }
2833          }
2834          if (method == ClpSolve::useBarrier || barrierStatus < 0) {
2835               if (maxIts && barrierStatus < 4 && !quadraticObj) {
2836                    //printf("***** crossover - needs more thought on difficult models\n");
2837#if SAVEIT==1
2838                    model2->ClpSimplex::saveModel("xx.save");
2839#endif
2840                    // make sure no status left
2841                    model2->createStatus();
2842                    // solve
2843                    if (!forceFixing)
2844                         model2->setPerturbation(100);
2845                    if (model2->factorizationFrequency() == 200) {
2846                         // User did not touch preset
2847                         model2->defaultFactorizationFrequency();
2848                    }
2849#if 1 //ndef ABC_INHERIT //#if 1
2850                    // throw some into basis
2851                    if(!forceFixing) {
2852                         int numberRows = model2->numberRows();
2853                         int numberColumns = model2->numberColumns();
2854                         double * dsort = new double[numberColumns];
2855                         int * sort = new int[numberColumns];
2856                         int n = 0;
2857                         const double * columnLower = model2->columnLower();
2858                         const double * columnUpper = model2->columnUpper();
2859                         double * primalSolution = model2->primalColumnSolution();
2860                         const double * dualSolution = model2->dualColumnSolution();
2861                         double tolerance = 10.0 * primalTolerance_;
2862                         int i;
2863                         for ( i = 0; i < numberRows; i++)
2864                              model2->setRowStatus(i, superBasic);
2865                         for ( i = 0; i < numberColumns; i++) {
2866                              double distance = CoinMin(columnUpper[i] - primalSolution[i],
2867                                                        primalSolution[i] - columnLower[i]);
2868                              if (distance > tolerance) {
2869                                   if (fabs(dualSolution[i]) < 1.0e-5)
2870                                        distance *= 100.0;
2871                                   dsort[n] = -distance;
2872                                   sort[n++] = i;
2873                                   model2->setStatus(i, superBasic);
2874                              } else if (distance > primalTolerance_) {
2875                                   model2->setStatus(i, superBasic);
2876                              } else if (primalSolution[i] <= columnLower[i] + primalTolerance_) {
2877                                   model2->setStatus(i, atLowerBound);
2878                                   primalSolution[i] = columnLower[i];
2879                              } else {
2880                                   model2->setStatus(i, atUpperBound);
2881                                   primalSolution[i] = columnUpper[i];
2882                              }
2883                         }
2884                         CoinSort_2(dsort, dsort + n, sort);
2885                         n = CoinMin(numberRows, n);
2886                         for ( i = 0; i < n; i++) {
2887                              int iColumn = sort[i];
2888                              model2->setStatus(iColumn, basic);
2889                         }
2890                         delete [] sort;
2891                         delete [] dsort;
2892                         // model2->allSlackBasis();
2893                         if (gap < 1.0e-3 * static_cast<double> (numberRows + numberColumns)) {
2894                              if (saveUpper) {
2895                                   int numberRows = model2->numberRows();
2896                                   int numberColumns = model2->numberColumns();
2897                                   CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
2898                                   CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
2899                                   CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
2900                                   CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
2901                                   delete [] saveLower;
2902                                   delete [] saveUpper;
2903                                   saveLower = NULL;
2904                                   saveUpper = NULL;
2905                              }
2906                              int numberRows = model2->numberRows();
2907                              int numberColumns = model2->numberColumns();
2908#ifdef ABC_INHERIT
2909                              model2->checkSolution(0);
2910                              printf("%d primal infeasibilities summing to %g\n",
2911                                     model2->numberPrimalInfeasibilities(),
2912                                     model2->sumPrimalInfeasibilities());
2913                              model2->dealWithAbc(1,1);
2914                         }
2915                    }
2916#else
2917                              // just primal values pass
2918                              double saveScale = model2->objectiveScale();
2919                              model2->setObjectiveScale(1.0e-3);
2920                              model2->primal(2);
2921                              model2->setObjectiveScale(saveScale);
2922                              // save primal solution and copy back dual
2923                              CoinMemcpyN(model2->primalRowSolution(),
2924                                          numberRows, rowPrimal);
2925                              CoinMemcpyN(rowDual,
2926                                          numberRows, model2->dualRowSolution());
2927                              CoinMemcpyN(model2->primalColumnSolution(),
2928                                          numberColumns, columnPrimal);
2929                              CoinMemcpyN(columnDual,
2930                                          numberColumns, model2->dualColumnSolution());
2931                              //model2->primal(1);
2932                              // clean up reduced costs and flag variables
2933                              {
2934                                   double * dj = model2->dualColumnSolution();
2935                                   double * cost = model2->objective();
2936                                   double * saveCost = new double[numberColumns];
2937                                   CoinMemcpyN(cost, numberColumns, saveCost);
2938                                   double * saveLower = new double[numberColumns];
2939                                   double * lower = model2->columnLower();
2940                                   CoinMemcpyN(lower, numberColumns, saveLower);
2941                                   double * saveUpper = new double[numberColumns];
2942                                   double * upper = model2->columnUpper();
2943                                   CoinMemcpyN(upper, numberColumns, saveUpper);
2944                                   int i;
2945                                   double tolerance = 10.0 * dualTolerance_;
2946                                   for ( i = 0; i < numberColumns; i++) {
2947                                        if (model2->getStatus(i) == basic) {
2948                                             dj[i] = 0.0;
2949                                        } else if (model2->getStatus(i) == atLowerBound) {
2950                                             if (optimizationDirection_ * dj[i] < tolerance) {
2951                                                  if (optimizationDirection_ * dj[i] < 0.0) {
2952                                                       //if (dj[i]<-1.0e-3)
2953                                                       //printf("bad dj at lb %d %g\n",i,dj[i]);
2954                                                       cost[i] -= dj[i];
2955                                                       dj[i] = 0.0;
2956                                                  }
2957                                             } else {
2958                                                  upper[i] = lower[i];
2959                                             }
2960                                        } else if (model2->getStatus(i) == atUpperBound) {
2961                                             if (optimizationDirection_ * dj[i] > tolerance) {
2962                                                  if (optimizationDirection_ * dj[i] > 0.0) {
2963                                                       //if (dj[i]>1.0e-3)
2964                                                       //printf("bad dj at ub %d %g\n",i,dj[i]);
2965                                                       cost[i] -= dj[i];
2966                                                       dj[i] = 0.0;
2967                                                  }
2968                                             } else {
2969                                                  lower[i] = upper[i];
2970                                             }
2971                                        }
2972                                   }
2973                                   // just dual values pass
2974                                   //model2->setLogLevel(63);
2975                                   //model2->setFactorizationFrequency(1);
2976                                   model2->dual(2);
2977                                   CoinMemcpyN(saveCost, numberColumns, cost);
2978                                   delete [] saveCost;
2979                                   CoinMemcpyN(saveLower, numberColumns, lower);
2980                                   delete [] saveLower;
2981                                   CoinMemcpyN(saveUpper, numberColumns, upper);
2982                                   delete [] saveUpper;
2983                              }
2984                         }
2985                         // and finish
2986                         // move solutions
2987                         CoinMemcpyN(rowPrimal,
2988                                     numberRows, model2->primalRowSolution());
2989                         CoinMemcpyN(columnPrimal,
2990                                     numberColumns, model2->primalColumnSolution());
2991                    }
2992                    double saveScale = model2->objectiveScale();
2993                    model2->setObjectiveScale(1.0e-3);
2994                    model2->primal(2);
2995                    model2->setObjectiveScale(saveScale);
2996                    model2->primal(1);
2997#endif
2998#else
2999                    // just primal
3000#ifdef ABC_INHERIT
3001                    model2->checkSolution(0);
3002                    printf("%d primal infeasibilities summing to %g\n",
3003                           model2->numberPrimalInfeasibilities(),
3004                           model2->sumPrimalInfeasibilities());
3005                    model2->dealWithAbc(1,1);
3006#else
3007                    model2->primal(1);
3008#endif
3009                    //model2->primal(1);
3010#endif
3011               } else if (barrierStatus == 4) {
3012                    // memory problems
3013                    model2->setPerturbation(savePerturbation);
3014                    model2->createStatus();
3015                    model2->dual();
3016               } else if (maxIts && quadraticObj) {
3017                    // make sure no status left
3018                    model2->createStatus();
3019                    // solve
3020                    model2->setPerturbation(100);
3021                    model2->reducedGradient(1);
3022               }
3023          }
3024
3025          //model2->setMaximumIterations(saveMaxIts);
3026#ifdef BORROW
3027          model2->setNumberIterations(model2->numberIterations()+saveNumberIterations);
3028          delete [] rowPrimal;
3029          delete [] columnPrimal;
3030          delete [] rowDual;
3031          delete [] columnDual;
3032#endif
3033          if (extraPresolve) {
3034               pinfo2.postsolve(true);
3035               delete model2;
3036               model2 = saveModel2;
3037          }
3038          if (saveUpper) {
3039               if (!forceFixing) {
3040                    int numberRows = model2->numberRows();
3041                    int numberColumns = model2->numberColumns();
3042                    CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3043                    CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3044                    CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3045                    CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3046               }
3047               delete [] saveLower;
3048               delete [] saveUpper;
3049               saveLower = NULL;
3050               saveUpper = NULL;
3051               if (method != ClpSolve::useBarrierNoCross)
3052                    model2->primal(1);
3053          }
3054          model2->setPerturbation(savePerturbation);
3055          time2 = CoinCpuTime();
3056          timeCore = time2 - timeX;
3057          handler_->message(CLP_INTERVAL_TIMING, messages_)
3058                    << "Crossover" << timeCore << time2 - time1
3059                    << CoinMessageEol;
3060          timeX = time2;
3061#else
3062          abort();
3063#endif
3064     } else {
3065          assert (method != ClpSolve::automatic); // later
3066          time2 = 0.0;
3067     }
3068     if (saveMatrix) {
3069          if (model2 == this) {
3070               // delete and replace
3071               delete model2->clpMatrix();
3072               model2->replaceMatrix(saveMatrix);
3073          } else {
3074               delete saveMatrix;
3075          }
3076     }
3077     numberIterations = model2->numberIterations();
3078     finalStatus = model2->status();
3079     int finalSecondaryStatus = model2->secondaryStatus();
3080     if (presolve == ClpSolve::presolveOn) {
3081          int saveLevel = logLevel();
3082          if ((specialOptions_ & 1024) == 0)
3083               setLogLevel(CoinMin(1, saveLevel));
3084          else
3085               setLogLevel(CoinMin(0, saveLevel));
3086          pinfo->postsolve(true);
3087          numberIterations_ = 0;
3088          delete pinfo;
3089          pinfo = NULL;
3090          factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
3091          time2 = CoinCpuTime();
3092          timePresolve += time2 - timeX;
3093          handler_->message(CLP_INTERVAL_TIMING, messages_)
3094                    << "Postsolve" << time2 - timeX << time2 - time1
3095                    << CoinMessageEol;
3096          timeX = time2;
3097          if (!presolveToFile) {
3098#if 1 //ndef ABC_INHERIT
3099               delete model2;
3100#else
3101               if (model2->abcSimplex())
3102                 delete model2->abcSimplex();
3103               else
3104                 delete model2;
3105#endif
3106          }
3107          if (interrupt)
3108               currentModel = this;
3109          // checkSolution(); already done by postSolve
3110          setLogLevel(saveLevel);
3111          int oldStatus=problemStatus_;
3112          setProblemStatus(finalStatus);
3113          setSecondaryStatus(finalSecondaryStatus);
3114          int rcode=eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
3115          if (finalStatus != 3 && rcode < 0 && (finalStatus || oldStatus == -1)) {
3116               int savePerturbation = perturbation();
3117               if (!finalStatus || (moreSpecialOptions_ & 2) == 0 ||
3118                   fabs(sumDualInfeasibilities_)+
3119                   fabs(sumPrimalInfeasibilities_)<1.0e-3) {
3120                    if (finalStatus == 2) {
3121                         // unbounded - get feasible first
3122                         double save = optimizationDirection_;
3123                         optimizationDirection_ = 0.0;
3124                         primal(1);
3125                         optimizationDirection_ = save;
3126                         primal(1);
3127                    } else if (finalStatus == 1) {
3128                         dual();
3129                    } else {
3130                        if ((moreSpecialOptions_&65536)==0) {
3131                          if (numberRows_<10000) 
3132                            setPerturbation(100); // probably better to perturb after n its
3133                          else if (savePerturbation<100)
3134                            setPerturbation(51); // probably better to perturb after n its
3135                        }
3136#ifndef ABC_INHERIT
3137                        primal(1);
3138#else
3139                        dealWithAbc(1,2,interrupt);
3140#endif
3141                    }
3142               } else {
3143                    // just set status
3144                    problemStatus_ = finalStatus;
3145               }
3146               setPerturbation(savePerturbation);
3147               numberIterations += numberIterations_;
3148               numberIterations_ = numberIterations;
3149               finalStatus = status();
3150               time2 = CoinCpuTime();
3151               handler_->message(CLP_INTERVAL_TIMING, messages_)
3152                         << "Cleanup" << time2 - timeX << time2 - time1
3153                         << CoinMessageEol;
3154               timeX = time2;
3155          } else if (rcode >= 0) {
3156#ifdef ABC_INHERIT
3157            dealWithAbc(1,2,true);
3158#else
3159            primal(1);
3160#endif
3161          } else {
3162               secondaryStatus_ = finalSecondaryStatus;
3163          }
3164     } else if (model2 != this) {
3165          // not presolved - but different model used (sprint probably)
3166          CoinMemcpyN(model2->primalRowSolution(),
3167                      numberRows_, this->primalRowSolution());
3168          CoinMemcpyN(model2->dualRowSolution(),
3169                      numberRows_, this->dualRowSolution());
3170          CoinMemcpyN(model2->primalColumnSolution(),
3171                      numberColumns_, this->primalColumnSolution());
3172          CoinMemcpyN(model2->dualColumnSolution(),
3173                      numberColumns_, this->dualColumnSolution());
3174          CoinMemcpyN(model2->statusArray(),
3175                      numberColumns_ + numberRows_, this->statusArray());
3176          objectiveValue_ = model2->objectiveValue_;
3177          numberIterations_ = model2->numberIterations_;
3178          problemStatus_ = model2->problemStatus_;
3179          secondaryStatus_ = model2->secondaryStatus_;
3180          delete model2;
3181     }
3182     if (method != ClpSolve::useBarrierNoCross &&
3183               method != ClpSolve::useBarrier)
3184          setMaximumIterations(saveMaxIterations);
3185     std::string statusMessage[] = {"Unknown", "Optimal", "PrimalInfeasible", "DualInfeasible", "Stopped",
3186                                    "Errors", "User stopped"
3187                                   };
3188     assert (finalStatus >= -1 && finalStatus <= 5);
3189     numberIterations_ = numberIterations;
3190     handler_->message(CLP_TIMING, messages_)
3191               << statusMessage[finalStatus+1] << objectiveValue() << numberIterations << time2 - time1;
3192     handler_->printing(presolve == ClpSolve::presolveOn)
3193               << timePresolve;
3194     handler_->printing(timeIdiot != 0.0)
3195               << timeIdiot;
3196     handler_->message() << CoinMessageEol;
3197     if (interrupt)
3198          signal(SIGINT, saveSignal);
3199     perturbation_ = savePerturbation;
3200     scalingFlag_ = saveScaling;
3201     // If faking objective - put back correct one
3202     if (savedObjective) {
3203          delete objective_;
3204          objective_ = savedObjective;
3205     }
3206     if (options.getSpecialOption(1) == 2 &&
3207               options.getExtraInfo(1) > 1000000) {
3208          ClpObjective * savedObjective = objective_;
3209          // make up zero objective
3210          double * obj = new double[numberColumns_];
3211          for (int i = 0; i < numberColumns_; i++)
3212               obj[i] = 0.0;
3213          objective_ = new ClpLinearObjective(obj, numberColumns_);
3214          delete [] obj;
3215          primal(1);
3216          delete objective_;
3217          objective_ = savedObjective;
3218          finalStatus = status();
3219     }
3220     eventHandler()->event(ClpEventHandler::presolveEnd);
3221     delete pinfo;
3222     return finalStatus;
3223}
3224// General solve
3225int
3226ClpSimplex::initialSolve()
3227{
3228     // Default so use dual
3229     ClpSolve options;
3230     return initialSolve(options);
3231}
3232// General dual solve
3233int
3234ClpSimplex::initialDualSolve()
3235{
3236     ClpSolve options;
3237     // Use dual
3238     options.setSolveType(ClpSolve::useDual);
3239     return initialSolve(options);
3240}
3241// General primal solve
3242int
3243ClpSimplex::initialPrimalSolve()
3244{
3245     ClpSolve options;
3246     // Use primal
3247     options.setSolveType(ClpSolve::usePrimal);
3248     return initialSolve(options);
3249}
3250// barrier solve, not to be followed by crossover
3251int
3252ClpSimplex::initialBarrierNoCrossSolve()
3253{
3254     ClpSolve options;
3255     // Use primal
3256     options.setSolveType(ClpSolve::useBarrierNoCross);
3257     return initialSolve(options);
3258}
3259
3260// General barrier solve
3261int
3262ClpSimplex::initialBarrierSolve()
3263{
3264     ClpSolve options;
3265     // Use primal
3266     options.setSolveType(ClpSolve::useBarrier);
3267     return initialSolve(options);
3268}
3269
3270// Default constructor
3271ClpSolve::ClpSolve (  )
3272{
3273     method_ = automatic;
3274     presolveType_ = presolveOn;
3275     numberPasses_ = 5;
3276     int i;
3277     for (i = 0; i < 7; i++)
3278          options_[i] = 0;
3279     // say no +-1 matrix
3280     options_[3] = 1;
3281     for (i = 0; i < 7; i++)
3282          extraInfo_[i] = -1;
3283     independentOptions_[0] = 0;
3284     // But switch off slacks
3285     independentOptions_[1] = 512;
3286     // Substitute up to 3
3287     independentOptions_[2] = 3;
3288
3289}
3290// Constructor when you really know what you are doing
3291ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
3292                     int numberPasses, int options[6],
3293                     int extraInfo[6], int independentOptions[3])
3294{
3295     method_ = method;
3296     presolveType_ = presolveType;
3297     numberPasses_ = numberPasses;
3298     int i;
3299     for (i = 0; i < 6; i++)
3300          options_[i] = options[i];
3301     options_[6] = 0;
3302     for (i = 0; i < 6; i++)
3303          extraInfo_[i] = extraInfo[i];
3304     extraInfo_[6] = 0;
3305     for (i = 0; i < 3; i++)
3306          independentOptions_[i] = independentOptions[i];
3307}
3308
3309// Copy constructor.
3310ClpSolve::ClpSolve(const ClpSolve & rhs)
3311{
3312     method_ = rhs.method_;
3313     presolveType_ = rhs.presolveType_;
3314     numberPasses_ = rhs.numberPasses_;
3315     int i;
3316     for ( i = 0; i < 7; i++)
3317          options_[i] = rhs.options_[i];
3318     for ( i = 0; i < 7; i++)
3319          extraInfo_[i] = rhs.extraInfo_[i];
3320     for (i = 0; i < 3; i++)
3321          independentOptions_[i] = rhs.independentOptions_[i];
3322}
3323// Assignment operator. This copies the data
3324ClpSolve &
3325ClpSolve::operator=(const ClpSolve & rhs)
3326{
3327     if (this != &rhs) {
3328          method_ = rhs.method_;
3329          presolveType_ = rhs.presolveType_;
3330          numberPasses_ = rhs.numberPasses_;
3331          int i;
3332          for (i = 0; i < 7; i++)
3333               options_[i] = rhs.options_[i];
3334          for (i = 0; i < 7; i++)
3335               extraInfo_[i] = rhs.extraInfo_[i];
3336          for (i = 0; i < 3; i++)
3337               independentOptions_[i] = rhs.independentOptions_[i];
3338     }
3339     return *this;
3340
3341}
3342// Destructor
3343ClpSolve::~ClpSolve (  )
3344{
3345}
3346// See header file for details
3347void
3348ClpSolve::setSpecialOption(int which, int value, int extraInfo)
3349{
3350     options_[which] = value;
3351     extraInfo_[which] = extraInfo;
3352}
3353int
3354ClpSolve::getSpecialOption(int which) const
3355{
3356     return options_[which];
3357}
3358
3359// Solve types
3360void
3361ClpSolve::setSolveType(SolveType method, int /*extraInfo*/)
3362{
3363     method_ = method;
3364}
3365
3366ClpSolve::SolveType
3367ClpSolve::getSolveType()
3368{
3369     return method_;
3370}
3371
3372// Presolve types
3373void
3374ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
3375{
3376     presolveType_ = amount;
3377     numberPasses_ = extraInfo;
3378}
3379ClpSolve::PresolveType
3380ClpSolve::getPresolveType()
3381{
3382     return presolveType_;
3383}
3384// Extra info for idiot (or sprint)
3385int
3386ClpSolve::getExtraInfo(int which) const
3387{
3388     return extraInfo_[which];
3389}
3390int
3391ClpSolve::getPresolvePasses() const
3392{
3393     return numberPasses_;
3394}
3395/* Say to return at once if infeasible,
3396   default is to solve */
3397void
3398ClpSolve::setInfeasibleReturn(bool trueFalse)
3399{
3400     independentOptions_[0] = trueFalse ? 1 : 0;
3401}
3402#include <string>
3403// Generates code for above constructor
3404void
3405ClpSolve::generateCpp(FILE * fp)
3406{
3407     std::string solveType[] = {
3408          "ClpSolve::useDual",
3409          "ClpSolve::usePrimal",
3410          "ClpSolve::usePrimalorSprint",
3411          "ClpSolve::useBarrier",
3412          "ClpSolve::useBarrierNoCross",
3413          "ClpSolve::automatic",
3414          "ClpSolve::notImplemented"
3415     };
3416     std::string presolveType[] =  {
3417          "ClpSolve::presolveOn",
3418          "ClpSolve::presolveOff",
3419          "ClpSolve::presolveNumber",
3420          "ClpSolve::presolveNumberCost"
3421     };
3422     fprintf(fp, "3  ClpSolve::SolveType method = %s;\n", solveType[method_].c_str());
3423     fprintf(fp, "3  ClpSolve::PresolveType presolveType = %s;\n",
3424             presolveType[presolveType_].c_str());
3425     fprintf(fp, "3  int numberPasses = %d;\n", numberPasses_);
3426     fprintf(fp, "3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
3427             options_[0], options_[1], options_[2],
3428             options_[3], options_[4], options_[5]);
3429     fprintf(fp, "3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
3430             extraInfo_[0], extraInfo_[1], extraInfo_[2],
3431             extraInfo_[3], extraInfo_[4], extraInfo_[5]);
3432     fprintf(fp, "3  int independentOptions[] = {%d,%d,%d};\n",
3433             independentOptions_[0], independentOptions_[1], independentOptions_[2]);
3434     fprintf(fp, "3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
3435     fprintf(fp, "3                    options,extraInfo,independentOptions);\n");
3436}
3437//#############################################################################
3438#include "ClpNonLinearCost.hpp"
3439
3440ClpSimplexProgress::ClpSimplexProgress ()
3441{
3442     int i;
3443     for (i = 0; i < CLP_PROGRESS; i++) {
3444          objective_[i] = COIN_DBL_MAX;
3445          infeasibility_[i] = -1.0; // set to an impossible value
3446          realInfeasibility_[i] = COIN_DBL_MAX;
3447          numberInfeasibilities_[i] = -1;
3448          iterationNumber_[i] = -1;
3449     }
3450#ifdef CLP_PROGRESS_WEIGHT
3451     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3452          objectiveWeight_[i] = COIN_DBL_MAX;
3453          infeasibilityWeight_[i] = -1.0; // set to an impossible value
3454          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3455          numberInfeasibilitiesWeight_[i] = -1;
3456          iterationNumberWeight_[i] = -1;
3457     }
3458     drop_ = 0.0;
3459     best_ = 0.0;
3460#endif
3461     initialWeight_ = 0.0;
3462     for (i = 0; i < CLP_CYCLE; i++) {
3463          //obj_[i]=COIN_DBL_MAX;
3464          in_[i] = -1;
3465          out_[i] = -1;
3466          way_[i] = 0;
3467     }
3468     numberTimes_ = 0;
3469     numberBadTimes_ = 0;
3470     numberReallyBadTimes_ = 0;
3471     numberTimesFlagged_ = 0;
3472     model_ = NULL;
3473     oddState_ = 0;
3474}
3475
3476
3477//-----------------------------------------------------------------------------
3478
3479ClpSimplexProgress::~ClpSimplexProgress ()
3480{
3481}
3482// Copy constructor.
3483ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
3484{
3485     int i;
3486     for (i = 0; i < CLP_PROGRESS; i++) {
3487          objective_[i] = rhs.objective_[i];
3488          infeasibility_[i] = rhs.infeasibility_[i];
3489          realInfeasibility_[i] = rhs.realInfeasibility_[i];
3490          numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
3491          iterationNumber_[i] = rhs.iterationNumber_[i];
3492     }
3493#ifdef CLP_PROGRESS_WEIGHT
3494     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3495          objectiveWeight_[i] = rhs.objectiveWeight_[i];
3496          infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
3497          realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
3498          numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
3499          iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
3500     }
3501     drop_ = rhs.drop_;
3502     best_ = rhs.best_;
3503#endif
3504     initialWeight_ = rhs.initialWeight_;
3505     for (i = 0; i < CLP_CYCLE; i++) {
3506          //obj_[i]=rhs.obj_[i];
3507          in_[i] = rhs.in_[i];
3508          out_[i] = rhs.out_[i];
3509          way_[i] = rhs.way_[i];
3510     }
3511     numberTimes_ = rhs.numberTimes_;
3512     numberBadTimes_ = rhs.numberBadTimes_;
3513     numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
3514     numberTimesFlagged_ = rhs.numberTimesFlagged_;
3515     model_ = rhs.model_;
3516     oddState_ = rhs.oddState_;
3517}
3518// Copy constructor.from model
3519ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model)
3520{
3521     model_ = model;
3522     reset();
3523     initialWeight_ = 0.0;
3524}
3525// Fill from model
3526void
3527ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
3528{
3529     model_ = model;
3530     reset();
3531     initialWeight_ = 0.0;
3532}
3533// Assignment operator. This copies the data
3534ClpSimplexProgress &
3535ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
3536{
3537     if (this != &rhs) {
3538          int i;
3539          for (i = 0; i < CLP_PROGRESS; i++) {
3540               objective_[i] = rhs.objective_[i];
3541               infeasibility_[i] = rhs.infeasibility_[i];
3542               realInfeasibility_[i] = rhs.realInfeasibility_[i];
3543               numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
3544               iterationNumber_[i] = rhs.iterationNumber_[i];
3545          }
3546#ifdef CLP_PROGRESS_WEIGHT
3547          for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3548               objectiveWeight_[i] = rhs.objectiveWeight_[i];
3549               infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
3550               realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
3551               numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
3552               iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
3553          }
3554          drop_ = rhs.drop_;
3555          best_ = rhs.best_;
3556#endif
3557          initialWeight_ = rhs.initialWeight_;
3558          for (i = 0; i < CLP_CYCLE; i++) {
3559               //obj_[i]=rhs.obj_[i];
3560               in_[i] = rhs.in_[i];
3561               out_[i] = rhs.out_[i];
3562               way_[i] = rhs.way_[i];
3563          }
3564          numberTimes_ = rhs.numberTimes_;
3565          numberBadTimes_ = rhs.numberBadTimes_;
3566          numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
3567          numberTimesFlagged_ = rhs.numberTimesFlagged_;
3568          model_ = rhs.model_;
3569          oddState_ = rhs.oddState_;
3570     }
3571     return *this;
3572}
3573// Seems to be something odd about exact comparison of doubles on linux
3574static bool equalDouble(double value1, double value2)
3575{
3576
3577     union {
3578          double d;
3579          int i[2];
3580     } v1, v2;
3581     v1.d = value1;
3582     v2.d = value2;
3583     if (sizeof(int) * 2 == sizeof(double))
3584          return (v1.i[0] == v2.i[0] && v1.i[1] == v2.i[1]);
3585     else
3586          return (v1.i[0] == v2.i[0]);
3587}
3588int
3589ClpSimplexProgress::looping()
3590{
3591     if (!model_)
3592          return -1;
3593     double objective;
3594     if (model_->algorithm() < 0) {
3595       objective = model_->rawObjectiveValue();
3596       objective -= model_->bestPossibleImprovement();
3597     } else {
3598       objective = model_->nonLinearCost()->feasibleReportCost();
3599     }
3600     double infeasibility;
3601     double realInfeasibility = 0.0;
3602     int numberInfeasibilities;
3603     int iterationNumber = model_->numberIterations();
3604     //numberTimesFlagged_ = 0;
3605     if (model_->algorithm() < 0) {
3606          // dual
3607          infeasibility = model_->sumPrimalInfeasibilities();
3608          numberInfeasibilities = model_->numberPrimalInfeasibilities();
3609     } else {
3610          //primal
3611          infeasibility = model_->sumDualInfeasibilities();
3612          realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
3613          numberInfeasibilities = model_->numberDualInfeasibilities();
3614     }
3615     int i;
3616     int numberMatched = 0;
3617     int matched = 0;
3618     int nsame = 0;
3619     for (i = 0; i < CLP_PROGRESS; i++) {
3620          bool matchedOnObjective = equalDouble(objective, objective_[i]);
3621          bool matchedOnInfeasibility = equalDouble(infeasibility, infeasibility_[i]);
3622          bool matchedOnInfeasibilities =
3623               (numberInfeasibilities == numberInfeasibilities_[i]);
3624
3625          if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
3626               matched |= (1 << i);
3627               // Check not same iteration
3628               if (iterationNumber != iterationNumber_[i]) {
3629                    numberMatched++;
3630                    // here mainly to get over compiler bug?
3631                    if (model_->messageHandler()->logLevel() > 10)
3632                         printf("%d %d %d %d %d loop check\n", i, numberMatched,
3633                                matchedOnObjective, matchedOnInfeasibility,
3634                                matchedOnInfeasibilities);
3635               } else {
3636                    // stuck but code should notice
3637                    nsame++;
3638               }
3639          }
3640          if (i) {
3641               objective_[i-1] = objective_[i];
3642               infeasibility_[i-1] = infeasibility_[i];
3643               realInfeasibility_[i-1] = realInfeasibility_[i];
3644               numberInfeasibilities_[i-1] = numberInfeasibilities_[i];
3645               iterationNumber_[i-1] = iterationNumber_[i];
3646          }
3647     }
3648     objective_[CLP_PROGRESS-1] = objective;
3649     infeasibility_[CLP_PROGRESS-1] = infeasibility;
3650     realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
3651     numberInfeasibilities_[CLP_PROGRESS-1] = numberInfeasibilities;
3652     iterationNumber_[CLP_PROGRESS-1] = iterationNumber;
3653     if (nsame == CLP_PROGRESS)
3654          numberMatched = CLP_PROGRESS; // really stuck
3655     if (model_->progressFlag())
3656          numberMatched = 0;
3657     numberTimes_++;
3658     if (numberTimes_ < 10)
3659          numberMatched = 0;
3660     // skip if just last time as may be checking something
3661     if (matched == (1 << (CLP_PROGRESS - 1)))
3662          numberMatched = 0;
3663     if (numberMatched && model_->clpMatrix()->type() < 15) {
3664          model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
3665                    << numberMatched
3666                    << matched
3667                    << numberTimes_
3668                    << CoinMessageEol;
3669          numberBadTimes_++;
3670          if (numberBadTimes_ < 10) {
3671               // make factorize every iteration
3672               model_->forceFactorization(1);
3673               if (numberBadTimes_ < 2) {
3674                    startCheck(); // clear other loop check
3675                    if (model_->algorithm() < 0) {
3676                         // dual - change tolerance
3677                         model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
3678                         // if infeasible increase dual bound
3679                         if (model_->dualBound() < 1.0e17) {
3680                              model_->setDualBound(model_->dualBound() * 1.1);
3681                              static_cast<ClpSimplexDual *>(model_)->resetFakeBounds(0);
3682                         }
3683                    } else {
3684                         // primal - change tolerance
3685                         if (numberBadTimes_ > 3)
3686                              model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
3687                         // if infeasible increase infeasibility cost
3688                         if (model_->nonLinearCost()->numberInfeasibilities() &&
3689                                   model_->infeasibilityCost() < 1.0e17) {
3690                              model_->setInfeasibilityCost(model_->infeasibilityCost() * 1.1);
3691                         }
3692                    }
3693               } else {
3694                    // flag
3695                    int iSequence;
3696                    if (model_->algorithm() < 0) {
3697                         // dual
3698                         if (model_->dualBound() > 1.0e14)
3699                              model_->setDualBound(1.0e14);
3700                         iSequence = in_[CLP_CYCLE-1];
3701                    } else {
3702                         // primal
3703                         if (model_->infeasibilityCost() > 1.0e14)
3704                              model_->setInfeasibilityCost(1.0e14);
3705                         iSequence = out_[CLP_CYCLE-1];
3706                    }
3707                    if (iSequence >= 0) {
3708                         char x = model_->isColumn(iSequence) ? 'C' : 'R';
3709                         if (model_->messageHandler()->logLevel() >= 63)
3710                              model_->messageHandler()->message(CLP_SIMPLEX_FLAG, model_->messages())
3711                                        << x << model_->sequenceWithin(iSequence)
3712                                        << CoinMessageEol;
3713                         // if Gub then needs to be sequenceIn_
3714                         int save = model_->sequenceIn();
3715                         model_->setSequenceIn(iSequence);
3716                         model_->setFlagged(iSequence);
3717                         model_->setSequenceIn(save);
3718                         //printf("flagging %d from loop\n",iSequence);
3719                         startCheck();
3720                    } else {
3721                         // Give up
3722                         if (model_->messageHandler()->logLevel() >= 63)
3723                              printf("***** All flagged?\n");
3724                         return 4;
3725                    }
3726                    // reset
3727                    numberBadTimes_ = 2;
3728               }
3729               return -2;
3730          } else {
3731               // look at solution and maybe declare victory
3732               if (infeasibility < 1.0e-4) {
3733                    return 0;
3734               } else {
3735                    model_->messageHandler()->message(CLP_LOOP, model_->messages())
3736                              << CoinMessageEol;
3737#ifndef NDEBUG
3738                    printf("debug loop ClpSimplex A\n");
3739                    abort();
3740#endif
3741                    return 3;
3742               }
3743          }
3744     }
3745     return -1;
3746}
3747// Resets as much as possible
3748void
3749ClpSimplexProgress::reset()
3750{
3751     int i;
3752     for (i = 0; i < CLP_PROGRESS; i++) {
3753          if (model_->algorithm() >= 0)
3754               objective_[i] = COIN_DBL_MAX;
3755          else
3756               objective_[i] = -COIN_DBL_MAX;
3757          infeasibility_[i] = -1.0; // set to an impossible value
3758          realInfeasibility_[i] = COIN_DBL_MAX;
3759          numberInfeasibilities_[i] = -1;
3760          iterationNumber_[i] = -1;
3761     }
3762#ifdef CLP_PROGRESS_WEIGHT
3763     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3764          objectiveWeight_[i] = COIN_DBL_MAX;
3765          infeasibilityWeight_[i] = -1.0; // set to an impossible value
3766          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3767          numberInfeasibilitiesWeight_[i] = -1;
3768          iterationNumberWeight_[i] = -1;
3769     }
3770     drop_ = 0.0;
3771     best_ = 0.0;
3772#endif
3773     for (i = 0; i < CLP_CYCLE; i++) {
3774          //obj_[i]=COIN_DBL_MAX;
3775          in_[i] = -1;
3776          out_[i] = -1;
3777          way_[i] = 0;
3778     }
3779     numberTimes_ = 0;
3780     numberBadTimes_ = 0;
3781     numberReallyBadTimes_ = 0;
3782     numberTimesFlagged_ = 0;
3783     oddState_ = 0;
3784}
3785// Returns previous objective (if -1) - current if (0)
3786double
3787ClpSimplexProgress::lastObjective(int back) const
3788{
3789     return objective_[CLP_PROGRESS-1-back];
3790}
3791// Returns previous infeasibility (if -1) - current if (0)
3792double
3793ClpSimplexProgress::lastInfeasibility(int back) const
3794{
3795     return realInfeasibility_[CLP_PROGRESS-1-back];
3796}
3797// Sets real primal infeasibility
3798void
3799ClpSimplexProgress::setInfeasibility(double value)
3800{
3801     for (int i = 1; i < CLP_PROGRESS; i++)
3802          realInfeasibility_[i-1] = realInfeasibility_[i];
3803     realInfeasibility_[CLP_PROGRESS-1] = value;
3804}
3805// Modify objective e.g. if dual infeasible in dual
3806void
3807ClpSimplexProgress::modifyObjective(double value)
3808{
3809     objective_[CLP_PROGRESS-1] = value;
3810}
3811// Returns previous iteration number (if -1) - current if (0)
3812int
3813ClpSimplexProgress::lastIterationNumber(int back) const
3814{
3815     return iterationNumber_[CLP_PROGRESS-1-back];
3816}
3817// clears iteration numbers (to switch off panic)
3818void
3819ClpSimplexProgress::clearIterationNumbers()
3820{
3821     for (int i = 0; i < CLP_PROGRESS; i++)
3822          iterationNumber_[i] = -1;
3823}
3824// Start check at beginning of whileIterating
3825void
3826ClpSimplexProgress::startCheck()
3827{
3828     int i;
3829     for (i = 0; i < CLP_CYCLE; i++) {
3830          //obj_[i]=COIN_DBL_MAX;
3831          in_[i] = -1;
3832          out_[i] = -1;
3833          way_[i] = 0;
3834     }
3835}
3836// Returns cycle length in whileIterating
3837int
3838ClpSimplexProgress::cycle(int in, int out, int wayIn, int wayOut)
3839{
3840     int i;
3841#if 0
3842     if (model_->numberIterations() > 206571) {
3843          printf("in %d out %d\n", in, out);
3844          for (i = 0; i < CLP_CYCLE; i++)
3845               printf("cy %d in %d out %d\n", i, in_[i], out_[i]);
3846     }
3847#endif
3848     int matched = 0;
3849     // first see if in matches any out
3850     for (i = 1; i < CLP_CYCLE; i++) {
3851          if (in == out_[i]) {
3852               // even if flip then suspicious
3853               matched = -1;
3854               break;
3855          }
3856     }
3857#if 0
3858     if (!matched || in_[0] < 0) {
3859          // can't be cycle
3860          for (i = 0; i < CLP_CYCLE - 1; i++) {
3861               //obj_[i]=obj_[i+1];
3862               in_[i] = in_[i+1];
3863               out_[i] = out_[i+1];
3864               way_[i] = way_[i+1];
3865          }
3866     } else {
3867          // possible cycle
3868          matched = 0;
3869          for (i = 0; i < CLP_CYCLE - 1; i++) {
3870               int k;
3871               char wayThis = way_[i];
3872               int inThis = in_[i];
3873               int outThis = out_[i];
3874               //double objThis = obj_[i];
3875               for(k = i + 1; k < CLP_CYCLE; k++) {
3876                    if (inThis == in_[k] && outThis == out_[k] && wayThis == way_[k]) {
3877                         int distance = k - i;
3878                         if (k + distance < CLP_CYCLE) {
3879                              // See if repeats
3880                              int j = k + distance;
3881                              if (inThis == in_[j] && outThis == out_[j] && wayThis == way_[j]) {
3882                                   matched = distance;
3883                                   break;
3884                              }
3885                         } else {
3886                              matched = distance;
3887                              break;
3888                         }
3889                    }
3890               }
3891               //obj_[i]=obj_[i+1];
3892               in_[i] = in_[i+1];
3893               out_[i] = out_[i+1];
3894               way_[i] = way_[i+1];
3895          }
3896     }
3897#else
3898     if (matched && in_[0] >= 0) {
3899          // possible cycle - only check [0] against all
3900          matched = 0;
3901          int nMatched = 0;
3902          char way0 = way_[0];
3903          int in0 = in_[0];
3904          int out0 = out_[0];
3905          //double obj0 = obj_[i];
3906          for(int k = 1; k < CLP_CYCLE - 4; k++) {
3907               if (in0 == in_[k] && out0 == out_[k] && way0 == way_[k]) {
3908                    nMatched++;
3909                    // See if repeats
3910                    int end = CLP_CYCLE - k;
3911                    int j;
3912                    for ( j = 1; j < end; j++) {
3913                         if (in_[j+k] != in_[j] || out_[j+k] != out_[j] || way_[j+k] != way_[j])
3914                              break;
3915                    }
3916                    if (j == end) {
3917                         matched = k;
3918                         break;
3919                    }
3920               }
3921          }
3922          // If three times then that is too much even if not regular
3923          if (matched <= 0 && nMatched > 1)
3924               matched = 100;
3925     }
3926     for (i = 0; i < CLP_CYCLE - 1; i++) {
3927          //obj_[i]=obj_[i+1];
3928          in_[i] = in_[i+1];
3929          out_[i] = out_[i+1];
3930          way_[i] = way_[i+1];
3931     }
3932#endif
3933     int way = 1 - wayIn + 4 * (1 - wayOut);
3934     //obj_[i]=model_->objectiveValue();
3935     in_[CLP_CYCLE-1] = in;
3936     out_[CLP_CYCLE-1] = out;
3937     way_[CLP_CYCLE-1] = static_cast<char>(way);
3938     return matched;
3939}
3940#include "CoinStructuredModel.hpp"
3941// Solve using structure of model and maybe in parallel
3942int
3943ClpSimplex::solve(CoinStructuredModel * model)
3944{
3945     // analyze structure
3946     int numberRowBlocks = model->numberRowBlocks();
3947     int numberColumnBlocks = model->numberColumnBlocks();
3948     int numberElementBlocks = model->numberElementBlocks();
3949     if (numberElementBlocks == 1) {
3950          loadProblem(*model, false);
3951          return dual();
3952     }
3953     // For now just get top level structure
3954     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3955     for (int i = 0; i < numberElementBlocks; i++) {
3956          CoinStructuredModel * subModel =
3957               dynamic_cast<CoinStructuredModel *>(model->block(i));
3958          CoinModel * thisBlock;
3959          if (subModel) {
3960               thisBlock = subModel->coinModelBlock(blockInfo[i]);
3961               model->setCoinModel(thisBlock, i);
3962          } else {
3963               thisBlock = dynamic_cast<CoinModel *>(model->block(i));
3964               assert (thisBlock);
3965               // just fill in info
3966               CoinModelBlockInfo info = CoinModelBlockInfo();
3967               int whatsSet = thisBlock->whatIsSet();
3968               info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
3969               info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
3970               info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
3971               info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
3972               info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
3973               info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
3974               // Which block
3975               int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
3976               info.rowBlock = iRowBlock;
3977               int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
3978               info.columnBlock = iColumnBlock;
3979               blockInfo[i] = info;
3980          }
3981     }
3982     int * rowCounts = new int [numberRowBlocks];
3983     CoinZeroN(rowCounts, numberRowBlocks);
3984     int * columnCounts = new int [numberColumnBlocks+1];
3985     CoinZeroN(columnCounts, numberColumnBlocks);
3986     int decomposeType = 0;
3987     for (int i = 0; i < numberElementBlocks; i++) {
3988          int iRowBlock = blockInfo[i].rowBlock;
3989          int iColumnBlock = blockInfo[i].columnBlock;
3990          rowCounts[iRowBlock]++;
3991          columnCounts[iColumnBlock]++;
3992     }
3993     if (numberRowBlocks == numberColumnBlocks ||
3994               numberRowBlocks == numberColumnBlocks + 1) {
3995          // could be Dantzig-Wolfe
3996          int numberG1 = 0;
3997          for (int i = 0; i < numberRowBlocks; i++) {
3998               if (rowCounts[i] > 1)
3999                    numberG1++;
4000          }
4001          bool masterColumns = (numberColumnBlocks == numberRowBlocks);
4002          if ((masterColumns && numberElementBlocks == 2 * numberRowBlocks - 1)
4003                    || (!masterColumns && numberElementBlocks == 2 * numberRowBlocks)) {
4004               if (numberG1 < 2)
4005                    decomposeType = 1;
4006          }
4007     }
4008     if (!decomposeType && (numberRowBlocks == numberColumnBlocks ||
4009                            numberRowBlocks == numberColumnBlocks - 1)) {
4010          // could be Benders
4011          int numberG1 = 0;
4012          for (int i = 0; i < numberColumnBlocks; i++) {
4013               if (columnCounts[i] > 1)
4014                    numberG1++;
4015          }
4016          bool masterRows = (numberColumnBlocks == numberRowBlocks);
4017          if ((masterRows && numberElementBlocks == 2 * numberColumnBlocks - 1)
4018                    || (!masterRows && numberElementBlocks == 2 * numberColumnBlocks)) {
4019               if (numberG1 < 2)
4020                    decomposeType = 2;
4021          }
4022     }
4023     delete [] rowCounts;
4024     delete [] columnCounts;
4025     delete [] blockInfo;
4026     // decide what to do
4027     switch (decomposeType) {
4028          // No good
4029     case 0:
4030          loadProblem(*model, false);
4031          return dual();
4032          // DW
4033     case 1:
4034          return solveDW(model);
4035          // Benders
4036     case 2:
4037          return solveBenders(model);
4038     }
4039     return 0; // to stop compiler warning
4040}
4041/* This loads a model from a CoinStructuredModel object - returns number of errors.
4042   If originalOrder then keep to order stored in blocks,
4043   otherwise first column/rows correspond to first block - etc.
4044   If keepSolution true and size is same as current then
4045   keeps current status and solution
4046*/
4047int
4048ClpSimplex::loadProblem (  CoinStructuredModel & coinModel,
4049                           bool originalOrder,
4050                           bool keepSolution)
4051{
4052     unsigned char * status = NULL;
4053     double * psol = NULL;
4054     double * dsol = NULL;
4055     int numberRows = coinModel.numberRows();
4056     int numberColumns = coinModel.numberColumns();
4057     int numberRowBlocks = coinModel.numberRowBlocks();
4058     int numberColumnBlocks = coinModel.numberColumnBlocks();
4059     int numberElementBlocks = coinModel.numberElementBlocks();
4060     if (status_ && numberRows_ && numberRows_ == numberRows &&
4061               numberColumns_ == numberColumns && keepSolution) {
4062          status = new unsigned char [numberRows_+numberColumns_];
4063          CoinMemcpyN(status_, numberRows_ + numberColumns_, status);
4064          psol = new double [numberRows_+numberColumns_];
4065          CoinMemcpyN(columnActivity_, numberColumns_, psol);
4066          CoinMemcpyN(rowActivity_, numberRows_, psol + numberColumns_);
4067          dsol = new double [numberRows_+numberColumns_];
4068          CoinMemcpyN(reducedCost_, numberColumns_, dsol);
4069          CoinMemcpyN(dual_, numberRows_, dsol + numberColumns_);
4070     }
4071     int returnCode = 0;
4072     double * rowLower = new double [numberRows];
4073     double * rowUpper = new double [numberRows];
4074     double * columnLower = new double [numberColumns];
4075     double * columnUpper = new double [numberColumns];
4076     double * objective = new double [numberColumns];
4077     int * integerType = new int [numberColumns];
4078     CoinBigIndex numberElements = 0;
4079     // Bases for blocks
4080     int * rowBase = new int[numberRowBlocks];
4081     CoinFillN(rowBase, numberRowBlocks, -1);
4082     // And row to put it
4083     int * whichRow = new int [numberRows];
4084     int * columnBase = new int[numberColumnBlocks];
4085     CoinFillN(columnBase, numberColumnBlocks, -1);
4086     // And column to put it
4087     int * whichColumn = new int [numberColumns];
4088     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4089          CoinModel * block = coinModel.coinBlock(iBlock);
4090          numberElements += block->numberElements();
4091          //and set up elements etc
4092          double * associated = block->associatedArray();
4093          // If strings then do copies
4094          if (block->stringsExist())
4095               returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
4096                                                 objective, integerType, associated);
4097          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4098          int iRowBlock = info.rowBlock;
4099          int iColumnBlock = info.columnBlock;
4100          if (rowBase[iRowBlock] < 0) {
4101               rowBase[iRowBlock] = block->numberRows();
4102               // Save block number
4103               whichRow[numberRows-numberRowBlocks+iRowBlock] = iBlock;
4104          } else {
4105               assert(rowBase[iRowBlock] == block->numberRows());
4106          }
4107          if (columnBase[iColumnBlock] < 0) {
4108               columnBase[iColumnBlock] = block->numberColumns();
4109               // Save block number
4110               whichColumn[numberColumns-numberColumnBlocks+iColumnBlock] = iBlock;
4111          } else {
4112               assert(columnBase[iColumnBlock] == block->numberColumns());
4113          }
4114     }
4115     // Fill arrays with defaults
4116     CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
4117     CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
4118     CoinFillN(columnLower, numberColumns, 0.0);
4119     CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
4120     CoinFillN(objective, numberColumns, 0.0);
4121     CoinFillN(integerType, numberColumns, 0);
4122     int n = 0;
4123     for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4124          int k = rowBase[iBlock];
4125          rowBase[iBlock] = n;
4126          assert (k >= 0);
4127          // block number
4128          int jBlock = whichRow[numberRows-numberRowBlocks+iBlock];
4129          if (originalOrder) {
4130               memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
4131          } else {
4132               CoinIotaN(whichRow + n, k, n);
4133          }
4134          n += k;
4135     }
4136     assert (n == numberRows);
4137     n = 0;
4138     for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4139          int k = columnBase[iBlock];
4140          columnBase[iBlock] = n;
4141          assert (k >= 0);
4142          if (k) {
4143               // block number
4144               int jBlock = whichColumn[numberColumns-numberColumnBlocks+iBlock];
4145               if (originalOrder) {
4146                    memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
4147                           k * sizeof(int));
4148               } else {
4149                    CoinIotaN(whichColumn + n, k, n);
4150               }
4151               n += k;
4152          }
4153     }
4154     assert (n == numberColumns);
4155     bool gotIntegers = false;
4156     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4157          CoinModel * block = coinModel.coinBlock(iBlock);
4158          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4159          int iRowBlock = info.rowBlock;
4160          int iRowBase = rowBase[iRowBlock];
4161          int iColumnBlock = info.columnBlock;
4162          int iColumnBase = columnBase[iColumnBlock];
4163          if (info.rhs) {
4164               int nRows = block->numberRows();
4165               const double * lower = block->rowLowerArray();
4166               const double * upper = block->rowUpperArray();
4167               for (int i = 0; i < nRows; i++) {
4168                    int put = whichRow[i+iRowBase];
4169                    rowLower[put] = lower[i];
4170                    rowUpper[put] = upper[i];
4171               }
4172          }
4173          if (info.bounds) {
4174               int nColumns = block->numberColumns();
4175               const double * lower = block->columnLowerArray();
4176               const double * upper = block->columnUpperArray();
4177               const double * obj = block->objectiveArray();
4178               for (int i = 0; i < nColumns; i++) {
4179                    int put = whichColumn[i+iColumnBase];
4180                    columnLower[put] = lower[i];
4181                    columnUpper[put] = upper[i];
4182                    objective[put] = obj[i];
4183               }
4184          }
4185          if (info.integer) {
4186               gotIntegers = true;
4187               int nColumns = block->numberColumns();
4188               const int * type = block->integerTypeArray();
4189               for (int i = 0; i < nColumns; i++) {
4190                    int put = whichColumn[i+iColumnBase];
4191                    integerType[put] = type[i];
4192               }
4193          }
4194     }
4195     gutsOfLoadModel(numberRows, numberColumns,
4196                     columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
4197     delete [] rowLower;
4198     delete [] rowUpper;
4199     delete [] columnLower;
4200     delete [] columnUpper;
4201     delete [] objective;
4202     // Do integers if wanted
4203     if (gotIntegers) {
4204          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4205               if (integerType[iColumn])
4206                    setInteger(iColumn);
4207          }
4208     }
4209     delete [] integerType;
4210     setObjectiveOffset(coinModel.objectiveOffset());
4211     // Space for elements
4212     int * row = new int[numberElements];
4213     int * column = new int[numberElements];
4214     double * element = new double[numberElements];
4215     numberElements = 0;
4216     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4217          CoinModel * block = coinModel.coinBlock(iBlock);
4218          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
4219          int iRowBlock = info.rowBlock;
4220          int iRowBase = rowBase[iRowBlock];
4221          int iColumnBlock = info.columnBlock;
4222          int iColumnBase = columnBase[iColumnBlock];
4223          if (info.rowName) {
4224               int numberItems = block->rowNames()->numberItems();
4225               assert( block->numberRows() >= numberItems);
4226               if (numberItems) {
4227                    const char *const * rowNames = block->rowNames()->names();
4228                    for (int i = 0; i < numberItems; i++) {
4229                         int put = whichRow[i+iRowBase];
4230                         std::string name = rowNames[i];
4231                         setRowName(put, name);
4232                    }
4233               }
4234          }
4235          if (info.columnName) {
4236               int numberItems = block->columnNames()->numberItems();
4237               assert( block->numberColumns() >= numberItems);
4238               if (numberItems) {
4239                    const char *const * columnNames = block->columnNames()->names();
4240                    for (int i = 0; i < numberItems; i++) {
4241                         int put = whichColumn[i+iColumnBase];
4242                         std::string name = columnNames[i];
4243                         setColumnName(put, name);
4244                    }
4245               }
4246          }
4247          if (info.matrix) {
4248               CoinPackedMatrix matrix2;
4249               const CoinPackedMatrix * matrix = block->packedMatrix();
4250               if (!matrix) {
4251                    double * associated = block->associatedArray();
4252                    block->createPackedMatrix(matrix2, associated);
4253                    matrix = &matrix2;
4254               }
4255               // get matrix data pointers
4256               const int * row2 = matrix->getIndices();
4257               const CoinBigIndex * columnStart = matrix->getVectorStarts();
4258               const double * elementByColumn = matrix->getElements();
4259               const int * columnLength = matrix->getVectorLengths();
4260               int n = matrix->getNumCols();
4261               assert (matrix->isColOrdered());
4262               for (int iColumn = 0; iColumn < n; iColumn++) {
4263                    CoinBigIndex j;
4264                    int jColumn = whichColumn[iColumn+iColumnBase];
4265                    for (j = columnStart[iColumn];
4266                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4267                         row[numberElements] = whichRow[row2[j] + iRowBase];
4268                         column[numberElements] = jColumn;
4269                         element[numberElements++] = elementByColumn[j];
4270                    }
4271               }
4272          }
4273     }
4274     delete [] whichRow;
4275     delete [] whichColumn;
4276     delete [] rowBase;
4277     delete [] columnBase;
4278     CoinPackedMatrix * matrix =
4279          new CoinPackedMatrix (true, row, column, element, numberElements);
4280     matrix_ = new ClpPackedMatrix(matrix);
4281     matrix_->setDimensions(numberRows, numberColumns);
4282     delete [] row;
4283     delete [] column;
4284     delete [] element;
4285     createStatus();
4286     if (status) {
4287          // copy back
4288          CoinMemcpyN(status, numberRows_ + numberColumns_, status_);
4289          CoinMemcpyN(psol, numberColumns_, columnActivity_);
4290          CoinMemcpyN(psol + numberColumns_, numberRows_, rowActivity_);
4291          CoinMemcpyN(dsol, numberColumns_, reducedCost_);
4292          CoinMemcpyN(dsol + numberColumns_, numberRows_, dual_);
4293          delete [] status;
4294          delete [] psol;
4295          delete [] dsol;
4296     }
4297     optimizationDirection_ = coinModel.optimizationDirection();
4298     return returnCode;
4299}
4300/*  If input negative scales objective so maximum <= -value
4301    and returns scale factor used.  If positive unscales and also
4302    redoes dual stuff
4303*/
4304double
4305ClpSimplex::scaleObjective(double value)
4306{
4307     double * obj = objective();
4308     double largest = 0.0;
4309     if (value < 0.0) {
4310          value = - value;
4311          for (int i = 0; i < numberColumns_; i++) {
4312               largest = CoinMax(largest, fabs(obj[i]));
4313          }
4314          if (largest > value) {
4315               double scaleFactor = value / largest;
4316               for (int i = 0; i < numberColumns_; i++) {
4317                    obj[i] *= scaleFactor;
4318                    reducedCost_[i] *= scaleFactor;
4319               }
4320               for (int i = 0; i < numberRows_; i++) {
4321                    dual_[i] *= scaleFactor;
4322               }
4323               largest /= value;
4324          } else {
4325               // no need
4326               largest = 1.0;
4327          }
4328     } else {
4329          // at end
4330          if (value != 1.0) {
4331               for (int i = 0; i < numberColumns_; i++) {
4332                    obj[i] *= value;
4333                    reducedCost_[i] *= value;
4334               }
4335               for (int i = 0; i < numberRows_; i++) {
4336                    dual_[i] *= value;
4337               }
4338               computeObjectiveValue();
4339          }
4340     }
4341     return largest;
4342}
4343// Solve using Dantzig-Wolfe decomposition and maybe in parallel
4344int
4345ClpSimplex::solveDW(CoinStructuredModel * model)
4346{
4347     double time1 = CoinCpuTime();
4348     int numberColumns = model->numberColumns();
4349     int numberRowBlocks = model->numberRowBlocks();
4350     int numberColumnBlocks = model->numberColumnBlocks();
4351     int numberElementBlocks = model->numberElementBlocks();
4352     // We already have top level structure
4353     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4354     for (int i = 0; i < numberElementBlocks; i++) {
4355          CoinModel * thisBlock = model->coinBlock(i);
4356          assert (thisBlock);
4357          // just fill in info
4358          CoinModelBlockInfo info = CoinModelBlockInfo();
4359          int whatsSet = thisBlock->whatIsSet();
4360          info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
4361          info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
4362          info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
4363          info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
4364          info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
4365          info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
4366          // Which block
4367          int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4368          info.rowBlock = iRowBlock;
4369          int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4370          info.columnBlock = iColumnBlock;
4371          blockInfo[i] = info;
4372     }
4373     // make up problems
4374     int numberBlocks = numberRowBlocks - 1;
4375     // Find master rows and columns
4376     int * rowCounts = new int [numberRowBlocks];
4377     CoinZeroN(rowCounts, numberRowBlocks);
4378     int * columnCounts = new int [numberColumnBlocks+1];
4379     CoinZeroN(columnCounts, numberColumnBlocks);
4380     int iBlock;
4381     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4382          int iRowBlock = blockInfo[iBlock].rowBlock;
4383          rowCounts[iRowBlock]++;
4384          int iColumnBlock = blockInfo[iBlock].columnBlock;
4385          columnCounts[iColumnBlock]++;
4386     }
4387     int * whichBlock = new int [numberElementBlocks];
4388     int masterRowBlock = -1;
4389     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4390          if (rowCounts[iBlock] > 1) {
4391               if (masterRowBlock == -1) {
4392                    masterRowBlock = iBlock;
4393               } else {
4394                    // Can't decode
4395                    masterRowBlock = -2;
4396                    break;
4397               }
4398          }
4399     }
4400     int masterColumnBlock = -1;
4401     int kBlock = 0;
4402     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4403          int count = columnCounts[iBlock];
4404          columnCounts[iBlock] = kBlock;
4405          kBlock += count;
4406     }
4407     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4408          int iColumnBlock = blockInfo[iBlock].columnBlock;
4409          whichBlock[columnCounts[iColumnBlock]] = iBlock;
4410          columnCounts[iColumnBlock]++;
4411     }
4412     for (iBlock = numberColumnBlocks - 1; iBlock >= 0; iBlock--)
4413          columnCounts[iBlock+1] = columnCounts[iBlock];
4414     columnCounts[0] = 0;
4415     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4416          int count = columnCounts[iBlock+1] - columnCounts[iBlock];
4417          if (count == 1) {
4418               int kBlock = whichBlock[columnCounts[iBlock]];
4419               int iRowBlock = blockInfo[kBlock].rowBlock;
4420               if (iRowBlock == masterRowBlock) {
4421                    if (masterColumnBlock == -1) {
4422                         masterColumnBlock = iBlock;
4423                    } else {
4424                         // Can't decode
4425                         masterColumnBlock = -2;
4426                         break;
4427                    }
4428               }
4429          }
4430     }
4431     if (masterRowBlock < 0 || masterColumnBlock == -2) {
4432          // What now
4433          abort();
4434     }
4435     delete [] rowCounts;
4436     // create all data
4437     const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
4438     ClpSimplex * sub = new ClpSimplex [numberBlocks];
4439     ClpSimplex master;
4440     // Set offset
4441     master.setObjectiveOffset(model->objectiveOffset());
4442     kBlock = 0;
4443     int masterBlock = -1;
4444     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4445          top[kBlock] = NULL;
4446          int start = columnCounts[iBlock];
4447          int end = columnCounts[iBlock+1];
4448          assert (end - start <= 2);
4449          for (int j = start; j < end; j++) {
4450               int jBlock = whichBlock[j];
4451               int iRowBlock = blockInfo[jBlock].rowBlock;
4452               int iColumnBlock = blockInfo[jBlock].columnBlock;
4453               assert (iColumnBlock == iBlock);
4454               if (iColumnBlock != masterColumnBlock && iRowBlock == masterRowBlock) {
4455                    // top matrix
4456                    top[kBlock] = model->coinBlock(jBlock)->packedMatrix();
4457               } else {
4458                    const CoinPackedMatrix * matrix
4459                    = model->coinBlock(jBlock)->packedMatrix();
4460                    // Get pointers to arrays
4461                    const double * rowLower;
4462                    const double * rowUpper;
4463                    const double * columnLower;
4464                    const double * columnUpper;
4465                    const double * objective;
4466                    model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
4467                                 columnLower, columnUpper, objective);
4468                    if (iColumnBlock != masterColumnBlock) {
4469                         // diagonal block
4470                         sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
4471                                                 objective, rowLower, rowUpper);
4472                         if (true) {
4473                              double * lower = sub[kBlock].columnLower();
4474                              double * upper = sub[kBlock].columnUpper();
4475                              int n = sub[kBlock].numberColumns();
4476                              for (int i = 0; i < n; i++) {
4477                                   lower[i] = CoinMax(-1.0e8, lower[i]);
4478                                   upper[i] = CoinMin(1.0e8, upper[i]);
4479                              }
4480                         }
4481                         if (optimizationDirection_ < 0.0) {
4482                              double * obj = sub[kBlock].objective();
4483                              int n = sub[kBlock].numberColumns();
4484                              for (int i = 0; i < n; i++)
4485                                   obj[i] = - obj[i];
4486                         }
4487                         if (this->factorizationFrequency() == 200) {
4488                              // User did not touch preset
4489                              sub[kBlock].defaultFactorizationFrequency();
4490                         } else {
4491                              // make sure model has correct value
4492                              sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
4493                         }
4494                         sub[kBlock].setPerturbation(50);
4495                         // Set columnCounts to be diagonal block index for cleanup
4496                         columnCounts[kBlock] = jBlock;
4497                    } else {
4498                         // master
4499                         masterBlock = jBlock;
4500                         master.loadProblem(*matrix, columnLower, columnUpper,
4501                                            objective, rowLower, rowUpper);
4502                         if (optimizationDirection_ < 0.0) {
4503                              double * obj = master.objective();
4504                              int n = master.numberColumns();
4505                              for (int i = 0; i < n; i++)
4506                                   obj[i] = - obj[i];
4507                         }
4508                    }
4509               }
4510          }
4511          if (iBlock != masterColumnBlock)
4512               kBlock++;
4513     }
4514     delete [] whichBlock;
4515     delete [] blockInfo;
4516     // For now master must have been defined (does not have to have columns)
4517     assert (master.numberRows());
4518     assert (masterBlock >= 0);
4519     int numberMasterRows = master.numberRows();
4520     // Overkill in terms of space
4521     int spaceNeeded = CoinMax(numberBlocks * (numberMasterRows + 1),
4522                               2 * numberMasterRows);
4523     int * rowAdd = new int[spaceNeeded];
4524     double * elementAdd = new double[spaceNeeded];
4525     spaceNeeded = numberBlocks;
4526     int * columnAdd = new int[spaceNeeded+1];
4527     double * objective = new double[spaceNeeded];
4528     // Add in costed slacks
4529     int firstArtificial = master.numberColumns();
4530     int lastArtificial = firstArtificial;
4531     if (true) {
4532          const double * lower = master.rowLower();
4533          const double * upper = master.rowUpper();
4534          int kCol = 0;
4535          for (int iRow = 0; iRow < numberMasterRows; iRow++) {
4536               if (lower[iRow] > -1.0e10) {
4537                    rowAdd[kCol] = iRow;
4538                    elementAdd[kCol++] = 1.0;
4539               }
4540               if (upper[iRow] < 1.0e10) {
4541                    rowAdd[kCol] = iRow;
4542                    elementAdd[kCol++] = -1.0;
4543               }
4544          }
4545          if (kCol > spaceNeeded) {
4546               spaceNeeded = kCol;
4547               delete [] columnAdd;
4548               delete [] objective;
4549               columnAdd = new int[spaceNeeded+1];
4550               objective = new double[spaceNeeded];
4551          }
4552          for (int i = 0; i < kCol; i++) {
4553               columnAdd[i] = i;
4554               objective[i] = 1.0e13;
4555          }
4556          columnAdd[kCol] = kCol;
4557          master.addColumns(kCol, NULL, NULL, objective,
4558                            columnAdd, rowAdd, elementAdd);
4559          lastArtificial = master.numberColumns();
4560     }
4561     int maxPass = 500;
4562     int iPass;
4563     double lastObjective = 1.0e31;
4564     // Create convexity rows for proposals
4565     int numberMasterColumns = master.numberColumns();
4566     master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
4567     if (this->factorizationFrequency() == 200) {
4568          // User did not touch preset
4569          master.defaultFactorizationFrequency();
4570     } else {
4571          // make sure model has correct value
4572          master.setFactorizationFrequency(this->factorizationFrequency());
4573     }
4574     master.setPerturbation(50);
4575     // Arrays to say which block and when created
4576     int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
4577     whichBlock = new int[maximumColumns];
4578     int * when = new int[maximumColumns];
4579     int numberColumnsGenerated = numberBlocks;
4580     // fill in rhs and add in artificials
4581     {
4582          double * rowLower = master.rowLower();
4583          double * rowUpper = master.rowUpper();
4584          int iBlock;
4585          columnAdd[0] = 0;
4586          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4587               int iRow = iBlock + numberMasterRows;;
4588               rowLower[iRow] = 1.0;
4589               rowUpper[iRow] = 1.0;
4590               rowAdd[iBlock] = iRow;
4591               elementAdd[iBlock] = 1.0;
4592               objective[iBlock] = 1.0e13;
4593               columnAdd[iBlock+1] = iBlock + 1;
4594               when[iBlock] = -1;
4595               whichBlock[iBlock] = iBlock;
4596          }
4597          master.addColumns(numberBlocks, NULL, NULL, objective,
4598                            columnAdd, rowAdd, elementAdd);
4599     }
4600     // and resize matrix to double check clp will be happy
4601     //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
4602     //                  numberMasterColumns+numberBlocks);
4603     std::cout << "Time to decompose " << CoinCpuTime() - time1 << " seconds" << std::endl;
4604     for (iPass = 0; iPass < maxPass; iPass++) {
4605          printf("Start of pass %d\n", iPass);
4606          // Solve master - may be infeasible
4607          //master.scaling(0);
4608          if (0) {
4609               master.writeMps("yy.mps");
4610          }
4611          // Correct artificials
4612          double sumArtificials = 0.0;
4613          if (iPass) {
4614               double * upper = master.columnUpper();
4615               double * solution = master.primalColumnSolution();
4616               double * obj = master.objective();
4617               sumArtificials = 0.0;
4618               for (int i = firstArtificial; i < lastArtificial; i++) {
4619                    sumArtificials += solution[i];
4620                    //assert (solution[i]>-1.0e-2);
4621                    if (solution[i] < 1.0e-6) {
4622#if 0
4623                         // Could take out
4624                         obj[i] = 0.0;
4625                         upper[i] = 0.0;
4626#else
4627                         obj[i] = 1.0e7;
4628                         upper[i] = 1.0e-1;
4629#endif
4630                         solution[i] = 0.0;
4631                         master.setColumnStatus(i, isFixed);
4632                    } else {
4633                         upper[i] = solution[i] + 1.0e-5 * (1.0 + solution[i]);
4634                    }
4635               }
4636               printf("Sum of artificials before solve is %g\n", sumArtificials);
4637          }
4638          // scale objective to be reasonable
4639          double scaleFactor = master.scaleObjective(-1.0e9);
4640          {
4641               double * dual = master.dualRowSolution();
4642               int n = master.numberRows();
4643               memset(dual, 0, n * sizeof(double));
4644               double * solution = master.primalColumnSolution();
4645               master.clpMatrix()->times(1.0, solution, dual);
4646               double sum = 0.0;
4647               double * lower = master.rowLower();
4648               double * upper = master.rowUpper();
4649               for (int iRow = 0; iRow < n; iRow++) {
4650                    double value = dual[iRow];
4651                    if (value > upper[iRow])
4652                         sum += value - upper[iRow];
4653                    else if (value < lower[iRow])
4654                         sum -= value - lower[iRow];
4655               }
4656               printf("suminf %g\n", sum);
4657               lower = master.columnLower();
4658               upper = master.columnUpper();
4659               n = master.numberColumns();
4660               for (int iColumn = 0; iColumn < n; iColumn++) {
4661                    double value = solution[iColumn];
4662                    if (value > upper[iColumn] + 1.0e-5)
4663                         sum += value - upper[iColumn];
4664                    else if (value < lower[iColumn] - 1.0e-5)
4665                         sum -= value - lower[iColumn];
4666               }
4667               printf("suminf %g\n", sum);
4668          }
4669          master.primal(1);
4670          // Correct artificials
4671          sumArtificials = 0.0;
4672          {
4673               double * solution = master.primalColumnSolution();
4674               for (int i = firstArtificial; i < lastArtificial; i++) {
4675                    sumArtificials += solution[i];
4676               }
4677               printf("Sum of artificials after solve is %g\n", sumArtificials);
4678          }
4679          master.scaleObjective(scaleFactor);
4680          int problemStatus = master.status(); // do here as can change (delcols)
4681          if (master.numberIterations() == 0 && iPass)
4682               break; // finished
4683          if (master.objectiveValue() > lastObjective - 1.0e-7 && iPass > 555)
4684               break; // finished
4685          lastObjective = master.objectiveValue();
4686          // mark basic ones and delete if necessary
4687          int iColumn;
4688          numberColumnsGenerated = master.numberColumns() - numberMasterColumns;
4689          for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
4690               if (master.getStatus(iColumn + numberMasterColumns) == ClpSimplex::basic)
4691                    when[iColumn] = iPass;
4692          }
4693          if (numberColumnsGenerated + numberBlocks > maximumColumns) {
4694               // delete
4695               int numberKeep = 0;
4696               int numberDelete = 0;
4697               int * whichDelete = new int[numberColumnsGenerated];
4698               for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
4699                    if (when[iColumn] > iPass - 7) {
4700                         // keep
4701                         when[numberKeep] = when[iColumn];
4702                         whichBlock[numberKeep++] = whichBlock[iColumn];
4703                    } else {
4704                         // delete
4705                         whichDelete[numberDelete++] = iColumn + numberMasterColumns;
4706                    }
4707               }
4708               numberColumnsGenerated -= numberDelete;
4709               master.deleteColumns(numberDelete, whichDelete);
4710               delete [] whichDelete;
4711          }
4712          const double * dual = NULL;
4713          bool deleteDual = false;
4714          if (problemStatus == 0) {
4715               dual = master.dualRowSolution();
4716          } else if (problemStatus == 1) {
4717               // could do composite objective
4718               dual = master.infeasibilityRay();
4719               deleteDual = true;
4720               printf("The sum of infeasibilities is %g\n",
4721                      master.sumPrimalInfeasibilities());
4722          } else if (!master.numberColumns()) {
4723               assert(!iPass);
4724               dual = master.dualRowSolution();
4725               memset(master.dualRowSolution(),
4726                      0, (numberMasterRows + numberBlocks)*sizeof(double));
4727          } else {
4728               abort();
4729          }
4730          // Scale back on first time
4731          if (!iPass) {
4732               double * dual2 = master.dualRowSolution();
4733               for (int i = 0; i < numberMasterRows + numberBlocks; i++) {
4734                    dual2[i] *= 1.0e-7;
4735               }
4736               dual = master.dualRowSolution();
4737          }
4738          // Create objective for sub problems and solve
4739          columnAdd[0] = 0;
4740          int numberProposals = 0;
4741          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4742               int numberColumns2 = sub[iBlock].numberColumns();
4743               double * saveObj = new double [numberColumns2];
4744               double * objective2 = sub[iBlock].objective();
4745               memcpy(saveObj, objective2, numberColumns2 * sizeof(double));
4746               // new objective
4747               top[iBlock]->transposeTimes(dual, objective2);
4748               int i;
4749               if (problemStatus == 0) {
4750                    for (i = 0; i < numberColumns2; i++)
4751                         objective2[i] = saveObj[i] - objective2[i];
4752               } else {
4753                    for (i = 0; i < numberColumns2; i++)
4754                         objective2[i] = -objective2[i];
4755               }
4756               // scale objective to be reasonable
4757               double scaleFactor =
4758                    sub[iBlock].scaleObjective((sumArtificials > 1.0e-5) ? -1.0e-4 : -1.0e9);
4759               if (iPass) {
4760                    sub[iBlock].primal();
4761               } else {
4762                    sub[iBlock].dual();
4763               }
4764               sub[iBlock].scaleObjective(scaleFactor);
4765               if (!sub[iBlock].isProvenOptimal() &&
4766                         !sub[iBlock].isProvenDualInfeasible()) {
4767                    memset(objective2, 0, numberColumns2 * sizeof(double));
4768                    sub[iBlock].primal();
4769                    if (problemStatus == 0) {
4770                         for (i = 0; i < numberColumns2; i++)
4771                              objective2[i] = saveObj[i] - objective2[i];
4772                    } else {
4773                         for (i = 0; i < numberColumns2; i++)
4774                              objective2[i] = -objective2[i];
4775                    }
4776                    double scaleFactor = sub[iBlock].scaleObjective(-1.0e9);
4777                    sub[iBlock].primal(1);
4778                    sub[iBlock].scaleObjective(scaleFactor);
4779               }
4780               memcpy(objective2, saveObj, numberColumns2 * sizeof(double));
4781               // get proposal
4782               if (sub[iBlock].numberIterations() || !iPass) {
4783                    double objValue = 0.0;
4784                    int start = columnAdd[numberProposals];
4785                    // proposal
4786                    if (sub[iBlock].isProvenOptimal()) {
4787                         const double * solution = sub[iBlock].primalColumnSolution();
4788                         top[iBlock]->times(solution, elementAdd + start);
4789                         for (i = 0; i < numberColumns2; i++)
4790                              objValue += solution[i] * saveObj[i];
4791                         // See if good dj and pack down
4792                         int number = start;
4793                         double dj = objValue;
4794                         if (problemStatus)
4795                              dj = 0.0;
4796                         double smallest = 1.0e100;
4797                         double largest = 0.0;
4798                         for (i = 0; i < numberMasterRows; i++) {
4799                              double value = elementAdd[start+i];
4800                              if (fabs(value) > 1.0e-15) {
4801                                   dj -= dual[i] * value;
4802                                   smallest = CoinMin(smallest, fabs(value));
4803                                   largest = CoinMax(largest, fabs(value));
4804                                   rowAdd[number] = i;
4805                                   elementAdd[number++] = value;
4806                              }
4807                         }
4808                         // and convexity
4809                         dj -= dual[numberMasterRows+iBlock];
4810                         rowAdd[number] = numberMasterRows + iBlock;
4811                         elementAdd[number++] = 1.0;
4812                         // if elements large then scale?
4813                         //if (largest>1.0e8||smallest<1.0e-8)
4814                         printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
4815                                iBlock, smallest, largest, dj);
4816                         if (dj < -1.0e-6 || !iPass) {
4817                              // take
4818                              objective[numberProposals] = objValue;
4819                              columnAdd[++numberProposals] = number;
4820                              when[numberColumnsGenerated] = iPass;
4821                              whichBlock[numberColumnsGenerated++] = iBlock;
4822                         }
4823                    } else if (sub[iBlock].isProvenDualInfeasible()) {
4824                         // use ray
4825                         const double * solution = sub[iBlock].unboundedRay();
4826                         top[iBlock]->times(solution, elementAdd + start);
4827                         for (i = 0; i < numberColumns2; i++)
4828                              objValue += solution[i] * saveObj[i];
4829                         // See if good dj and pack down
4830                         int number = start;
4831                         double dj = objValue;
4832                         double smallest = 1.0e100;
4833                         double largest = 0.0;
4834                         for (i = 0; i < numberMasterRows; i++) {
4835                              double value = elementAdd[start+i];
4836                              if (fabs(value) > 1.0e-15) {
4837                                   dj -= dual[i] * value;
4838                                   smallest = CoinMin(smallest, fabs(value));
4839                                   largest = CoinMax(largest, fabs(value));
4840                                   rowAdd[number] = i;
4841                                   elementAdd[number++] = value;
4842                              }
4843                         }
4844                         // if elements large or small then scale?
4845                         //if (largest>1.0e8||smallest<1.0e-8)
4846                         printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
4847                                iBlock, smallest, largest, dj);
4848                         if (dj < -1.0e-6) {
4849                              // take
4850                              objective[numberProposals] = objValue;
4851                              columnAdd[++numberProposals] = number;
4852                              when[numberColumnsGenerated] = iPass;
4853                              whichBlock[numberColumnsGenerated++] = iBlock;
4854                         }
4855                    } else {
4856                         abort();
4857                    }
4858               }
4859               delete [] saveObj;
4860          }
4861          if (deleteDual)
4862               delete [] dual;
4863          if (numberProposals)
4864               master.addColumns(numberProposals, NULL, NULL, objective,
4865                                 columnAdd, rowAdd, elementAdd);
4866     }
4867     std::cout << "Time at end of D-W " << CoinCpuTime() - time1 << " seconds" << std::endl;
4868     //master.scaling(0);
4869     //master.primal(1);
4870     loadProblem(*model);
4871     // now put back a good solution
4872     double * lower = new double[numberMasterRows+numberBlocks];
4873     double * upper = new double[numberMasterRows+numberBlocks];
4874     numberColumnsGenerated  += numberMasterColumns;
4875     double * sol = new double[numberColumnsGenerated];
4876     const double * solution = master.primalColumnSolution();
4877     const double * masterLower = master.rowLower();
4878     const double * masterUpper = master.rowUpper();
4879     double * fullSolution = primalColumnSolution();
4880     const double * fullLower = columnLower();
4881     const double * fullUpper = columnUpper();
4882     const double * rowSolution = master.primalRowSolution();
4883     double * fullRowSolution = primalRowSolution();
4884     const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4885     int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4886     const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4887     int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4888     for (int iRow = 0; iRow < numberRows2; iRow++) {
4889          int kRow = rowBack[iRow];
4890          setRowStatus(kRow, master.getRowStatus(iRow));
4891          fullRowSolution[kRow] = rowSolution[iRow];
4892     }
4893     for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4894          int kColumn = columnBack[iColumn];
4895          setStatus(kColumn, master.getStatus(iColumn));
4896          fullSolution[kColumn] = solution[iColumn];
4897     }
4898     for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4899          // move basis
4900          int kBlock = columnCounts[iBlock];
4901          const int * rowBack = model->coinBlock(kBlock)->originalRows();
4902          int numberRows2 = model->coinBlock(kBlock)->numberRows();
4903          const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4904          int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4905          for (int iRow = 0; iRow < numberRows2; iRow++) {
4906               int kRow = rowBack[iRow];
4907               setRowStatus(kRow, sub[iBlock].getRowStatus(iRow));
4908          }
4909          for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4910               int kColumn = columnBack[iColumn];
4911               setStatus(kColumn, sub[iBlock].getStatus(iColumn));
4912          }
4913          // convert top bit to by rows
4914          CoinPackedMatrix topMatrix = *top[iBlock];
4915          topMatrix.reverseOrdering();
4916          // zero solution
4917          memset(sol, 0, numberColumnsGenerated * sizeof(double));
4918
4919          for (int i = numberMasterColumns; i < numberColumnsGenerated; i++) {
4920               if (whichBlock[i-numberMasterColumns] == iBlock)
4921                    sol[i] = solution[i];
4922          }
4923          memset(lower, 0, (numberMasterRows + numberBlocks)*sizeof(double));
4924          master.clpMatrix()->times(1.0, sol, lower);
4925          for (int iRow = 0; iRow < numberMasterRows; iRow++) {
4926               double value = lower[iRow];
4927               if (masterUpper[iRow] < 1.0e20)
4928                    upper[iRow] = value;
4929               else
4930                    upper[iRow] = COIN_DBL_MAX;
4931               if (masterLower[iRow] > -1.0e20)
4932                    lower[iRow] = value;
4933               else
4934                    lower[iRow] = -COIN_DBL_MAX;
4935          }
4936          sub[iBlock].addRows(numberMasterRows, lower, upper,
4937                              topMatrix.getVectorStarts(),
4938                              topMatrix.getVectorLengths(),
4939                              topMatrix.getIndices(),
4940                              topMatrix.getElements());
4941          sub[iBlock].primal(1);
4942          const double * subSolution = sub[iBlock].primalColumnSolution();
4943          const double * subRowSolution = sub[iBlock].primalRowSolution();
4944          // move solution
4945          for (int iRow = 0; iRow < numberRows2; iRow++) {
4946               int kRow = rowBack[iRow];
4947               fullRowSolution[kRow] = subRowSolution[iRow];
4948          }
4949          for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4950               int kColumn = columnBack[iColumn];
4951               fullSolution[kColumn] = subSolution[iColumn];
4952          }
4953     }
4954     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4955          if (fullSolution[iColumn] < fullUpper[iColumn] - 1.0e-8 &&
4956                    fullSolution[iColumn] > fullLower[iColumn] + 1.0e-8) {
4957               if (getStatus(iColumn) != ClpSimplex::basic) {
4958                    if (columnLower_[iColumn] > -1.0e30 ||
4959                              columnUpper_[iColumn] < 1.0e30)
4960                         setStatus(iColumn, ClpSimplex::superBasic);
4961                    else
4962                         setStatus(iColumn, ClpSimplex::isFree);
4963               }
4964          } else if (fullSolution[iColumn] >= fullUpper[iColumn] - 1.0e-8) {
4965               // may help to make rest non basic
4966               if (getStatus(iColumn) != ClpSimplex::basic)
4967                    setStatus(iColumn, ClpSimplex::atUpperBound);
4968          } else if (fullSolution[iColumn] <= fullLower[iColumn] + 1.0e-8) {
4969               // may help to make rest non basic
4970               if (getStatus(iColumn) != ClpSimplex::basic)
4971                    setStatus(iColumn, ClpSimplex::atLowerBound);
4972          }
4973     }
4974     //int numberRows=model->numberRows();
4975     //for (int iRow=0;iRow<numberRows;iRow++)
4976     //setRowStatus(iRow,ClpSimplex::superBasic);
4977     std::cout << "Time before cleanup of full model " << CoinCpuTime() - time1 << " seconds" << std::endl;
4978     primal(1);
4979     std::cout << "Total time " << CoinCpuTime() - time1 << " seconds" << std::endl;
4980     delete [] columnCounts;
4981     delete [] sol;
4982     delete [] lower;
4983     delete [] upper;
4984     delete [] whichBlock;
4985     delete [] when;
4986     delete [] columnAdd;
4987     delete [] rowAdd;
4988     delete [] elementAdd;
4989     delete [] objective;
4990     delete [] top;
4991     delete [] sub;
4992     return 0;
4993}
4994// Solve using Benders decomposition and maybe in parallel
4995int
4996ClpSimplex::solveBenders(CoinStructuredModel * model)
4997{
4998     double time1 = CoinCpuTime();
4999     //int numberColumns = model->numberColumns();
5000     int numberRowBlocks = model->numberRowBlocks();
5001     int numberColumnBlocks = model->numberColumnBlocks();
5002     int numberElementBlocks = model->numberElementBlocks();
5003     // We already have top level structure
5004     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
5005     for (int i = 0; i < numberElementBlocks; i++) {
5006          CoinModel * thisBlock = model->coinBlock(i);
5007          assert (thisBlock);
5008          // just fill in info
5009          CoinModelBlockInfo info = CoinModelBlockInfo();
5010          int whatsSet = thisBlock->whatIsSet();
5011          info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
5012          info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
5013          info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
5014          info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
5015          info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
5016          info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
5017          // Which block
5018          int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
5019          info.rowBlock = iRowBlock;
5020          int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
5021          info.columnBlock = iColumnBlock;
5022          blockInfo[i] = info;
5023     }
5024     // make up problems
5025     int numberBlocks = numberColumnBlocks - 1;
5026     // Find master columns and rows
5027     int * columnCounts = new int [numberColumnBlocks];
5028     CoinZeroN(columnCounts, numberColumnBlocks);
5029     int * rowCounts = new int [numberRowBlocks+1];
5030     CoinZeroN(rowCounts, numberRowBlocks);
5031     int iBlock;
5032     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5033          int iColumnBlock = blockInfo[iBlock].columnBlock;
5034          columnCounts[iColumnBlock]++;
5035          int iRowBlock = blockInfo[iBlock].rowBlock;
5036          rowCounts[iRowBlock]++;
5037     }
5038     int * whichBlock = new int [numberElementBlocks];
5039     int masterColumnBlock = -1;
5040     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5041          if (columnCounts[iBlock] > 1) {
5042               if (masterColumnBlock == -1) {
5043                    masterColumnBlock = iBlock;
5044               } else {
5045                    // Can't decode
5046                    masterColumnBlock = -2;
5047                    break;
5048               }
5049          }
5050     }
5051     int masterRowBlock = -1;
5052     int kBlock = 0;
5053     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5054          int count = rowCounts[iBlock];
5055          rowCounts[iBlock] = kBlock;
5056          kBlock += count;
5057     }
5058     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5059          int iRowBlock = blockInfo[iBlock].rowBlock;
5060          whichBlock[rowCounts[iRowBlock]] = iBlock;
5061          rowCounts[iRowBlock]++;
5062     }
5063     for (iBlock = numberRowBlocks - 1; iBlock >= 0; iBlock--)
5064          rowCounts[iBlock+1] = rowCounts[iBlock];
5065     rowCounts[0] = 0;
5066     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5067          int count = rowCounts[iBlock+1] - rowCounts[iBlock];
5068          if (count == 1) {
5069               int kBlock = whichBlock[rowCounts[iBlock]];
5070               int iColumnBlock = blockInfo[kBlock].columnBlock;
5071               if (iColumnBlock == masterColumnBlock) {
5072                    if (masterRowBlock == -1) {
5073                         masterRowBlock = iBlock;
5074                    } else {
5075                         // Can't decode
5076                         masterRowBlock = -2;
5077                         break;
5078                    }
5079               }
5080          }
5081     }
5082     if (masterColumnBlock < 0 || masterRowBlock == -2) {
5083          // What now
5084          abort();
5085     }
5086     delete [] columnCounts;
5087     // create all data
5088     const CoinPackedMatrix ** first = new const CoinPackedMatrix * [numberRowBlocks];
5089     ClpSimplex * sub = new ClpSimplex [numberBlocks];
5090     ClpSimplex master;
5091     // Set offset
5092     master.setObjectiveOffset(model->objectiveOffset());
5093     kBlock = 0;
5094     int masterBlock = -1;
5095     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5096          first[kBlock] = NULL;
5097          int start = rowCounts[iBlock];
5098          int end = rowCounts[iBlock+1];
5099          assert (end - start <= 2);
5100          for (int j = start; j < end; j++) {
5101               int jBlock = whichBlock[j];
5102               int iColumnBlock = blockInfo[jBlock].columnBlock;
5103               int iRowBlock = blockInfo[jBlock].rowBlock;
5104               assert (iRowBlock == iBlock);
5105               if (iRowBlock != masterRowBlock && iColumnBlock == masterColumnBlock) {
5106                    // first matrix
5107                    first[kBlock] = model->coinBlock(jBlock)->packedMatrix();
5108               } else {
5109                    const CoinPackedMatrix * matrix
5110                    = model->coinBlock(jBlock)->packedMatrix();
5111                    // Get pointers to arrays
5112                    const double * columnLower;
5113                    const double * columnUpper;
5114                    const double * rowLower;
5115                    const double * rowUpper;
5116                    const double * objective;
5117                    model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
5118                                 columnLower, columnUpper, objective);
5119                    if (iRowBlock != masterRowBlock) {
5120                         // diagonal block
5121                         sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
5122                                                 objective, rowLower, rowUpper);
5123                         if (optimizationDirection_ < 0.0) {
5124                              double * obj = sub[kBlock].objective();
5125                              int n = sub[kBlock].numberColumns();
5126