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

Last change on this file since 1141 was 1141, checked in by forrest, 12 years ago

for threaded random numbers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 78.3 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4// This file has higher level solve functions
5
6#include "ClpConfig.h"
7#include "CoinPragma.hpp"
8
9#include <math.h>
10
11#include "CoinHelperFunctions.hpp"
12#include "ClpHelperFunctions.hpp"
13#include "CoinSort.hpp"
14#include "ClpFactorization.hpp"
15#include "ClpSimplex.hpp"
16#include "ClpSimplexOther.hpp"
17#ifndef SLIM_CLP
18#include "ClpQuadraticObjective.hpp"
19#include "ClpInterior.hpp"
20#include "ClpCholeskyDense.hpp"
21#include "ClpCholeskyBase.hpp"
22#include "ClpPlusMinusOneMatrix.hpp"
23#include "ClpNetworkMatrix.hpp"
24#endif
25#include "ClpLinearObjective.hpp"
26#include "ClpSolve.hpp"
27#include "ClpPackedMatrix.hpp"
28#include "ClpMessage.hpp"
29#include "CoinTime.hpp"
30
31#include "ClpPresolve.hpp"
32#ifndef SLIM_CLP
33#include "Idiot.hpp"
34#ifdef WSSMP_BARRIER
35#include "ClpCholeskyWssmp.hpp"
36#include "ClpCholeskyWssmpKKT.hpp"
37#define FAST_BARRIER
38#endif
39#ifdef UFL_BARRIER
40#include "ClpCholeskyUfl.hpp"
41#define FAST_BARRIER
42#endif
43#ifdef TAUCS_BARRIER
44#include "ClpCholeskyTaucs.hpp"
45#define FAST_BARRIER
46#endif
47#ifdef COIN_DEVELOP
48#ifndef FAST_BARRIER
49static int numberBarrier=0;
50#endif
51#endif
52#ifdef COIN_HAS_VOL
53#include "VolVolume.hpp"
54#include "CoinHelperFunctions.hpp"
55#include "CoinPackedMatrix.hpp"
56#include "CoinMpsIO.hpp"
57
58//#############################################################################
59
60class lpHook : public VOL_user_hooks {
61private:
62   lpHook(const lpHook&);
63   lpHook& operator=(const lpHook&);
64private:
65   /// Pointer to dense vector of structural variable upper bounds
66   double  *colupper_;
67   /// Pointer to dense vector of structural variable lower bounds
68   double  *collower_;
69   /// Pointer to dense vector of objective coefficients
70   double  *objcoeffs_;
71   /// Pointer to dense vector of right hand sides
72   double  *rhs_;
73   /// Pointer to dense vector of senses
74   char    *sense_;
75
76   /// The problem matrix in a row ordered form
77   CoinPackedMatrix rowMatrix_;
78   /// The problem matrix in a column ordered form
79   CoinPackedMatrix colMatrix_;
80
81public:
82   lpHook(double* clb, double* cub, double* obj,
83          double* rhs, char* sense, const CoinPackedMatrix& mat);
84   virtual ~lpHook();
85   
86public:
87   // for all hooks: return value of -1 means that volume should quit
88   /** compute reduced costs   
89       @param u (IN) the dual variables
90       @param rc (OUT) the reduced cost with respect to the dual values
91   */
92   virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc);
93
94   /** Solve the subproblem for the subgradient step.
95       @param dual (IN) the dual variables
96       @param rc (IN) the reduced cost with respect to the dual values
97       @param lcost (OUT) the lagrangean cost with respect to the dual values
98       @param x (OUT) the primal result of solving the subproblem
99       @param v (OUT) b-Ax for the relaxed constraints
100       @param pcost (OUT) the primal objective value of <code>x</code>
101   */
102   virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
103                                double& lcost, VOL_dvector& x, VOL_dvector& v,
104                                double& pcost);
105   /** Starting from the primal vector x, run a heuristic to produce
106       an integer solution 
107       @param x (IN) the primal vector
108       @param heur_val (OUT) the value of the integer solution (return
109       <code>DBL_MAX</code> here if no feas sol was found
110   */
111   virtual int heuristics(const VOL_problem& p, 
112                          const VOL_dvector& x, double& heur_val) {
113      return 0;
114   }
115};
116 
117//#############################################################################
118
119lpHook::lpHook(double* clb, double* cub, double* obj,
120               double* rhs, char* sense,
121               const CoinPackedMatrix& mat)
122{
123   colupper_ = cub;
124   collower_ = clb;
125   objcoeffs_ = obj;
126   rhs_ = rhs;
127   sense_ = sense;
128   assert (mat.isColOrdered());
129   colMatrix_.copyOf(mat);
130   rowMatrix_.reverseOrderedCopyOf(mat);
131}
132
133//-----------------------------------------------------------------------------
134
135lpHook::~lpHook()
136{
137}
138
139//#############################################################################
140
141int
142lpHook::compute_rc(const VOL_dvector& u, VOL_dvector& rc)
143{
144   rowMatrix_.transposeTimes(u.v, rc.v);
145   const int psize = rowMatrix_.getNumCols();
146
147   for (int i = 0; i < psize; ++i)
148      rc[i] = objcoeffs_[i] - rc[i];
149   return 0;
150}
151
152//-----------------------------------------------------------------------------
153
154int
155lpHook::solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
156                         double& lcost, VOL_dvector& x, VOL_dvector& v,
157                         double& pcost)
158{
159   int i;
160   const int psize = x.size();
161   const int dsize = v.size();
162
163   // compute the lagrangean solution corresponding to the reduced costs
164   for (i = 0; i < psize; ++i) 
165      x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
166
167   // compute the lagrangean value (rhs*dual + primal*rc)
168   lcost = 0;
169   for (i = 0; i < dsize; ++i)
170      lcost += rhs_[i] * dual[i];
171   for (i = 0; i < psize; ++i)
172      lcost += x[i] * rc[i];
173
174   // compute the rhs - lhs
175   colMatrix_.times(x.v, v.v);
176   for (i = 0; i < dsize; ++i)
177      v[i] = rhs_[i] - v[i];
178
179   // compute the lagrangean primal objective
180   pcost = 0;
181   for (i = 0; i < psize; ++i)
182      pcost += x[i] * objcoeffs_[i];
183
184   return 0;
185}
186
187//#############################################################################
188/** A quick inlined function to convert from lb/ub style constraint
189    definition to sense/rhs/range style */
190inline void
191convertBoundToSense(const double lower, const double upper,
192                                        char& sense, double& right,
193                                        double& range) 
194{
195  range = 0.0;
196  if (lower > -1.0e20) {
197    if (upper < 1.0e20) {
198      right = upper;
199      if (upper==lower) {
200        sense = 'E';
201      } else {
202        sense = 'R';
203        range = upper - lower;
204      }
205    } else {
206      sense = 'G';
207      right = lower;
208    }
209  } else {
210    if (upper < 1.0e20) {
211      sense = 'L';
212      right = upper;
213    } else {
214      sense = 'N';
215      right = 0.0;
216    }
217  }
218}
219
220static int
221solveWithVolume(ClpSimplex * model, int numberPasses, int doIdiot)
222{
223   VOL_problem volprob;
224   volprob.parm.gap_rel_precision=0.00001;
225   volprob.parm.maxsgriters=3000;
226   if(numberPasses>3000) {
227     volprob.parm.maxsgriters=numberPasses;
228     volprob.parm.primal_abs_precision=0.0;
229     volprob.parm.minimum_rel_ascent=0.00001;
230   } else if (doIdiot>0) {
231     volprob.parm.maxsgriters=doIdiot;
232   }
233   if (model->logLevel()<2) 
234     volprob.parm.printflag=0;
235   else
236     volprob.parm.printflag=3;
237   const CoinPackedMatrix* mat = model->matrix();
238   int psize = model->numberColumns();
239   int dsize = model->numberRows();
240   char * sense = new char[dsize];
241   double * rhs = new double [dsize];
242
243   // Set the lb/ub on the duals
244   volprob.dsize = dsize;
245   volprob.psize = psize;
246   volprob.dual_lb.allocate(dsize);
247   volprob.dual_ub.allocate(dsize);
248   int i;
249   const double * rowLower = model->rowLower();
250   const double * rowUpper = model->rowUpper();
251   for (i = 0; i < dsize; ++i) {
252     double range;
253     convertBoundToSense(rowLower[i],rowUpper[i],
254                         sense[i],rhs[i],range);
255      switch (sense[i]) {
256       case 'E':
257         volprob.dual_lb[i] = -1.0e31;
258         volprob.dual_ub[i] = 1.0e31;
259         break;
260       case 'L':
261         volprob.dual_lb[i] = -1.0e31;
262         volprob.dual_ub[i] = 0.0;
263         break;
264       case 'G':
265         volprob.dual_lb[i] = 0.0;
266         volprob.dual_ub[i] = 1.0e31;
267         break;
268       default:
269         printf("Volume Algorithm can't work if there is a non ELG row\n");
270         return 1;
271      }
272   }
273   // Check bounds
274   double * saveLower = model->columnLower();
275   double * saveUpper = model->columnUpper();
276   bool good=true;
277   for (i=0;i<psize;i++) {
278     if (saveLower[i]<-1.0e20||saveUpper[i]>1.0e20) {
279       good=false;
280       break;
281     }
282   }
283   if (!good) {
284     saveLower = CoinCopyOfArray(model->columnLower(),psize);
285     saveUpper = CoinCopyOfArray(model->columnUpper(),psize);
286     for (i=0;i<psize;i++) {
287       if (saveLower[i]<-1.0e20)
288         saveLower[i]=-1.0e20;
289       if(saveUpper[i]>1.0e20) 
290         saveUpper[i]=1.0e20;
291     }
292   }
293   lpHook myHook(saveLower, saveUpper,
294                 model->objective(),
295                 rhs, sense, *mat);
296
297   volprob.solve(myHook, false /* no warmstart */);
298
299   if (saveLower!=model->columnLower()) {
300     delete [] saveLower;
301     delete [] saveUpper;
302   }
303   //------------- extract the solution ---------------------------
304
305   //printf("Best lagrangean value: %f\n", volprob.value);
306
307   double avg = 0;
308   for (i = 0; i < dsize; ++i) {
309      switch (sense[i]) {
310       case 'E':
311         avg += CoinAbs(volprob.viol[i]);
312         break;
313       case 'L':
314         if (volprob.viol[i] < 0)
315            avg +=  (-volprob.viol[i]);
316         break;
317       case 'G':
318         if (volprob.viol[i] > 0)
319            avg +=  volprob.viol[i];
320         break;
321      }
322   }
323     
324   //printf("Average primal constraint violation: %f\n", avg/dsize);
325
326   // volprob.dsol contains the dual solution (dual feasible)
327   // volprob.psol contains the primal solution
328   //              (NOT necessarily primal feasible)
329   CoinMemcpyN(volprob.dsol.v,dsize,model->dualRowSolution());
330   CoinMemcpyN(volprob.psol.v,psize,model->primalColumnSolution());
331   return 0;
332}
333#endif
334static ClpInterior * currentModel2 = NULL;
335#endif
336//#############################################################################
337// Allow for interrupts
338// But is this threadsafe ? (so switched off by option)
339
340#include "CoinSignal.hpp"
341static ClpSimplex * currentModel = NULL;
342
343extern "C" {
344   static void signal_handler(int whichSignal)
345   {
346      if (currentModel!=NULL) 
347         currentModel->setMaximumIterations(0); // stop at next iterations
348#ifndef SLIM_CLP
349      if (currentModel2!=NULL) 
350         currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
351#endif
352      return;
353   }
354}
355
356/** General solve algorithm which can do presolve
357    special options (bits)
358    1 - do not perturb
359    2 - do not scale
360    4 - use crash (default allslack in dual, idiot in primal)
361    8 - all slack basis in primal
362    16 - switch off interrupt handling
363    32 - do not try and make plus minus one matrix
364    64 - do not use sprint even if problem looks good
365 */
366int 
367ClpSimplex::initialSolve(ClpSolve & options)
368{
369  ClpSolve::SolveType method=options.getSolveType();
370  //ClpSolve::SolveType originalMethod=method;
371  ClpSolve::PresolveType presolve = options.getPresolveType();
372  int saveMaxIterations = maximumIterations();
373  int finalStatus=-1;
374  int numberIterations=0;
375  double time1 = CoinCpuTime();
376  double timeX = time1;
377  double time2;
378  ClpMatrixBase * saveMatrix=NULL;
379  ClpObjective * savedObjective=NULL;
380  if (!objective_||!matrix_) {
381    // totally empty
382    handler_->message(CLP_EMPTY_PROBLEM,messages_)
383      <<0
384      <<0
385      <<0
386      <<CoinMessageEol;
387    return -1;
388  } else if (!numberRows_||!numberColumns_||!getNumElements()) {
389    presolve = ClpSolve::presolveOff;
390  }
391  if (objective_->type()>=2&&optimizationDirection_==0) {
392    // pretend linear
393    savedObjective=objective_;
394    // make up objective
395    double * obj = new double[numberColumns_];
396    for (int i=0;i<numberColumns_;i++) {
397      double l = fabs(columnLower_[i]);
398      double u = fabs(columnUpper_[i]);
399      obj[i]=0.0;
400      if (CoinMin(l,u)<1.0e20) {
401        if (l<u) 
402          obj[i]=1.0+randomNumberGenerator_.randomDouble()*1.0e-2;
403        else
404          obj[i]=-1.0-randomNumberGenerator_.randomDouble()*1.0e-2;
405      }
406    }
407    objective_= new ClpLinearObjective(obj,numberColumns_);
408    delete [] obj;
409  }
410  ClpSimplex * model2 = this;
411  bool interrupt = (options.getSpecialOption(2)==0);
412  CoinSighandler_t saveSignal=SIG_DFL;
413  if (interrupt) {
414    currentModel = model2;
415    // register signal handler
416    saveSignal = signal(SIGINT,signal_handler);
417  }
418  // If no status array - set up basis
419  if (!status_)
420    allSlackBasis();
421  ClpPresolve pinfo;
422  pinfo.setSubstitution(options.substitution());
423  int presolveOptions = options.presolveActions();
424  bool presolveToFile = (presolveOptions&0x40000000)!=0;
425  presolveOptions &= ~0x40000000;
426  if ((presolveOptions&0xffff)!=0)
427    pinfo.setPresolveActions(presolveOptions);
428  // switch off singletons to slacks
429  //pinfo.setDoSingletonColumn(false); // done by bits
430  int printOptions = options.getSpecialOption(5);
431  if ((printOptions&1)!=0)
432    pinfo.statistics();
433  double timePresolve=0.0;
434  double timeIdiot=0.0;
435  double timeCore=0.0;
436  int savePerturbation=perturbation_;
437  int saveScaling = scalingFlag_;
438#ifndef SLIM_CLP
439#ifndef NO_RTTI
440  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
441    // network - switch off stuff
442    presolve = ClpSolve::presolveOff;
443  }
444#else
445  if (matrix_->type()==11) {
446    // network - switch off stuff
447    presolve = ClpSolve::presolveOff;
448  }
449#endif
450#endif
451  // For below >0 overrides
452  // 0 means no, -1 means maybe
453  int doIdiot=0;
454  int doCrash=0;
455  int doSprint=0;
456  int doSlp=0;
457  int primalStartup=1;
458  if (method!=ClpSolve::useDual&&method!=ClpSolve::useBarrier
459      &&method!=ClpSolve::useBarrierNoCross) {
460    switch (options.getSpecialOption(1)) {
461    case 0:
462      doIdiot=-1;
463      doCrash=-1;
464      doSprint=-1;
465      break;
466    case 1:
467      doIdiot=0;
468      doCrash=1;
469      if (options.getExtraInfo(1)>0)
470        doCrash = options.getExtraInfo(1);
471      doSprint=0;
472      break;
473    case 2:
474      doIdiot=1;
475      if (options.getExtraInfo(1)>0)
476        doIdiot = options.getExtraInfo(1);
477      doCrash=0;
478      doSprint=0;
479      break;
480    case 3:
481      doIdiot=0;
482      doCrash=0;
483      doSprint=1;
484      break;
485    case 4:
486      doIdiot=0;
487      doCrash=0;
488      doSprint=0;
489      break;
490    case 5:
491      doIdiot=0;
492      doCrash=-1;
493      doSprint=-1;
494      break;
495    case 6:
496      doIdiot=-1;
497      doCrash=-1;
498      doSprint=0;
499      break;
500    case 7:
501      doIdiot=-1;
502      doCrash=0;
503      doSprint=-1;
504      break;
505    case 8:
506      doIdiot=-1;
507      doCrash=0;
508      doSprint=0;
509      break;
510    case 9:
511      doIdiot=0;
512      doCrash=0;
513      doSprint=-1;
514      break;
515    case 10:
516      doIdiot=0;
517      doCrash=0;
518      doSprint=0;
519      if (options.getExtraInfo(1)>0)
520        doSlp = options.getExtraInfo(1);
521      break;
522    case 11:
523      doIdiot=0;
524      doCrash=0;
525      doSprint=0;
526      primalStartup=0;
527      break;
528    default:
529      abort();
530    }
531  } else {
532    // Dual
533    switch (options.getSpecialOption(0)) {
534    case 0:
535      doIdiot=0;
536      doCrash=0;
537      doSprint=0;
538      break;
539    case 1:
540      doIdiot=0;
541      doCrash=1;
542      if (options.getExtraInfo(0)>0)
543        doCrash = options.getExtraInfo(0);
544      doSprint=0;
545      break;
546    case 2:
547      doIdiot=-1;
548      if (options.getExtraInfo(0)>0)
549        doIdiot = options.getExtraInfo(0);
550      doCrash=0;
551      doSprint=0;
552      break;
553    default:
554      abort();
555    }
556  }
557#ifndef NO_RTTI
558  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objectiveAsObject()));
559#else
560  ClpQuadraticObjective * quadraticObj = NULL;
561  if (objective_->type()==2)
562    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
563#endif
564  // If quadratic then primal or barrier or slp
565  if (quadraticObj) {
566    doSprint=0;
567    doIdiot=0;
568    // off
569    if (method==ClpSolve::useBarrier)
570      method=ClpSolve::useBarrierNoCross;
571    else if (method!=ClpSolve::useBarrierNoCross)
572      method=ClpSolve::usePrimal;
573  }
574#ifdef COIN_HAS_VOL
575  // Save number of idiot
576  int saveDoIdiot=doIdiot;
577#endif
578  // Just do this number of passes in Sprint
579  int maxSprintPass=100;
580  bool costedSlacks=false;
581  if (presolve!=ClpSolve::presolveOff) {
582    int numberPasses=5;
583    if (presolve==ClpSolve::presolveNumber) {
584      numberPasses=options.getPresolvePasses();
585      presolve = ClpSolve::presolveOn;
586    } else if (presolve==ClpSolve::presolveNumberCost) {
587      numberPasses=options.getPresolvePasses();
588      presolve = ClpSolve::presolveOn;
589      costedSlacks=true;
590      // switch on singletons to slacks
591      pinfo.setDoSingletonColumn(true);
592    }
593#ifndef CLP_NO_STD
594    if (presolveToFile) {
595      // PreSolve to file - not fully tested
596      printf("Presolving to file - presolve.save\n");
597      pinfo.presolvedModelToFile(*this,"presolve.save",dblParam_[ClpPresolveTolerance],
598                           false,numberPasses);
599      model2=this;
600    } else {
601#endif
602      model2 = pinfo.presolvedModel(*this,dblParam_[ClpPresolveTolerance],
603                                    false,numberPasses,true,costedSlacks);
604#ifndef CLP_NO_STD
605    }
606#endif
607    time2 = CoinCpuTime();
608    timePresolve = time2-timeX;
609    handler_->message(CLP_INTERVAL_TIMING,messages_)
610      <<"Presolve"<<timePresolve<<time2-time1
611      <<CoinMessageEol;
612    timeX=time2;
613    if (!model2) {
614      handler_->message(CLP_INFEASIBLE,messages_)
615        <<CoinMessageEol;
616      model2 = this;
617      problemStatus_=1; // may be unbounded but who cares
618      if (options.infeasibleReturn()||(moreSpecialOptions_&1)!=0) {
619        return -1;
620      }
621      presolve=ClpSolve::presolveOff;
622    }
623    // We may be better off using original (but if dual leave because of bounds)
624    if (presolve!=ClpSolve::presolveOff&&
625        numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_
626        &&model2!=this) {
627      if(method!=ClpSolve::useDual||
628         (numberRows_==model2->numberRows_&&numberColumns_==model2->numberColumns_)) {
629        delete model2;
630        model2 = this;
631        presolve=ClpSolve::presolveOff;
632      }
633    }
634  }
635  if (interrupt)
636    currentModel = model2;
637  // See if worth trying +- one matrix
638  bool plusMinus=false;
639  int numberElements=model2->getNumElements();
640#ifndef SLIM_CLP
641#ifndef NO_RTTI
642  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
643    // network - switch off stuff
644    doIdiot=0;
645    if (doSprint<0)
646      doSprint=0;
647  }
648#else
649  if (matrix_->type()==11) {
650    // network - switch off stuff
651    doIdiot=0;
652    //doSprint=0;
653  }
654#endif
655#endif
656  int numberColumns = model2->numberColumns();
657  int numberRows = model2->numberRows();
658  // If not all slack basis - switch off all except sprint
659  int numberRowsBasic=0;
660  int iRow;
661  for (iRow=0;iRow<numberRows;iRow++)
662    if (model2->getRowStatus(iRow)==basic)
663      numberRowsBasic++;
664  if (numberRowsBasic<numberRows) {
665    doIdiot=0;
666    doCrash=0;
667    //doSprint=0;
668  }
669  if (options.getSpecialOption(3)==0) {
670    if(numberElements>100000)
671      plusMinus=true;
672    if(numberElements>10000&&(doIdiot||doSprint)) 
673      plusMinus=true;
674  }
675  //plusMinus=true;
676#ifndef SLIM_CLP
677  // Statistics (+1,-1, other) - used to decide on strategy if not +-1
678  CoinBigIndex statistics[3]={-1,0,0};
679  if (plusMinus) {
680    saveMatrix = model2->clpMatrix();
681#ifndef NO_RTTI
682    ClpPackedMatrix* clpMatrix =
683      dynamic_cast< ClpPackedMatrix*>(saveMatrix);
684#else
685    ClpPackedMatrix* clpMatrix = NULL;
686    if (saveMatrix->type()==1)
687      clpMatrix =
688        static_cast< ClpPackedMatrix*>(saveMatrix);
689#endif
690    if (clpMatrix) {
691      ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
692      if (newMatrix->getIndices()) {
693        model2->replaceMatrix(newMatrix);
694      } else {
695        handler_->message(CLP_MATRIX_CHANGE,messages_)
696          <<"+- 1"
697          <<CoinMessageEol;
698        CoinMemcpyN(newMatrix->startPositive(),3,statistics);
699        saveMatrix=NULL;
700        plusMinus=false;
701        delete newMatrix;
702      }
703    } else {
704      saveMatrix=NULL;
705      plusMinus=false;
706    }
707  }
708#endif
709  if (this->factorizationFrequency()==200) {
710    // User did not touch preset
711    model2->defaultFactorizationFrequency();
712  } else if (model2!=this) {
713    // make sure model2 has correct value
714    model2->setFactorizationFrequency(this->factorizationFrequency());
715  }
716  if (method==ClpSolve::automatic) {
717    if (doSprint==0&&doIdiot==0) {
718      // off
719      method=ClpSolve::useDual;
720    } else {
721      // only do primal if sprint or idiot
722      if (doSprint>0) {
723        method=ClpSolve::usePrimalorSprint;
724      } else if (doIdiot>0) {
725        method=ClpSolve::usePrimal;
726      } else {
727        if (numberElements<500000) {
728          // Small problem
729          if(numberRows*10>numberColumns||numberColumns<6000
730             ||(numberRows*20>numberColumns&&!plusMinus))
731            doSprint=0; // switch off sprint
732        } else {
733          // larger problem
734          if(numberRows*8>numberColumns)
735            doSprint=0; // switch off sprint
736        }
737        // switch off sprint or idiot if any free variable
738        int iColumn;
739        double * columnLower = model2->columnLower();
740        double * columnUpper = model2->columnUpper();
741        for (iColumn=0;iColumn<numberColumns;iColumn++) {
742          if (columnLower[iColumn]<-1.0e10&&columnUpper[iColumn]>1.0e10) {
743            doSprint=0;
744            doIdiot=0;
745            break;
746          }
747        }
748        int nPasses=0;
749        // look at rhs
750        int iRow;
751        double largest=0.0;
752        double smallest = 1.0e30;
753        double largestGap=0.0;
754        int numberNotE=0;
755        bool notInteger=false;
756        for (iRow=0;iRow<numberRows;iRow++) {
757          double value1 = model2->rowLower_[iRow];
758          if (value1&&value1>-1.0e31) {
759            largest = CoinMax(largest,fabs(value1));
760            smallest=CoinMin(smallest,fabs(value1));
761            if (fabs(value1-floor(value1+0.5))>1.0e-8) {
762              notInteger=true;
763              break;
764            }
765          }
766          double value2 = model2->rowUpper_[iRow];
767          if (value2&&value2<1.0e31) {
768            largest = CoinMax(largest,fabs(value2));
769            smallest=CoinMin(smallest,fabs(value2));
770            if (fabs(value2-floor(value2+0.5))>1.0e-8) {
771              notInteger=true;
772              break;
773            }
774          }
775          if (value2>value1) {
776            numberNotE++;
777            if (value2>1.0e31||value1<-1.0e31)
778              largestGap = COIN_DBL_MAX;
779            else
780              largestGap = value2-value1;
781          }
782        }
783        bool tryIt= numberRows>200&&numberColumns>2000&&numberColumns>2*numberRows;
784        if (numberRows<1000&&numberColumns<3000)
785          tryIt=false;
786        if (notInteger)
787          tryIt=false;
788        if (largest/smallest>10||(largest/smallest>2.0&&largest>50))
789          tryIt=false;
790        if (tryIt) {
791          if (largest/smallest>2.0) {
792            nPasses = 10+numberColumns/100000;
793            nPasses = CoinMin(nPasses,50);
794            nPasses = CoinMax(nPasses,15);
795            if (numberRows>20000&&nPasses>5) {
796              // Might as well go for it
797              nPasses = CoinMax(nPasses,71);
798            } else if (numberRows>2000&&nPasses>5) {
799              nPasses = CoinMax(nPasses,50);
800            } else if (numberElements<3*numberColumns) {
801              nPasses=CoinMin(nPasses,10); // probably not worh it
802            }
803          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
804            nPasses = 10+numberColumns/1000;
805            nPasses = CoinMin(nPasses,100);
806            nPasses = CoinMax(nPasses,30);
807            if (numberRows>25000) {
808              // Might as well go for it
809              nPasses = CoinMax(nPasses,71);
810            }
811            if (!largestGap)
812              nPasses *= 2;
813          } else {
814            nPasses = 10+numberColumns/1000;
815            nPasses = CoinMin(nPasses,200);
816            nPasses = CoinMax(nPasses,100);
817            if (!largestGap)
818              nPasses *= 2;
819          }
820        }
821        //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
822        //     numberRows,numberColumns,plusMinus ? 'Y' : 'N',
823        //     tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
824        //exit(0);
825        if (!tryIt||nPasses<=5)
826          doIdiot=0;
827        if (doSprint) {
828          method = ClpSolve::usePrimalorSprint;
829        } else if (doIdiot) {
830          method = ClpSolve::usePrimal;
831        } else {
832          method = ClpSolve::useDual;
833        }
834      }
835    }
836  }
837  if (method==ClpSolve::usePrimalorSprint) {
838    if (doSprint<0) { 
839      if (numberElements<500000) {
840        // Small problem
841        if(numberRows*10>numberColumns||numberColumns<6000
842           ||(numberRows*20>numberColumns&&!plusMinus))
843          method=ClpSolve::usePrimal; // switch off sprint
844      } else {
845        // larger problem
846        if(numberRows*8>numberColumns)
847          method=ClpSolve::usePrimal; // switch off sprint
848        // but make lightweight
849        if(numberRows*10>numberColumns||numberColumns<6000
850           ||(numberRows*20>numberColumns&&!plusMinus))
851          maxSprintPass=10;
852      }
853    } else if (doSprint==0) {
854      method=ClpSolve::usePrimal; // switch off sprint
855    }
856  }
857  if (method==ClpSolve::useDual) {
858    double * saveLower=NULL;
859    double * saveUpper=NULL;
860    if (presolve==ClpSolve::presolveOn) {
861      int numberInfeasibilities = model2->tightenPrimalBounds(0.0,0);
862      if (numberInfeasibilities) {
863        handler_->message(CLP_INFEASIBLE,messages_)
864          <<CoinMessageEol;
865        model2 = this;
866        presolve=ClpSolve::presolveOff;
867      }
868    } else if (numberRows_+numberColumns_>5000) {
869      // do anyway
870      saveLower = new double[numberRows_+numberColumns_];
871      CoinMemcpyN(model2->columnLower(),numberColumns_,saveLower);
872      CoinMemcpyN(model2->rowLower(),numberRows_,saveLower+numberColumns_);
873      saveUpper = new double[numberRows_+numberColumns_];
874      CoinMemcpyN(model2->columnUpper(),numberColumns_,saveUpper);
875      CoinMemcpyN(model2->rowUpper(),numberRows_,saveUpper+numberColumns_);
876      int numberInfeasibilities = model2->tightenPrimalBounds();
877      if (numberInfeasibilities) {
878        handler_->message(CLP_INFEASIBLE,messages_)
879          <<CoinMessageEol;
880        CoinMemcpyN(saveLower,numberColumns_,model2->columnLower());
881        CoinMemcpyN(saveLower+numberColumns_,numberRows_,model2->rowLower());
882        delete [] saveLower;
883        saveLower=NULL;
884        CoinMemcpyN(saveUpper,numberColumns_,model2->columnUpper());
885        CoinMemcpyN(saveUpper+numberColumns_,numberRows_,model2->rowUpper());
886        delete [] saveUpper;
887        saveUpper=NULL;
888      }
889    }
890#ifndef COIN_HAS_VOL
891    // switch off idiot and volume for now
892    doIdiot=0; 
893#endif
894    // pick up number passes
895    int nPasses=0;
896    int numberNotE=0;
897#ifndef SLIM_CLP
898    if ((doIdiot<0&&plusMinus)||doIdiot>0) {
899      // See if candidate for idiot
900      nPasses=0;
901      Idiot info(*model2);
902      // Get average number of elements per column
903      double ratio  = ((double) numberElements/(double) numberColumns);
904      // look at rhs
905      int iRow;
906      double largest=0.0;
907      double smallest = 1.0e30;
908      double largestGap=0.0;
909      for (iRow=0;iRow<numberRows;iRow++) {
910        double value1 = model2->rowLower_[iRow];
911        if (value1&&value1>-1.0e31) {
912          largest = CoinMax(largest,fabs(value1));
913          smallest=CoinMin(smallest,fabs(value1));
914        }
915        double value2 = model2->rowUpper_[iRow];
916        if (value2&&value2<1.0e31) {
917          largest = CoinMax(largest,fabs(value2));
918          smallest=CoinMin(smallest,fabs(value2));
919        }
920        if (value2>value1) {
921          numberNotE++;
922          if (value2>1.0e31||value1<-1.0e31)
923            largestGap = COIN_DBL_MAX;
924          else
925            largestGap = value2-value1;
926        }
927      }
928      if (doIdiot<0) {
929        if (numberRows>200&&numberColumns>5000&&ratio>=3.0&&
930            largest/smallest<1.1&&!numberNotE) {
931          nPasses = 71;
932        }
933      } 
934      if (doIdiot>0) {
935        nPasses=CoinMax(nPasses,doIdiot);
936        if (nPasses>70) {
937          info.setStartingWeight(1.0e3);
938          info.setDropEnoughFeasibility(0.01);
939        }
940      }
941      if (nPasses>20) {
942#ifdef COIN_HAS_VOL
943        int returnCode = solveWithVolume(model2,nPasses,saveDoIdiot);
944        if (!returnCode) {
945          time2 = CoinCpuTime();
946          timeIdiot = time2-timeX;
947          handler_->message(CLP_INTERVAL_TIMING,messages_)
948            <<"Idiot Crash"<<timeIdiot<<time2-time1
949            <<CoinMessageEol;
950          timeX=time2;
951        } else {
952          nPasses=0;
953        }
954#else
955        nPasses=0;
956#endif
957      } else {
958        nPasses=0;
959      }
960    }
961#endif
962    if (doCrash) {
963      switch(doCrash) {
964        // standard
965      case 1:
966        model2->crash(1000,1);
967        break;
968        // As in paper by Solow and Halim (approx)
969      case 2:
970      case 3:
971        model2->crash(model2->dualBound(),0);
972        break;
973        // Just put free in basis
974      case 4:
975        model2->crash(0.0,3);
976        break;
977      }
978    }
979    if (!nPasses) {
980      int saveOptions = model2->specialOptions();
981      if (model2->numberRows()>100)
982        model2->setSpecialOptions(saveOptions|64); // go as far as possible
983      //int numberRows = model2->numberRows();
984      //int numberColumns = model2->numberColumns();
985      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
986        // See if original wanted vector
987        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
988        ClpMatrixBase * matrix = model2->clpMatrix();
989        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
990          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
991          clpMatrix->makeSpecialColumnCopy();
992          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
993          model2->dual(0);
994          clpMatrix->releaseSpecialColumnCopy();
995        } else {
996          model2->dual(0);
997        }
998      } else {
999        model2->dual(0);
1000      }
1001    } else if (!numberNotE&&0) {
1002      // E so we can do in another way
1003      double * pi = model2->dualRowSolution();
1004      int i;
1005      int numberColumns = model2->numberColumns();
1006      int numberRows = model2->numberRows();
1007      double * saveObj = new double[numberColumns];
1008      CoinMemcpyN(model2->objective(),numberColumns,saveObj);
1009      CoinMemcpyN(model2->objective(),
1010             numberColumns,model2->dualColumnSolution());
1011      model2->clpMatrix()->transposeTimes(-1.0,pi,model2->dualColumnSolution());
1012      CoinMemcpyN(model2->dualColumnSolution(),
1013             numberColumns,model2->objective());
1014      const double * rowsol = model2->primalRowSolution();
1015      double offset=0.0;
1016      for (i=0;i<numberRows;i++) {
1017        offset += pi[i]*rowsol[i];
1018      }
1019      double value2;
1020      model2->getDblParam(ClpObjOffset,value2);
1021      //printf("Offset %g %g\n",offset,value2);
1022      model2->setDblParam(ClpObjOffset,value2-offset);
1023      model2->setPerturbation(51);
1024      //model2->setRowObjective(pi);
1025      // zero out pi
1026      //memset(pi,0,numberRows*sizeof(double));
1027      // Could put some in basis - only partially tested
1028      model2->allSlackBasis(); 
1029      //model2->factorization()->maximumPivots(200);
1030      //model2->setLogLevel(63);
1031      // solve
1032      model2->dual(0);
1033      model2->setDblParam(ClpObjOffset,value2);
1034      CoinMemcpyN(saveObj,numberColumns,model2->objective());
1035      // zero out pi
1036      //memset(pi,0,numberRows*sizeof(double));
1037      //model2->setRowObjective(pi);
1038      delete [] saveObj;
1039      //model2->dual(0);
1040      model2->setPerturbation(50);
1041      model2->primal();
1042    } else {
1043      // solve
1044      model2->setPerturbation(100);
1045      model2->dual(2);
1046      model2->setPerturbation(50);
1047      model2->dual(0);
1048    }
1049    if (saveLower) {
1050      CoinMemcpyN(saveLower,numberColumns_,model2->columnLower());
1051      CoinMemcpyN(saveLower+numberColumns_,numberRows_,model2->rowLower());
1052      delete [] saveLower;
1053      saveLower=NULL;
1054      CoinMemcpyN(saveUpper,numberColumns_,model2->columnUpper());
1055      CoinMemcpyN(saveUpper+numberColumns_,numberRows_,model2->rowUpper());
1056      delete [] saveUpper;
1057      saveUpper=NULL;
1058    }
1059    time2 = CoinCpuTime();
1060    timeCore = time2-timeX;
1061    handler_->message(CLP_INTERVAL_TIMING,messages_)
1062      <<"Dual"<<timeCore<<time2-time1
1063      <<CoinMessageEol;
1064    timeX=time2;
1065  } else if (method==ClpSolve::usePrimal) {
1066#ifndef SLIM_CLP
1067    if (doIdiot) {
1068      int nPasses=0;
1069      Idiot info(*model2);
1070      // Get average number of elements per column
1071      double ratio  = ((double) numberElements/(double) numberColumns);
1072      // look at rhs
1073      int iRow;
1074      double largest=0.0;
1075      double smallest = 1.0e30;
1076      double largestGap=0.0;
1077      int numberNotE=0;
1078      for (iRow=0;iRow<numberRows;iRow++) {
1079        double value1 = model2->rowLower_[iRow];
1080        if (value1&&value1>-1.0e31) {
1081          largest = CoinMax(largest,fabs(value1));
1082          smallest=CoinMin(smallest,fabs(value1));
1083        }
1084        double value2 = model2->rowUpper_[iRow];
1085        if (value2&&value2<1.0e31) {
1086          largest = CoinMax(largest,fabs(value2));
1087          smallest=CoinMin(smallest,fabs(value2));
1088        }
1089        if (value2>value1) {
1090          numberNotE++;
1091          if (value2>1.0e31||value1<-1.0e31)
1092            largestGap = COIN_DBL_MAX;
1093          else
1094            largestGap = value2-value1;
1095        }
1096      }
1097      bool increaseSprint=plusMinus;
1098      if (!plusMinus) {
1099        // If 90% +- 1 then go for sprint
1100        if (statistics[0]>=0&&10*statistics[2]<statistics[0]+statistics[1])
1101          increaseSprint=true;
1102      }
1103      bool tryIt= numberRows>200&&numberColumns>2000&&numberColumns>2*numberRows;
1104      if (numberRows<1000&&numberColumns<3000)
1105        tryIt=false;
1106      if (tryIt) {
1107        if (increaseSprint) {
1108          info.setStartingWeight(1.0e3);
1109          info.setReduceIterations(6);
1110          // also be more lenient on infeasibilities
1111          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
1112          info.setDropEnoughWeighted(-2.0);
1113          if (largest/smallest>2.0) {
1114            nPasses = 10+numberColumns/100000;
1115            nPasses = CoinMin(nPasses,50);
1116            nPasses = CoinMax(nPasses,15);
1117            if (numberRows>20000&&nPasses>5) {
1118              // Might as well go for it
1119              nPasses = CoinMax(nPasses,71);
1120            } else if (numberRows>2000&&nPasses>5) {
1121              nPasses = CoinMax(nPasses,50);
1122            } else if (numberElements<3*numberColumns) {
1123              nPasses=CoinMin(nPasses,10); // probably not worh it
1124              if (doIdiot<0)
1125                info.setLightweight(1); // say lightweight idiot
1126            } else {
1127              if (doIdiot<0)
1128                info.setLightweight(1); // say lightweight idiot
1129            }
1130          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
1131            nPasses = 10+numberColumns/1000;
1132            nPasses = CoinMin(nPasses,100);
1133            nPasses = CoinMax(nPasses,30);
1134            if (numberRows>25000) {
1135              // Might as well go for it
1136              nPasses = CoinMax(nPasses,71);
1137            }
1138            if (!largestGap)
1139              nPasses *= 2;
1140          } else {
1141            nPasses = 10+numberColumns/1000;
1142            nPasses = CoinMin(nPasses,200);
1143            nPasses = CoinMax(nPasses,100);
1144            info.setStartingWeight(1.0e-1);
1145            info.setReduceIterations(6);
1146            if (!largestGap)
1147              nPasses *= 2;
1148            //info.setFeasibilityTolerance(1.0e-7);
1149          }
1150          // If few passes - don't bother
1151          if (nPasses<=5)
1152            nPasses=0;
1153        } else {
1154          if (doIdiot<0)
1155            info.setLightweight(1); // say lightweight idiot
1156          if (largest/smallest>1.01||numberNotE||statistics[2]>statistics[0]+statistics[1]) {
1157            if (numberRows>25000||numberColumns>5*numberRows) {
1158              nPasses = 50;
1159            } else if (numberColumns>4*numberRows) {
1160              nPasses = 20;
1161            } else {
1162              nPasses=5;
1163            }
1164          } else {
1165            if (numberRows>25000||numberColumns>5*numberRows) {
1166              nPasses = 50;
1167              info.setLightweight(0); // say not lightweight idiot
1168            } else if (numberColumns>4*numberRows) {
1169              nPasses = 20;
1170            } else {
1171              nPasses=15;
1172            }
1173          }
1174          if (ratio<3.0) { 
1175            nPasses=(int) ((ratio*(double) nPasses)/4.0); // probably not worh it
1176          } else {
1177            nPasses = CoinMax(nPasses,5);
1178          }
1179          if (numberRows>25000&&nPasses>5) {
1180            // Might as well go for it
1181            nPasses = CoinMax(nPasses,71);
1182          } else if (increaseSprint) {
1183            nPasses *= 2;
1184            nPasses=CoinMin(nPasses,71);
1185          } else if (nPasses==5&&ratio>5.0) {
1186            nPasses = (int) (((double) nPasses)*(ratio/5.0)); // increase if lots of elements per column
1187          }
1188          if (nPasses<=5)
1189            nPasses=0;
1190          //info.setStartingWeight(1.0e-1);
1191        }
1192      }
1193      if (doIdiot>0) {
1194        // pick up number passes
1195        nPasses=options.getExtraInfo(1);
1196        if (nPasses>70) {
1197          info.setStartingWeight(1.0e3);
1198          info.setReduceIterations(6);
1199          if (nPasses>=5000) {
1200            int k= nPasses&100;
1201            nPasses /= 100;
1202            info.setReduceIterations(3);
1203            if (k)
1204              info.setStartingWeight(1.0e2);
1205          }
1206          // also be more lenient on infeasibilities
1207          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
1208          info.setDropEnoughWeighted(-2.0);
1209        } else if (nPasses>=50) {
1210          info.setStartingWeight(1.0e3);
1211          //info.setReduceIterations(6);
1212        } 
1213        // For experimenting
1214        if (nPasses<70&&(nPasses%10)>0&&(nPasses%10)<4) {
1215          info.setStartingWeight(1.0e3);
1216          info.setLightweight(nPasses%10); // special testing
1217#ifdef COIN_DEVELOP
1218          printf("warning - odd lightweight %d\n",nPasses%10);
1219          //info.setReduceIterations(6);
1220#endif
1221        }
1222      }
1223      if (nPasses) {
1224        doCrash=0;
1225#if 0
1226        double * solution = model2->primalColumnSolution();
1227        int iColumn;
1228        double * saveLower = new double[numberColumns];
1229        CoinMemcpyN(model2->columnLower(),numberColumns,saveLower);
1230        double * saveUpper = new double[numberColumns];
1231        CoinMemcpyN(model2->columnUpper(),numberColumns,saveUpper);
1232        printf("doing tighten before idiot\n");
1233        model2->tightenPrimalBounds();
1234        // Move solution
1235        double * columnLower = model2->columnLower();
1236        double * columnUpper = model2->columnUpper();
1237        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1238          if (columnLower[iColumn]>0.0)
1239            solution[iColumn]=columnLower[iColumn];
1240          else if (columnUpper[iColumn]<0.0)
1241            solution[iColumn]=columnUpper[iColumn];
1242          else
1243            solution[iColumn]=0.0;
1244        }
1245        CoinMemcpyN(saveLower,numberColumns,columnLower);
1246        CoinMemcpyN(saveUpper,numberColumns,columnUpper);
1247        delete [] saveLower;
1248        delete [] saveUpper;
1249#else
1250        // Allow for crossover
1251        //if (doIdiot>0)
1252          info.setStrategy(512|info.getStrategy());
1253        // Allow for scaling
1254        info.setStrategy(32|info.getStrategy());
1255        info.crash(nPasses,model2->messageHandler(),model2->messagesPointer());
1256#endif
1257        time2 = CoinCpuTime();
1258        timeIdiot = time2-timeX;
1259        handler_->message(CLP_INTERVAL_TIMING,messages_)
1260          <<"Idiot Crash"<<timeIdiot<<time2-time1
1261          <<CoinMessageEol;
1262        timeX=time2;
1263      }
1264    }
1265#endif
1266    // ?
1267    if (doCrash) {
1268      switch(doCrash) {
1269        // standard
1270      case 1:
1271        model2->crash(1000,1);
1272        break;
1273        // As in paper by Solow and Halim (approx)
1274      case 2:
1275        model2->crash(model2->dualBound(),0);
1276        break;
1277        // My take on it
1278      case 3:
1279        model2->crash(model2->dualBound(),-1);
1280        break;
1281        // Just put free in basis
1282      case 4:
1283        model2->crash(0.0,3);
1284        break;
1285      }
1286    }
1287#ifndef SLIM_CLP
1288    if (doSlp>0&&objective_->type()==2) {
1289      model2->nonlinearSLP(doSlp,1.0e-5);
1290    }
1291#endif
1292    if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1293      // See if original wanted vector
1294      ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1295      ClpMatrixBase * matrix = model2->clpMatrix();
1296      if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
1297        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1298        clpMatrix->makeSpecialColumnCopy();
1299        //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
1300        model2->primal(primalStartup);
1301        clpMatrix->releaseSpecialColumnCopy();
1302      } else {
1303        model2->primal(primalStartup);
1304      }
1305    } else {
1306      model2->primal(primalStartup);
1307    }
1308    time2 = CoinCpuTime();
1309    timeCore = time2-timeX;
1310    handler_->message(CLP_INTERVAL_TIMING,messages_)
1311      <<"Primal"<<timeCore<<time2-time1
1312      <<CoinMessageEol;
1313    timeX=time2;
1314  } else if (method==ClpSolve::usePrimalorSprint) {
1315    // Sprint
1316    /*
1317      This driver implements what I called Sprint when I introduced the idea
1318      many years ago.  Cplex calls it "sifting" which I think is just as silly.
1319      When I thought of this trivial idea
1320      it reminded me of an LP code of the 60's called sprint which after
1321      every factorization took a subset of the matrix into memory (all
1322      64K words!) and then iterated very fast on that subset.  On the
1323      problems of those days it did not work very well, but it worked very
1324      well on aircrew scheduling problems where there were very large numbers
1325      of columns all with the same flavor.
1326    */
1327   
1328    /* The idea works best if you can get feasible easily.  To make it
1329       more general we can add in costed slacks */
1330   
1331    int originalNumberColumns = model2->numberColumns();
1332    int numberRows = model2->numberRows();
1333    ClpSimplex * originalModel2 = model2;
1334   
1335    // We will need arrays to choose variables.  These are too big but ..
1336    double * weight = new double [numberRows+originalNumberColumns];
1337    int * sort = new int [numberRows+originalNumberColumns];
1338    int numberSort=0;
1339    // We are going to add slacks to get feasible.
1340    // initial list will just be artificials
1341    int iColumn;
1342    const double * columnLower = model2->columnLower();
1343    const double * columnUpper = model2->columnUpper();
1344    double * columnSolution = model2->primalColumnSolution();
1345
1346    // See if we have costed slacks
1347    int * negSlack = new int[numberRows];
1348    int * posSlack = new int[numberRows];
1349    int iRow;
1350    for (iRow=0;iRow<numberRows;iRow++) {
1351      negSlack[iRow]=-1;
1352      posSlack[iRow]=-1;
1353    }
1354    const double * element = model2->matrix()->getElements();
1355    const int * row = model2->matrix()->getIndices();
1356    const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
1357    const int * columnLength = model2->matrix()->getVectorLengths();
1358    //bool allSlack = (numberRowsBasic==numberRows);
1359    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1360      if (!columnSolution[iColumn]||fabs(columnSolution[iColumn])>1.0e20) {
1361        double value =0.0;
1362        if (columnLower[iColumn]>0.0)
1363          value = columnLower[iColumn];
1364        else if (columnUpper[iColumn]<0.0)
1365          value = columnUpper[iColumn];
1366        columnSolution[iColumn]=value;
1367      }
1368      if (columnLength[iColumn]==1) {
1369        int jRow=row[columnStart[iColumn]];
1370        if (!columnLower[iColumn]) {
1371          if (element[columnStart[iColumn]]>0.0&&posSlack[jRow]<0)
1372            posSlack[jRow]=iColumn;
1373          else if (element[columnStart[iColumn]]<0.0&&negSlack[jRow]<0)
1374            negSlack[jRow]=iColumn;
1375        } else if (!columnUpper[iColumn]) {
1376          if (element[columnStart[iColumn]]<0.0&&posSlack[jRow]<0)
1377            posSlack[jRow]=iColumn;
1378          else if (element[columnStart[iColumn]]>0.0&&negSlack[jRow]<0)
1379            negSlack[jRow]=iColumn;
1380        }
1381      }
1382    }
1383    // now see what that does to row solution
1384    double * rowSolution = model2->primalRowSolution();
1385    CoinZeroN (rowSolution,numberRows);
1386    model2->clpMatrix()->times(1.0,columnSolution,rowSolution);
1387    // See if we can adjust using costed slacks
1388    double penalty=CoinMax(1.0e5,CoinMin(infeasibilityCost_*0.01,1.0e10))*optimizationDirection_;
1389    const double * lower = model2->rowLower();
1390    const double * upper = model2->rowUpper();
1391    for (iRow=0;iRow<numberRows;iRow++) {
1392      if (lower[iRow]>rowSolution[iRow]+1.0e-8) {
1393        int jColumn = posSlack[iRow];
1394        if (jColumn>=0) {
1395          if (columnSolution[jColumn])
1396            continue;
1397          double difference = lower[iRow]-rowSolution[iRow];
1398          double elementValue = element[columnStart[jColumn]];
1399          if (elementValue>0.0) {
1400            double movement = CoinMin(difference/elementValue,columnUpper[jColumn]);
1401            columnSolution[jColumn] = movement;
1402            rowSolution[iRow] += movement*elementValue;
1403          } else {
1404            double movement = CoinMax(difference/elementValue,columnLower[jColumn]);
1405            columnSolution[jColumn] = movement;
1406            rowSolution[iRow] += movement*elementValue;
1407          }
1408        }
1409      } else if (upper[iRow]<rowSolution[iRow]-1.0e-8) {
1410        int jColumn = negSlack[iRow];
1411        if (jColumn>=0) {
1412          if (columnSolution[jColumn])
1413            continue;
1414          double difference = upper[iRow]-rowSolution[iRow];
1415          double elementValue = element[columnStart[jColumn]];
1416          if (elementValue<0.0) {
1417            double movement = CoinMin(difference/elementValue,columnUpper[jColumn]);
1418            columnSolution[jColumn] = movement;
1419            rowSolution[iRow] += movement*elementValue;
1420          } else {
1421            double movement = CoinMax(difference/elementValue,columnLower[jColumn]);
1422            columnSolution[jColumn] = movement;
1423            rowSolution[iRow] += movement*elementValue;
1424          }
1425        }
1426      }
1427    }
1428    delete [] negSlack;
1429    delete [] posSlack;
1430    int nRow=numberRows;
1431    bool network=false;
1432    if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
1433      network=true;
1434      nRow *= 2;
1435    }
1436    int * addStarts = new int [nRow+1];
1437    int * addRow = new int[nRow];
1438    double * addElement = new double[nRow];
1439    addStarts[0]=0;
1440    int numberArtificials=0;
1441    int numberAdd=0;
1442    double * addCost = new double [numberRows];
1443    for (iRow=0;iRow<numberRows;iRow++) {
1444      if (lower[iRow]>rowSolution[iRow]+1.0e-8) {
1445        addRow[numberAdd]=iRow;
1446        addElement[numberAdd++]=1.0;
1447        if (network) {
1448          addRow[numberAdd]=numberRows;
1449          addElement[numberAdd++]=-1.0;
1450        }
1451        addCost[numberArtificials]=penalty;
1452        numberArtificials++;
1453        addStarts[numberArtificials]=numberAdd;
1454      } else if (upper[iRow]<rowSolution[iRow]-1.0e-8) {
1455        addRow[numberAdd]=iRow;
1456        addElement[numberAdd++]=-1.0;
1457        if (network) {
1458          addRow[numberAdd]=numberRows;
1459          addElement[numberAdd++]=1.0;
1460        }
1461        addCost[numberArtificials]=penalty;
1462        numberArtificials++;
1463        addStarts[numberArtificials]=numberAdd;
1464      }
1465    }
1466    if (numberArtificials) {
1467      // need copy so as not to disturb original
1468      model2 = new ClpSimplex(*model2);
1469      if (network) {
1470        // network - add a null row
1471        model2->addRow(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX);
1472        numberRows++;
1473      }
1474      model2->addColumns(numberArtificials,NULL,NULL,addCost,
1475                         addStarts,addRow,addElement);
1476    }
1477    delete [] addStarts;
1478    delete [] addRow;
1479    delete [] addElement;
1480    delete [] addCost;
1481    // look at rhs to see if to perturb
1482    double largest=0.0;
1483    double smallest = 1.0e30;
1484    for (iRow=0;iRow<numberRows;iRow++) {
1485      double value;
1486      value = fabs(model2->rowLower_[iRow]);
1487      if (value&&value<1.0e30) {
1488        largest = CoinMax(largest,value);
1489        smallest=CoinMin(smallest,value);
1490      }
1491      value = fabs(model2->rowUpper_[iRow]);
1492      if (value&&value<1.0e30) {
1493        largest = CoinMax(largest,value);
1494        smallest=CoinMin(smallest,value);
1495      }
1496    }
1497    double * saveLower = NULL;
1498    double * saveUpper = NULL;
1499    if (largest<2.01*smallest) {
1500      // perturb - so switch off standard
1501      model2->setPerturbation(100);
1502      saveLower = new double[numberRows];
1503      CoinMemcpyN(model2->rowLower_,numberRows,saveLower);
1504      saveUpper = new double[numberRows];
1505      CoinMemcpyN(model2->rowUpper_,numberRows,saveUpper);
1506      double * lower = model2->rowLower();
1507      double * upper = model2->rowUpper();
1508      for (iRow=0;iRow<numberRows;iRow++) {
1509        double lowerValue=lower[iRow], upperValue=upper[iRow];
1510        double value = randomNumberGenerator_.randomDouble();
1511        if (upperValue>lowerValue+primalTolerance_) {
1512          if (lowerValue>-1.0e20&&lowerValue)
1513            lowerValue -= value * 1.0e-4*fabs(lowerValue); 
1514          if (upperValue<1.0e20&&upperValue)
1515            upperValue += value * 1.0e-4*fabs(upperValue); 
1516        } else if (upperValue>0.0) {
1517          upperValue -= value * 1.0e-4*fabs(lowerValue); 
1518          lowerValue -= value * 1.0e-4*fabs(lowerValue); 
1519        } else if (upperValue<0.0) {
1520          upperValue += value * 1.0e-4*fabs(lowerValue); 
1521          lowerValue += value * 1.0e-4*fabs(lowerValue); 
1522        } else {
1523        }
1524        lower[iRow]=lowerValue;
1525        upper[iRow]=upperValue;
1526      }
1527    }
1528    int i;
1529    // Just do this number of passes in Sprint
1530    if (doSprint>0)
1531      maxSprintPass=options.getExtraInfo(1);
1532    // but if big use to get ratio
1533    double ratio=3;
1534    if (maxSprintPass>1000) {
1535      ratio = ((double) maxSprintPass)*0.0001;
1536      ratio = CoinMax(ratio,1.1);
1537      maxSprintPass= maxSprintPass %1000;
1538#ifdef COIN_DEVELOP
1539      printf("%d passes wanted with ratio of %g\n",maxSprintPass,ratio);
1540#endif
1541    }
1542    // Just take this number of columns in small problem
1543    int smallNumberColumns = (int) CoinMin(ratio*numberRows,(double) numberColumns);
1544    smallNumberColumns = CoinMax(smallNumberColumns,3000);
1545    smallNumberColumns = CoinMin(smallNumberColumns,numberColumns);
1546    //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
1547    //smallNumberColumns = CoinMax(smallNumberColumns,3000);
1548    //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
1549    // redo as may have changed
1550    columnLower = model2->columnLower();
1551    columnUpper = model2->columnUpper();
1552    columnSolution = model2->primalColumnSolution();
1553    // Set up initial list
1554    numberSort=0;
1555    if (numberArtificials) {
1556      numberSort=numberArtificials;
1557      for (i=0;i<numberSort;i++)
1558        sort[i] = i+originalNumberColumns;
1559    } 
1560    // maybe a solution there already
1561    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1562      if (model2->getColumnStatus(iColumn)==basic)
1563        sort[numberSort++]=iColumn;
1564    }
1565    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1566      if (model2->getColumnStatus(iColumn)!=basic) {
1567        if (columnSolution[iColumn]>columnLower[iColumn]&&
1568            columnSolution[iColumn]<columnUpper[iColumn]&&
1569            columnSolution[iColumn])
1570          sort[numberSort++]=iColumn;
1571      }
1572    }
1573    numberSort = CoinMin(numberSort,smallNumberColumns);
1574   
1575    int numberColumns = model2->numberColumns();
1576    double * fullSolution = model2->primalColumnSolution();
1577   
1578   
1579    int iPass;
1580    double lastObjective=1.0e31;
1581    // It will be safe to allow dense
1582    model2->setInitialDenseFactorization(true);
1583   
1584    // We will be using all rows
1585    int * whichRows = new int [numberRows];
1586    for (iRow=0;iRow<numberRows;iRow++)
1587      whichRows[iRow]=iRow;
1588    double originalOffset;
1589    model2->getDblParam(ClpObjOffset,originalOffset);
1590    int totalIterations=0;
1591    double lastSumArtificials=COIN_DBL_MAX;
1592    int originalMaxSprintPass=maxSprintPass;
1593    maxSprintPass=20; // so we do that many if infeasible
1594    for (iPass=0;iPass<maxSprintPass;iPass++) {
1595      //printf("Bug until submodel new version\n");
1596      //CoinSort_2(sort,sort+numberSort,weight);
1597      // Create small problem
1598      ClpSimplex small(model2,numberRows,whichRows,numberSort,sort);
1599      small.setPerturbation(model2->perturbation());
1600      small.setInfeasibilityCost(model2->infeasibilityCost());
1601      if (model2->factorizationFrequency()==200) {
1602        // User did not touch preset
1603        small.defaultFactorizationFrequency();
1604      }
1605      // now see what variables left out do to row solution
1606      double * rowSolution = model2->primalRowSolution();
1607      double * sumFixed = new double[numberRows];
1608      CoinZeroN (sumFixed,numberRows);
1609      int iRow,iColumn;
1610      // zero out ones in small problem
1611      for (iColumn=0;iColumn<numberSort;iColumn++) {
1612        int kColumn = sort[iColumn];
1613        fullSolution[kColumn]=0.0;
1614      }
1615      // Get objective offset
1616      double offset=0.0;
1617      const double * objective = model2->objective();
1618      for (iColumn=0;iColumn<numberColumns;iColumn++) 
1619        offset += fullSolution[iColumn]*objective[iColumn];
1620      small.setDblParam(ClpObjOffset,originalOffset-offset);
1621      model2->clpMatrix()->times(1.0,fullSolution,sumFixed);
1622     
1623      double * lower = small.rowLower();
1624      double * upper = small.rowUpper();
1625      for (iRow=0;iRow<numberRows;iRow++) {
1626        if (lower[iRow]>-1.0e50) 
1627          lower[iRow] -= sumFixed[iRow];
1628        if (upper[iRow]<1.0e50)
1629          upper[iRow] -= sumFixed[iRow];
1630        rowSolution[iRow] -= sumFixed[iRow];
1631      }
1632      delete [] sumFixed;
1633      // Solve
1634      if (interrupt)
1635        currentModel = &small;
1636      small.defaultFactorizationFrequency();
1637      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1638        // See if original wanted vector
1639        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1640        ClpMatrixBase * matrix = small.clpMatrix();
1641        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
1642          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1643          clpMatrix->makeSpecialColumnCopy();
1644          small.primal(1);
1645          clpMatrix->releaseSpecialColumnCopy();
1646        } else {
1647#if 1
1648          small.primal(1);
1649#else
1650          int numberColumns = small.numberColumns();
1651          int numberRows = small.numberRows();
1652          // Use dual region
1653          double * rhs = small.dualRowSolution();
1654          int * whichRow = new int[3*numberRows];
1655          int * whichColumn = new int[2*numberColumns];
1656          int nBound;
1657          ClpSimplex * small2 = ((ClpSimplexOther *) (&small))->crunch(rhs,whichRow,whichColumn,
1658                                                                        nBound,false,false);
1659          if (small2) {
1660            small2->primal(1);
1661            if (small2->problemStatus()==0) {
1662              small.setProblemStatus(0);
1663              ((ClpSimplexOther *) (&small))->afterCrunch(*small2,whichRow,whichColumn,nBound);
1664            } else {
1665              small2->primal(1);
1666              if (small2->problemStatus())
1667                small.primal(1);
1668            }
1669            delete small2;
1670          } else {
1671            small.primal(1);
1672          }
1673          delete [] whichRow;
1674          delete [] whichColumn;
1675#endif
1676        }
1677      } else {
1678        small.primal(1);
1679      }
1680      totalIterations += small.numberIterations();
1681      // move solution back
1682      const double * solution = small.primalColumnSolution();
1683      for (iColumn=0;iColumn<numberSort;iColumn++) {
1684        int kColumn = sort[iColumn];
1685        model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
1686        fullSolution[kColumn]=solution[iColumn];
1687      }
1688      for (iRow=0;iRow<numberRows;iRow++) 
1689        model2->setRowStatus(iRow,small.getRowStatus(iRow));
1690      CoinMemcpyN(small.primalRowSolution(),
1691             numberRows,model2->primalRowSolution());
1692      double sumArtificials = 0.0;
1693      for (i=0;i<numberArtificials;i++)
1694        sumArtificials += fullSolution[i + originalNumberColumns];
1695      if (sumArtificials&&iPass>5&&sumArtificials>=lastSumArtificials) {
1696        // increase costs
1697        double * cost = model2->objective()+originalNumberColumns;
1698        double newCost = CoinMin(1.0e10,cost[0]*1.5);
1699        for (i=0;i<numberArtificials;i++)
1700          cost[i]=newCost;
1701      }
1702      lastSumArtificials = sumArtificials;
1703      // get reduced cost for large problem
1704      double * djs = model2->dualColumnSolution();
1705      CoinMemcpyN(model2->objective(),numberColumns,djs);
1706      model2->clpMatrix()->transposeTimes(-1.0,small.dualRowSolution(),djs);
1707      int numberNegative=0;
1708      double sumNegative = 0.0;
1709      // now massage weight so all basic in plus good djs
1710      // first count and do basic
1711      numberSort=0;
1712      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1713        double dj = djs[iColumn]*optimizationDirection_;
1714        double value = fullSolution[iColumn];
1715        if (model2->getColumnStatus(iColumn)==ClpSimplex::basic) {
1716          sort[numberSort++] = iColumn;
1717        } else if (dj<-dualTolerance_&&value<columnUpper[iColumn]) {
1718          numberNegative++;
1719          sumNegative -= dj;
1720        } else if (dj>dualTolerance_&&value>columnLower[iColumn]) {
1721          numberNegative++;
1722          sumNegative += dj;
1723        }
1724      }
1725      handler_->message(CLP_SPRINT,messages_)
1726        <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative
1727        <<numberNegative
1728        <<CoinMessageEol;
1729      if (sumArtificials<1.0e-8&&originalMaxSprintPass>=0) {
1730        maxSprintPass = iPass+originalMaxSprintPass;
1731        originalMaxSprintPass=-1;
1732      }
1733      if (iPass>20)
1734        sumArtificials=0.0;
1735      if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5&&sumArtificials<1.0e-8)||
1736          (!small.numberIterations()&&iPass)||
1737          iPass==maxSprintPass-1||small.status()==3) {
1738       
1739        break; // finished
1740      } else {
1741        lastObjective = small.objectiveValue()*optimizationDirection_;
1742        double tolerance;
1743        double averageNegDj = sumNegative/((double) (numberNegative+1));
1744        if (numberNegative+numberSort>smallNumberColumns)
1745          tolerance = -dualTolerance_;
1746        else 
1747          tolerance = 10.0*averageNegDj;
1748        int saveN = numberSort;
1749        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1750          double dj = djs[iColumn]*optimizationDirection_;
1751          double value = fullSolution[iColumn];
1752          if (model2->getColumnStatus(iColumn)!=ClpSimplex::basic) {
1753            if (dj<-dualTolerance_&&value<columnUpper[iColumn])
1754              dj = dj;
1755            else if (dj>dualTolerance_&&value>columnLower[iColumn])
1756              dj = -dj;
1757            else if (columnUpper[iColumn]>columnLower[iColumn])
1758              dj = fabs(dj);
1759            else
1760              dj = 1.0e50;
1761            if (dj<tolerance) {
1762              weight[numberSort] = dj;
1763              sort[numberSort++] = iColumn;
1764            }
1765          }
1766        }
1767        // sort
1768        CoinSort_2(weight+saveN,weight+numberSort,sort+saveN);
1769        numberSort = CoinMin(smallNumberColumns,numberSort);
1770      }
1771    }
1772    if (interrupt) 
1773      currentModel = model2;
1774    for (i=0;i<numberArtificials;i++)
1775      sort[i] = i + originalNumberColumns;
1776    model2->deleteColumns(numberArtificials,sort);
1777    if (network) {
1778      int iRow=numberRows-1;
1779      model2->deleteRows(1,&iRow);
1780    }
1781    delete [] weight;
1782    delete [] sort;
1783    delete [] whichRows;
1784    if (saveLower) {
1785      // unperturb and clean
1786      for (iRow=0;iRow<numberRows;iRow++) {
1787        double diffLower = saveLower[iRow]-model2->rowLower_[iRow];
1788        double diffUpper = saveUpper[iRow]-model2->rowUpper_[iRow];
1789        model2->rowLower_[iRow]=saveLower[iRow];
1790        model2->rowUpper_[iRow]=saveUpper[iRow];
1791        if (diffLower) 
1792          assert (!diffUpper||fabs(diffLower-diffUpper)<1.0e-5);
1793        else
1794          diffLower = diffUpper;
1795        model2->rowActivity_[iRow] += diffLower;
1796      }
1797      delete [] saveLower;
1798      delete [] saveUpper;
1799    }
1800    model2->primal(1);
1801    if (model2!=originalModel2)
1802      originalModel2->moveInfo(*model2);
1803    model2->setPerturbation(savePerturbation);
1804    time2 = CoinCpuTime();
1805    timeCore = time2-timeX;
1806    handler_->message(CLP_INTERVAL_TIMING,messages_)
1807      <<"Sprint"<<timeCore<<time2-time1
1808      <<CoinMessageEol;
1809    timeX=time2;
1810    model2->setNumberIterations(model2->numberIterations()+totalIterations);
1811  } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
1812#ifndef SLIM_CLP
1813    //printf("***** experimental pretty crude barrier\n");
1814    //#define SAVEIT 2
1815#ifndef SAVEIT
1816#define BORROW
1817#endif
1818#ifdef BORROW
1819    ClpInterior barrier;
1820    barrier.borrowModel(*model2);
1821#else
1822    ClpInterior barrier(*model2);
1823#endif
1824    if (interrupt) 
1825      currentModel2 = &barrier;
1826    int barrierOptions = options.getSpecialOption(4);
1827    int aggressiveGamma=0;
1828    bool presolveInCrossover=false;
1829    bool scale=false;
1830    bool doKKT=false;
1831    if (barrierOptions&16) {
1832      barrierOptions &= ~16;
1833      doKKT=true;
1834    }
1835    if (barrierOptions&(32+64+128)) {
1836      aggressiveGamma=(barrierOptions&(32+64+128))>>5;
1837      barrierOptions &= ~(32+64+128);
1838    }
1839    if (barrierOptions&256) {
1840      barrierOptions &= ~256;
1841      presolveInCrossover=true;
1842    }
1843    if (barrierOptions&8) {
1844      barrierOptions &= ~8;
1845      scale=true;
1846    }
1847#ifdef COIN_DEVELOP
1848#ifndef FAST_BARRIER
1849    if (!numberBarrier)
1850      std::cout<<"Warning - the default ordering is just on row counts! "
1851               <<"The factorization is being improved"<<std::endl;
1852    numberBarrier++;
1853#endif
1854#endif
1855    // If quadratic force KKT
1856    if (quadraticObj) {
1857      doKKT=true;
1858    }
1859    switch (barrierOptions) {
1860    case 0:
1861    default:
1862      if (!doKKT) {
1863        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1864        barrier.setCholesky(cholesky);
1865      } else {
1866        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1867        cholesky->setKKT(true);
1868        barrier.setCholesky(cholesky);
1869      }
1870      break;
1871    case 1:
1872      if (!doKKT) {
1873        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1874        barrier.setCholesky(cholesky);
1875      } else {
1876        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1877        cholesky->setKKT(true);
1878        barrier.setCholesky(cholesky);
1879      }
1880      break;
1881#ifdef WSSMP_BARRIER
1882    case 2:
1883      {
1884        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
1885        barrier.setCholesky(cholesky);
1886        assert (!doKKT);
1887      }
1888      break;
1889    case 3:
1890      if (!doKKT) {
1891        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
1892        barrier.setCholesky(cholesky);
1893      } else {
1894        ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
1895        barrier.setCholesky(cholesky);
1896      }
1897      break;
1898#endif
1899#ifdef UFL_BARRIER
1900    case 4:
1901      if (!doKKT) {
1902        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1903        barrier.setCholesky(cholesky);
1904      } else {
1905        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1906        cholesky->setKKT(true);
1907        barrier.setCholesky(cholesky);
1908      }
1909      break;
1910#endif
1911#ifdef TAUCS_BARRIER
1912    case 5:
1913      {
1914        ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
1915        barrier.setCholesky(cholesky);
1916        assert (!doKKT);
1917      }
1918      break;
1919#endif
1920    }
1921    int numberRows = model2->numberRows();
1922    int numberColumns = model2->numberColumns();
1923    int saveMaxIts = model2->maximumIterations();
1924    if (saveMaxIts<1000) {
1925      barrier.setMaximumBarrierIterations(saveMaxIts);
1926      model2->setMaximumIterations(1000000);
1927    }
1928#ifndef SAVEIT
1929    //barrier.setDiagonalPerturbation(1.0e-25);
1930    if (aggressiveGamma) {
1931      switch (aggressiveGamma) {
1932      case 1:
1933        barrier.setGamma(1.0e-5);
1934        barrier.setDelta(1.0e-5);
1935        break;
1936      case 2:
1937        barrier.setGamma(1.0e-5);
1938        break;
1939      case 3:
1940        barrier.setDelta(1.0e-5);
1941        break;
1942      case 4:
1943        barrier.setGamma(1.0e-3);
1944        barrier.setDelta(1.0e-3);
1945        break;
1946      case 5:
1947        barrier.setGamma(1.0e-3);
1948        break;
1949      case 6:
1950        barrier.setDelta(1.0e-3);
1951        break;
1952      }
1953    }
1954    if (scale)
1955      barrier.scaling(1);
1956    else
1957      barrier.scaling(0);
1958    barrier.primalDual();
1959#elif SAVEIT==1
1960    barrier.primalDual();
1961#else
1962    model2->restoreModel("xx.save");
1963    // move solutions
1964    CoinMemcpyN(model2->primalRowSolution(),
1965                numberRows,barrier.primalRowSolution());
1966    CoinMemcpyN(model2->dualRowSolution(),
1967                numberRows,barrier.dualRowSolution());
1968    CoinMemcpyN(model2->primalColumnSolution(),
1969                numberColumns,barrier.primalColumnSolution());
1970    CoinMemcpyN(model2->dualColumnSolution(),
1971                numberColumns,barrier.dualColumnSolution());
1972#endif
1973    time2 = CoinCpuTime();
1974    timeCore = time2-timeX;
1975    handler_->message(CLP_INTERVAL_TIMING,messages_)
1976      <<"Barrier"<<timeCore<<time2-time1
1977      <<CoinMessageEol;
1978    timeX=time2;
1979    int maxIts = barrier.maximumBarrierIterations();
1980    int barrierStatus=barrier.status();
1981    double gap = barrier.complementarityGap();
1982    // get which variables are fixed
1983    double * saveLower=NULL;
1984    double * saveUpper=NULL;
1985    ClpPresolve pinfo2;
1986    ClpSimplex * saveModel2=NULL;
1987    bool extraPresolve=false;
1988    int numberFixed = barrier.numberFixed();
1989    if (numberFixed) {
1990      int numberRows = barrier.numberRows();
1991      int numberColumns = barrier.numberColumns();
1992      int numberTotal = numberRows+numberColumns;
1993      saveLower = new double [numberTotal];
1994      saveUpper = new double [numberTotal];
1995      CoinMemcpyN(barrier.columnLower(),numberColumns,saveLower);
1996      CoinMemcpyN(barrier.rowLower(),numberRows,saveLower+numberColumns);
1997      CoinMemcpyN(barrier.columnUpper(),numberColumns,saveUpper);
1998      CoinMemcpyN(barrier.rowUpper(),numberRows,saveUpper+numberColumns);
1999    }
2000    if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&
2001        presolveInCrossover) {
2002      // may as well do presolve
2003      barrier.fixFixed();
2004      saveModel2=model2;
2005      extraPresolve=true;
2006    } else if (numberFixed) {
2007      // Set fixed to bounds (may have restored earlier solution)
2008      barrier.fixFixed(false);
2009    }
2010#ifdef BORROW   
2011    barrier.returnModel(*model2);
2012    double * rowPrimal = new double [numberRows];
2013    double * columnPrimal = new double [numberColumns];
2014    double * rowDual = new double [numberRows];
2015    double * columnDual = new double [numberColumns];
2016    // move solutions other way
2017    CoinMemcpyN(model2->primalRowSolution(),
2018                numberRows,rowPrimal);
2019    CoinMemcpyN(model2->dualRowSolution(),
2020                numberRows,rowDual);
2021    CoinMemcpyN(model2->primalColumnSolution(),
2022                numberColumns,columnPrimal);
2023    CoinMemcpyN(model2->dualColumnSolution(),
2024                  numberColumns,columnDual);
2025#else
2026    double * rowPrimal = barrier.primalRowSolution();
2027    double * columnPrimal = barrier.primalColumnSolution();
2028    double * rowDual = barrier.dualRowSolution();
2029    double * columnDual = barrier.dualColumnSolution();
2030    // move solutions
2031    CoinMemcpyN(rowPrimal,
2032                numberRows,model2->primalRowSolution());
2033    CoinMemcpyN(rowDual,
2034                numberRows,model2->dualRowSolution());
2035    CoinMemcpyN(columnPrimal,
2036                numberColumns,model2->primalColumnSolution());
2037    CoinMemcpyN(columnDual,
2038                  numberColumns,model2->dualColumnSolution());
2039#endif
2040    if (saveModel2) {
2041      // do presolve
2042      model2 = pinfo2.presolvedModel(*model2,dblParam_[ClpPresolveTolerance],
2043                                    false,5,true);
2044      if (!model2) {
2045        model2=saveModel2;
2046        saveModel2=NULL;
2047        int numberRows = model2->numberRows();
2048        int numberColumns = model2->numberColumns();
2049        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2050        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2051        delete [] saveLower;
2052        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2053        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2054        delete [] saveUpper;
2055        saveLower=NULL;
2056        saveUpper=NULL;
2057      }
2058    }
2059    if (method==ClpSolve::useBarrier) {
2060      if (maxIts&&barrierStatus<4&&!quadraticObj) {
2061        //printf("***** crossover - needs more thought on difficult models\n");
2062#if SAVEIT==1
2063        model2->ClpSimplex::saveModel("xx.save");
2064#endif
2065        // make sure no status left
2066        model2->createStatus();
2067        // solve
2068        model2->setPerturbation(100);
2069        if (model2->factorizationFrequency()==200) {
2070          // User did not touch preset
2071          model2->defaultFactorizationFrequency();
2072        }
2073#if 1
2074        // throw some into basis
2075        {
2076          int numberRows = model2->numberRows();
2077          int numberColumns = model2->numberColumns();
2078          double * dsort = new double[numberColumns];
2079          int * sort = new int[numberColumns];
2080          int n=0;
2081          const double * columnLower = model2->columnLower();
2082          const double * columnUpper = model2->columnUpper();
2083          double * primalSolution = model2->primalColumnSolution();
2084          const double * dualSolution = model2->dualColumnSolution();
2085          double tolerance = 10.0*primalTolerance_;
2086          int i;
2087          for ( i=0;i<numberRows;i++) 
2088            model2->setRowStatus(i,superBasic);
2089          for ( i=0;i<numberColumns;i++) {
2090            double distance = CoinMin(columnUpper[i]-primalSolution[i],
2091                                      primalSolution[i]-columnLower[i]);
2092            if (distance>tolerance) {
2093              if (fabs(dualSolution[i])<1.0e-5)
2094                distance *= 100.0;
2095              dsort[n]=-distance;
2096              sort[n++]=i;
2097              model2->setStatus(i,superBasic);
2098            } else if (distance>primalTolerance_) {
2099              model2->setStatus(i,superBasic);
2100            } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
2101              model2->setStatus(i,atLowerBound);
2102              primalSolution[i]=columnLower[i];
2103            } else {
2104              model2->setStatus(i,atUpperBound);
2105              primalSolution[i]=columnUpper[i];
2106            }
2107          }
2108          CoinSort_2(dsort,dsort+n,sort);
2109          n = CoinMin(numberRows,n);
2110          for ( i=0;i<n;i++) {
2111            int iColumn = sort[i];
2112            model2->setStatus(iColumn,basic);
2113          }
2114          delete [] sort;
2115          delete [] dsort;
2116        }
2117        // model2->allSlackBasis();
2118        if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
2119          if (saveUpper) {
2120            int numberRows = model2->numberRows();
2121            int numberColumns = model2->numberColumns();
2122            CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2123            CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2124            delete [] saveLower;
2125            CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2126            CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2127            delete [] saveUpper;
2128            saveLower=NULL;
2129            saveUpper=NULL;
2130          }
2131          int numberRows = model2->numberRows();
2132          int numberColumns = model2->numberColumns();
2133          // just primal values pass
2134          double saveScale = model2->objectiveScale();
2135          model2->setObjectiveScale(1.0e-3);
2136          model2->primal(2);
2137          model2->setObjectiveScale(saveScale);
2138          // save primal solution and copy back dual
2139          CoinMemcpyN(model2->primalRowSolution(),
2140                      numberRows,rowPrimal);
2141          CoinMemcpyN(rowDual,
2142                      numberRows,model2->dualRowSolution());
2143          CoinMemcpyN(model2->primalColumnSolution(),
2144                      numberColumns,columnPrimal);
2145          CoinMemcpyN(columnDual,
2146                      numberColumns,model2->dualColumnSolution());
2147          //model2->primal(1);
2148          // clean up reduced costs and flag variables
2149          {
2150            double * dj = model2->dualColumnSolution();
2151            double * cost = model2->objective();
2152            double * saveCost = new double[numberColumns];
2153            CoinMemcpyN(cost,numberColumns,saveCost);
2154            double * saveLower = new double[numberColumns];
2155            double * lower = model2->columnLower();
2156            CoinMemcpyN(lower,numberColumns,saveLower);
2157            double * saveUpper = new double[numberColumns];
2158            double * upper = model2->columnUpper();
2159            CoinMemcpyN(upper,numberColumns,saveUpper);
2160            int i;
2161            double tolerance = 10.0*dualTolerance_;
2162            for ( i=0;i<numberColumns;i++) {
2163              if (model2->getStatus(i)==basic) {
2164                dj[i]=0.0;
2165              } else if (model2->getStatus(i)==atLowerBound) {
2166                if (optimizationDirection_*dj[i]<tolerance) {
2167                  if (optimizationDirection_*dj[i]<0.0) {
2168                    //if (dj[i]<-1.0e-3)
2169                    //printf("bad dj at lb %d %g\n",i,dj[i]);
2170                    cost[i] -= dj[i];
2171                    dj[i]=0.0;
2172                  }
2173                } else {
2174                  upper[i]=lower[i];
2175                }
2176              } else if (model2->getStatus(i)==atUpperBound) {
2177                if (optimizationDirection_*dj[i]>tolerance) {
2178                  if (optimizationDirection_*dj[i]>0.0) {
2179                    //if (dj[i]>1.0e-3)
2180                    //printf("bad dj at ub %d %g\n",i,dj[i]);
2181                    cost[i] -= dj[i];
2182                    dj[i]=0.0;
2183                  }
2184                } else {
2185                  lower[i]=upper[i];
2186                }
2187              }
2188            }
2189            // just dual values pass
2190            //model2->setLogLevel(63);
2191            //model2->setFactorizationFrequency(1);
2192            model2->dual(2);
2193            CoinMemcpyN(saveCost,numberColumns,cost);
2194            delete [] saveCost;
2195            CoinMemcpyN(saveLower,numberColumns,lower);
2196            delete [] saveLower;
2197            CoinMemcpyN(saveUpper,numberColumns,upper);
2198            delete [] saveUpper;
2199          }
2200          // and finish
2201          // move solutions
2202          CoinMemcpyN(rowPrimal,
2203                      numberRows,model2->primalRowSolution());
2204          CoinMemcpyN(columnPrimal,
2205                      numberColumns,model2->primalColumnSolution());
2206        }
2207        double saveScale = model2->objectiveScale();
2208        model2->setObjectiveScale(1.0e-3);
2209        model2->primal(2);
2210        model2->setObjectiveScale(saveScale);
2211        model2->primal(1);
2212#else
2213        // just primal
2214        model2->primal(1);
2215#endif
2216      } else if (barrierStatus==4) {
2217        // memory problems
2218        model2->setPerturbation(savePerturbation);
2219        model2->createStatus();
2220        model2->dual();
2221      } else if (maxIts&&quadraticObj) {
2222        // make sure no status left
2223        model2->createStatus();
2224        // solve
2225        model2->setPerturbation(100);
2226        model2->reducedGradient(1);
2227      }
2228    }
2229    model2->setMaximumIterations(saveMaxIts);
2230#ifdef BORROW
2231    delete [] rowPrimal;
2232    delete [] columnPrimal;
2233    delete [] rowDual;
2234    delete [] columnDual;
2235#endif
2236    if (extraPresolve) {
2237      pinfo2.postsolve(true);
2238      delete model2;
2239      model2=saveModel2;
2240    }
2241    if (saveUpper) {
2242      int numberRows = model2->numberRows();
2243      int numberColumns = model2->numberColumns();
2244      CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2245      CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2246      delete [] saveLower;
2247      CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2248      CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2249      delete [] saveUpper;
2250      saveLower=NULL;
2251      saveUpper=NULL;
2252      if (method!=ClpSolve::useBarrierNoCross) 
2253        model2->primal(1);
2254    }
2255    model2->setPerturbation(savePerturbation);
2256    time2 = CoinCpuTime();
2257    timeCore = time2-timeX;
2258    handler_->message(CLP_INTERVAL_TIMING,messages_)
2259      <<"Crossover"<<timeCore<<time2-time1
2260      <<CoinMessageEol;
2261    timeX=time2;
2262#else
2263    abort();
2264#endif
2265  } else {
2266    assert (method!=ClpSolve::automatic); // later
2267    time2=0.0;
2268  }
2269  if (saveMatrix) {
2270    if (model2==this) {
2271      // delete and replace
2272      delete model2->clpMatrix();
2273      model2->replaceMatrix(saveMatrix);
2274    } else {
2275      delete saveMatrix;
2276    }
2277  }
2278  numberIterations = model2->numberIterations();
2279  finalStatus=model2->status();
2280  int finalSecondaryStatus = model2->secondaryStatus();
2281  if (presolve==ClpSolve::presolveOn) {
2282    int saveLevel = logLevel();
2283    if ((specialOptions_&1024)==0)
2284      setLogLevel(CoinMin(1,saveLevel));
2285    else
2286      setLogLevel(CoinMin(0,saveLevel));
2287    pinfo.postsolve(true);
2288    factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
2289    time2 = CoinCpuTime();
2290    timePresolve += time2-timeX;
2291    handler_->message(CLP_INTERVAL_TIMING,messages_)
2292      <<"Postsolve"<<time2-timeX<<time2-time1
2293      <<CoinMessageEol;
2294    timeX=time2;
2295    if (!presolveToFile)
2296      delete model2;
2297    if (interrupt)
2298      currentModel = this;
2299    // checkSolution(); already done by postSolve
2300    setLogLevel(saveLevel);
2301    if (finalStatus!=3&&(finalStatus||status()==-1)) {
2302      int savePerturbation = perturbation();
2303      if (!finalStatus||(moreSpecialOptions_&2)==0) {
2304        if (finalStatus==2) {
2305          // unbounded - get feasible first
2306          double save = optimizationDirection_;
2307          optimizationDirection_=0.0;
2308          primal(1);
2309          optimizationDirection_=save;
2310          primal(1);
2311        } else if (finalStatus==1) {
2312          dual();
2313        } else {
2314          setPerturbation(100);
2315          primal(1);
2316        }
2317      } else {
2318        // just set status
2319        problemStatus_=finalStatus;
2320      }
2321      setPerturbation(savePerturbation);
2322      numberIterations += numberIterations_;
2323      numberIterations_ = numberIterations;
2324      finalStatus=status();
2325      time2 = CoinCpuTime();
2326      handler_->message(CLP_INTERVAL_TIMING,messages_)
2327      <<"Cleanup"<<time2-timeX<<time2-time1
2328      <<CoinMessageEol;
2329      timeX=time2;
2330    } else {
2331      secondaryStatus_=finalSecondaryStatus;
2332    }
2333  } else if (model2!=this) {
2334    // not presolved - but different model used (sprint probably)
2335    CoinMemcpyN(model2->primalRowSolution(),
2336                numberRows_,this->primalRowSolution());
2337    CoinMemcpyN(model2->dualRowSolution(),
2338                numberRows_,this->dualRowSolution());
2339    CoinMemcpyN(model2->primalColumnSolution(),
2340                numberColumns_,this->primalColumnSolution());
2341    CoinMemcpyN(model2->dualColumnSolution(),
2342                numberColumns_,this->dualColumnSolution());
2343    CoinMemcpyN(model2->statusArray(),
2344                numberColumns_+numberRows_,this->statusArray());
2345    objectiveValue_=model2->objectiveValue_;
2346    numberIterations_ =model2->numberIterations_;
2347    problemStatus_ =model2->problemStatus_;
2348    secondaryStatus_ =model2->secondaryStatus_;
2349    delete model2;
2350  }
2351  setMaximumIterations(saveMaxIterations);
2352  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
2353                               "Errors","User stopped"};
2354  assert (finalStatus>=-1&&finalStatus<=5);
2355  handler_->message(CLP_TIMING,messages_)
2356    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
2357  handler_->printing(presolve==ClpSolve::presolveOn)
2358    <<timePresolve;
2359  handler_->printing(timeIdiot!=0.0)
2360    <<timeIdiot;
2361  handler_->message()<<CoinMessageEol;
2362  if (interrupt) 
2363    signal(SIGINT,saveSignal);
2364  perturbation_=savePerturbation;
2365  scalingFlag_=saveScaling;
2366  // If faking objective - put back correct one
2367  if (savedObjective) {
2368    delete objective_;
2369    objective_=savedObjective;
2370  }
2371  return finalStatus;
2372}
2373// General solve
2374int 
2375ClpSimplex::initialSolve()
2376{
2377  // Default so use dual
2378  ClpSolve options;
2379  return initialSolve(options);
2380}
2381// General dual solve
2382int 
2383ClpSimplex::initialDualSolve()
2384{
2385  ClpSolve options;
2386  // Use dual
2387  options.setSolveType(ClpSolve::useDual);
2388  return initialSolve(options);
2389}
2390// General dual solve
2391int 
2392ClpSimplex::initialPrimalSolve()
2393{
2394  ClpSolve options;
2395  // Use primal
2396  options.setSolveType(ClpSolve::usePrimal);
2397  return initialSolve(options);
2398}
2399// barrier solve, not to be followed by crossover
2400int 
2401ClpSimplex::initialBarrierNoCrossSolve()
2402{
2403  ClpSolve options;
2404  // Use primal
2405  options.setSolveType(ClpSolve::useBarrierNoCross);
2406  return initialSolve(options);
2407}
2408
2409// General barrier solve
2410int 
2411ClpSimplex::initialBarrierSolve()
2412{
2413  ClpSolve options;
2414  // Use primal
2415  options.setSolveType(ClpSolve::useBarrier);
2416  return initialSolve(options);
2417}
2418
2419// Default constructor
2420ClpSolve::ClpSolve (  )
2421{
2422  method_ = automatic;
2423  presolveType_=presolveOn;
2424  numberPasses_=5;
2425  int i;
2426  for (i=0;i<7;i++)
2427    options_[i]=0;
2428  // say no +-1 matrix
2429  options_[3]=1;
2430  for (i=0;i<7;i++)
2431    extraInfo_[i]=-1;
2432  independentOptions_[0]=0;
2433  // But switch off slacks
2434  independentOptions_[1]=512;
2435  // Substitute up to 3
2436  independentOptions_[2]=3;
2437 
2438}
2439// Constructor when you really know what you are doing
2440ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
2441             int numberPasses, int options[6],
2442             int extraInfo[6], int independentOptions[3])
2443{
2444  method_ = method;
2445  presolveType_=presolveType;
2446  numberPasses_=numberPasses;
2447  int i;
2448  for (i=0;i<6;i++)
2449    options_[i]=options[i];
2450  options_[6]=0;
2451  for (i=0;i<6;i++)
2452    extraInfo_[i]=extraInfo[i];
2453  extraInfo_[6]=0;
2454  for (i=0;i<3;i++)
2455    independentOptions_[i]=independentOptions[i];
2456}
2457
2458// Copy constructor.
2459ClpSolve::ClpSolve(const ClpSolve & rhs)
2460{
2461  method_ = rhs.method_;
2462  presolveType_=rhs.presolveType_;
2463  numberPasses_=rhs.numberPasses_;
2464  int i;
2465  for ( i=0;i<7;i++)
2466    options_[i]=rhs.options_[i];
2467  for ( i=0;i<7;i++)
2468    extraInfo_[i]=rhs.extraInfo_[i];
2469  for (i=0;i<3;i++)
2470    independentOptions_[i]=rhs.independentOptions_[i];
2471}
2472// Assignment operator. This copies the data
2473ClpSolve & 
2474ClpSolve::operator=(const ClpSolve & rhs)
2475{
2476  if (this != &rhs) {
2477    method_ = rhs.method_;
2478    presolveType_=rhs.presolveType_;
2479    numberPasses_=rhs.numberPasses_;
2480    int i;
2481    for (i=0;i<7;i++)
2482      options_[i]=rhs.options_[i];
2483    for (i=0;i<7;i++)
2484      extraInfo_[i]=rhs.extraInfo_[i];
2485    for (i=0;i<3;i++)
2486      independentOptions_[i]=rhs.independentOptions_[i];
2487  }
2488  return *this;
2489
2490}
2491// Destructor
2492ClpSolve::~ClpSolve (  )
2493{
2494}
2495// See header file for deatils
2496void 
2497ClpSolve::setSpecialOption(int which,int value,int extraInfo)
2498{
2499  options_[which]=value;
2500  extraInfo_[which]=extraInfo;
2501}
2502int 
2503ClpSolve::getSpecialOption(int which) const
2504{
2505  return options_[which];
2506}
2507
2508// Solve types
2509void 
2510ClpSolve::setSolveType(SolveType method, int extraInfo)
2511{
2512  method_=method;
2513}
2514
2515ClpSolve::SolveType
2516ClpSolve::getSolveType()
2517{
2518  return method_;
2519}
2520
2521// Presolve types
2522void 
2523ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
2524{
2525  presolveType_=amount;
2526  numberPasses_=extraInfo;
2527}
2528ClpSolve::PresolveType
2529ClpSolve::getPresolveType()
2530{
2531  return presolveType_;
2532}
2533// Extra info for idiot (or sprint)
2534int 
2535ClpSolve::getExtraInfo(int which) const
2536{
2537  return extraInfo_[which];
2538}
2539int 
2540ClpSolve::getPresolvePasses() const
2541{
2542  return numberPasses_;
2543}
2544/* Say to return at once if infeasible,
2545   default is to solve */
2546void 
2547ClpSolve::setInfeasibleReturn(bool trueFalse)
2548{
2549  independentOptions_[0]= trueFalse ? 1 : 0;
2550}
2551#include <string>
2552// Generates code for above constructor
2553void 
2554ClpSolve::generateCpp(FILE * fp)
2555{
2556  std::string solveType[] = {
2557    "ClpSolve::useDual",
2558    "ClpSolve::usePrimal",
2559    "ClpSolve::usePrimalorSprint",
2560    "ClpSolve::useBarrier",
2561    "ClpSolve::useBarrierNoCross",
2562    "ClpSolve::automatic",
2563    "ClpSolve::notImplemented"
2564  };
2565  std::string presolveType[] =  {
2566    "ClpSolve::presolveOn",
2567    "ClpSolve::presolveOff",
2568    "ClpSolve::presolveNumber",
2569    "ClpSolve::presolveNumberCost"
2570  };
2571  fprintf(fp,"3  ClpSolve::SolveType method = %s;\n",solveType[method_].c_str());
2572  fprintf(fp,"3  ClpSolve::PresolveType presolveType = %s;\n",
2573    presolveType[presolveType_].c_str());
2574  fprintf(fp,"3  int numberPasses = %d;\n",numberPasses_);
2575  fprintf(fp,"3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
2576    options_[0],options_[1],options_[2],
2577    options_[3],options_[4],options_[5]);
2578  fprintf(fp,"3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
2579    extraInfo_[0],extraInfo_[1],extraInfo_[2],
2580    extraInfo_[3],extraInfo_[4],extraInfo_[5]);
2581  fprintf(fp,"3  int independentOptions[] = {%d,%d,%d};\n",
2582    independentOptions_[0],independentOptions_[1],independentOptions_[2]);
2583  fprintf(fp,"3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
2584  fprintf(fp,"3                    options,extraInfo,independentOptions);\n");
2585}
Note: See TracBrowser for help on using the repository browser.