source: branches/BSP/trunk/Clp/src/ClpSolve.cpp @ 1078

Last change on this file since 1078 was 1078, checked in by jpfasano, 12 years ago

changed to remove MS compiler message
d:\coin\coinall-trunk\clp\src\clpsolve.cpp(1718) : warning C4244: '=' : conversion from 'double' to 'float', possible loss of data

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 76.5 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+CoinDrand48()*1.0e-2;
403        else
404          obj[i]=-1.0-CoinDrand48()*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      if (options.infeasibleReturn()) {
618        return -1;
619      }
620      presolve=ClpSolve::presolveOff;
621    }
622    // We may be better off using original (but if dual leave because of bounds)
623    if (presolve!=ClpSolve::presolveOff&&
624        numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_
625        &&model2!=this) {
626      if(method!=ClpSolve::useDual||
627         (numberRows_==model2->numberRows_&&numberColumns_==model2->numberColumns_)) {
628        delete model2;
629        model2 = this;
630        presolve=ClpSolve::presolveOff;
631      }
632    }
633  }
634  if (interrupt)
635    currentModel = model2;
636  // See if worth trying +- one matrix
637  bool plusMinus=false;
638  int numberElements=model2->getNumElements();
639#ifndef SLIM_CLP
640#ifndef NO_RTTI
641  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
642    // network - switch off stuff
643    doIdiot=0;
644    doSprint=0;
645  }
646#else
647  if (matrix_->type()==11) {
648    // network - switch off stuff
649    doIdiot=0;
650    doSprint=0;
651  }
652#endif
653#endif
654  int numberColumns = model2->numberColumns();
655  int numberRows = model2->numberRows();
656  // If not all slack basis - switch off all except sprint
657  int numberRowsBasic=0;
658  int iRow;
659  for (iRow=0;iRow<numberRows;iRow++)
660    if (model2->getRowStatus(iRow)==basic)
661      numberRowsBasic++;
662  if (numberRowsBasic<numberRows) {
663    doIdiot=0;
664    doCrash=0;
665    //doSprint=0;
666  }
667  if (options.getSpecialOption(3)==0) {
668    if(numberElements>100000)
669      plusMinus=true;
670    if(numberElements>10000&&(doIdiot||doSprint)) 
671      plusMinus=true;
672  }
673#ifndef SLIM_CLP
674  // Statistics (+1,-1, other) - used to decide on strategy if not +-1
675  CoinBigIndex statistics[3]={-1,0,0};
676  if (plusMinus) {
677    saveMatrix = model2->clpMatrix();
678#ifndef NO_RTTI
679    ClpPackedMatrix* clpMatrix =
680      dynamic_cast< ClpPackedMatrix*>(saveMatrix);
681#else
682    ClpPackedMatrix* clpMatrix = NULL;
683    if (saveMatrix->type()==1)
684      clpMatrix =
685        static_cast< ClpPackedMatrix*>(saveMatrix);
686#endif
687    if (clpMatrix) {
688      ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
689      if (newMatrix->getIndices()) {
690        model2->replaceMatrix(newMatrix);
691      } else {
692        handler_->message(CLP_MATRIX_CHANGE,messages_)
693          <<"+- 1"
694          <<CoinMessageEol;
695        CoinMemcpyN(newMatrix->startPositive(),3,statistics);
696        saveMatrix=NULL;
697        plusMinus=false;
698        delete newMatrix;
699      }
700    } else {
701      saveMatrix=NULL;
702      plusMinus=false;
703    }
704  }
705#endif
706  if (this->factorizationFrequency()==200) {
707    // User did not touch preset
708    model2->defaultFactorizationFrequency();
709  } else if (model2!=this) {
710    // make sure model2 has correct value
711    model2->setFactorizationFrequency(this->factorizationFrequency());
712  }
713  if (method==ClpSolve::automatic) {
714    if (doSprint==0&&doIdiot==0) {
715      // off
716      method=ClpSolve::useDual;
717    } else {
718      // only do primal if sprint or idiot
719      if (doSprint>0) {
720        method=ClpSolve::usePrimalorSprint;
721      } else if (doIdiot>0) {
722        method=ClpSolve::usePrimal;
723      } else {
724        if (numberElements<500000) {
725          // Small problem
726          if(numberRows*10>numberColumns||numberColumns<6000
727             ||(numberRows*20>numberColumns&&!plusMinus))
728            doSprint=0; // switch off sprint
729        } else {
730          // larger problem
731          if(numberRows*8>numberColumns)
732            doSprint=0; // switch off sprint
733        }
734        // switch off sprint or idiot if any free variable
735        int iColumn;
736        double * columnLower = model2->columnLower();
737        double * columnUpper = model2->columnUpper();
738        for (iColumn=0;iColumn<numberColumns;iColumn++) {
739          if (columnLower[iColumn]<-1.0e10&&columnUpper[iColumn]>1.0e10) {
740            doSprint=0;
741            doIdiot=0;
742            break;
743          }
744        }
745        int nPasses=0;
746        // look at rhs
747        int iRow;
748        double largest=0.0;
749        double smallest = 1.0e30;
750        double largestGap=0.0;
751        int numberNotE=0;
752        bool notInteger=false;
753        for (iRow=0;iRow<numberRows;iRow++) {
754          double value1 = model2->rowLower_[iRow];
755          if (value1&&value1>-1.0e31) {
756            largest = CoinMax(largest,fabs(value1));
757            smallest=CoinMin(smallest,fabs(value1));
758            if (fabs(value1-floor(value1+0.5))>1.0e-8) {
759              notInteger=true;
760              break;
761            }
762          }
763          double value2 = model2->rowUpper_[iRow];
764          if (value2&&value2<1.0e31) {
765            largest = CoinMax(largest,fabs(value2));
766            smallest=CoinMin(smallest,fabs(value2));
767            if (fabs(value2-floor(value2+0.5))>1.0e-8) {
768              notInteger=true;
769              break;
770            }
771          }
772          if (value2>value1) {
773            numberNotE++;
774            if (value2>1.0e31||value1<-1.0e31)
775              largestGap = COIN_DBL_MAX;
776            else
777              largestGap = value2-value1;
778          }
779        }
780        bool tryIt= numberRows>200&&numberColumns>2000&&numberColumns>2*numberRows;
781        if (numberRows<1000&&numberColumns<3000)
782          tryIt=false;
783        if (notInteger)
784          tryIt=false;
785        if (largest/smallest>10||(largest/smallest>2.0&&largest>50))
786          tryIt=false;
787        if (tryIt) {
788          if (largest/smallest>2.0) {
789            nPasses = 10+numberColumns/100000;
790            nPasses = CoinMin(nPasses,50);
791            nPasses = CoinMax(nPasses,15);
792            if (numberRows>20000&&nPasses>5) {
793              // Might as well go for it
794              nPasses = CoinMax(nPasses,71);
795            } else if (numberRows>2000&&nPasses>5) {
796              nPasses = CoinMax(nPasses,50);
797            } else if (numberElements<3*numberColumns) {
798              nPasses=CoinMin(nPasses,10); // probably not worh it
799            }
800          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
801            nPasses = 10+numberColumns/1000;
802            nPasses = CoinMin(nPasses,100);
803            nPasses = CoinMax(nPasses,30);
804            if (numberRows>25000) {
805              // Might as well go for it
806              nPasses = CoinMax(nPasses,71);
807            }
808            if (!largestGap)
809              nPasses *= 2;
810          } else {
811            nPasses = 10+numberColumns/1000;
812            nPasses = CoinMin(nPasses,200);
813            nPasses = CoinMax(nPasses,100);
814            if (!largestGap)
815              nPasses *= 2;
816          }
817        }
818        //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
819        //     numberRows,numberColumns,plusMinus ? 'Y' : 'N',
820        //     tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
821        //exit(0);
822        if (!tryIt||nPasses<=5)
823          doIdiot=0;
824        if (doSprint) {
825          method = ClpSolve::usePrimalorSprint;
826        } else if (doIdiot) {
827          method = ClpSolve::usePrimal;
828        } else {
829          method = ClpSolve::useDual;
830        }
831      }
832    }
833  }
834  if (method==ClpSolve::usePrimalorSprint) {
835    if (doSprint<0) { 
836      if (numberElements<500000) {
837        // Small problem
838        if(numberRows*10>numberColumns||numberColumns<6000
839           ||(numberRows*20>numberColumns&&!plusMinus))
840          method=ClpSolve::usePrimal; // switch off sprint
841      } else {
842        // larger problem
843        if(numberRows*8>numberColumns)
844          method=ClpSolve::usePrimal; // switch off sprint
845        // but make lightweight
846        if(numberRows*10>numberColumns||numberColumns<6000
847           ||(numberRows*20>numberColumns&&!plusMinus))
848          maxSprintPass=5;
849      }
850    } else if (doSprint==0) {
851      method=ClpSolve::usePrimal; // switch off sprint
852    }
853  }
854  if (method==ClpSolve::useDual) {
855    double * saveLower=NULL;
856    double * saveUpper=NULL;
857    if (presolve==ClpSolve::presolveOn) {
858      int numberInfeasibilities = model2->tightenPrimalBounds(0.0,0);
859      if (numberInfeasibilities) {
860        handler_->message(CLP_INFEASIBLE,messages_)
861          <<CoinMessageEol;
862        model2 = this;
863        presolve=ClpSolve::presolveOff;
864      }
865    } else if (numberRows_+numberColumns_>5000) {
866      // do anyway
867      saveLower = new double[numberRows_+numberColumns_];
868      CoinMemcpyN(model2->columnLower(),numberColumns_,saveLower);
869      CoinMemcpyN(model2->rowLower(),numberRows_,saveLower+numberColumns_);
870      saveUpper = new double[numberRows_+numberColumns_];
871      CoinMemcpyN(model2->columnUpper(),numberColumns_,saveUpper);
872      CoinMemcpyN(model2->rowUpper(),numberRows_,saveUpper+numberColumns_);
873      int numberInfeasibilities = model2->tightenPrimalBounds();
874      if (numberInfeasibilities) {
875        handler_->message(CLP_INFEASIBLE,messages_)
876          <<CoinMessageEol;
877        CoinMemcpyN(saveLower,numberColumns_,model2->columnLower());
878        CoinMemcpyN(saveLower+numberColumns_,numberRows_,model2->rowLower());
879        delete [] saveLower;
880        saveLower=NULL;
881        CoinMemcpyN(saveUpper,numberColumns_,model2->columnUpper());
882        CoinMemcpyN(saveUpper+numberColumns_,numberRows_,model2->rowUpper());
883        delete [] saveUpper;
884        saveUpper=NULL;
885      }
886    }
887#ifndef COIN_HAS_VOL
888    // switch off idiot and volume for now
889    doIdiot=0; 
890#endif
891    // pick up number passes
892    int nPasses=0;
893    int numberNotE=0;
894#ifndef SLIM_CLP
895    if ((doIdiot<0&&plusMinus)||doIdiot>0) {
896      // See if candidate for idiot
897      nPasses=0;
898      Idiot info(*model2);
899      // Get average number of elements per column
900      double ratio  = ((double) numberElements/(double) numberColumns);
901      // look at rhs
902      int iRow;
903      double largest=0.0;
904      double smallest = 1.0e30;
905      double largestGap=0.0;
906      for (iRow=0;iRow<numberRows;iRow++) {
907        double value1 = model2->rowLower_[iRow];
908        if (value1&&value1>-1.0e31) {
909          largest = CoinMax(largest,fabs(value1));
910          smallest=CoinMin(smallest,fabs(value1));
911        }
912        double value2 = model2->rowUpper_[iRow];
913        if (value2&&value2<1.0e31) {
914          largest = CoinMax(largest,fabs(value2));
915          smallest=CoinMin(smallest,fabs(value2));
916        }
917        if (value2>value1) {
918          numberNotE++;
919          if (value2>1.0e31||value1<-1.0e31)
920            largestGap = COIN_DBL_MAX;
921          else
922            largestGap = value2-value1;
923        }
924      }
925      if (doIdiot<0) {
926        if (numberRows>200&&numberColumns>5000&&ratio>=3.0&&
927            largest/smallest<1.1&&!numberNotE) {
928          nPasses = 71;
929        }
930      } 
931      if (doIdiot>0) {
932        nPasses=CoinMax(nPasses,doIdiot);
933        if (nPasses>70) {
934          info.setStartingWeight(1.0e3);
935          info.setDropEnoughFeasibility(0.01);
936        }
937      }
938      if (nPasses>20) {
939#ifdef COIN_HAS_VOL
940        int returnCode = solveWithVolume(model2,nPasses,saveDoIdiot);
941        if (!returnCode) {
942          time2 = CoinCpuTime();
943          timeIdiot = time2-timeX;
944          handler_->message(CLP_INTERVAL_TIMING,messages_)
945            <<"Idiot Crash"<<timeIdiot<<time2-time1
946            <<CoinMessageEol;
947          timeX=time2;
948        } else {
949          nPasses=0;
950        }
951#else
952        nPasses=0;
953#endif
954      } else {
955        nPasses=0;
956      }
957    }
958#endif
959    if (doCrash) {
960      switch(doCrash) {
961        // standard
962      case 1:
963        model2->crash(1000,1);
964        break;
965        // As in paper by Solow and Halim (approx)
966      case 2:
967      case 3:
968        model2->crash(model2->dualBound(),0);
969        break;
970        // Just put free in basis
971      case 4:
972        model2->crash(0.0,3);
973        break;
974      }
975    }
976    if (!nPasses) {
977      int saveOptions = model2->specialOptions();
978      if (model2->numberRows()>100)
979        model2->setSpecialOptions(saveOptions|64); // go as far as possible
980      //int numberRows = model2->numberRows();
981      //int numberColumns = model2->numberColumns();
982      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
983        // See if original wanted vector
984        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
985        ClpMatrixBase * matrix = model2->clpMatrix();
986        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
987          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
988          clpMatrix->makeSpecialColumnCopy();
989          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
990          model2->dual(0);
991          clpMatrix->releaseSpecialColumnCopy();
992        } else {
993          model2->dual(0);
994        }
995      } else {
996        model2->dual(0);
997      }
998    } else if (!numberNotE&&0) {
999      // E so we can do in another way
1000      double * pi = model2->dualRowSolution();
1001      int i;
1002      int numberColumns = model2->numberColumns();
1003      int numberRows = model2->numberRows();
1004      double * saveObj = new double[numberColumns];
1005      CoinMemcpyN(model2->objective(),numberColumns,saveObj);
1006      CoinMemcpyN(model2->objective(),
1007             numberColumns,model2->dualColumnSolution());
1008      model2->clpMatrix()->transposeTimes(-1.0,pi,model2->dualColumnSolution());
1009      CoinMemcpyN(model2->dualColumnSolution(),
1010             numberColumns,model2->objective());
1011      const double * rowsol = model2->primalRowSolution();
1012      double offset=0.0;
1013      for (i=0;i<numberRows;i++) {
1014        offset += pi[i]*rowsol[i];
1015      }
1016      double value2;
1017      model2->getDblParam(ClpObjOffset,value2);
1018      //printf("Offset %g %g\n",offset,value2);
1019      model2->setDblParam(ClpObjOffset,value2-offset);
1020      model2->setPerturbation(51);
1021      //model2->setRowObjective(pi);
1022      // zero out pi
1023      //memset(pi,0,numberRows*sizeof(double));
1024      // Could put some in basis - only partially tested
1025      model2->allSlackBasis(); 
1026      //model2->factorization()->maximumPivots(200);
1027      //model2->setLogLevel(63);
1028      // solve
1029      model2->dual(0);
1030      model2->setDblParam(ClpObjOffset,value2);
1031      CoinMemcpyN(saveObj,numberColumns,model2->objective());
1032      // zero out pi
1033      //memset(pi,0,numberRows*sizeof(double));
1034      //model2->setRowObjective(pi);
1035      delete [] saveObj;
1036      //model2->dual(0);
1037      model2->setPerturbation(50);
1038      model2->primal();
1039    } else {
1040      // solve
1041      model2->setPerturbation(100);
1042      model2->dual(2);
1043      model2->setPerturbation(50);
1044      model2->dual(0);
1045    }
1046    if (saveLower) {
1047      CoinMemcpyN(saveLower,numberColumns_,model2->columnLower());
1048      CoinMemcpyN(saveLower+numberColumns_,numberRows_,model2->rowLower());
1049      delete [] saveLower;
1050      saveLower=NULL;
1051      CoinMemcpyN(saveUpper,numberColumns_,model2->columnUpper());
1052      CoinMemcpyN(saveUpper+numberColumns_,numberRows_,model2->rowUpper());
1053      delete [] saveUpper;
1054      saveUpper=NULL;
1055    }
1056    time2 = CoinCpuTime();
1057    timeCore = time2-timeX;
1058    handler_->message(CLP_INTERVAL_TIMING,messages_)
1059      <<"Dual"<<timeCore<<time2-time1
1060      <<CoinMessageEol;
1061    timeX=time2;
1062  } else if (method==ClpSolve::usePrimal) {
1063#ifndef SLIM_CLP
1064    if (doIdiot) {
1065      int nPasses=0;
1066      Idiot info(*model2);
1067      // Get average number of elements per column
1068      double ratio  = ((double) numberElements/(double) numberColumns);
1069      // look at rhs
1070      int iRow;
1071      double largest=0.0;
1072      double smallest = 1.0e30;
1073      double largestGap=0.0;
1074      int numberNotE=0;
1075      for (iRow=0;iRow<numberRows;iRow++) {
1076        double value1 = model2->rowLower_[iRow];
1077        if (value1&&value1>-1.0e31) {
1078          largest = CoinMax(largest,fabs(value1));
1079          smallest=CoinMin(smallest,fabs(value1));
1080        }
1081        double value2 = model2->rowUpper_[iRow];
1082        if (value2&&value2<1.0e31) {
1083          largest = CoinMax(largest,fabs(value2));
1084          smallest=CoinMin(smallest,fabs(value2));
1085        }
1086        if (value2>value1) {
1087          numberNotE++;
1088          if (value2>1.0e31||value1<-1.0e31)
1089            largestGap = COIN_DBL_MAX;
1090          else
1091            largestGap = value2-value1;
1092        }
1093      }
1094      bool increaseSprint=plusMinus;
1095      if (!plusMinus) {
1096        // If 90% +- 1 then go for sprint
1097        if (statistics[0]>=0&&10*statistics[2]<statistics[0]+statistics[1])
1098          increaseSprint=true;
1099      }
1100      bool tryIt= numberRows>200&&numberColumns>2000&&numberColumns>2*numberRows;
1101      if (numberRows<1000&&numberColumns<3000)
1102        tryIt=false;
1103      if (tryIt) {
1104        if (increaseSprint) {
1105          info.setStartingWeight(1.0e3);
1106          info.setReduceIterations(6);
1107          // also be more lenient on infeasibilities
1108          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
1109          info.setDropEnoughWeighted(-2.0);
1110          if (largest/smallest>2.0) {
1111            nPasses = 10+numberColumns/100000;
1112            nPasses = CoinMin(nPasses,50);
1113            nPasses = CoinMax(nPasses,15);
1114            if (numberRows>20000&&nPasses>5) {
1115              // Might as well go for it
1116              nPasses = CoinMax(nPasses,71);
1117            } else if (numberRows>2000&&nPasses>5) {
1118              nPasses = CoinMax(nPasses,50);
1119            } else if (numberElements<3*numberColumns) {
1120              nPasses=CoinMin(nPasses,10); // probably not worh it
1121              if (doIdiot<0)
1122                info.setLightweight(1); // say lightweight idiot
1123            } else {
1124              if (doIdiot<0)
1125                info.setLightweight(1); // say lightweight idiot
1126            }
1127          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
1128            nPasses = 10+numberColumns/1000;
1129            nPasses = CoinMin(nPasses,100);
1130            nPasses = CoinMax(nPasses,30);
1131            if (numberRows>25000) {
1132              // Might as well go for it
1133              nPasses = CoinMax(nPasses,71);
1134            }
1135            if (!largestGap)
1136              nPasses *= 2;
1137          } else {
1138            nPasses = 10+numberColumns/1000;
1139            nPasses = CoinMin(nPasses,200);
1140            nPasses = CoinMax(nPasses,100);
1141            info.setStartingWeight(1.0e-1);
1142            info.setReduceIterations(6);
1143            if (!largestGap)
1144              nPasses *= 2;
1145            //info.setFeasibilityTolerance(1.0e-7);
1146          }
1147          // If few passes - don't bother
1148          if (nPasses<=5)
1149            nPasses=0;
1150        } else {
1151          if (doIdiot<0)
1152            info.setLightweight(1); // say lightweight idiot
1153          if (largest/smallest>1.01||numberNotE||statistics[2]>statistics[0]+statistics[1]) {
1154            if (numberRows>25000||numberColumns>5*numberRows) {
1155              nPasses = 50;
1156            } else if (numberColumns>4*numberRows) {
1157              nPasses = 20;
1158            } else {
1159              nPasses=5;
1160            }
1161          } else {
1162            if (numberRows>25000||numberColumns>5*numberRows) {
1163              nPasses = 50;
1164              info.setLightweight(0); // say not lightweight idiot
1165            } else if (numberColumns>4*numberRows) {
1166              nPasses = 20;
1167            } else {
1168              nPasses=15;
1169            }
1170          }
1171          if (ratio<3.0) { 
1172            nPasses=(int) ((ratio*(double) nPasses)/4.0); // probably not worh it
1173          } else {
1174            nPasses = CoinMax(nPasses,5);
1175          }
1176          if (numberRows>25000&&nPasses>5) {
1177            // Might as well go for it
1178            nPasses = CoinMax(nPasses,71);
1179          } else if (increaseSprint) {
1180            nPasses *= 2;
1181            nPasses=CoinMin(nPasses,71);
1182          } else if (nPasses==5&&ratio>5.0) {
1183            nPasses = (int) (((double) nPasses)*(ratio/5.0)); // increase if lots of elements per column
1184          }
1185          if (nPasses<=5)
1186            nPasses=0;
1187          //info.setStartingWeight(1.0e-1);
1188        }
1189      }
1190      if (doIdiot>0) {
1191        // pick up number passes
1192        nPasses=options.getExtraInfo(1);
1193        if (nPasses>70) {
1194          info.setStartingWeight(1.0e3);
1195          info.setReduceIterations(6);
1196          if (nPasses>=5000) {
1197            int k= nPasses&100;
1198            nPasses /= 100;
1199            info.setReduceIterations(3);
1200            if (k)
1201              info.setStartingWeight(1.0e2);
1202          }
1203          // also be more lenient on infeasibilities
1204          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
1205          info.setDropEnoughWeighted(-2.0);
1206        } else if (nPasses>=50) {
1207          info.setStartingWeight(1.0e3);
1208          //info.setReduceIterations(6);
1209        } 
1210        // For experimenting
1211        if (nPasses<70&&(nPasses%10)>0&&(nPasses%10)<4) {
1212          info.setStartingWeight(1.0e3);
1213          info.setLightweight(nPasses%10); // special testing
1214#ifdef COIN_DEVELOP
1215          printf("warning - odd lightweight %d\n",nPasses%10);
1216          //info.setReduceIterations(6);
1217#endif
1218        }
1219      }
1220      if (nPasses) {
1221        doCrash=0;
1222#if 0
1223        double * solution = model2->primalColumnSolution();
1224        int iColumn;
1225        double * saveLower = new double[numberColumns];
1226        CoinMemcpyN(model2->columnLower(),numberColumns,saveLower);
1227        double * saveUpper = new double[numberColumns];
1228        CoinMemcpyN(model2->columnUpper(),numberColumns,saveUpper);
1229        printf("doing tighten before idiot\n");
1230        model2->tightenPrimalBounds();
1231        // Move solution
1232        double * columnLower = model2->columnLower();
1233        double * columnUpper = model2->columnUpper();
1234        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1235          if (columnLower[iColumn]>0.0)
1236            solution[iColumn]=columnLower[iColumn];
1237          else if (columnUpper[iColumn]<0.0)
1238            solution[iColumn]=columnUpper[iColumn];
1239          else
1240            solution[iColumn]=0.0;
1241        }
1242        CoinMemcpyN(saveLower,numberColumns,columnLower);
1243        CoinMemcpyN(saveUpper,numberColumns,columnUpper);
1244        delete [] saveLower;
1245        delete [] saveUpper;
1246#else
1247        // Allow for crossover
1248        //if (doIdiot>0)
1249          info.setStrategy(512|info.getStrategy());
1250        // Allow for scaling
1251        info.setStrategy(32|info.getStrategy());
1252        info.crash(nPasses,model2->messageHandler(),model2->messagesPointer());
1253#endif
1254        time2 = CoinCpuTime();
1255        timeIdiot = time2-timeX;
1256        handler_->message(CLP_INTERVAL_TIMING,messages_)
1257          <<"Idiot Crash"<<timeIdiot<<time2-time1
1258          <<CoinMessageEol;
1259        timeX=time2;
1260      }
1261    }
1262#endif
1263    // ?
1264    if (doCrash) {
1265      switch(doCrash) {
1266        // standard
1267      case 1:
1268        model2->crash(1000,1);
1269        break;
1270        // As in paper by Solow and Halim (approx)
1271      case 2:
1272        model2->crash(model2->dualBound(),0);
1273        break;
1274        // My take on it
1275      case 3:
1276        model2->crash(model2->dualBound(),-1);
1277        break;
1278        // Just put free in basis
1279      case 4:
1280        model2->crash(0.0,3);
1281        break;
1282      }
1283    }
1284#ifndef SLIM_CLP
1285    if (doSlp>0&&objective_->type()==2) {
1286      model2->nonlinearSLP(doSlp,1.0e-5);
1287    }
1288#endif
1289    if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1290      // See if original wanted vector
1291      ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1292      ClpMatrixBase * matrix = model2->clpMatrix();
1293      if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
1294        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1295        clpMatrix->makeSpecialColumnCopy();
1296        //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
1297        model2->primal(primalStartup);
1298        clpMatrix->releaseSpecialColumnCopy();
1299      } else {
1300        model2->primal(primalStartup);
1301      }
1302    } else {
1303      model2->primal(primalStartup);
1304    }
1305    time2 = CoinCpuTime();
1306    timeCore = time2-timeX;
1307    handler_->message(CLP_INTERVAL_TIMING,messages_)
1308      <<"Primal"<<timeCore<<time2-time1
1309      <<CoinMessageEol;
1310    timeX=time2;
1311  } else if (method==ClpSolve::usePrimalorSprint) {
1312    // Sprint
1313    /*
1314      This driver implements what I called Sprint when I introduced the idea
1315      many years ago.  Cplex calls it "sifting" which I think is just as silly.
1316      When I thought of this trivial idea
1317      it reminded me of an LP code of the 60's called sprint which after
1318      every factorization took a subset of the matrix into memory (all
1319      64K words!) and then iterated very fast on that subset.  On the
1320      problems of those days it did not work very well, but it worked very
1321      well on aircrew scheduling problems where there were very large numbers
1322      of columns all with the same flavor.
1323    */
1324   
1325    /* The idea works best if you can get feasible easily.  To make it
1326       more general we can add in costed slacks */
1327   
1328    int originalNumberColumns = model2->numberColumns();
1329    int numberRows = model2->numberRows();
1330   
1331    // We will need arrays to choose variables.  These are too big but ..
1332    double * weight = new double [numberRows+originalNumberColumns];
1333    int * sort = new int [numberRows+originalNumberColumns];
1334    int numberSort=0;
1335    // We are going to add slacks to get feasible.
1336    // initial list will just be artificials
1337    int iColumn;
1338    const double * columnLower = model2->columnLower();
1339    const double * columnUpper = model2->columnUpper();
1340    double * columnSolution = model2->primalColumnSolution();
1341
1342    // See if we have costed slacks
1343    int * negSlack = new int[numberRows];
1344    int * posSlack = new int[numberRows];
1345    int iRow;
1346    for (iRow=0;iRow<numberRows;iRow++) {
1347      negSlack[iRow]=-1;
1348      posSlack[iRow]=-1;
1349    }
1350    const double * element = model2->matrix()->getElements();
1351    const int * row = model2->matrix()->getIndices();
1352    const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
1353    const int * columnLength = model2->matrix()->getVectorLengths();
1354    //bool allSlack = (numberRowsBasic==numberRows);
1355    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1356      if (!columnSolution[iColumn]||fabs(columnSolution[iColumn])>1.0e20) {
1357        double value =0.0;
1358        if (columnLower[iColumn]>0.0)
1359          value = columnLower[iColumn];
1360        else if (columnUpper[iColumn]<0.0)
1361          value = columnUpper[iColumn];
1362        columnSolution[iColumn]=value;
1363      }
1364      if (columnLength[iColumn]==1) {
1365        int jRow=row[columnStart[iColumn]];
1366        if (!columnLower[iColumn]) {
1367          if (element[columnStart[iColumn]]>0.0&&posSlack[jRow]<0)
1368            posSlack[jRow]=iColumn;
1369          else if (element[columnStart[iColumn]]<0.0&&negSlack[jRow]<0)
1370            negSlack[jRow]=iColumn;
1371        } else if (!columnUpper[iColumn]) {
1372          if (element[columnStart[iColumn]]<0.0&&posSlack[jRow]<0)
1373            posSlack[jRow]=iColumn;
1374          else if (element[columnStart[iColumn]]>0.0&&negSlack[jRow]<0)
1375            negSlack[jRow]=iColumn;
1376        }
1377      }
1378    }
1379    // now see what that does to row solution
1380    double * rowSolution = model2->primalRowSolution();
1381    CoinZeroN (rowSolution,numberRows);
1382    model2->clpMatrix()->times(1.0,columnSolution,rowSolution);
1383    // See if we can adjust using costed slacks
1384    double penalty=CoinMin(infeasibilityCost_,1.0e10)*optimizationDirection_;
1385    const double * lower = model2->rowLower();
1386    const double * upper = model2->rowUpper();
1387    for (iRow=0;iRow<numberRows;iRow++) {
1388      if (lower[iRow]>rowSolution[iRow]+1.0e-8) {
1389        int jColumn = posSlack[iRow];
1390        if (jColumn>=0) {
1391          if (columnSolution[jColumn])
1392            continue;
1393          double difference = lower[iRow]-rowSolution[iRow];
1394          double elementValue = element[columnStart[jColumn]];
1395          if (elementValue>0.0) {
1396            double movement = CoinMin(difference/elementValue,columnUpper[jColumn]);
1397            columnSolution[jColumn] = movement;
1398            rowSolution[iRow] += movement*elementValue;
1399          } else {
1400            double movement = CoinMax(difference/elementValue,columnLower[jColumn]);
1401            columnSolution[jColumn] = movement;
1402            rowSolution[iRow] += movement*elementValue;
1403          }
1404        }
1405      } else if (upper[iRow]<rowSolution[iRow]-1.0e-8) {
1406        int jColumn = negSlack[iRow];
1407        if (jColumn>=0) {
1408          if (columnSolution[jColumn])
1409            continue;
1410          double difference = upper[iRow]-rowSolution[iRow];
1411          double elementValue = element[columnStart[jColumn]];
1412          if (elementValue<0.0) {
1413            double movement = CoinMin(difference/elementValue,columnUpper[jColumn]);
1414            columnSolution[jColumn] = movement;
1415            rowSolution[iRow] += movement*elementValue;
1416          } else {
1417            double movement = CoinMax(difference/elementValue,columnLower[jColumn]);
1418            columnSolution[jColumn] = movement;
1419            rowSolution[iRow] += movement*elementValue;
1420          }
1421        }
1422      }
1423    }
1424    delete [] negSlack;
1425    delete [] posSlack;
1426    int * addStarts = new int [numberRows+1];
1427    int * addRow = new int[numberRows];
1428    double * addElement = new double[numberRows];
1429    addStarts[0]=0;
1430    int numberArtificials=0;
1431    double * addCost = new double [numberRows];
1432    for (iRow=0;iRow<numberRows;iRow++) {
1433      if (lower[iRow]>rowSolution[iRow]+1.0e-8) {
1434        addRow[numberArtificials]=iRow;
1435        addElement[numberArtificials]=1.0;
1436        addCost[numberArtificials]=penalty;
1437        numberArtificials++;
1438        addStarts[numberArtificials]=numberArtificials;
1439      } else if (upper[iRow]<rowSolution[iRow]-1.0e-8) {
1440        addRow[numberArtificials]=iRow;
1441        addElement[numberArtificials]=-1.0;
1442        addCost[numberArtificials]=penalty;
1443        numberArtificials++;
1444        addStarts[numberArtificials]=numberArtificials;
1445      }
1446    }
1447    if (numberArtificials) {
1448      // need copy so as not to disturb original
1449      model2 = new ClpSimplex(*model2);
1450    }
1451    model2->addColumns(numberArtificials,NULL,NULL,addCost,
1452                       addStarts,addRow,addElement);
1453    delete [] addStarts;
1454    delete [] addRow;
1455    delete [] addElement;
1456    delete [] addCost;
1457    // look at rhs to see if to perturb
1458    double largest=0.0;
1459    double smallest = 1.0e30;
1460    for (iRow=0;iRow<numberRows;iRow++) {
1461      double value;
1462      value = fabs(model2->rowLower_[iRow]);
1463      if (value&&value<1.0e30) {
1464        largest = CoinMax(largest,value);
1465        smallest=CoinMin(smallest,value);
1466      }
1467      value = fabs(model2->rowUpper_[iRow]);
1468      if (value&&value<1.0e30) {
1469        largest = CoinMax(largest,value);
1470        smallest=CoinMin(smallest,value);
1471      }
1472    }
1473    double * saveLower = NULL;
1474    double * saveUpper = NULL;
1475    if (largest<2.01*smallest) {
1476      // perturb - so switch off standard
1477      model2->setPerturbation(100);
1478      saveLower = new double[numberRows];
1479      CoinMemcpyN(model2->rowLower_,numberRows,saveLower);
1480      saveUpper = new double[numberRows];
1481      CoinMemcpyN(model2->rowUpper_,numberRows,saveUpper);
1482      double * lower = model2->rowLower();
1483      double * upper = model2->rowUpper();
1484      for (iRow=0;iRow<numberRows;iRow++) {
1485        double lowerValue=lower[iRow], upperValue=upper[iRow];
1486        double value = CoinDrand48();
1487        if (upperValue>lowerValue+primalTolerance_) {
1488          if (lowerValue>-1.0e20&&lowerValue)
1489            lowerValue -= value * 1.0e-4*fabs(lowerValue); 
1490          if (upperValue<1.0e20&&upperValue)
1491            upperValue += value * 1.0e-4*fabs(upperValue); 
1492        } else if (upperValue>0.0) {
1493          upperValue -= value * 1.0e-4*fabs(lowerValue); 
1494          lowerValue -= value * 1.0e-4*fabs(lowerValue); 
1495        } else if (upperValue<0.0) {
1496          upperValue += value * 1.0e-4*fabs(lowerValue); 
1497          lowerValue += value * 1.0e-4*fabs(lowerValue); 
1498        } else {
1499        }
1500        lower[iRow]=lowerValue;
1501        upper[iRow]=upperValue;
1502      }
1503    }
1504    int i;
1505    // Just do this number of passes in Sprint
1506    if (doSprint>0)
1507      maxSprintPass=options.getExtraInfo(1);
1508    // but if big use to get ratio
1509    double ratio=3;
1510    if (maxSprintPass>1000) {
1511      ratio = ((double) maxSprintPass)*0.0001;
1512      ratio = CoinMax(ratio,1.1);
1513      maxSprintPass= maxSprintPass %1000;
1514#ifdef COIN_DEVELOP
1515      printf("%d passes wanted with ratio of %g\n",maxSprintPass,ratio);
1516#endif
1517    }
1518    // Just take this number of columns in small problem
1519    int smallNumberColumns = (int) CoinMin(ratio*numberRows,(double) numberColumns);
1520    smallNumberColumns = CoinMax(smallNumberColumns,3000);
1521    smallNumberColumns = CoinMin(smallNumberColumns,numberColumns);
1522    //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
1523    //smallNumberColumns = CoinMax(smallNumberColumns,3000);
1524    //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
1525    // redo as may have changed
1526    columnLower = model2->columnLower();
1527    columnUpper = model2->columnUpper();
1528    columnSolution = model2->primalColumnSolution();
1529    // Set up initial list
1530    numberSort=0;
1531    if (numberArtificials) {
1532      numberSort=numberArtificials;
1533      for (i=0;i<numberSort;i++)
1534        sort[i] = i+originalNumberColumns;
1535    } 
1536    // maybe a solution there already
1537    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1538      if (model2->getColumnStatus(iColumn)==basic)
1539        sort[numberSort++]=iColumn;
1540    }
1541    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1542      if (model2->getColumnStatus(iColumn)!=basic) {
1543        if (columnSolution[iColumn]>columnLower[iColumn]&&
1544            columnSolution[iColumn]<columnUpper[iColumn]&&
1545            columnSolution[iColumn])
1546          sort[numberSort++]=iColumn;
1547      }
1548    }
1549    numberSort = CoinMin(numberSort,smallNumberColumns);
1550   
1551    int numberColumns = model2->numberColumns();
1552    double * fullSolution = model2->primalColumnSolution();
1553   
1554   
1555    int iPass;
1556    double lastObjective=1.0e31;
1557    // It will be safe to allow dense
1558    model2->setInitialDenseFactorization(true);
1559   
1560    // We will be using all rows
1561    int * whichRows = new int [numberRows];
1562    for (iRow=0;iRow<numberRows;iRow++)
1563      whichRows[iRow]=iRow;
1564    double originalOffset;
1565    model2->getDblParam(ClpObjOffset,originalOffset);
1566    int totalIterations=0;
1567    for (iPass=0;iPass<maxSprintPass;iPass++) {
1568      //printf("Bug until submodel new version\n");
1569      //CoinSort_2(sort,sort+numberSort,weight);
1570      // Create small problem
1571      ClpSimplex small(model2,numberRows,whichRows,numberSort,sort);
1572      small.setPerturbation(model2->perturbation());
1573      small.setInfeasibilityCost(model2->infeasibilityCost());
1574      if (model2->factorizationFrequency()==200) {
1575        // User did not touch preset
1576        small.defaultFactorizationFrequency();
1577      }
1578      // now see what variables left out do to row solution
1579      double * rowSolution = model2->primalRowSolution();
1580      double * sumFixed = new double[numberRows];
1581      CoinZeroN (sumFixed,numberRows);
1582      int iRow,iColumn;
1583      // zero out ones in small problem
1584      for (iColumn=0;iColumn<numberSort;iColumn++) {
1585        int kColumn = sort[iColumn];
1586        fullSolution[kColumn]=0.0;
1587      }
1588      // Get objective offset
1589      double offset=0.0;
1590      const double * objective = model2->objective();
1591      for (iColumn=0;iColumn<numberColumns;iColumn++) 
1592        offset += fullSolution[iColumn]*objective[iColumn];
1593      small.setDblParam(ClpObjOffset,originalOffset-offset);
1594      model2->clpMatrix()->times(1.0,fullSolution,sumFixed);
1595     
1596      double * lower = small.rowLower();
1597      double * upper = small.rowUpper();
1598      for (iRow=0;iRow<numberRows;iRow++) {
1599        if (lower[iRow]>-1.0e50) 
1600          lower[iRow] -= sumFixed[iRow];
1601        if (upper[iRow]<1.0e50)
1602          upper[iRow] -= sumFixed[iRow];
1603        rowSolution[iRow] -= sumFixed[iRow];
1604      }
1605      delete [] sumFixed;
1606      // Solve
1607      if (interrupt)
1608        currentModel = &small;
1609      small.defaultFactorizationFrequency();
1610      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1611        // See if original wanted vector
1612        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1613        ClpMatrixBase * matrix = small.clpMatrix();
1614        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
1615          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1616          clpMatrix->makeSpecialColumnCopy();
1617          small.primal(1);
1618          clpMatrix->releaseSpecialColumnCopy();
1619        } else {
1620#if 1
1621          small.primal(1);
1622#else
1623          int numberColumns = small.numberColumns();
1624          int numberRows = small.numberRows();
1625          // Use dual region
1626          double * rhs = small.dualRowSolution();
1627          int * whichRow = new int[3*numberRows];
1628          int * whichColumn = new int[2*numberColumns];
1629          int nBound;
1630          ClpSimplex * small2 = ((ClpSimplexOther *) (&small))->crunch(rhs,whichRow,whichColumn,
1631                                                                        nBound,false,false);
1632          if (small2) {
1633            small2->primal(1);
1634            if (small2->problemStatus()==0) {
1635              small.setProblemStatus(0);
1636              ((ClpSimplexOther *) (&small))->afterCrunch(*small2,whichRow,whichColumn,nBound);
1637            } else {
1638              small2->primal(1);
1639              if (small2->problemStatus())
1640                small.primal(1);
1641            }
1642            delete small2;
1643          } else {
1644            small.primal(1);
1645          }
1646          delete [] whichRow;
1647          delete [] whichColumn;
1648#endif
1649        }
1650      } else {
1651        small.primal(1);
1652      }
1653      totalIterations += small.numberIterations();
1654      // move solution back
1655      const double * solution = small.primalColumnSolution();
1656      for (iColumn=0;iColumn<numberSort;iColumn++) {
1657        int kColumn = sort[iColumn];
1658        model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
1659        fullSolution[kColumn]=solution[iColumn];
1660      }
1661      for (iRow=0;iRow<numberRows;iRow++) 
1662        model2->setRowStatus(iRow,small.getRowStatus(iRow));
1663      CoinMemcpyN(small.primalRowSolution(),
1664             numberRows,model2->primalRowSolution());
1665      // get reduced cost for large problem
1666      double * djs = model2->dualColumnSolution();
1667      CoinMemcpyN(model2->objective(),numberColumns,djs);
1668      model2->clpMatrix()->transposeTimes(-1.0,small.dualRowSolution(),djs);
1669      int numberNegative=0;
1670      double sumNegative = 0.0;
1671      // now massage weight so all basic in plus good djs
1672      // first count and do basic
1673      numberSort=0;
1674      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1675        double dj = djs[iColumn]*optimizationDirection_;
1676        double value = fullSolution[iColumn];
1677        if (model2->getColumnStatus(iColumn)==ClpSimplex::basic) {
1678          sort[numberSort++] = iColumn;
1679        } else if (dj<-dualTolerance_&&value<columnUpper[iColumn]) {
1680          numberNegative++;
1681          sumNegative -= dj;
1682        } else if (dj>dualTolerance_&&value>columnLower[iColumn]) {
1683          numberNegative++;
1684          sumNegative += dj;
1685        }
1686      }
1687      handler_->message(CLP_SPRINT,messages_)
1688        <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative
1689        <<numberNegative
1690        <<CoinMessageEol;
1691      if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5)||
1692          (!small.numberIterations()&&iPass)||
1693          iPass==maxSprintPass-1||small.status()==3) {
1694       
1695        break; // finished
1696      } else {
1697        lastObjective = small.objectiveValue()*optimizationDirection_;
1698        double tolerance;
1699        double averageNegDj = sumNegative/((double) (numberNegative+1));
1700        if (numberNegative+numberSort>smallNumberColumns)
1701          tolerance = -dualTolerance_;
1702        else 
1703          tolerance = 10.0*averageNegDj;
1704        int saveN = numberSort;
1705        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1706          double dj = djs[iColumn]*optimizationDirection_;
1707          double value = fullSolution[iColumn];
1708          if (model2->getColumnStatus(iColumn)!=ClpSimplex::basic) {
1709            if (dj<-dualTolerance_&&value<columnUpper[iColumn])
1710              dj = dj;
1711            else if (dj>dualTolerance_&&value>columnLower[iColumn])
1712              dj = -dj;
1713            else if (columnUpper[iColumn]>columnLower[iColumn])
1714              dj = fabs(dj);
1715            else
1716              dj = 1.0e50;
1717            if (dj<tolerance) {
1718              weight[numberSort] = dj;
1719              sort[numberSort++] = iColumn;
1720            }
1721          }
1722        }
1723        // sort
1724        CoinSort_2(weight+saveN,weight+numberSort,sort+saveN);
1725        numberSort = CoinMin(smallNumberColumns,numberSort);
1726      }
1727    }
1728    if (interrupt) 
1729      currentModel = model2;
1730    for (i=0;i<numberArtificials;i++)
1731      sort[i] = i + originalNumberColumns;
1732    model2->deleteColumns(numberArtificials,sort);
1733    delete [] weight;
1734    delete [] sort;
1735    delete [] whichRows;
1736    if (saveLower) {
1737      // unperturb and clean
1738      for (iRow=0;iRow<numberRows;iRow++) {
1739        double diffLower = saveLower[iRow]-model2->rowLower_[iRow];
1740        double diffUpper = saveUpper[iRow]-model2->rowUpper_[iRow];
1741        model2->rowLower_[iRow]=saveLower[iRow];
1742        model2->rowUpper_[iRow]=saveUpper[iRow];
1743        if (diffLower) 
1744          assert (!diffUpper||fabs(diffLower-diffUpper)<1.0e-5);
1745        else
1746          diffLower = diffUpper;
1747        model2->rowActivity_[iRow] += diffLower;
1748      }
1749      delete [] saveLower;
1750      delete [] saveUpper;
1751    }
1752    model2->primal(1);
1753    model2->setPerturbation(savePerturbation);
1754    time2 = CoinCpuTime();
1755    timeCore = time2-timeX;
1756    handler_->message(CLP_INTERVAL_TIMING,messages_)
1757      <<"Sprint"<<timeCore<<time2-time1
1758      <<CoinMessageEol;
1759    timeX=time2;
1760    model2->setNumberIterations(model2->numberIterations()+totalIterations);
1761  } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
1762#ifndef SLIM_CLP
1763    //printf("***** experimental pretty crude barrier\n");
1764    //#define SAVEIT 2
1765#ifndef SAVEIT
1766#define BORROW
1767#endif
1768#ifdef BORROW
1769    ClpInterior barrier;
1770    barrier.borrowModel(*model2);
1771#else
1772    ClpInterior barrier(*model2);
1773#endif
1774    if (interrupt) 
1775      currentModel2 = &barrier;
1776    int barrierOptions = options.getSpecialOption(4);
1777    int aggressiveGamma=0;
1778    bool presolveInCrossover=false;
1779    bool scale=false;
1780    bool doKKT=false;
1781    if (barrierOptions&16) {
1782      barrierOptions &= ~16;
1783      doKKT=true;
1784    }
1785    if (barrierOptions&(32+64+128)) {
1786      aggressiveGamma=(barrierOptions&(32+64+128))>>5;
1787      barrierOptions &= ~(32+64+128);
1788    }
1789    if (barrierOptions&256) {
1790      barrierOptions &= ~256;
1791      presolveInCrossover=true;
1792    }
1793    if (barrierOptions&8) {
1794      barrierOptions &= ~8;
1795      scale=true;
1796    }
1797#ifdef COIN_DEVELOP
1798#ifndef FAST_BARRIER
1799    if (!numberBarrier)
1800      std::cout<<"Warning - the default ordering is just on row counts! "
1801               <<"The factorization is being improved"<<std::endl;
1802    numberBarrier++;
1803#endif
1804#endif
1805    // If quadratic force KKT
1806    if (quadraticObj) {
1807      doKKT=true;
1808    }
1809    switch (barrierOptions) {
1810    case 0:
1811    default:
1812      if (!doKKT) {
1813        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1814        barrier.setCholesky(cholesky);
1815      } else {
1816        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1817        cholesky->setKKT(true);
1818        barrier.setCholesky(cholesky);
1819      }
1820      break;
1821    case 1:
1822      if (!doKKT) {
1823        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1824        barrier.setCholesky(cholesky);
1825      } else {
1826        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1827        cholesky->setKKT(true);
1828        barrier.setCholesky(cholesky);
1829      }
1830      break;
1831#ifdef WSSMP_BARRIER
1832    case 2:
1833      {
1834        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
1835        barrier.setCholesky(cholesky);
1836        assert (!doKKT);
1837      }
1838      break;
1839    case 3:
1840      if (!doKKT) {
1841        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
1842        barrier.setCholesky(cholesky);
1843      } else {
1844        ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
1845        barrier.setCholesky(cholesky);
1846      }
1847      break;
1848#endif
1849#ifdef UFL_BARRIER
1850    case 4:
1851      if (!doKKT) {
1852        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1853        barrier.setCholesky(cholesky);
1854      } else {
1855        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1856        cholesky->setKKT(true);
1857        barrier.setCholesky(cholesky);
1858      }
1859      break;
1860#endif
1861#ifdef TAUCS_BARRIER
1862    case 5:
1863      {
1864        ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
1865        barrier.setCholesky(cholesky);
1866        assert (!doKKT);
1867      }
1868      break;
1869#endif
1870    }
1871    int numberRows = model2->numberRows();
1872    int numberColumns = model2->numberColumns();
1873    int saveMaxIts = model2->maximumIterations();
1874    if (saveMaxIts<1000) {
1875      barrier.setMaximumBarrierIterations(saveMaxIts);
1876      model2->setMaximumIterations(1000000);
1877    }
1878#ifndef SAVEIT
1879    //barrier.setDiagonalPerturbation(1.0e-25);
1880    if (aggressiveGamma) {
1881      switch (aggressiveGamma) {
1882      case 1:
1883        barrier.setGamma(1.0e-5);
1884        barrier.setDelta(1.0e-5);
1885        break;
1886      case 2:
1887        barrier.setGamma(1.0e-5);
1888        break;
1889      case 3:
1890        barrier.setDelta(1.0e-5);
1891        break;
1892      case 4:
1893        barrier.setGamma(1.0e-3);
1894        barrier.setDelta(1.0e-3);
1895        break;
1896      case 5:
1897        barrier.setGamma(1.0e-3);
1898        break;
1899      case 6:
1900        barrier.setDelta(1.0e-3);
1901        break;
1902      }
1903    }
1904    if (scale)
1905      barrier.scaling(1);
1906    else
1907      barrier.scaling(0);
1908    barrier.primalDual();
1909#elif SAVEIT==1
1910    barrier.primalDual();
1911#else
1912    model2->restoreModel("xx.save");
1913    // move solutions
1914    CoinMemcpyN(model2->primalRowSolution(),
1915                numberRows,barrier.primalRowSolution());
1916    CoinMemcpyN(model2->dualRowSolution(),
1917                numberRows,barrier.dualRowSolution());
1918    CoinMemcpyN(model2->primalColumnSolution(),
1919                numberColumns,barrier.primalColumnSolution());
1920    CoinMemcpyN(model2->dualColumnSolution(),
1921                numberColumns,barrier.dualColumnSolution());
1922#endif
1923    time2 = CoinCpuTime();
1924    timeCore = time2-timeX;
1925    handler_->message(CLP_INTERVAL_TIMING,messages_)
1926      <<"Barrier"<<timeCore<<time2-time1
1927      <<CoinMessageEol;
1928    timeX=time2;
1929    int maxIts = barrier.maximumBarrierIterations();
1930    int barrierStatus=barrier.status();
1931    double gap = barrier.complementarityGap();
1932    // get which variables are fixed
1933    double * saveLower=NULL;
1934    double * saveUpper=NULL;
1935    ClpPresolve pinfo2;
1936    ClpSimplex * saveModel2=NULL;
1937    bool extraPresolve=false;
1938    int numberFixed = barrier.numberFixed();
1939    if (numberFixed) {
1940      int numberRows = barrier.numberRows();
1941      int numberColumns = barrier.numberColumns();
1942      int numberTotal = numberRows+numberColumns;
1943      saveLower = new double [numberTotal];
1944      saveUpper = new double [numberTotal];
1945      CoinMemcpyN(barrier.columnLower(),numberColumns,saveLower);
1946      CoinMemcpyN(barrier.rowLower(),numberRows,saveLower+numberColumns);
1947      CoinMemcpyN(barrier.columnUpper(),numberColumns,saveUpper);
1948      CoinMemcpyN(barrier.rowUpper(),numberRows,saveUpper+numberColumns);
1949    }
1950    if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&
1951        presolveInCrossover) {
1952      // may as well do presolve
1953      barrier.fixFixed();
1954      saveModel2=model2;
1955      extraPresolve=true;
1956    } else if (numberFixed) {
1957      // Set fixed to bounds (may have restored earlier solution)
1958      barrier.fixFixed(false);
1959    }
1960#ifdef BORROW   
1961    barrier.returnModel(*model2);
1962    double * rowPrimal = new double [numberRows];
1963    double * columnPrimal = new double [numberColumns];
1964    double * rowDual = new double [numberRows];
1965    double * columnDual = new double [numberColumns];
1966    // move solutions other way
1967    CoinMemcpyN(model2->primalRowSolution(),
1968                numberRows,rowPrimal);
1969    CoinMemcpyN(model2->dualRowSolution(),
1970                numberRows,rowDual);
1971    CoinMemcpyN(model2->primalColumnSolution(),
1972                numberColumns,columnPrimal);
1973    CoinMemcpyN(model2->dualColumnSolution(),
1974                  numberColumns,columnDual);
1975#else
1976    double * rowPrimal = barrier.primalRowSolution();
1977    double * columnPrimal = barrier.primalColumnSolution();
1978    double * rowDual = barrier.dualRowSolution();
1979    double * columnDual = barrier.dualColumnSolution();
1980    // move solutions
1981    CoinMemcpyN(rowPrimal,
1982                numberRows,model2->primalRowSolution());
1983    CoinMemcpyN(rowDual,
1984                numberRows,model2->dualRowSolution());
1985    CoinMemcpyN(columnPrimal,
1986                numberColumns,model2->primalColumnSolution());
1987    CoinMemcpyN(columnDual,
1988                  numberColumns,model2->dualColumnSolution());
1989#endif
1990    if (saveModel2) {
1991      // do presolve
1992      model2 = pinfo2.presolvedModel(*model2,dblParam_[ClpPresolveTolerance],
1993                                    false,5,true);
1994      if (!model2) {
1995        model2=saveModel2;
1996        saveModel2=NULL;
1997        int numberRows = model2->numberRows();
1998        int numberColumns = model2->numberColumns();
1999        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2000        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2001        delete [] saveLower;
2002        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2003        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2004        delete [] saveUpper;
2005        saveLower=NULL;
2006        saveUpper=NULL;
2007      }
2008    }
2009    if (method==ClpSolve::useBarrier) {
2010      if (maxIts&&barrierStatus<4&&!quadraticObj) {
2011        //printf("***** crossover - needs more thought on difficult models\n");
2012#if SAVEIT==1
2013        model2->ClpSimplex::saveModel("xx.save");
2014#endif
2015        // make sure no status left
2016        model2->createStatus();
2017        // solve
2018        model2->setPerturbation(100);
2019        if (model2->factorizationFrequency()==200) {
2020          // User did not touch preset
2021          model2->defaultFactorizationFrequency();
2022        }
2023#if 1
2024        // throw some into basis
2025        {
2026          int numberRows = model2->numberRows();
2027          int numberColumns = model2->numberColumns();
2028          double * dsort = new double[numberColumns];
2029          int * sort = new int[numberColumns];
2030          int n=0;
2031          const double * columnLower = model2->columnLower();
2032          const double * columnUpper = model2->columnUpper();
2033          double * primalSolution = model2->primalColumnSolution();
2034          const double * dualSolution = model2->dualColumnSolution();
2035          double tolerance = 10.0*primalTolerance_;
2036          int i;
2037          for ( i=0;i<numberRows;i++) 
2038            model2->setRowStatus(i,superBasic);
2039          for ( i=0;i<numberColumns;i++) {
2040            double distance = CoinMin(columnUpper[i]-primalSolution[i],
2041                                      primalSolution[i]-columnLower[i]);
2042            if (distance>tolerance) {
2043              if (fabs(dualSolution[i])<1.0e-5)
2044                distance *= 100.0;
2045              dsort[n]=-distance;
2046              sort[n++]=i;
2047              model2->setStatus(i,superBasic);
2048            } else if (distance>primalTolerance_) {
2049              model2->setStatus(i,superBasic);
2050            } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
2051              model2->setStatus(i,atLowerBound);
2052              primalSolution[i]=columnLower[i];
2053            } else {
2054              model2->setStatus(i,atUpperBound);
2055              primalSolution[i]=columnUpper[i];
2056            }
2057          }
2058          CoinSort_2(dsort,dsort+n,sort);
2059          n = CoinMin(numberRows,n);
2060          for ( i=0;i<n;i++) {
2061            int iColumn = sort[i];
2062            model2->setStatus(iColumn,basic);
2063          }
2064          delete [] sort;
2065          delete [] dsort;
2066        }
2067        // model2->allSlackBasis();
2068        if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
2069          if (saveUpper) {
2070            int numberRows = model2->numberRows();
2071            int numberColumns = model2->numberColumns();
2072            CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2073            CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2074            delete [] saveLower;
2075            CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2076            CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2077            delete [] saveUpper;
2078            saveLower=NULL;
2079            saveUpper=NULL;
2080          }
2081          int numberRows = model2->numberRows();
2082          int numberColumns = model2->numberColumns();
2083          // just primal values pass
2084          double saveScale = model2->objectiveScale();
2085          model2->setObjectiveScale(1.0e-3);
2086          model2->primal(2);
2087          model2->setObjectiveScale(saveScale);
2088          // save primal solution and copy back dual
2089          CoinMemcpyN(model2->primalRowSolution(),
2090                      numberRows,rowPrimal);
2091          CoinMemcpyN(rowDual,
2092                      numberRows,model2->dualRowSolution());
2093          CoinMemcpyN(model2->primalColumnSolution(),
2094                      numberColumns,columnPrimal);
2095          CoinMemcpyN(columnDual,
2096                      numberColumns,model2->dualColumnSolution());
2097          //model2->primal(1);
2098          // clean up reduced costs and flag variables
2099          {
2100            double * dj = model2->dualColumnSolution();
2101            double * cost = model2->objective();
2102            double * saveCost = new double[numberColumns];
2103            CoinMemcpyN(cost,numberColumns,saveCost);
2104            double * saveLower = new double[numberColumns];
2105            double * lower = model2->columnLower();
2106            CoinMemcpyN(lower,numberColumns,saveLower);
2107            double * saveUpper = new double[numberColumns];
2108            double * upper = model2->columnUpper();
2109            CoinMemcpyN(upper,numberColumns,saveUpper);
2110            int i;
2111            double tolerance = 10.0*dualTolerance_;
2112            for ( i=0;i<numberColumns;i++) {
2113              if (model2->getStatus(i)==basic) {
2114                dj[i]=0.0;
2115              } else if (model2->getStatus(i)==atLowerBound) {
2116                if (optimizationDirection_*dj[i]<tolerance) {
2117                  if (optimizationDirection_*dj[i]<0.0) {
2118                    //if (dj[i]<-1.0e-3)
2119                    //printf("bad dj at lb %d %g\n",i,dj[i]);
2120                    cost[i] -= dj[i];
2121                    dj[i]=0.0;
2122                  }
2123                } else {
2124                  upper[i]=lower[i];
2125                }
2126              } else if (model2->getStatus(i)==atUpperBound) {
2127                if (optimizationDirection_*dj[i]>tolerance) {
2128                  if (optimizationDirection_*dj[i]>0.0) {
2129                    //if (dj[i]>1.0e-3)
2130                    //printf("bad dj at ub %d %g\n",i,dj[i]);
2131                    cost[i] -= dj[i];
2132                    dj[i]=0.0;
2133                  }
2134                } else {
2135                  lower[i]=upper[i];
2136                }
2137              }
2138            }
2139            // just dual values pass
2140            //model2->setLogLevel(63);
2141            //model2->setFactorizationFrequency(1);
2142            model2->dual(2);
2143            CoinMemcpyN(saveCost,numberColumns,cost);
2144            delete [] saveCost;
2145            CoinMemcpyN(saveLower,numberColumns,lower);
2146            delete [] saveLower;
2147            CoinMemcpyN(saveUpper,numberColumns,upper);
2148            delete [] saveUpper;
2149          }
2150          // and finish
2151          // move solutions
2152          CoinMemcpyN(rowPrimal,
2153                      numberRows,model2->primalRowSolution());
2154          CoinMemcpyN(columnPrimal,
2155                      numberColumns,model2->primalColumnSolution());
2156        }
2157        double saveScale = model2->objectiveScale();
2158        model2->setObjectiveScale(1.0e-3);
2159        model2->primal(2);
2160        model2->setObjectiveScale(saveScale);
2161        model2->primal(1);
2162#else
2163        // just primal
2164        model2->primal(1);
2165#endif
2166      } else if (barrierStatus==4) {
2167        // memory problems
2168        model2->setPerturbation(savePerturbation);
2169        model2->createStatus();
2170        model2->dual();
2171      } else if (maxIts&&quadraticObj) {
2172        // make sure no status left
2173        model2->createStatus();
2174        // solve
2175        model2->setPerturbation(100);
2176        model2->reducedGradient(1);
2177      }
2178    }
2179    model2->setMaximumIterations(saveMaxIts);
2180#ifdef BORROW
2181    delete [] rowPrimal;
2182    delete [] columnPrimal;
2183    delete [] rowDual;
2184    delete [] columnDual;
2185#endif
2186    if (extraPresolve) {
2187      pinfo2.postsolve(true);
2188      delete model2;
2189      model2=saveModel2;
2190    }
2191    if (saveUpper) {
2192      int numberRows = model2->numberRows();
2193      int numberColumns = model2->numberColumns();
2194      CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2195      CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2196      delete [] saveLower;
2197      CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2198      CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2199      delete [] saveUpper;
2200      saveLower=NULL;
2201      saveUpper=NULL;
2202      if (method!=ClpSolve::useBarrierNoCross) 
2203        model2->primal(1);
2204    }
2205    model2->setPerturbation(savePerturbation);
2206    time2 = CoinCpuTime();
2207    timeCore = time2-timeX;
2208    handler_->message(CLP_INTERVAL_TIMING,messages_)
2209      <<"Crossover"<<timeCore<<time2-time1
2210      <<CoinMessageEol;
2211    timeX=time2;
2212#else
2213    abort();
2214#endif
2215  } else {
2216    assert (method!=ClpSolve::automatic); // later
2217    time2=0.0;
2218  }
2219  if (saveMatrix&&model2==this) {
2220    // delete and replace
2221    delete model2->clpMatrix();
2222    model2->replaceMatrix(saveMatrix);
2223  }
2224  numberIterations = model2->numberIterations();
2225  finalStatus=model2->status();
2226  int finalSecondaryStatus = model2->secondaryStatus();
2227  if (presolve==ClpSolve::presolveOn) {
2228    int saveLevel = logLevel();
2229    if ((specialOptions_&1024)==0)
2230      setLogLevel(CoinMin(1,saveLevel));
2231    else
2232      setLogLevel(CoinMin(0,saveLevel));
2233    pinfo.postsolve(true);
2234    factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
2235    time2 = CoinCpuTime();
2236    timePresolve += time2-timeX;
2237    handler_->message(CLP_INTERVAL_TIMING,messages_)
2238      <<"Postsolve"<<time2-timeX<<time2-time1
2239      <<CoinMessageEol;
2240    timeX=time2;
2241    if (!presolveToFile)
2242      delete model2;
2243    if (interrupt)
2244      currentModel = this;
2245    // checkSolution(); already done by postSolve
2246    setLogLevel(saveLevel);
2247    if (finalStatus!=3&&(finalStatus||status()==-1)) {
2248      int savePerturbation = perturbation();
2249      setPerturbation(100);
2250      if (finalStatus==2) {
2251        // unbounded - get feasible first
2252        double save = optimizationDirection_;
2253        optimizationDirection_=0.0;
2254        primal(1);
2255        optimizationDirection_=save;
2256      }
2257      primal(1);
2258      setPerturbation(savePerturbation);
2259      numberIterations += numberIterations_;
2260      numberIterations_ = numberIterations;
2261      finalStatus=status();
2262      time2 = CoinCpuTime();
2263      handler_->message(CLP_INTERVAL_TIMING,messages_)
2264      <<"Cleanup"<<time2-timeX<<time2-time1
2265      <<CoinMessageEol;
2266      timeX=time2;
2267    } else {
2268      secondaryStatus_=finalSecondaryStatus;
2269    }
2270  } else if (model2!=this) {
2271    // not presolved - but different model used (sprint probably)
2272    CoinMemcpyN(model2->primalRowSolution(),
2273                numberRows_,this->primalRowSolution());
2274    CoinMemcpyN(model2->dualRowSolution(),
2275                numberRows_,this->dualRowSolution());
2276    CoinMemcpyN(model2->primalColumnSolution(),
2277                numberColumns_,this->primalColumnSolution());
2278    CoinMemcpyN(model2->dualColumnSolution(),
2279                numberColumns_,this->dualColumnSolution());
2280    CoinMemcpyN(model2->statusArray(),
2281                numberColumns_+numberRows_,this->statusArray());
2282    objectiveValue_=model2->objectiveValue_;
2283    numberIterations_ =model2->numberIterations_;
2284    problemStatus_ =model2->problemStatus_;
2285    secondaryStatus_ =model2->secondaryStatus_;
2286    delete model2;
2287  }
2288  setMaximumIterations(saveMaxIterations);
2289  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
2290                               "Errors","User stopped"};
2291  assert (finalStatus>=-1&&finalStatus<=5);
2292  handler_->message(CLP_TIMING,messages_)
2293    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
2294  handler_->printing(presolve==ClpSolve::presolveOn)
2295    <<timePresolve;
2296  handler_->printing(timeIdiot!=0.0)
2297    <<timeIdiot;
2298  handler_->message()<<CoinMessageEol;
2299  if (interrupt) 
2300    signal(SIGINT,saveSignal);
2301  perturbation_=savePerturbation;
2302  scalingFlag_=saveScaling;
2303  // If faking objective - put back correct one
2304  if (savedObjective) {
2305    delete objective_;
2306    objective_=savedObjective;
2307  }
2308  return finalStatus;
2309}
2310// General solve
2311int 
2312ClpSimplex::initialSolve()
2313{
2314  // Default so use dual
2315  ClpSolve options;
2316  return initialSolve(options);
2317}
2318// General dual solve
2319int 
2320ClpSimplex::initialDualSolve()
2321{
2322  ClpSolve options;
2323  // Use dual
2324  options.setSolveType(ClpSolve::useDual);
2325  return initialSolve(options);
2326}
2327// General dual solve
2328int 
2329ClpSimplex::initialPrimalSolve()
2330{
2331  ClpSolve options;
2332  // Use primal
2333  options.setSolveType(ClpSolve::usePrimal);
2334  return initialSolve(options);
2335}
2336// barrier solve, not to be followed by crossover
2337int 
2338ClpSimplex::initialBarrierNoCrossSolve()
2339{
2340  ClpSolve options;
2341  // Use primal
2342  options.setSolveType(ClpSolve::useBarrierNoCross);
2343  return initialSolve(options);
2344}
2345
2346// General barrier solve
2347int 
2348ClpSimplex::initialBarrierSolve()
2349{
2350  ClpSolve options;
2351  // Use primal
2352  options.setSolveType(ClpSolve::useBarrier);
2353  return initialSolve(options);
2354}
2355
2356// Default constructor
2357ClpSolve::ClpSolve (  )
2358{
2359  method_ = automatic;
2360  presolveType_=presolveOn;
2361  numberPasses_=5;
2362  int i;
2363  for (i=0;i<7;i++)
2364    options_[i]=0;
2365  // say no +-1 matrix
2366  options_[3]=1;
2367  for (i=0;i<7;i++)
2368    extraInfo_[i]=-1;
2369  independentOptions_[0]=0;
2370  // But switch off slacks
2371  independentOptions_[1]=512;
2372  // Substitute up to 3
2373  independentOptions_[2]=3;
2374 
2375}
2376// Constructor when you really know what you are doing
2377ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
2378             int numberPasses, int options[6],
2379             int extraInfo[6], int independentOptions[3])
2380{
2381  method_ = method;
2382  presolveType_=presolveType;
2383  numberPasses_=numberPasses;
2384  int i;
2385  for (i=0;i<6;i++)
2386    options_[i]=options[i];
2387  options_[6]=0;
2388  for (i=0;i<6;i++)
2389    extraInfo_[i]=extraInfo[i];
2390  extraInfo_[6]=0;
2391  for (i=0;i<3;i++)
2392    independentOptions_[i]=independentOptions[i];
2393}
2394
2395// Copy constructor.
2396ClpSolve::ClpSolve(const ClpSolve & rhs)
2397{
2398  method_ = rhs.method_;
2399  presolveType_=rhs.presolveType_;
2400  numberPasses_=rhs.numberPasses_;
2401  int i;
2402  for ( i=0;i<7;i++)
2403    options_[i]=rhs.options_[i];
2404  for ( i=0;i<7;i++)
2405    extraInfo_[i]=rhs.extraInfo_[i];
2406  for (i=0;i<3;i++)
2407    independentOptions_[i]=rhs.independentOptions_[i];
2408}
2409// Assignment operator. This copies the data
2410ClpSolve & 
2411ClpSolve::operator=(const ClpSolve & rhs)
2412{
2413  if (this != &rhs) {
2414    method_ = rhs.method_;
2415    presolveType_=rhs.presolveType_;
2416    numberPasses_=rhs.numberPasses_;
2417    int i;
2418    for (i=0;i<7;i++)
2419      options_[i]=rhs.options_[i];
2420    for (i=0;i<7;i++)
2421      extraInfo_[i]=rhs.extraInfo_[i];
2422    for (i=0;i<3;i++)
2423      independentOptions_[i]=rhs.independentOptions_[i];
2424  }
2425  return *this;
2426
2427}
2428// Destructor
2429ClpSolve::~ClpSolve (  )
2430{
2431}
2432// See header file for deatils
2433void 
2434ClpSolve::setSpecialOption(int which,int value,int extraInfo)
2435{
2436  options_[which]=value;
2437  extraInfo_[which]=extraInfo;
2438}
2439int 
2440ClpSolve::getSpecialOption(int which) const
2441{
2442  return options_[which];
2443}
2444
2445// Solve types
2446void 
2447ClpSolve::setSolveType(SolveType method, int extraInfo)
2448{
2449  method_=method;
2450}
2451
2452ClpSolve::SolveType
2453ClpSolve::getSolveType()
2454{
2455  return method_;
2456}
2457
2458// Presolve types
2459void 
2460ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
2461{
2462  presolveType_=amount;
2463  numberPasses_=extraInfo;
2464}
2465ClpSolve::PresolveType
2466ClpSolve::getPresolveType()
2467{
2468  return presolveType_;
2469}
2470// Extra info for idiot (or sprint)
2471int 
2472ClpSolve::getExtraInfo(int which) const
2473{
2474  return extraInfo_[which];
2475}
2476int 
2477ClpSolve::getPresolvePasses() const
2478{
2479  return numberPasses_;
2480}
2481/* Say to return at once if infeasible,
2482   default is to solve */
2483void 
2484ClpSolve::setInfeasibleReturn(bool trueFalse)
2485{
2486  independentOptions_[0]= trueFalse ? 1 : 0;
2487}
2488#include <string>
2489// Generates code for above constructor
2490void 
2491ClpSolve::generateCpp(FILE * fp)
2492{
2493  std::string solveType[] = {
2494    "ClpSolve::useDual",
2495    "ClpSolve::usePrimal",
2496    "ClpSolve::usePrimalorSprint",
2497    "ClpSolve::useBarrier",
2498    "ClpSolve::useBarrierNoCross",
2499    "ClpSolve::automatic",
2500    "ClpSolve::notImplemented"
2501  };
2502  std::string presolveType[] =  {
2503    "ClpSolve::presolveOn",
2504    "ClpSolve::presolveOff",
2505    "ClpSolve::presolveNumber",
2506    "ClpSolve::presolveNumberCost"
2507  };
2508  fprintf(fp,"3  ClpSolve::SolveType method = %s;\n",solveType[method_].c_str());
2509  fprintf(fp,"3  ClpSolve::PresolveType presolveType = %s;\n",
2510    presolveType[presolveType_].c_str());
2511  fprintf(fp,"3  int numberPasses = %d;\n",numberPasses_);
2512  fprintf(fp,"3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
2513    options_[0],options_[1],options_[2],
2514    options_[3],options_[4],options_[5]);
2515  fprintf(fp,"3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
2516    extraInfo_[0],extraInfo_[1],extraInfo_[2],
2517    extraInfo_[3],extraInfo_[4],extraInfo_[5]);
2518  fprintf(fp,"3  int independentOptions[] = {%d,%d,%d};\n",
2519    independentOptions_[0],independentOptions_[1],independentOptions_[2]);
2520  fprintf(fp,"3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
2521  fprintf(fp,"3                    options,extraInfo,independentOptions);\n");
2522}
Note: See TracBrowser for help on using the repository browser.