source: stable/1.16/Clp/src/ClpSolve.cpp @ 2188

Last change on this file since 2188 was 2188, checked in by forrest, 5 years ago

fix misplaced ;

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