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

Last change on this file since 798 was 798, checked in by ladanyi, 14 years ago

finished switch to subversion

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