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

Last change on this file since 1304 was 1277, checked in by forrest, 12 years ago

out printout and bad presolve

  • 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    if (model2!=originalModel2)
1840      originalModel2->moveInfo(*model2);
1841    model2->setPerturbation(savePerturbation);
1842    time2 = CoinCpuTime();
1843    timeCore = time2-timeX;
1844    handler_->message(CLP_INTERVAL_TIMING,messages_)
1845      <<"Sprint"<<timeCore<<time2-time1
1846      <<CoinMessageEol;
1847    timeX=time2;
1848    model2->setNumberIterations(model2->numberIterations()+totalIterations);
1849  } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
1850#ifndef SLIM_CLP
1851    //printf("***** experimental pretty crude barrier\n");
1852    //#define SAVEIT 2
1853#ifndef SAVEIT
1854#define BORROW
1855#endif
1856#ifdef BORROW
1857    ClpInterior barrier;
1858    barrier.borrowModel(*model2);
1859#else
1860    ClpInterior barrier(*model2);
1861#endif
1862    if (interrupt) 
1863      currentModel2 = &barrier;
1864    int barrierOptions = options.getSpecialOption(4);
1865    int aggressiveGamma=0;
1866    bool presolveInCrossover=false;
1867    bool scale=false;
1868    bool doKKT=false;
1869    if (barrierOptions&16) {
1870      barrierOptions &= ~16;
1871      doKKT=true;
1872    }
1873    if (barrierOptions&(32+64+128)) {
1874      aggressiveGamma=(barrierOptions&(32+64+128))>>5;
1875      barrierOptions &= ~(32+64+128);
1876    }
1877    if (barrierOptions&256) {
1878      barrierOptions &= ~256;
1879      presolveInCrossover=true;
1880    }
1881    if (barrierOptions&8) {
1882      barrierOptions &= ~8;
1883      scale=true;
1884    }
1885#ifdef COIN_DEVELOP
1886#ifndef FAST_BARRIER
1887    if (!numberBarrier)
1888      std::cout<<"Warning - the default ordering is just on row counts! "
1889               <<"The factorization is being improved"<<std::endl;
1890    numberBarrier++;
1891#endif
1892#endif
1893    // If quadratic force KKT
1894    if (quadraticObj) {
1895      doKKT=true;
1896    }
1897    switch (barrierOptions) {
1898    case 0:
1899    default:
1900      if (!doKKT) {
1901        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1902        barrier.setCholesky(cholesky);
1903      } else {
1904        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1905        cholesky->setKKT(true);
1906        barrier.setCholesky(cholesky);
1907      }
1908      break;
1909    case 1:
1910      if (!doKKT) {
1911        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1912        barrier.setCholesky(cholesky);
1913      } else {
1914        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1915        cholesky->setKKT(true);
1916        barrier.setCholesky(cholesky);
1917      }
1918      break;
1919#ifdef WSSMP_BARRIER
1920    case 2:
1921      {
1922        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
1923        barrier.setCholesky(cholesky);
1924        assert (!doKKT);
1925      }
1926      break;
1927    case 3:
1928      if (!doKKT) {
1929        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
1930        barrier.setCholesky(cholesky);
1931      } else {
1932        ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
1933        barrier.setCholesky(cholesky);
1934      }
1935      break;
1936#endif
1937#ifdef UFL_BARRIER
1938    case 4:
1939      if (!doKKT) {
1940        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1941        barrier.setCholesky(cholesky);
1942      } else {
1943        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1944        cholesky->setKKT(true);
1945        barrier.setCholesky(cholesky);
1946      }
1947      break;
1948#endif
1949#ifdef TAUCS_BARRIER
1950    case 5:
1951      {
1952        ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
1953        barrier.setCholesky(cholesky);
1954        assert (!doKKT);
1955      }
1956      break;
1957#endif
1958    }
1959    int numberRows = model2->numberRows();
1960    int numberColumns = model2->numberColumns();
1961    int saveMaxIts = model2->maximumIterations();
1962    if (saveMaxIts<1000) {
1963      barrier.setMaximumBarrierIterations(saveMaxIts);
1964      model2->setMaximumIterations(1000000);
1965    }
1966#ifndef SAVEIT
1967    //barrier.setDiagonalPerturbation(1.0e-25);
1968    if (aggressiveGamma) {
1969      switch (aggressiveGamma) {
1970      case 1:
1971        barrier.setGamma(1.0e-5);
1972        barrier.setDelta(1.0e-5);
1973        break;
1974      case 2:
1975        barrier.setGamma(1.0e-5);
1976        break;
1977      case 3:
1978        barrier.setDelta(1.0e-5);
1979        break;
1980      case 4:
1981        barrier.setGamma(1.0e-3);
1982        barrier.setDelta(1.0e-3);
1983        break;
1984      case 5:
1985        barrier.setGamma(1.0e-3);
1986        break;
1987      case 6:
1988        barrier.setDelta(1.0e-3);
1989        break;
1990      }
1991    }
1992    if (scale)
1993      barrier.scaling(1);
1994    else
1995      barrier.scaling(0);
1996    barrier.primalDual();
1997#elif SAVEIT==1
1998    barrier.primalDual();
1999#else
2000    model2->restoreModel("xx.save");
2001    // move solutions
2002    CoinMemcpyN(model2->primalRowSolution(),
2003                numberRows,barrier.primalRowSolution());
2004    CoinMemcpyN(model2->dualRowSolution(),
2005                numberRows,barrier.dualRowSolution());
2006    CoinMemcpyN(model2->primalColumnSolution(),
2007                numberColumns,barrier.primalColumnSolution());
2008    CoinMemcpyN(model2->dualColumnSolution(),
2009                numberColumns,barrier.dualColumnSolution());
2010#endif
2011    time2 = CoinCpuTime();
2012    timeCore = time2-timeX;
2013    handler_->message(CLP_INTERVAL_TIMING,messages_)
2014      <<"Barrier"<<timeCore<<time2-time1
2015      <<CoinMessageEol;
2016    timeX=time2;
2017    int maxIts = barrier.maximumBarrierIterations();
2018    int barrierStatus=barrier.status();
2019    double gap = barrier.complementarityGap();
2020    // get which variables are fixed
2021    double * saveLower=NULL;
2022    double * saveUpper=NULL;
2023    ClpPresolve pinfo2;
2024    ClpSimplex * saveModel2=NULL;
2025    bool extraPresolve=false;
2026    int numberFixed = barrier.numberFixed();
2027    if (numberFixed) {
2028      int numberRows = barrier.numberRows();
2029      int numberColumns = barrier.numberColumns();
2030      int numberTotal = numberRows+numberColumns;
2031      saveLower = new double [numberTotal];
2032      saveUpper = new double [numberTotal];
2033      CoinMemcpyN(barrier.columnLower(),numberColumns,saveLower);
2034      CoinMemcpyN(barrier.rowLower(),numberRows,saveLower+numberColumns);
2035      CoinMemcpyN(barrier.columnUpper(),numberColumns,saveUpper);
2036      CoinMemcpyN(barrier.rowUpper(),numberRows,saveUpper+numberColumns);
2037    }
2038    if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&
2039        presolveInCrossover) {
2040      // may as well do presolve
2041      barrier.fixFixed();
2042      saveModel2=model2;
2043      extraPresolve=true;
2044    } else if (numberFixed) {
2045      // Set fixed to bounds (may have restored earlier solution)
2046      barrier.fixFixed(false);
2047    }
2048#ifdef BORROW   
2049    barrier.returnModel(*model2);
2050    double * rowPrimal = new double [numberRows];
2051    double * columnPrimal = new double [numberColumns];
2052    double * rowDual = new double [numberRows];
2053    double * columnDual = new double [numberColumns];
2054    // move solutions other way
2055    CoinMemcpyN(model2->primalRowSolution(),
2056                numberRows,rowPrimal);
2057    CoinMemcpyN(model2->dualRowSolution(),
2058                numberRows,rowDual);
2059    CoinMemcpyN(model2->primalColumnSolution(),
2060                numberColumns,columnPrimal);
2061    CoinMemcpyN(model2->dualColumnSolution(),
2062                  numberColumns,columnDual);
2063#else
2064    double * rowPrimal = barrier.primalRowSolution();
2065    double * columnPrimal = barrier.primalColumnSolution();
2066    double * rowDual = barrier.dualRowSolution();
2067    double * columnDual = barrier.dualColumnSolution();
2068    // move solutions
2069    CoinMemcpyN(rowPrimal,
2070                numberRows,model2->primalRowSolution());
2071    CoinMemcpyN(rowDual,
2072                numberRows,model2->dualRowSolution());
2073    CoinMemcpyN(columnPrimal,
2074                numberColumns,model2->primalColumnSolution());
2075    CoinMemcpyN(columnDual,
2076                  numberColumns,model2->dualColumnSolution());
2077#endif
2078    if (saveModel2) {
2079      // do presolve
2080      model2 = pinfo2.presolvedModel(*model2,dblParam_[ClpPresolveTolerance],
2081                                    false,5,true);
2082      if (!model2) {
2083        model2=saveModel2;
2084        saveModel2=NULL;
2085        int numberRows = model2->numberRows();
2086        int numberColumns = model2->numberColumns();
2087        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2088        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2089        delete [] saveLower;
2090        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2091        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2092        delete [] saveUpper;
2093        saveLower=NULL;
2094        saveUpper=NULL;
2095      }
2096    }
2097    if (method==ClpSolve::useBarrier) {
2098      if (maxIts&&barrierStatus<4&&!quadraticObj) {
2099        //printf("***** crossover - needs more thought on difficult models\n");
2100#if SAVEIT==1
2101        model2->ClpSimplex::saveModel("xx.save");
2102#endif
2103        // make sure no status left
2104        model2->createStatus();
2105        // solve
2106        model2->setPerturbation(100);
2107        if (model2->factorizationFrequency()==200) {
2108          // User did not touch preset
2109          model2->defaultFactorizationFrequency();
2110        }
2111#if 1
2112        // throw some into basis
2113        {
2114          int numberRows = model2->numberRows();
2115          int numberColumns = model2->numberColumns();
2116          double * dsort = new double[numberColumns];
2117          int * sort = new int[numberColumns];
2118          int n=0;
2119          const double * columnLower = model2->columnLower();
2120          const double * columnUpper = model2->columnUpper();
2121          double * primalSolution = model2->primalColumnSolution();
2122          const double * dualSolution = model2->dualColumnSolution();
2123          double tolerance = 10.0*primalTolerance_;
2124          int i;
2125          for ( i=0;i<numberRows;i++) 
2126            model2->setRowStatus(i,superBasic);
2127          for ( i=0;i<numberColumns;i++) {
2128            double distance = CoinMin(columnUpper[i]-primalSolution[i],
2129                                      primalSolution[i]-columnLower[i]);
2130            if (distance>tolerance) {
2131              if (fabs(dualSolution[i])<1.0e-5)
2132                distance *= 100.0;
2133              dsort[n]=-distance;
2134              sort[n++]=i;
2135              model2->setStatus(i,superBasic);
2136            } else if (distance>primalTolerance_) {
2137              model2->setStatus(i,superBasic);
2138            } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
2139              model2->setStatus(i,atLowerBound);
2140              primalSolution[i]=columnLower[i];
2141            } else {
2142              model2->setStatus(i,atUpperBound);
2143              primalSolution[i]=columnUpper[i];
2144            }
2145          }
2146          CoinSort_2(dsort,dsort+n,sort);
2147          n = CoinMin(numberRows,n);
2148          for ( i=0;i<n;i++) {
2149            int iColumn = sort[i];
2150            model2->setStatus(iColumn,basic);
2151          }
2152          delete [] sort;
2153          delete [] dsort;
2154        }
2155        // model2->allSlackBasis();
2156        if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
2157          if (saveUpper) {
2158            int numberRows = model2->numberRows();
2159            int numberColumns = model2->numberColumns();
2160            CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2161            CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2162            delete [] saveLower;
2163            CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2164            CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2165            delete [] saveUpper;
2166            saveLower=NULL;
2167            saveUpper=NULL;
2168          }
2169          int numberRows = model2->numberRows();
2170          int numberColumns = model2->numberColumns();
2171          // just primal values pass
2172          double saveScale = model2->objectiveScale();
2173          model2->setObjectiveScale(1.0e-3);
2174          model2->primal(2);
2175          model2->setObjectiveScale(saveScale);
2176          // save primal solution and copy back dual
2177          CoinMemcpyN(model2->primalRowSolution(),
2178                      numberRows,rowPrimal);
2179          CoinMemcpyN(rowDual,
2180                      numberRows,model2->dualRowSolution());
2181          CoinMemcpyN(model2->primalColumnSolution(),
2182                      numberColumns,columnPrimal);
2183          CoinMemcpyN(columnDual,
2184                      numberColumns,model2->dualColumnSolution());
2185          //model2->primal(1);
2186          // clean up reduced costs and flag variables
2187          {
2188            double * dj = model2->dualColumnSolution();
2189            double * cost = model2->objective();
2190            double * saveCost = new double[numberColumns];
2191            CoinMemcpyN(cost,numberColumns,saveCost);
2192            double * saveLower = new double[numberColumns];
2193            double * lower = model2->columnLower();
2194            CoinMemcpyN(lower,numberColumns,saveLower);
2195            double * saveUpper = new double[numberColumns];
2196            double * upper = model2->columnUpper();
2197            CoinMemcpyN(upper,numberColumns,saveUpper);
2198            int i;
2199            double tolerance = 10.0*dualTolerance_;
2200            for ( i=0;i<numberColumns;i++) {
2201              if (model2->getStatus(i)==basic) {
2202                dj[i]=0.0;
2203              } else if (model2->getStatus(i)==atLowerBound) {
2204                if (optimizationDirection_*dj[i]<tolerance) {
2205                  if (optimizationDirection_*dj[i]<0.0) {
2206                    //if (dj[i]<-1.0e-3)
2207                    //printf("bad dj at lb %d %g\n",i,dj[i]);
2208                    cost[i] -= dj[i];
2209                    dj[i]=0.0;
2210                  }
2211                } else {
2212                  upper[i]=lower[i];
2213                }
2214              } else if (model2->getStatus(i)==atUpperBound) {
2215                if (optimizationDirection_*dj[i]>tolerance) {
2216                  if (optimizationDirection_*dj[i]>0.0) {
2217                    //if (dj[i]>1.0e-3)
2218                    //printf("bad dj at ub %d %g\n",i,dj[i]);
2219                    cost[i] -= dj[i];
2220                    dj[i]=0.0;
2221                  }
2222                } else {
2223                  lower[i]=upper[i];
2224                }
2225              }
2226            }
2227            // just dual values pass
2228            //model2->setLogLevel(63);
2229            //model2->setFactorizationFrequency(1);
2230            model2->dual(2);
2231            CoinMemcpyN(saveCost,numberColumns,cost);
2232            delete [] saveCost;
2233            CoinMemcpyN(saveLower,numberColumns,lower);
2234            delete [] saveLower;
2235            CoinMemcpyN(saveUpper,numberColumns,upper);
2236            delete [] saveUpper;
2237          }
2238          // and finish
2239          // move solutions
2240          CoinMemcpyN(rowPrimal,
2241                      numberRows,model2->primalRowSolution());
2242          CoinMemcpyN(columnPrimal,
2243                      numberColumns,model2->primalColumnSolution());
2244        }
2245        double saveScale = model2->objectiveScale();
2246        model2->setObjectiveScale(1.0e-3);
2247        model2->primal(2);
2248        model2->setObjectiveScale(saveScale);
2249        model2->primal(1);
2250#else
2251        // just primal
2252        model2->primal(1);
2253#endif
2254      } else if (barrierStatus==4) {
2255        // memory problems
2256        model2->setPerturbation(savePerturbation);
2257        model2->createStatus();
2258        model2->dual();
2259      } else if (maxIts&&quadraticObj) {
2260        // make sure no status left
2261        model2->createStatus();
2262        // solve
2263        model2->setPerturbation(100);
2264        model2->reducedGradient(1);
2265      }
2266    }
2267    model2->setMaximumIterations(saveMaxIts);
2268#ifdef BORROW
2269    delete [] rowPrimal;
2270    delete [] columnPrimal;
2271    delete [] rowDual;
2272    delete [] columnDual;
2273#endif
2274    if (extraPresolve) {
2275      pinfo2.postsolve(true);
2276      delete model2;
2277      model2=saveModel2;
2278    }
2279    if (saveUpper) {
2280      int numberRows = model2->numberRows();
2281      int numberColumns = model2->numberColumns();
2282      CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2283      CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2284      delete [] saveLower;
2285      CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2286      CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2287      delete [] saveUpper;
2288      saveLower=NULL;
2289      saveUpper=NULL;
2290      if (method!=ClpSolve::useBarrierNoCross) 
2291        model2->primal(1);
2292    }
2293    model2->setPerturbation(savePerturbation);
2294    time2 = CoinCpuTime();
2295    timeCore = time2-timeX;
2296    handler_->message(CLP_INTERVAL_TIMING,messages_)
2297      <<"Crossover"<<timeCore<<time2-time1
2298      <<CoinMessageEol;
2299    timeX=time2;
2300#else
2301    abort();
2302#endif
2303  } else {
2304    assert (method!=ClpSolve::automatic); // later
2305    time2=0.0;
2306  }
2307  if (saveMatrix) {
2308    if (model2==this) {
2309      // delete and replace
2310      delete model2->clpMatrix();
2311      model2->replaceMatrix(saveMatrix);
2312    } else {
2313      delete saveMatrix;
2314    }
2315  }
2316  numberIterations = model2->numberIterations();
2317  finalStatus=model2->status();
2318  int finalSecondaryStatus = model2->secondaryStatus();
2319  if (presolve==ClpSolve::presolveOn) {
2320    int saveLevel = logLevel();
2321    if ((specialOptions_&1024)==0)
2322      setLogLevel(CoinMin(1,saveLevel));
2323    else
2324      setLogLevel(CoinMin(0,saveLevel));
2325    pinfo.postsolve(true);
2326    factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
2327    time2 = CoinCpuTime();
2328    timePresolve += time2-timeX;
2329    handler_->message(CLP_INTERVAL_TIMING,messages_)
2330      <<"Postsolve"<<time2-timeX<<time2-time1
2331      <<CoinMessageEol;
2332    timeX=time2;
2333    if (!presolveToFile)
2334      delete model2;
2335    if (interrupt)
2336      currentModel = this;
2337    // checkSolution(); already done by postSolve
2338    setLogLevel(saveLevel);
2339    if (finalStatus!=3&&(finalStatus||status()==-1)) {
2340      int savePerturbation = perturbation();
2341      if (!finalStatus||(moreSpecialOptions_&2)==0) {
2342        if (finalStatus==2) {
2343          // unbounded - get feasible first
2344          double save = optimizationDirection_;
2345          optimizationDirection_=0.0;
2346          primal(1);
2347          optimizationDirection_=save;
2348          primal(1);
2349        } else if (finalStatus==1) {
2350          dual();
2351        } else {
2352          setPerturbation(100);
2353          primal(1);
2354        }
2355      } else {
2356        // just set status
2357        problemStatus_=finalStatus;
2358      }
2359      setPerturbation(savePerturbation);
2360      numberIterations += numberIterations_;
2361      numberIterations_ = numberIterations;
2362      finalStatus=status();
2363      time2 = CoinCpuTime();
2364      handler_->message(CLP_INTERVAL_TIMING,messages_)
2365      <<"Cleanup"<<time2-timeX<<time2-time1
2366      <<CoinMessageEol;
2367      timeX=time2;
2368    } else {
2369      secondaryStatus_=finalSecondaryStatus;
2370    }
2371  } else if (model2!=this) {
2372    // not presolved - but different model used (sprint probably)
2373    CoinMemcpyN(model2->primalRowSolution(),
2374                numberRows_,this->primalRowSolution());
2375    CoinMemcpyN(model2->dualRowSolution(),
2376                numberRows_,this->dualRowSolution());
2377    CoinMemcpyN(model2->primalColumnSolution(),
2378                numberColumns_,this->primalColumnSolution());
2379    CoinMemcpyN(model2->dualColumnSolution(),
2380                numberColumns_,this->dualColumnSolution());
2381    CoinMemcpyN(model2->statusArray(),
2382                numberColumns_+numberRows_,this->statusArray());
2383    objectiveValue_=model2->objectiveValue_;
2384    numberIterations_ =model2->numberIterations_;
2385    problemStatus_ =model2->problemStatus_;
2386    secondaryStatus_ =model2->secondaryStatus_;
2387    delete model2;
2388  }
2389  setMaximumIterations(saveMaxIterations);
2390  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
2391                               "Errors","User stopped"};
2392  assert (finalStatus>=-1&&finalStatus<=5);
2393  handler_->message(CLP_TIMING,messages_)
2394    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
2395  handler_->printing(presolve==ClpSolve::presolveOn)
2396    <<timePresolve;
2397  handler_->printing(timeIdiot!=0.0)
2398    <<timeIdiot;
2399  handler_->message()<<CoinMessageEol;
2400  if (interrupt) 
2401    signal(SIGINT,saveSignal);
2402  perturbation_=savePerturbation;
2403  scalingFlag_=saveScaling;
2404  // If faking objective - put back correct one
2405  if (savedObjective) {
2406    delete objective_;
2407    objective_=savedObjective;
2408  }
2409  return finalStatus;
2410}
2411// General solve
2412int 
2413ClpSimplex::initialSolve()
2414{
2415  // Default so use dual
2416  ClpSolve options;
2417  return initialSolve(options);
2418}
2419// General dual solve
2420int 
2421ClpSimplex::initialDualSolve()
2422{
2423  ClpSolve options;
2424  // Use dual
2425  options.setSolveType(ClpSolve::useDual);
2426  return initialSolve(options);
2427}
2428// General dual solve
2429int 
2430ClpSimplex::initialPrimalSolve()
2431{
2432  ClpSolve options;
2433  // Use primal
2434  options.setSolveType(ClpSolve::usePrimal);
2435  return initialSolve(options);
2436}
2437// barrier solve, not to be followed by crossover
2438int 
2439ClpSimplex::initialBarrierNoCrossSolve()
2440{
2441  ClpSolve options;
2442  // Use primal
2443  options.setSolveType(ClpSolve::useBarrierNoCross);
2444  return initialSolve(options);
2445}
2446
2447// General barrier solve
2448int 
2449ClpSimplex::initialBarrierSolve()
2450{
2451  ClpSolve options;
2452  // Use primal
2453  options.setSolveType(ClpSolve::useBarrier);
2454  return initialSolve(options);
2455}
2456
2457// Default constructor
2458ClpSolve::ClpSolve (  )
2459{
2460  method_ = automatic;
2461  presolveType_=presolveOn;
2462  numberPasses_=5;
2463  int i;
2464  for (i=0;i<7;i++)
2465    options_[i]=0;
2466  // say no +-1 matrix
2467  options_[3]=1;
2468  for (i=0;i<7;i++)
2469    extraInfo_[i]=-1;
2470  independentOptions_[0]=0;
2471  // But switch off slacks
2472  independentOptions_[1]=512;
2473  // Substitute up to 3
2474  independentOptions_[2]=3;
2475 
2476}
2477// Constructor when you really know what you are doing
2478ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
2479             int numberPasses, int options[6],
2480             int extraInfo[6], int independentOptions[3])
2481{
2482  method_ = method;
2483  presolveType_=presolveType;
2484  numberPasses_=numberPasses;
2485  int i;
2486  for (i=0;i<6;i++)
2487    options_[i]=options[i];
2488  options_[6]=0;
2489  for (i=0;i<6;i++)
2490    extraInfo_[i]=extraInfo[i];
2491  extraInfo_[6]=0;
2492  for (i=0;i<3;i++)
2493    independentOptions_[i]=independentOptions[i];
2494}
2495
2496// Copy constructor.
2497ClpSolve::ClpSolve(const ClpSolve & rhs)
2498{
2499  method_ = rhs.method_;
2500  presolveType_=rhs.presolveType_;
2501  numberPasses_=rhs.numberPasses_;
2502  int i;
2503  for ( i=0;i<7;i++)
2504    options_[i]=rhs.options_[i];
2505  for ( i=0;i<7;i++)
2506    extraInfo_[i]=rhs.extraInfo_[i];
2507  for (i=0;i<3;i++)
2508    independentOptions_[i]=rhs.independentOptions_[i];
2509}
2510// Assignment operator. This copies the data
2511ClpSolve & 
2512ClpSolve::operator=(const ClpSolve & rhs)
2513{
2514  if (this != &rhs) {
2515    method_ = rhs.method_;
2516    presolveType_=rhs.presolveType_;
2517    numberPasses_=rhs.numberPasses_;
2518    int i;
2519    for (i=0;i<7;i++)
2520      options_[i]=rhs.options_[i];
2521    for (i=0;i<7;i++)
2522      extraInfo_[i]=rhs.extraInfo_[i];
2523    for (i=0;i<3;i++)
2524      independentOptions_[i]=rhs.independentOptions_[i];
2525  }
2526  return *this;
2527
2528}
2529// Destructor
2530ClpSolve::~ClpSolve (  )
2531{
2532}
2533// See header file for deatils
2534void 
2535ClpSolve::setSpecialOption(int which,int value,int extraInfo)
2536{
2537  options_[which]=value;
2538  extraInfo_[which]=extraInfo;
2539}
2540int 
2541ClpSolve::getSpecialOption(int which) const
2542{
2543  return options_[which];
2544}
2545
2546// Solve types
2547void 
2548ClpSolve::setSolveType(SolveType method, int extraInfo)
2549{
2550  method_=method;
2551}
2552
2553ClpSolve::SolveType
2554ClpSolve::getSolveType()
2555{
2556  return method_;
2557}
2558
2559// Presolve types
2560void 
2561ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
2562{
2563  presolveType_=amount;
2564  numberPasses_=extraInfo;
2565}
2566ClpSolve::PresolveType
2567ClpSolve::getPresolveType()
2568{
2569  return presolveType_;
2570}
2571// Extra info for idiot (or sprint)
2572int 
2573ClpSolve::getExtraInfo(int which) const
2574{
2575  return extraInfo_[which];
2576}
2577int 
2578ClpSolve::getPresolvePasses() const
2579{
2580  return numberPasses_;
2581}
2582/* Say to return at once if infeasible,
2583   default is to solve */
2584void 
2585ClpSolve::setInfeasibleReturn(bool trueFalse)
2586{
2587  independentOptions_[0]= trueFalse ? 1 : 0;
2588}
2589#include <string>
2590// Generates code for above constructor
2591void 
2592ClpSolve::generateCpp(FILE * fp)
2593{
2594  std::string solveType[] = {
2595    "ClpSolve::useDual",
2596    "ClpSolve::usePrimal",
2597    "ClpSolve::usePrimalorSprint",
2598    "ClpSolve::useBarrier",
2599    "ClpSolve::useBarrierNoCross",
2600    "ClpSolve::automatic",
2601    "ClpSolve::notImplemented"
2602  };
2603  std::string presolveType[] =  {
2604    "ClpSolve::presolveOn",
2605    "ClpSolve::presolveOff",
2606    "ClpSolve::presolveNumber",
2607    "ClpSolve::presolveNumberCost"
2608  };
2609  fprintf(fp,"3  ClpSolve::SolveType method = %s;\n",solveType[method_].c_str());
2610  fprintf(fp,"3  ClpSolve::PresolveType presolveType = %s;\n",
2611    presolveType[presolveType_].c_str());
2612  fprintf(fp,"3  int numberPasses = %d;\n",numberPasses_);
2613  fprintf(fp,"3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
2614    options_[0],options_[1],options_[2],
2615    options_[3],options_[4],options_[5]);
2616  fprintf(fp,"3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
2617    extraInfo_[0],extraInfo_[1],extraInfo_[2],
2618    extraInfo_[3],extraInfo_[4],extraInfo_[5]);
2619  fprintf(fp,"3  int independentOptions[] = {%d,%d,%d};\n",
2620    independentOptions_[0],independentOptions_[1],independentOptions_[2]);
2621  fprintf(fp,"3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
2622  fprintf(fp,"3                    options,extraInfo,independentOptions);\n");
2623}
2624//#############################################################################
2625#include "ClpNonLinearCost.hpp"
2626
2627ClpSimplexProgress::ClpSimplexProgress () 
2628{
2629  int i;
2630  for (i=0;i<CLP_PROGRESS;i++) {
2631    objective_[i] = COIN_DBL_MAX;
2632    infeasibility_[i] = -1.0; // set to an impossible value
2633    realInfeasibility_[i] = COIN_DBL_MAX;
2634    numberInfeasibilities_[i]=-1; 
2635    iterationNumber_[i]=-1;
2636  }
2637#ifdef CLP_PROGRESS_WEIGHT
2638  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2639    objectiveWeight_[i] = COIN_DBL_MAX;
2640    infeasibilityWeight_[i] = -1.0; // set to an impossible value
2641    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
2642    numberInfeasibilitiesWeight_[i]=-1; 
2643    iterationNumberWeight_[i]=-1;
2644  }
2645  drop_ =0.0;
2646  best_ =0.0;
2647#endif
2648  initialWeight_=0.0;
2649  for (i=0;i<CLP_CYCLE;i++) {
2650    //obj_[i]=COIN_DBL_MAX;
2651    in_[i]=-1;
2652    out_[i]=-1;
2653    way_[i]=0;
2654  }
2655  numberTimes_ = 0;
2656  numberBadTimes_ = 0;
2657  model_ = NULL;
2658  oddState_=0;
2659}
2660
2661
2662//-----------------------------------------------------------------------------
2663
2664ClpSimplexProgress::~ClpSimplexProgress ()
2665{
2666}
2667// Copy constructor.
2668ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
2669{
2670  int i;
2671  for (i=0;i<CLP_PROGRESS;i++) {
2672    objective_[i] = rhs.objective_[i];
2673    infeasibility_[i] = rhs.infeasibility_[i];
2674    realInfeasibility_[i] = rhs.realInfeasibility_[i];
2675    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2676    iterationNumber_[i]=rhs.iterationNumber_[i];
2677  }
2678#ifdef CLP_PROGRESS_WEIGHT
2679  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2680    objectiveWeight_[i] = rhs.objectiveWeight_[i];
2681    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2682    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2683    numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i]; 
2684    iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
2685  }
2686  drop_ = rhs.drop_;
2687  best_ = rhs.best_;
2688#endif
2689  initialWeight_ = rhs.initialWeight_;
2690  for (i=0;i<CLP_CYCLE;i++) {
2691    //obj_[i]=rhs.obj_[i];
2692    in_[i]=rhs.in_[i];
2693    out_[i]=rhs.out_[i];
2694    way_[i]=rhs.way_[i];
2695  }
2696  numberTimes_ = rhs.numberTimes_;
2697  numberBadTimes_ = rhs.numberBadTimes_;
2698  model_ = rhs.model_;
2699  oddState_=rhs.oddState_;
2700}
2701// Copy constructor.from model
2702ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
2703{
2704  model_ = model;
2705  reset();
2706  initialWeight_=0.0;
2707}
2708// Fill from model
2709void 
2710ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
2711{
2712  model_ = model;
2713  reset();
2714  initialWeight_=0.0;
2715}
2716// Assignment operator. This copies the data
2717ClpSimplexProgress & 
2718ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2719{
2720  if (this != &rhs) {
2721    int i;
2722    for (i=0;i<CLP_PROGRESS;i++) {
2723      objective_[i] = rhs.objective_[i];
2724      infeasibility_[i] = rhs.infeasibility_[i];
2725      realInfeasibility_[i] = rhs.realInfeasibility_[i];
2726      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2727      iterationNumber_[i]=rhs.iterationNumber_[i];
2728    }
2729#ifdef CLP_PROGRESS_WEIGHT
2730    for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2731      objectiveWeight_[i] = rhs.objectiveWeight_[i];
2732      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2733      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2734      numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i]; 
2735      iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
2736    }
2737    drop_ = rhs.drop_;
2738    best_ = rhs.best_;
2739#endif
2740    initialWeight_ = rhs.initialWeight_;
2741    for (i=0;i<CLP_CYCLE;i++) {
2742      //obj_[i]=rhs.obj_[i];
2743      in_[i]=rhs.in_[i];
2744      out_[i]=rhs.out_[i];
2745      way_[i]=rhs.way_[i];
2746    }
2747    numberTimes_ = rhs.numberTimes_;
2748    numberBadTimes_ = rhs.numberBadTimes_;
2749    model_ = rhs.model_;
2750    oddState_=rhs.oddState_;
2751  }
2752  return *this;
2753}
2754// Seems to be something odd about exact comparison of doubles on linux
2755static bool equalDouble(double value1, double value2)
2756{
2757
2758  union { double d; int i[2]; } v1,v2;
2759  v1.d = value1;
2760  v2.d = value2;
2761  if (sizeof(int)*2==sizeof(double)) 
2762    return (v1.i[0]==v2.i[0]&&v1.i[1]==v2.i[1]);
2763  else
2764    return (v1.i[0]==v2.i[0]);
2765}
2766int
2767ClpSimplexProgress::looping()
2768{
2769  if (!model_)
2770    return -1;
2771  double objective = model_->rawObjectiveValue();
2772  double infeasibility;
2773  double realInfeasibility=0.0;
2774  int numberInfeasibilities;
2775  int iterationNumber = model_->numberIterations();
2776  if (model_->algorithm()<0) {
2777    // dual
2778    infeasibility = model_->sumPrimalInfeasibilities();
2779    numberInfeasibilities = model_->numberPrimalInfeasibilities();
2780  } else {
2781    //primal
2782    infeasibility = model_->sumDualInfeasibilities();
2783    realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
2784    numberInfeasibilities = model_->numberDualInfeasibilities();
2785  }
2786  int i;
2787  int numberMatched=0;
2788  int matched=0;
2789  int nsame=0;
2790  for (i=0;i<CLP_PROGRESS;i++) {
2791    bool matchedOnObjective = equalDouble(objective,objective_[i]);
2792    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
2793    bool matchedOnInfeasibilities = 
2794      (numberInfeasibilities==numberInfeasibilities_[i]);
2795   
2796    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
2797      matched |= (1<<i);
2798      // Check not same iteration
2799      if (iterationNumber!=iterationNumber_[i]) {
2800        numberMatched++;
2801        // here mainly to get over compiler bug?
2802        if (model_->messageHandler()->logLevel()>10)
2803          printf("%d %d %d %d %d loop check\n",i,numberMatched,
2804                 matchedOnObjective, matchedOnInfeasibility, 
2805                 matchedOnInfeasibilities);
2806      } else {
2807        // stuck but code should notice
2808        nsame++;
2809      }
2810    }
2811    if (i) {
2812      objective_[i-1] = objective_[i];
2813      infeasibility_[i-1] = infeasibility_[i];
2814      realInfeasibility_[i-1] = realInfeasibility_[i];
2815      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
2816      iterationNumber_[i-1]=iterationNumber_[i];
2817    }
2818  }
2819  objective_[CLP_PROGRESS-1] = objective;
2820  infeasibility_[CLP_PROGRESS-1] = infeasibility;
2821  realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
2822  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
2823  iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
2824  if (nsame==CLP_PROGRESS)
2825    numberMatched=CLP_PROGRESS; // really stuck
2826  if (model_->progressFlag())
2827    numberMatched=0;
2828  numberTimes_++;
2829  if (numberTimes_<10)
2830    numberMatched=0;
2831  // skip if just last time as may be checking something
2832  if (matched==(1<<(CLP_PROGRESS-1)))
2833    numberMatched=0;
2834  if (numberMatched) {
2835    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
2836      <<numberMatched
2837      <<matched
2838      <<numberTimes_
2839      <<CoinMessageEol;
2840    numberBadTimes_++;
2841    if (numberBadTimes_<10) {
2842      // make factorize every iteration
2843      model_->forceFactorization(1);
2844      if (numberBadTimes_<2) {
2845        startCheck(); // clear other loop check
2846        if (model_->algorithm()<0) {
2847          // dual - change tolerance
2848          model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
2849          // if infeasible increase dual bound
2850          if (model_->dualBound()<1.0e17) {
2851            model_->setDualBound(model_->dualBound()*1.1);
2852          }
2853        } else {
2854          // primal - change tolerance
2855          if (numberBadTimes_>3)
2856            model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
2857          // if infeasible increase infeasibility cost
2858          if (model_->nonLinearCost()->numberInfeasibilities()&&
2859              model_->infeasibilityCost()<1.0e17) {
2860            model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
2861          }
2862        }
2863      } else {
2864        // flag
2865        int iSequence;
2866        if (model_->algorithm()<0) {
2867          // dual
2868          if (model_->dualBound()>1.0e14) 
2869            model_->setDualBound(1.0e14);
2870          iSequence=in_[CLP_CYCLE-1];
2871        } else {
2872          // primal
2873          if (model_->infeasibilityCost()>1.0e14) 
2874            model_->setInfeasibilityCost(1.0e14);
2875          iSequence=out_[CLP_CYCLE-1];
2876        }
2877        if (iSequence>=0) {
2878          char x = model_->isColumn(iSequence) ? 'C' :'R';
2879          if (model_->messageHandler()->logLevel()>=63)
2880            model_->messageHandler()->message(CLP_SIMPLEX_FLAG,model_->messages())
2881              <<x<<model_->sequenceWithin(iSequence)
2882              <<CoinMessageEol;
2883          // if Gub then needs to be sequenceIn_
2884          int save=model_->sequenceIn();
2885          model_->setSequenceIn(iSequence);
2886          model_->setFlagged(iSequence);
2887          model_->setSequenceIn(save);
2888          //printf("flagging %d from loop\n",iSequence);
2889          startCheck();
2890        } else {
2891          // Give up
2892          if (model_->messageHandler()->logLevel()>=63)
2893            printf("***** All flagged?\n");
2894          return 4;
2895        }
2896        // reset
2897        numberBadTimes_=2;
2898      }
2899      return -2;
2900    } else {
2901      // look at solution and maybe declare victory
2902      if (infeasibility<1.0e-4) {
2903        return 0;
2904      } else {
2905        model_->messageHandler()->message(CLP_LOOP,model_->messages())
2906          <<CoinMessageEol;
2907#ifndef NDEBUG
2908        printf("debug loop ClpSimplex A\n");
2909        abort();
2910#endif
2911        return 3;
2912      }
2913    }
2914  }
2915  return -1;
2916}
2917// Resets as much as possible
2918void 
2919ClpSimplexProgress::reset()
2920{
2921  int i;
2922  for (i=0;i<CLP_PROGRESS;i++) {
2923    if (model_->algorithm()>=0)
2924      objective_[i] = COIN_DBL_MAX;
2925    else
2926      objective_[i] = -COIN_DBL_MAX;
2927    infeasibility_[i] = -1.0; // set to an impossible value
2928    realInfeasibility_[i] = COIN_DBL_MAX;
2929    numberInfeasibilities_[i]=-1; 
2930    iterationNumber_[i]=-1;
2931  }
2932#ifdef CLP_PROGRESS_WEIGHT
2933  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2934    objectiveWeight_[i] = COIN_DBL_MAX;
2935    infeasibilityWeight_[i] = -1.0; // set to an impossible value
2936    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
2937    numberInfeasibilitiesWeight_[i]=-1; 
2938    iterationNumberWeight_[i]=-1;
2939  }
2940  drop_ =0.0;
2941  best_ =0.0;
2942#endif
2943  for (i=0;i<CLP_CYCLE;i++) {
2944    //obj_[i]=COIN_DBL_MAX;
2945    in_[i]=-1;
2946    out_[i]=-1;
2947    way_[i]=0;
2948  }
2949  numberTimes_ = 0;
2950  numberBadTimes_ = 0;
2951  oddState_=0;
2952}
2953// Returns previous objective (if -1) - current if (0)
2954double 
2955ClpSimplexProgress::lastObjective(int back) const
2956{
2957  return objective_[CLP_PROGRESS-1-back];
2958}
2959// Returns previous infeasibility (if -1) - current if (0)
2960double 
2961ClpSimplexProgress::lastInfeasibility(int back) const
2962{
2963  return realInfeasibility_[CLP_PROGRESS-1-back];
2964}
2965// Sets real primal infeasibility
2966void
2967ClpSimplexProgress::setInfeasibility(double value)
2968{
2969  for (int i=1;i<CLP_PROGRESS;i++) 
2970    realInfeasibility_[i-1] = realInfeasibility_[i];
2971  realInfeasibility_[CLP_PROGRESS-1]=value;
2972}
2973// Modify objective e.g. if dual infeasible in dual
2974void 
2975ClpSimplexProgress::modifyObjective(double value)
2976{
2977  objective_[CLP_PROGRESS-1]=value;
2978}
2979// Returns previous iteration number (if -1) - current if (0)
2980int 
2981ClpSimplexProgress::lastIterationNumber(int back) const
2982{
2983  return iterationNumber_[CLP_PROGRESS-1-back];
2984}
2985// clears iteration numbers (to switch off panic)
2986void 
2987ClpSimplexProgress::clearIterationNumbers()
2988{
2989  for (int i=0;i<CLP_PROGRESS;i++) 
2990    iterationNumber_[i]=-1;
2991}
2992// Start check at beginning of whileIterating
2993void 
2994ClpSimplexProgress::startCheck()
2995{
2996  int i;
2997  for (i=0;i<CLP_CYCLE;i++) {
2998    //obj_[i]=COIN_DBL_MAX;
2999    in_[i]=-1;
3000    out_[i]=-1;
3001    way_[i]=0;
3002  }
3003}
3004// Returns cycle length in whileIterating
3005int 
3006ClpSimplexProgress::cycle(int in, int out,int wayIn,int wayOut)
3007{
3008  int i;
3009#if 0
3010  if (model_->numberIterations()>206571) {
3011    printf("in %d out %d\n",in,out);
3012    for (i=0;i<CLP_CYCLE;i++)
3013      printf("cy %d in %d out %d\n",i,in_[i],out_[i]);
3014  }
3015#endif
3016  int matched=0;
3017  // first see if in matches any out
3018  for (i=1;i<CLP_CYCLE;i++) {
3019    if (in==out_[i]) {
3020      // even if flip then suspicious
3021      matched=-1;
3022      break;
3023    }
3024  }
3025#if 0
3026  if (!matched||in_[0]<0) {
3027    // can't be cycle
3028    for (i=0;i<CLP_CYCLE-1;i++) {
3029      //obj_[i]=obj_[i+1];
3030      in_[i]=in_[i+1];
3031      out_[i]=out_[i+1];
3032      way_[i]=way_[i+1];
3033    }
3034  } else {
3035    // possible cycle
3036    matched=0;
3037    for (i=0;i<CLP_CYCLE-1;i++) {
3038      int k;
3039      char wayThis = way_[i];
3040      int inThis = in_[i];
3041      int outThis = out_[i];
3042      //double objThis = obj_[i];
3043      for(k=i+1;k<CLP_CYCLE;k++) {
3044        if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
3045          int distance = k-i;
3046          if (k+distance<CLP_CYCLE) {
3047            // See if repeats
3048            int j=k+distance;
3049            if (inThis==in_[j]&&outThis==out_[j]&&wayThis==way_[j]) {
3050              matched=distance;
3051              break;
3052            }
3053          } else {
3054            matched=distance;
3055            break;
3056          }
3057        }
3058      }
3059      //obj_[i]=obj_[i+1];
3060      in_[i]=in_[i+1];
3061      out_[i]=out_[i+1];
3062      way_[i]=way_[i+1];
3063    }
3064  }
3065#else
3066  if (matched&&in_[0]>=0) {
3067    // possible cycle - only check [0] against all
3068    matched=0;
3069    int nMatched=0;
3070    char way0 = way_[0];
3071    int in0 = in_[0];
3072    int out0 = out_[0];
3073    //double obj0 = obj_[i];
3074    for(int k=1;k<CLP_CYCLE-4;k++) {
3075      if (in0==in_[k]&&out0==out_[k]&&way0==way_[k]) {
3076        nMatched++;
3077        // See if repeats
3078        int end = CLP_CYCLE-k;
3079        int j;
3080        for ( j=1;j<end;j++) {
3081          if (in_[j+k]!=in_[j]||out_[j+k]!=out_[j]||way_[j+k]!=way_[j])
3082            break;
3083        }
3084        if (j==end) {
3085          matched=k;
3086          break;
3087        }
3088      }
3089    }
3090    // If three times then that is too much even if not regular
3091    if (matched<=0&&nMatched>1)
3092      matched=100;
3093  }
3094  for (i=0;i<CLP_CYCLE-1;i++) {
3095    //obj_[i]=obj_[i+1];
3096    in_[i]=in_[i+1];
3097    out_[i]=out_[i+1];
3098    way_[i]=way_[i+1];
3099  }
3100#endif
3101  int way = 1-wayIn+4*(1-wayOut);
3102  //obj_[i]=model_->objectiveValue();
3103  in_[CLP_CYCLE-1]=in;
3104  out_[CLP_CYCLE-1]=out;
3105  way_[CLP_CYCLE-1]=static_cast<char>(way);
3106  return matched;
3107}
Note: See TracBrowser for help on using the repository browser.