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

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

Minor stability fixes and an option to allow perturbation after presolve

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