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

Last change on this file since 1197 was 1197, checked in by forrest, 13 years ago

many changes to try and improve performance

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