source: tags/move-to-subversion/ClpSolve.cpp @ 753

Last change on this file since 753 was 719, checked in by forrest, 14 years ago

minor stuff

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