source: stable/1.8/Clp/src/ClpSolve.cpp @ 1236

Last change on this file since 1236 was 1236, checked in by ladanyi, 13 years ago

force cdecl on signal_handler() and main() if MS compiler is used (ticket #20)

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