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

Last change on this file since 1307 was 1307, checked in by forrest, 11 years ago

fix memory leak

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