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

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

fix a few problems with Aboca

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