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

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

changes for cholesky including code from Anshul Gupta

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 146.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#include "ClpSimplexDual.hpp"
18#ifndef SLIM_CLP
19#include "ClpQuadraticObjective.hpp"
20#include "ClpInterior.hpp"
21#include "ClpCholeskyDense.hpp"
22#include "ClpCholeskyBase.hpp"
23#include "ClpPlusMinusOneMatrix.hpp"
24#include "ClpNetworkMatrix.hpp"
25#endif
26#include "ClpLinearObjective.hpp"
27#include "ClpSolve.hpp"
28#include "ClpPackedMatrix.hpp"
29#include "ClpMessage.hpp"
30#include "CoinTime.hpp"
31
32#include "ClpPresolve.hpp"
33#ifndef SLIM_CLP
34#include "Idiot.hpp"
35#ifdef WSSMP_BARRIER
36#include "ClpCholeskyWssmp.hpp"
37#include "ClpCholeskyWssmpKKT.hpp"
38#endif
39#ifdef UFL_BARRIER
40#include "ClpCholeskyUfl.hpp"
41#endif
42#ifdef TAUCS_BARRIER
43#include "ClpCholeskyTaucs.hpp"
44#endif
45#ifdef MUMPS_BARRIER
46#include "ClpCholeskyMumps.hpp"
47#endif
48#ifdef COIN_HAS_VOL
49#include "VolVolume.hpp"
50#include "CoinHelperFunctions.hpp"
51#include "CoinPackedMatrix.hpp"
52#include "CoinMpsIO.hpp"
53
54//#############################################################################
55
56class lpHook : public VOL_user_hooks {
57private:
58   lpHook(const lpHook&);
59   lpHook& operator=(const lpHook&);
60private:
61   /// Pointer to dense vector of structural variable upper bounds
62   double  *colupper_;
63   /// Pointer to dense vector of structural variable lower bounds
64   double  *collower_;
65   /// Pointer to dense vector of objective coefficients
66   double  *objcoeffs_;
67   /// Pointer to dense vector of right hand sides
68   double  *rhs_;
69   /// Pointer to dense vector of senses
70   char    *sense_;
71
72   /// The problem matrix in a row ordered form
73   CoinPackedMatrix rowMatrix_;
74   /// The problem matrix in a column ordered form
75   CoinPackedMatrix colMatrix_;
76
77public:
78   lpHook(double* clb, double* cub, double* obj,
79          double* rhs, char* sense, const CoinPackedMatrix& mat);
80   virtual ~lpHook();
81   
82public:
83   // for all hooks: return value of -1 means that volume should quit
84   /** compute reduced costs   
85       @param u (IN) the dual variables
86       @param rc (OUT) the reduced cost with respect to the dual values
87   */
88   virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc);
89
90   /** Solve the subproblem for the subgradient step.
91       @param dual (IN) the dual variables
92       @param rc (IN) the reduced cost with respect to the dual values
93       @param lcost (OUT) the lagrangean cost with respect to the dual values
94       @param x (OUT) the primal result of solving the subproblem
95       @param v (OUT) b-Ax for the relaxed constraints
96       @param pcost (OUT) the primal objective value of <code>x</code>
97   */
98   virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
99                                double& lcost, VOL_dvector& x, VOL_dvector& v,
100                                double& pcost);
101   /** Starting from the primal vector x, run a heuristic to produce
102       an integer solution 
103       @param x (IN) the primal vector
104       @param heur_val (OUT) the value of the integer solution (return
105       <code>DBL_MAX</code> here if no feas sol was found
106   */
107   virtual int heuristics(const VOL_problem& p, 
108                          const VOL_dvector& x, double& heur_val) {
109      return 0;
110   }
111};
112 
113//#############################################################################
114
115lpHook::lpHook(double* clb, double* cub, double* obj,
116               double* rhs, char* sense,
117               const CoinPackedMatrix& mat)
118{
119   colupper_ = cub;
120   collower_ = clb;
121   objcoeffs_ = obj;
122   rhs_ = rhs;
123   sense_ = sense;
124   assert (mat.isColOrdered());
125   colMatrix_.copyOf(mat);
126   rowMatrix_.reverseOrderedCopyOf(mat);
127}
128
129//-----------------------------------------------------------------------------
130
131lpHook::~lpHook()
132{
133}
134
135//#############################################################################
136
137int
138lpHook::compute_rc(const VOL_dvector& u, VOL_dvector& rc)
139{
140   rowMatrix_.transposeTimes(u.v, rc.v);
141   const int psize = rowMatrix_.getNumCols();
142
143   for (int i = 0; i < psize; ++i)
144      rc[i] = objcoeffs_[i] - rc[i];
145   return 0;
146}
147
148//-----------------------------------------------------------------------------
149
150int
151lpHook::solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
152                         double& lcost, VOL_dvector& x, VOL_dvector& v,
153                         double& pcost)
154{
155   int i;
156   const int psize = x.size();
157   const int dsize = v.size();
158
159   // compute the lagrangean solution corresponding to the reduced costs
160   for (i = 0; i < psize; ++i) 
161      x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
162
163   // compute the lagrangean value (rhs*dual + primal*rc)
164   lcost = 0;
165   for (i = 0; i < dsize; ++i)
166      lcost += rhs_[i] * dual[i];
167   for (i = 0; i < psize; ++i)
168      lcost += x[i] * rc[i];
169
170   // compute the rhs - lhs
171   colMatrix_.times(x.v, v.v);
172   for (i = 0; i < dsize; ++i)
173      v[i] = rhs_[i] - v[i];
174
175   // compute the lagrangean primal objective
176   pcost = 0;
177   for (i = 0; i < psize; ++i)
178      pcost += x[i] * objcoeffs_[i];
179
180   return 0;
181}
182
183//#############################################################################
184/** A quick inlined function to convert from lb/ub style constraint
185    definition to sense/rhs/range style */
186inline void
187convertBoundToSense(const double lower, const double upper,
188                                        char& sense, double& right,
189                                        double& range) 
190{
191  range = 0.0;
192  if (lower > -1.0e20) {
193    if (upper < 1.0e20) {
194      right = upper;
195      if (upper==lower) {
196        sense = 'E';
197      } else {
198        sense = 'R';
199        range = upper - lower;
200      }
201    } else {
202      sense = 'G';
203      right = lower;
204    }
205  } else {
206    if (upper < 1.0e20) {
207      sense = 'L';
208      right = upper;
209    } else {
210      sense = 'N';
211      right = 0.0;
212    }
213  }
214}
215
216static int
217solveWithVolume(ClpSimplex * model, int numberPasses, int doIdiot)
218{
219   VOL_problem volprob;
220   volprob.parm.gap_rel_precision=0.00001;
221   volprob.parm.maxsgriters=3000;
222   if(numberPasses>3000) {
223     volprob.parm.maxsgriters=numberPasses;
224     volprob.parm.primal_abs_precision=0.0;
225     volprob.parm.minimum_rel_ascent=0.00001;
226   } else if (doIdiot>0) {
227     volprob.parm.maxsgriters=doIdiot;
228   }
229   if (model->logLevel()<2) 
230     volprob.parm.printflag=0;
231   else
232     volprob.parm.printflag=3;
233   const CoinPackedMatrix* mat = model->matrix();
234   int psize = model->numberColumns();
235   int dsize = model->numberRows();
236   char * sense = new char[dsize];
237   double * rhs = new double [dsize];
238
239   // Set the lb/ub on the duals
240   volprob.dsize = dsize;
241   volprob.psize = psize;
242   volprob.dual_lb.allocate(dsize);
243   volprob.dual_ub.allocate(dsize);
244   int i;
245   const double * rowLower = model->rowLower();
246   const double * rowUpper = model->rowUpper();
247   for (i = 0; i < dsize; ++i) {
248     double range;
249     convertBoundToSense(rowLower[i],rowUpper[i],
250                         sense[i],rhs[i],range);
251      switch (sense[i]) {
252       case 'E':
253         volprob.dual_lb[i] = -1.0e31;
254         volprob.dual_ub[i] = 1.0e31;
255         break;
256       case 'L':
257         volprob.dual_lb[i] = -1.0e31;
258         volprob.dual_ub[i] = 0.0;
259         break;
260       case 'G':
261         volprob.dual_lb[i] = 0.0;
262         volprob.dual_ub[i] = 1.0e31;
263         break;
264       default:
265         printf("Volume Algorithm can't work if there is a non ELG row\n");
266         return 1;
267      }
268   }
269   // Check bounds
270   double * saveLower = model->columnLower();
271   double * saveUpper = model->columnUpper();
272   bool good=true;
273   for (i=0;i<psize;i++) {
274     if (saveLower[i]<-1.0e20||saveUpper[i]>1.0e20) {
275       good=false;
276       break;
277     }
278   }
279   if (!good) {
280     saveLower = CoinCopyOfArray(model->columnLower(),psize);
281     saveUpper = CoinCopyOfArray(model->columnUpper(),psize);
282     for (i=0;i<psize;i++) {
283       if (saveLower[i]<-1.0e20)
284         saveLower[i]=-1.0e20;
285       if(saveUpper[i]>1.0e20) 
286         saveUpper[i]=1.0e20;
287     }
288   }
289   lpHook myHook(saveLower, saveUpper,
290                 model->objective(),
291                 rhs, sense, *mat);
292
293   volprob.solve(myHook, false /* no warmstart */);
294
295   if (saveLower!=model->columnLower()) {
296     delete [] saveLower;
297     delete [] saveUpper;
298   }
299   //------------- extract the solution ---------------------------
300
301   //printf("Best lagrangean value: %f\n", volprob.value);
302
303   double avg = 0;
304   for (i = 0; i < dsize; ++i) {
305      switch (sense[i]) {
306       case 'E':
307         avg += CoinAbs(volprob.viol[i]);
308         break;
309       case 'L':
310         if (volprob.viol[i] < 0)
311            avg +=  (-volprob.viol[i]);
312         break;
313       case 'G':
314         if (volprob.viol[i] > 0)
315            avg +=  volprob.viol[i];
316         break;
317      }
318   }
319     
320   //printf("Average primal constraint violation: %f\n", avg/dsize);
321
322   // volprob.dsol contains the dual solution (dual feasible)
323   // volprob.psol contains the primal solution
324   //              (NOT necessarily primal feasible)
325   CoinMemcpyN(volprob.dsol.v,dsize,model->dualRowSolution());
326   CoinMemcpyN(volprob.psol.v,psize,model->primalColumnSolution());
327   return 0;
328}
329#endif
330static ClpInterior * currentModel2 = NULL;
331#endif
332//#############################################################################
333// Allow for interrupts
334// But is this threadsafe ? (so switched off by option)
335
336#include "CoinSignal.hpp"
337static ClpSimplex * currentModel = NULL;
338
339extern "C" {
340   static void 
341#if defined(_MSC_VER)
342   __cdecl
343#endif // _MSC_VER
344   signal_handler(int whichSignal)
345   {
346      if (currentModel!=NULL) 
347         currentModel->setMaximumIterations(0); // stop at next iterations
348#ifndef SLIM_CLP
349      if (currentModel2!=NULL) 
350         currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
351#endif
352      return;
353   }
354}
355
356/** General solve algorithm which can do presolve
357    special options (bits)
358    1 - do not perturb
359    2 - do not scale
360    4 - use crash (default allslack in dual, idiot in primal)
361    8 - all slack basis in primal
362    16 - switch off interrupt handling
363    32 - do not try and make plus minus one matrix
364    64 - do not use sprint even if problem looks good
365 */
366int 
367ClpSimplex::initialSolve(ClpSolve & options)
368{
369  ClpSolve::SolveType method=options.getSolveType();
370  //ClpSolve::SolveType originalMethod=method;
371  ClpSolve::PresolveType presolve = options.getPresolveType();
372  int saveMaxIterations = maximumIterations();
373  int finalStatus=-1;
374  int numberIterations=0;
375  double time1 = CoinCpuTime();
376  double timeX = time1;
377  double time2;
378  ClpMatrixBase * saveMatrix=NULL;
379  ClpObjective * savedObjective=NULL;
380  if (!objective_||!matrix_) {
381    // totally empty
382    handler_->message(CLP_EMPTY_PROBLEM,messages_)
383      <<0
384      <<0
385      <<0
386      <<CoinMessageEol;
387    return -1;
388  } else if (!numberRows_||!numberColumns_||!getNumElements()) {
389    presolve = ClpSolve::presolveOff;
390  }
391  if (objective_->type()>=2&&optimizationDirection_==0) {
392    // pretend linear
393    savedObjective=objective_;
394    // make up objective
395    double * obj = new double[numberColumns_];
396    for (int i=0;i<numberColumns_;i++) {
397      double l = fabs(columnLower_[i]);
398      double u = fabs(columnUpper_[i]);
399      obj[i]=0.0;
400      if (CoinMin(l,u)<1.0e20) {
401        if (l<u) 
402          obj[i]=1.0+randomNumberGenerator_.randomDouble()*1.0e-2;
403        else
404          obj[i]=-1.0-randomNumberGenerator_.randomDouble()*1.0e-2;
405      }
406    }
407    objective_= new ClpLinearObjective(obj,numberColumns_);
408    delete [] obj;
409  }
410  ClpSimplex * model2 = this;
411  bool interrupt = (options.getSpecialOption(2)==0);
412  CoinSighandler_t saveSignal=static_cast<CoinSighandler_t> (0);
413  if (interrupt) {
414    currentModel = model2;
415    // register signal handler
416    saveSignal = signal(SIGINT,signal_handler);
417  }
418  // If no status array - set up basis
419  if (!status_)
420    allSlackBasis();
421  ClpPresolve pinfo;
422  pinfo.setSubstitution(options.substitution());
423  int presolveOptions = options.presolveActions();
424  bool presolveToFile = (presolveOptions&0x40000000)!=0;
425  presolveOptions &= ~0x40000000;
426  if ((presolveOptions&0xffff)!=0)
427    pinfo.setPresolveActions(presolveOptions);
428  // switch off singletons to slacks
429  //pinfo.setDoSingletonColumn(false); // done by bits
430  int printOptions = options.getSpecialOption(5);
431  if ((printOptions&1)!=0)
432    pinfo.statistics();
433  double timePresolve=0.0;
434  double timeIdiot=0.0;
435  double timeCore=0.0;
436  int savePerturbation=perturbation_;
437  int saveScaling = scalingFlag_;
438#ifndef SLIM_CLP
439#ifndef NO_RTTI
440  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
441    // network - switch off stuff
442    presolve = ClpSolve::presolveOff;
443  }
444#else
445  if (matrix_->type()==11) {
446    // network - switch off stuff
447    presolve = ClpSolve::presolveOff;
448  }
449#endif
450#endif
451  if (presolve!=ClpSolve::presolveOff) {
452    bool costedSlacks=false;
453    int numberPasses=5;
454    if (presolve==ClpSolve::presolveNumber) {
455      numberPasses=options.getPresolvePasses();
456      presolve = ClpSolve::presolveOn;
457    } else if (presolve==ClpSolve::presolveNumberCost) {
458      numberPasses=options.getPresolvePasses();
459      presolve = ClpSolve::presolveOn;
460      costedSlacks=true;
461      // switch on singletons to slacks
462      pinfo.setDoSingletonColumn(true);
463      // gub stuff for testing
464      //pinfo.setDoGubrow(true);
465    }
466#ifndef CLP_NO_STD
467    if (presolveToFile) {
468      // PreSolve to file - not fully tested
469      printf("Presolving to file - presolve.save\n");
470      pinfo.presolvedModelToFile(*this,"presolve.save",dblParam_[ClpPresolveTolerance],
471                           false,numberPasses);
472      model2=this;
473    } else {
474#endif
475      model2 = pinfo.presolvedModel(*this,dblParam_[ClpPresolveTolerance],
476                                    false,numberPasses,true,costedSlacks);
477#ifndef CLP_NO_STD
478    }
479#endif
480    time2 = CoinCpuTime();
481    timePresolve = time2-timeX;
482    handler_->message(CLP_INTERVAL_TIMING,messages_)
483      <<"Presolve"<<timePresolve<<time2-time1
484      <<CoinMessageEol;
485    timeX=time2;
486    if (!model2) {
487      handler_->message(CLP_INFEASIBLE,messages_)
488        <<CoinMessageEol;
489      model2 = this;
490      problemStatus_=1; // may be unbounded but who cares
491      if (options.infeasibleReturn()||(moreSpecialOptions_&1)!=0) {
492        return -1;
493      }
494      presolve=ClpSolve::presolveOff;
495    }
496    // We may be better off using original (but if dual leave because of bounds)
497    if (presolve!=ClpSolve::presolveOff&&
498        numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_
499        &&model2!=this) {
500      if(method!=ClpSolve::useDual||
501         (numberRows_==model2->numberRows_&&numberColumns_==model2->numberColumns_)) {
502        delete model2;
503        model2 = this;
504        presolve=ClpSolve::presolveOff;
505      }
506    }
507  }
508  if (interrupt)
509    currentModel = model2;
510  // For below >0 overrides
511  // 0 means no, -1 means maybe
512  int doIdiot=0;
513  int doCrash=0;
514  int doSprint=0;
515  int doSlp=0;
516  int primalStartup=1;
517  bool tryItSave = false;
518  // switch to primal from automatic if just one cost entry
519  if (method==ClpSolve::automatic&&model2->numberColumns()>5000&&
520      (specialOptions_&1024)!=0) {
521    int numberColumns = model2->numberColumns();
522    int numberRows = model2->numberRows();
523    const double * obj = model2->objective();
524    int nNon=0;
525    for (int i=0;i<numberColumns;i++) {
526      if (obj[i])
527        nNon++;
528    }
529    if (nNon==1) {
530#ifdef COIN_DEVELOP
531      printf("Forcing primal\n");
532#endif
533      method=ClpSolve::usePrimal;
534      tryItSave= numberRows>200&&numberColumns>2000&&
535        (numberColumns>2*numberRows||(specialOptions_&1024)!=0);
536    }
537  }
538  if (method!=ClpSolve::useDual&&method!=ClpSolve::useBarrier
539      &&method!=ClpSolve::useBarrierNoCross) {
540    switch (options.getSpecialOption(1)) {
541    case 0:
542      doIdiot=-1;
543      doCrash=-1;
544      doSprint=-1;
545      break;
546    case 1:
547      doIdiot=0;
548      doCrash=1;
549      if (options.getExtraInfo(1)>0)
550        doCrash = options.getExtraInfo(1);
551      doSprint=0;
552      break;
553    case 2:
554      doIdiot=1;
555      if (options.getExtraInfo(1)>0)
556        doIdiot = options.getExtraInfo(1);
557      doCrash=0;
558      doSprint=0;
559      break;
560    case 3:
561      doIdiot=0;
562      doCrash=0;
563      doSprint=1;
564      break;
565    case 4:
566      doIdiot=0;
567      doCrash=0;
568      doSprint=0;
569      break;
570    case 5:
571      doIdiot=0;
572      doCrash=-1;
573      doSprint=-1;
574      break;
575    case 6:
576      doIdiot=-1;
577      doCrash=-1;
578      doSprint=0;
579      break;
580    case 7:
581      doIdiot=-1;
582      doCrash=0;
583      doSprint=-1;
584      break;
585    case 8:
586      doIdiot=-1;
587      doCrash=0;
588      doSprint=0;
589      break;
590    case 9:
591      doIdiot=0;
592      doCrash=0;
593      doSprint=-1;
594      break;
595    case 10:
596      doIdiot=0;
597      doCrash=0;
598      doSprint=0;
599      if (options.getExtraInfo(1)>0)
600        doSlp = options.getExtraInfo(1);
601      break;
602    case 11:
603      doIdiot=0;
604      doCrash=0;
605      doSprint=0;
606      primalStartup=0;
607      break;
608    default:
609      abort();
610    }
611  } else {
612    // Dual
613    switch (options.getSpecialOption(0)) {
614    case 0:
615      doIdiot=0;
616      doCrash=0;
617      doSprint=0;
618      break;
619    case 1:
620      doIdiot=0;
621      doCrash=1;
622      if (options.getExtraInfo(0)>0)
623        doCrash = options.getExtraInfo(0);
624      doSprint=0;
625      break;
626    case 2:
627      doIdiot=-1;
628      if (options.getExtraInfo(0)>0)
629        doIdiot = options.getExtraInfo(0);
630      doCrash=0;
631      doSprint=0;
632      break;
633    default:
634      abort();
635    }
636  }
637#ifndef NO_RTTI
638  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objectiveAsObject()));
639#else
640  ClpQuadraticObjective * quadraticObj = NULL;
641  if (objective_->type()==2)
642    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
643#endif
644  // If quadratic then primal or barrier or slp
645  if (quadraticObj) {
646    doSprint=0;
647    doIdiot=0;
648    // off
649    if (method==ClpSolve::useBarrier)
650      method=ClpSolve::useBarrierNoCross;
651    else if (method!=ClpSolve::useBarrierNoCross)
652      method=ClpSolve::usePrimal;
653  }
654#ifdef COIN_HAS_VOL
655  // Save number of idiot
656  int saveDoIdiot=doIdiot;
657#endif
658  // Just do this number of passes in Sprint
659  int maxSprintPass=100;
660  // See if worth trying +- one matrix
661  bool plusMinus=false;
662  int numberElements=model2->getNumElements();
663#ifndef SLIM_CLP
664#ifndef NO_RTTI
665  if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
666    // network - switch off stuff
667    doIdiot=0;
668    if (doSprint<0)
669      doSprint=0;
670  }
671#else
672  if (matrix_->type()==11) {
673    // network - switch off stuff
674    doIdiot=0;
675    //doSprint=0;
676  }
677#endif
678#endif
679  int numberColumns = model2->numberColumns();
680  int numberRows = model2->numberRows();
681  // If not all slack basis - switch off all except sprint
682  int numberRowsBasic=0;
683  int iRow;
684  for (iRow=0;iRow<numberRows;iRow++)
685    if (model2->getRowStatus(iRow)==basic)
686      numberRowsBasic++;
687  if (numberRowsBasic<numberRows) {
688    doIdiot=0;
689    doCrash=0;
690    //doSprint=0;
691  }
692  if (options.getSpecialOption(3)==0) {
693    if(numberElements>100000)
694      plusMinus=true;
695    if(numberElements>10000&&(doIdiot||doSprint)) 
696      plusMinus=true;
697  } else if ((specialOptions_&1024)!=0) {
698    plusMinus=true;
699  }
700#ifndef SLIM_CLP
701  // Statistics (+1,-1, other) - used to decide on strategy if not +-1
702  CoinBigIndex statistics[3]={-1,0,0};
703  if (plusMinus) {
704    saveMatrix = model2->clpMatrix();
705#ifndef NO_RTTI
706    ClpPackedMatrix* clpMatrix =
707      dynamic_cast< ClpPackedMatrix*>(saveMatrix);
708#else
709    ClpPackedMatrix* clpMatrix = NULL;
710    if (saveMatrix->type()==1)
711      clpMatrix =
712        static_cast< ClpPackedMatrix*>(saveMatrix);
713#endif
714    if (clpMatrix) {
715      ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
716      if (newMatrix->getIndices()) {
717        if ((specialOptions_&1024)==0) {
718          model2->replaceMatrix(newMatrix);
719        } else {
720          // in integer - just use for sprint/idiot
721          saveMatrix=NULL;
722          delete newMatrix;
723        }
724      } else {
725        handler_->message(CLP_MATRIX_CHANGE,messages_)
726          <<"+- 1"
727          <<CoinMessageEol;
728        CoinMemcpyN(newMatrix->startPositive(),3,statistics);
729        saveMatrix=NULL;
730        plusMinus=false;
731        delete newMatrix;
732      }
733    } else {
734      saveMatrix=NULL;
735      plusMinus=false;
736    }
737  }
738#endif
739  if (this->factorizationFrequency()==200) {
740    // User did not touch preset
741    model2->defaultFactorizationFrequency();
742  } else if (model2!=this) {
743    // make sure model2 has correct value
744    model2->setFactorizationFrequency(this->factorizationFrequency());
745  }
746  if (method==ClpSolve::automatic) {
747    if (doSprint==0&&doIdiot==0) {
748      // off
749      method=ClpSolve::useDual;
750    } else {
751      // only do primal if sprint or idiot
752      if (doSprint>0) {
753        method=ClpSolve::usePrimalorSprint;
754      } else if (doIdiot>0) {
755        method=ClpSolve::usePrimal;
756      } else {
757        if (numberElements<500000) {
758          // Small problem
759          if(numberRows*10>numberColumns||numberColumns<6000
760             ||(numberRows*20>numberColumns&&!plusMinus))
761            doSprint=0; // switch off sprint
762        } else {
763          // larger problem
764          if(numberRows*8>numberColumns)
765            doSprint=0; // switch off sprint
766        }
767        // switch off sprint or idiot if any free variable
768        int iColumn;
769        double * columnLower = model2->columnLower();
770        double * columnUpper = model2->columnUpper();
771        for (iColumn=0;iColumn<numberColumns;iColumn++) {
772          if (columnLower[iColumn]<-1.0e10&&columnUpper[iColumn]>1.0e10) {
773            doSprint=0;
774            doIdiot=0;
775            break;
776          }
777        }
778        int nPasses=0;
779        // look at rhs
780        int iRow;
781        double largest=0.0;
782        double smallest = 1.0e30;
783        double largestGap=0.0;
784        int numberNotE=0;
785        bool notInteger=false;
786        for (iRow=0;iRow<numberRows;iRow++) {
787          double value1 = model2->rowLower_[iRow];
788          if (value1&&value1>-1.0e31) {
789            largest = CoinMax(largest,fabs(value1));
790            smallest=CoinMin(smallest,fabs(value1));
791            if (fabs(value1-floor(value1+0.5))>1.0e-8) {
792              notInteger=true;
793              break;
794            }
795          }
796          double value2 = model2->rowUpper_[iRow];
797          if (value2&&value2<1.0e31) {
798            largest = CoinMax(largest,fabs(value2));
799            smallest=CoinMin(smallest,fabs(value2));
800            if (fabs(value2-floor(value2+0.5))>1.0e-8) {
801              notInteger=true;
802              break;
803            }
804          }
805          if (value2>value1) {
806            numberNotE++;
807            if (value2>1.0e31||value1<-1.0e31)
808              largestGap = COIN_DBL_MAX;
809            else
810              largestGap = value2-value1;
811          }
812        }
813        bool tryIt= numberRows>200&&numberColumns>2000&&
814          (numberColumns>2*numberRows||(method!=ClpSolve::useDual&&(specialOptions_&1024)!=0));
815        tryItSave = tryIt;
816        if (numberRows<1000&&numberColumns<3000)
817          tryIt=false;
818        if (notInteger)
819          tryIt=false;
820        if (largest/smallest>10||(largest/smallest>2.0&&largest>50))
821          tryIt=false;
822        if (tryIt) {
823          if (largest/smallest>2.0) {
824            nPasses = 10+numberColumns/100000;
825            nPasses = CoinMin(nPasses,50);
826            nPasses = CoinMax(nPasses,15);
827            if (numberRows>20000&&nPasses>5) {
828              // Might as well go for it
829              nPasses = CoinMax(nPasses,71);
830            } else if (numberRows>2000&&nPasses>5) {
831              nPasses = CoinMax(nPasses,50);
832            } else if (numberElements<3*numberColumns) {
833              nPasses=CoinMin(nPasses,10); // probably not worh it
834            }
835          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
836            nPasses = 10+numberColumns/1000;
837            nPasses = CoinMin(nPasses,100);
838            nPasses = CoinMax(nPasses,30);
839            if (numberRows>25000) {
840              // Might as well go for it
841              nPasses = CoinMax(nPasses,71);
842            }
843            if (!largestGap)
844              nPasses *= 2;
845          } else {
846            nPasses = 10+numberColumns/1000;
847            nPasses = CoinMin(nPasses,200);
848            nPasses = CoinMax(nPasses,100);
849            if (!largestGap)
850              nPasses *= 2;
851          }
852        }
853        //printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
854        //     numberRows,numberColumns,plusMinus ? 'Y' : 'N',
855        //     tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
856        //exit(0);
857        if (!tryIt||nPasses<=5)
858          doIdiot=0;
859        if (doSprint) {
860          method = ClpSolve::usePrimalorSprint;
861        } else if (doIdiot) {
862          method = ClpSolve::usePrimal;
863        } else {
864          method = ClpSolve::useDual;
865        }
866      }
867    }
868  }
869  if (method==ClpSolve::usePrimalorSprint) {
870    if (doSprint<0) { 
871      if (numberElements<500000) {
872        // Small problem
873        if(numberRows*10>numberColumns||numberColumns<6000
874           ||(numberRows*20>numberColumns&&!plusMinus))
875          method=ClpSolve::usePrimal; // switch off sprint
876      } else {
877        // larger problem
878        if(numberRows*8>numberColumns)
879          method=ClpSolve::usePrimal; // switch off sprint
880        // but make lightweight
881        if(numberRows*10>numberColumns||numberColumns<6000
882           ||(numberRows*20>numberColumns&&!plusMinus))
883          maxSprintPass=10;
884      }
885    } else if (doSprint==0) {
886      method=ClpSolve::usePrimal; // switch off sprint
887    }
888  }
889  if (method==ClpSolve::useDual) {
890    double * saveLower=NULL;
891    double * saveUpper=NULL;
892    if (presolve==ClpSolve::presolveOn) {
893      int numberInfeasibilities = model2->tightenPrimalBounds(0.0,0);
894      if (numberInfeasibilities) {
895        handler_->message(CLP_INFEASIBLE,messages_)
896          <<CoinMessageEol;
897        model2 = this;
898        presolve=ClpSolve::presolveOff;
899      }
900    } else if (numberRows_+numberColumns_>5000) {
901      // do anyway
902      saveLower = new double[numberRows_+numberColumns_];
903      CoinMemcpyN(model2->columnLower(),numberColumns_,saveLower);
904      CoinMemcpyN(model2->rowLower(),numberRows_,saveLower+numberColumns_);
905      saveUpper = new double[numberRows_+numberColumns_];
906      CoinMemcpyN(model2->columnUpper(),numberColumns_,saveUpper);
907      CoinMemcpyN(model2->rowUpper(),numberRows_,saveUpper+numberColumns_);
908      int numberInfeasibilities = model2->tightenPrimalBounds();
909      if (numberInfeasibilities) {
910        handler_->message(CLP_INFEASIBLE,messages_)
911          <<CoinMessageEol;
912        CoinMemcpyN(saveLower,numberColumns_,model2->columnLower());
913        CoinMemcpyN(saveLower+numberColumns_,numberRows_,model2->rowLower());
914        delete [] saveLower;
915        saveLower=NULL;
916        CoinMemcpyN(saveUpper,numberColumns_,model2->columnUpper());
917        CoinMemcpyN(saveUpper+numberColumns_,numberRows_,model2->rowUpper());
918        delete [] saveUpper;
919        saveUpper=NULL;
920      }
921    }
922#ifndef COIN_HAS_VOL
923    // switch off idiot and volume for now
924    doIdiot=0; 
925#endif
926    // pick up number passes
927    int nPasses=0;
928    int numberNotE=0;
929#ifndef SLIM_CLP
930    if ((doIdiot<0&&plusMinus)||doIdiot>0) {
931      // See if candidate for idiot
932      nPasses=0;
933      Idiot info(*model2);
934      // Get average number of elements per column
935      double ratio  = static_cast<double> (numberElements)/static_cast<double> (numberColumns);
936      // look at rhs
937      int iRow;
938      double largest=0.0;
939      double smallest = 1.0e30;
940      double largestGap=0.0;
941      for (iRow=0;iRow<numberRows;iRow++) {
942        double value1 = model2->rowLower_[iRow];
943        if (value1&&value1>-1.0e31) {
944          largest = CoinMax(largest,fabs(value1));
945          smallest=CoinMin(smallest,fabs(value1));
946        }
947        double value2 = model2->rowUpper_[iRow];
948        if (value2&&value2<1.0e31) {
949          largest = CoinMax(largest,fabs(value2));
950          smallest=CoinMin(smallest,fabs(value2));
951        }
952        if (value2>value1) {
953          numberNotE++;
954          if (value2>1.0e31||value1<-1.0e31)
955            largestGap = COIN_DBL_MAX;
956          else
957            largestGap = value2-value1;
958        }
959      }
960      if (doIdiot<0) {
961        if (numberRows>200&&numberColumns>5000&&ratio>=3.0&&
962            largest/smallest<1.1&&!numberNotE) {
963          nPasses = 71;
964        }
965      } 
966      if (doIdiot>0) {
967        nPasses=CoinMax(nPasses,doIdiot);
968        if (nPasses>70) {
969          info.setStartingWeight(1.0e3);
970          info.setDropEnoughFeasibility(0.01);
971        }
972      }
973      if (nPasses>20) {
974#ifdef COIN_HAS_VOL
975        int returnCode = solveWithVolume(model2,nPasses,saveDoIdiot);
976        if (!returnCode) {
977          time2 = CoinCpuTime();
978          timeIdiot = time2-timeX;
979          handler_->message(CLP_INTERVAL_TIMING,messages_)
980            <<"Idiot Crash"<<timeIdiot<<time2-time1
981            <<CoinMessageEol;
982          timeX=time2;
983        } else {
984          nPasses=0;
985        }
986#else
987        nPasses=0;
988#endif
989      } else {
990        nPasses=0;
991      }
992    }
993#endif
994    if (doCrash) {
995      switch(doCrash) {
996        // standard
997      case 1:
998        model2->crash(1000,1);
999        break;
1000        // As in paper by Solow and Halim (approx)
1001      case 2:
1002      case 3:
1003        model2->crash(model2->dualBound(),0);
1004        break;
1005        // Just put free in basis
1006      case 4:
1007        model2->crash(0.0,3);
1008        break;
1009      }
1010    }
1011    if (!nPasses) {
1012      int saveOptions = model2->specialOptions();
1013      if (model2->numberRows()>100)
1014        model2->setSpecialOptions(saveOptions|64); // go as far as possible
1015      //int numberRows = model2->numberRows();
1016      //int numberColumns = model2->numberColumns();
1017      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1018        // See if original wanted vector
1019        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1020        ClpMatrixBase * matrix = model2->clpMatrix();
1021        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
1022          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1023          clpMatrix->makeSpecialColumnCopy();
1024          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
1025          model2->dual(0);
1026          clpMatrix->releaseSpecialColumnCopy();
1027        } else {
1028          model2->dual(0);
1029        }
1030      } else {
1031        model2->dual(0);
1032      }
1033    } else if (!numberNotE&&0) {
1034      // E so we can do in another way
1035      double * pi = model2->dualRowSolution();
1036      int i;
1037      int numberColumns = model2->numberColumns();
1038      int numberRows = model2->numberRows();
1039      double * saveObj = new double[numberColumns];
1040      CoinMemcpyN(model2->objective(),numberColumns,saveObj);
1041      CoinMemcpyN(model2->objective(),
1042             numberColumns,model2->dualColumnSolution());
1043      model2->clpMatrix()->transposeTimes(-1.0,pi,model2->dualColumnSolution());
1044      CoinMemcpyN(model2->dualColumnSolution(),
1045             numberColumns,model2->objective());
1046      const double * rowsol = model2->primalRowSolution();
1047      double offset=0.0;
1048      for (i=0;i<numberRows;i++) {
1049        offset += pi[i]*rowsol[i];
1050      }
1051      double value2;
1052      model2->getDblParam(ClpObjOffset,value2);
1053      //printf("Offset %g %g\n",offset,value2);
1054      model2->setDblParam(ClpObjOffset,value2-offset);
1055      model2->setPerturbation(51);
1056      //model2->setRowObjective(pi);
1057      // zero out pi
1058      //memset(pi,0,numberRows*sizeof(double));
1059      // Could put some in basis - only partially tested
1060      model2->allSlackBasis(); 
1061      //model2->factorization()->maximumPivots(200);
1062      //model2->setLogLevel(63);
1063      // solve
1064      model2->dual(0);
1065      model2->setDblParam(ClpObjOffset,value2);
1066      CoinMemcpyN(saveObj,numberColumns,model2->objective());
1067      // zero out pi
1068      //memset(pi,0,numberRows*sizeof(double));
1069      //model2->setRowObjective(pi);
1070      delete [] saveObj;
1071      //model2->dual(0);
1072      model2->setPerturbation(50);
1073      model2->primal();
1074    } else {
1075      // solve
1076      model2->setPerturbation(100);
1077      model2->dual(2);
1078      model2->setPerturbation(50);
1079      model2->dual(0);
1080    }
1081    if (saveLower) {
1082      CoinMemcpyN(saveLower,numberColumns_,model2->columnLower());
1083      CoinMemcpyN(saveLower+numberColumns_,numberRows_,model2->rowLower());
1084      delete [] saveLower;
1085      saveLower=NULL;
1086      CoinMemcpyN(saveUpper,numberColumns_,model2->columnUpper());
1087      CoinMemcpyN(saveUpper+numberColumns_,numberRows_,model2->rowUpper());
1088      delete [] saveUpper;
1089      saveUpper=NULL;
1090    }
1091    time2 = CoinCpuTime();
1092    timeCore = time2-timeX;
1093    handler_->message(CLP_INTERVAL_TIMING,messages_)
1094      <<"Dual"<<timeCore<<time2-time1
1095      <<CoinMessageEol;
1096    timeX=time2;
1097  } else if (method==ClpSolve::usePrimal) {
1098#ifndef SLIM_CLP
1099    if (doIdiot) {
1100      int nPasses=0;
1101      Idiot info(*model2);
1102      // Get average number of elements per column
1103      double ratio  = static_cast<double> (numberElements)/static_cast<double> (numberColumns);
1104      // look at rhs
1105      int iRow;
1106      double largest=0.0;
1107      double smallest = 1.0e30;
1108      double largestGap=0.0;
1109      int numberNotE=0;
1110      for (iRow=0;iRow<numberRows;iRow++) {
1111        double value1 = model2->rowLower_[iRow];
1112        if (value1&&value1>-1.0e31) {
1113          largest = CoinMax(largest,fabs(value1));
1114          smallest=CoinMin(smallest,fabs(value1));
1115        }
1116        double value2 = model2->rowUpper_[iRow];
1117        if (value2&&value2<1.0e31) {
1118          largest = CoinMax(largest,fabs(value2));
1119          smallest=CoinMin(smallest,fabs(value2));
1120        }
1121        if (value2>value1) {
1122          numberNotE++;
1123          if (value2>1.0e31||value1<-1.0e31)
1124            largestGap = COIN_DBL_MAX;
1125          else
1126            largestGap = value2-value1;
1127        }
1128      }
1129      bool increaseSprint=plusMinus;
1130      if ((specialOptions_&1024)!=0)
1131        increaseSprint=false;
1132      if (!plusMinus) {
1133        // If 90% +- 1 then go for sprint
1134        if (statistics[0]>=0&&10*statistics[2]<statistics[0]+statistics[1])
1135          increaseSprint=true;
1136      }
1137      bool tryIt= tryItSave;
1138      if (numberRows<1000&&numberColumns<3000)
1139        tryIt=false;
1140      if (tryIt) {
1141        if (increaseSprint) {
1142          info.setStartingWeight(1.0e3);
1143          info.setReduceIterations(6);
1144          // also be more lenient on infeasibilities
1145          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
1146          info.setDropEnoughWeighted(-2.0);
1147          if (largest/smallest>2.0) {
1148            nPasses = 10+numberColumns/100000;
1149            nPasses = CoinMin(nPasses,50);
1150            nPasses = CoinMax(nPasses,15);
1151            if (numberRows>20000&&nPasses>5) {
1152              // Might as well go for it
1153              nPasses = CoinMax(nPasses,71);
1154            } else if (numberRows>2000&&nPasses>5) {
1155              nPasses = CoinMax(nPasses,50);
1156            } else if (numberElements<3*numberColumns) {
1157              nPasses=CoinMin(nPasses,10); // probably not worh it
1158              if (doIdiot<0)
1159                info.setLightweight(1); // say lightweight idiot
1160            } else {
1161              if (doIdiot<0)
1162                info.setLightweight(1); // say lightweight idiot
1163            }
1164          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
1165            nPasses = 10+numberColumns/1000;
1166            nPasses = CoinMin(nPasses,100);
1167            nPasses = CoinMax(nPasses,30);
1168            if (numberRows>25000) {
1169              // Might as well go for it
1170              nPasses = CoinMax(nPasses,71);
1171            }
1172            if (!largestGap)
1173              nPasses *= 2;
1174          } else {
1175            nPasses = 10+numberColumns/1000;
1176            nPasses = CoinMin(nPasses,200);
1177            nPasses = CoinMax(nPasses,100);
1178            info.setStartingWeight(1.0e-1);
1179            info.setReduceIterations(6);
1180            if (!largestGap)
1181              nPasses *= 2;
1182            //info.setFeasibilityTolerance(1.0e-7);
1183          }
1184          // If few passes - don't bother
1185          if (nPasses<=5&&!plusMinus)
1186            nPasses=0;
1187        } else {
1188          if (doIdiot<0)
1189            info.setLightweight(1); // say lightweight idiot
1190          if (largest/smallest>1.01||numberNotE||statistics[2]>statistics[0]+statistics[1]) {
1191            if (numberRows>25000||numberColumns>5*numberRows) {
1192              nPasses = 50;
1193            } else if (numberColumns>4*numberRows) {
1194              nPasses = 20;
1195            } else {
1196              nPasses=5;
1197            }
1198          } else {
1199            if (numberRows>25000||numberColumns>5*numberRows) {
1200              nPasses = 50;
1201              info.setLightweight(0); // say not lightweight idiot
1202            } else if (numberColumns>4*numberRows) {
1203              nPasses = 20;
1204            } else {
1205              nPasses=15;
1206            }
1207          }
1208          if (ratio<3.0) { 
1209            nPasses=static_cast<int> (ratio*static_cast<double> (nPasses)/4.0); // probably not worth it
1210          } else {
1211            nPasses = CoinMax(nPasses,5);
1212          }
1213          if (numberRows>25000&&nPasses>5) {
1214            // Might as well go for it
1215            nPasses = CoinMax(nPasses,71);
1216          } else if (increaseSprint) {
1217            nPasses *= 2;
1218            nPasses=CoinMin(nPasses,71);
1219          } else if (nPasses==5&&ratio>5.0) {
1220            nPasses = static_cast<int> (static_cast<double> (nPasses)*(ratio/5.0)); // increase if lots of elements per column
1221          }
1222          if (nPasses<=5&&!plusMinus)
1223            nPasses=0;
1224          //info.setStartingWeight(1.0e-1);
1225        }
1226      }
1227      if (doIdiot>0) {
1228        // pick up number passes
1229        nPasses=options.getExtraInfo(1);
1230        if (nPasses>70) {
1231          info.setStartingWeight(1.0e3);
1232          info.setReduceIterations(6);
1233          if (nPasses>=5000) {
1234            int k= nPasses&100;
1235            nPasses /= 100;
1236            info.setReduceIterations(3);
1237            if (k)
1238              info.setStartingWeight(1.0e2);
1239          }
1240          // also be more lenient on infeasibilities
1241          info.setDropEnoughFeasibility(0.5*info.getDropEnoughFeasibility());
1242          info.setDropEnoughWeighted(-2.0);
1243        } else if (nPasses>=50) {
1244          info.setStartingWeight(1.0e3);
1245          //info.setReduceIterations(6);
1246        } 
1247        // For experimenting
1248        if (nPasses<70&&(nPasses%10)>0&&(nPasses%10)<4) {
1249          info.setStartingWeight(1.0e3);
1250          info.setLightweight(nPasses%10); // special testing
1251#ifdef COIN_DEVELOP
1252          printf("warning - odd lightweight %d\n",nPasses%10);
1253          //info.setReduceIterations(6);
1254#endif
1255        }
1256      }
1257      if (nPasses) {
1258        doCrash=0;
1259#if 0
1260        double * solution = model2->primalColumnSolution();
1261        int iColumn;
1262        double * saveLower = new double[numberColumns];
1263        CoinMemcpyN(model2->columnLower(),numberColumns,saveLower);
1264        double * saveUpper = new double[numberColumns];
1265        CoinMemcpyN(model2->columnUpper(),numberColumns,saveUpper);
1266        printf("doing tighten before idiot\n");
1267        model2->tightenPrimalBounds();
1268        // Move solution
1269        double * columnLower = model2->columnLower();
1270        double * columnUpper = model2->columnUpper();
1271        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1272          if (columnLower[iColumn]>0.0)
1273            solution[iColumn]=columnLower[iColumn];
1274          else if (columnUpper[iColumn]<0.0)
1275            solution[iColumn]=columnUpper[iColumn];
1276          else
1277            solution[iColumn]=0.0;
1278        }
1279        CoinMemcpyN(saveLower,numberColumns,columnLower);
1280        CoinMemcpyN(saveUpper,numberColumns,columnUpper);
1281        delete [] saveLower;
1282        delete [] saveUpper;
1283#else
1284        // Allow for crossover
1285        //if (doIdiot>0)
1286          info.setStrategy(512|info.getStrategy());
1287        // Allow for scaling
1288        info.setStrategy(32|info.getStrategy());
1289        info.crash(nPasses,model2->messageHandler(),model2->messagesPointer());
1290#endif
1291        time2 = CoinCpuTime();
1292        timeIdiot = time2-timeX;
1293        handler_->message(CLP_INTERVAL_TIMING,messages_)
1294          <<"Idiot Crash"<<timeIdiot<<time2-time1
1295          <<CoinMessageEol;
1296        timeX=time2;
1297      }
1298    }
1299#endif
1300    // ?
1301    if (doCrash) {
1302      switch(doCrash) {
1303        // standard
1304      case 1:
1305        model2->crash(1000,1);
1306        break;
1307        // As in paper by Solow and Halim (approx)
1308      case 2:
1309        model2->crash(model2->dualBound(),0);
1310        break;
1311        // My take on it
1312      case 3:
1313        model2->crash(model2->dualBound(),-1);
1314        break;
1315        // Just put free in basis
1316      case 4:
1317        model2->crash(0.0,3);
1318        break;
1319      }
1320    }
1321#ifndef SLIM_CLP
1322    if (doSlp>0&&objective_->type()==2) {
1323      model2->nonlinearSLP(doSlp,1.0e-5);
1324    }
1325#endif
1326    if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1327      // See if original wanted vector
1328      ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1329      ClpMatrixBase * matrix = model2->clpMatrix();
1330      if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
1331        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1332        clpMatrix->makeSpecialColumnCopy();
1333        //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
1334        model2->primal(primalStartup);
1335        clpMatrix->releaseSpecialColumnCopy();
1336      } else {
1337        model2->primal(primalStartup);
1338      }
1339    } else {
1340      model2->primal(primalStartup);
1341    }
1342    time2 = CoinCpuTime();
1343    timeCore = time2-timeX;
1344    handler_->message(CLP_INTERVAL_TIMING,messages_)
1345      <<"Primal"<<timeCore<<time2-time1
1346      <<CoinMessageEol;
1347    timeX=time2;
1348  } else if (method==ClpSolve::usePrimalorSprint) {
1349    // Sprint
1350    /*
1351      This driver implements what I called Sprint when I introduced the idea
1352      many years ago.  Cplex calls it "sifting" which I think is just as silly.
1353      When I thought of this trivial idea
1354      it reminded me of an LP code of the 60's called sprint which after
1355      every factorization took a subset of the matrix into memory (all
1356      64K words!) and then iterated very fast on that subset.  On the
1357      problems of those days it did not work very well, but it worked very
1358      well on aircrew scheduling problems where there were very large numbers
1359      of columns all with the same flavor.
1360    */
1361   
1362    /* The idea works best if you can get feasible easily.  To make it
1363       more general we can add in costed slacks */
1364   
1365    int originalNumberColumns = model2->numberColumns();
1366    int numberRows = model2->numberRows();
1367    ClpSimplex * originalModel2 = model2;
1368   
1369    // We will need arrays to choose variables.  These are too big but ..
1370    double * weight = new double [numberRows+originalNumberColumns];
1371    int * sort = new int [numberRows+originalNumberColumns];
1372    int numberSort=0;
1373    // We are going to add slacks to get feasible.
1374    // initial list will just be artificials
1375    int iColumn;
1376    const double * columnLower = model2->columnLower();
1377    const double * columnUpper = model2->columnUpper();
1378    double * columnSolution = model2->primalColumnSolution();
1379
1380    // See if we have costed slacks
1381    int * negSlack = new int[numberRows];
1382    int * posSlack = new int[numberRows];
1383    int iRow;
1384    for (iRow=0;iRow<numberRows;iRow++) {
1385      negSlack[iRow]=-1;
1386      posSlack[iRow]=-1;
1387    }
1388    const double * element = model2->matrix()->getElements();
1389    const int * row = model2->matrix()->getIndices();
1390    const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
1391    const int * columnLength = model2->matrix()->getVectorLengths();
1392    //bool allSlack = (numberRowsBasic==numberRows);
1393    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1394      if (!columnSolution[iColumn]||fabs(columnSolution[iColumn])>1.0e20) {
1395        double value =0.0;
1396        if (columnLower[iColumn]>0.0)
1397          value = columnLower[iColumn];
1398        else if (columnUpper[iColumn]<0.0)
1399          value = columnUpper[iColumn];
1400        columnSolution[iColumn]=value;
1401      }
1402      if (columnLength[iColumn]==1) {
1403        int jRow=row[columnStart[iColumn]];
1404        if (!columnLower[iColumn]) {
1405          if (element[columnStart[iColumn]]>0.0&&posSlack[jRow]<0)
1406            posSlack[jRow]=iColumn;
1407          else if (element[columnStart[iColumn]]<0.0&&negSlack[jRow]<0)
1408            negSlack[jRow]=iColumn;
1409        } else if (!columnUpper[iColumn]) {
1410          if (element[columnStart[iColumn]]<0.0&&posSlack[jRow]<0)
1411            posSlack[jRow]=iColumn;
1412          else if (element[columnStart[iColumn]]>0.0&&negSlack[jRow]<0)
1413            negSlack[jRow]=iColumn;
1414        }
1415      }
1416    }
1417    // now see what that does to row solution
1418    double * rowSolution = model2->primalRowSolution();
1419    CoinZeroN (rowSolution,numberRows);
1420    model2->clpMatrix()->times(1.0,columnSolution,rowSolution);
1421    // See if we can adjust using costed slacks
1422    double penalty=CoinMax(1.0e5,CoinMin(infeasibilityCost_*0.01,1.0e10))*optimizationDirection_;
1423    const double * lower = model2->rowLower();
1424    const double * upper = model2->rowUpper();
1425    for (iRow=0;iRow<numberRows;iRow++) {
1426      if (lower[iRow]>rowSolution[iRow]+1.0e-8) {
1427        int jColumn = posSlack[iRow];
1428        if (jColumn>=0) {
1429          if (columnSolution[jColumn])
1430            continue;
1431          double difference = lower[iRow]-rowSolution[iRow];
1432          double elementValue = element[columnStart[jColumn]];
1433          if (elementValue>0.0) {
1434            double movement = CoinMin(difference/elementValue,columnUpper[jColumn]);
1435            columnSolution[jColumn] = movement;
1436            rowSolution[iRow] += movement*elementValue;
1437          } else {
1438            double movement = CoinMax(difference/elementValue,columnLower[jColumn]);
1439            columnSolution[jColumn] = movement;
1440            rowSolution[iRow] += movement*elementValue;
1441          }
1442        }
1443      } else if (upper[iRow]<rowSolution[iRow]-1.0e-8) {
1444        int jColumn = negSlack[iRow];
1445        if (jColumn>=0) {
1446          if (columnSolution[jColumn])
1447            continue;
1448          double difference = upper[iRow]-rowSolution[iRow];
1449          double elementValue = element[columnStart[jColumn]];
1450          if (elementValue<0.0) {
1451            double movement = CoinMin(difference/elementValue,columnUpper[jColumn]);
1452            columnSolution[jColumn] = movement;
1453            rowSolution[iRow] += movement*elementValue;
1454          } else {
1455            double movement = CoinMax(difference/elementValue,columnLower[jColumn]);
1456            columnSolution[jColumn] = movement;
1457            rowSolution[iRow] += movement*elementValue;
1458          }
1459        }
1460      }
1461    }
1462    delete [] negSlack;
1463    delete [] posSlack;
1464    int nRow=numberRows;
1465    bool network=false;
1466    if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
1467      network=true;
1468      nRow *= 2;
1469    }
1470    int * addStarts = new int [nRow+1];
1471    int * addRow = new int[nRow];
1472    double * addElement = new double[nRow];
1473    addStarts[0]=0;
1474    int numberArtificials=0;
1475    int numberAdd=0;
1476    double * addCost = new double [numberRows];
1477    for (iRow=0;iRow<numberRows;iRow++) {
1478      if (lower[iRow]>rowSolution[iRow]+1.0e-8) {
1479        addRow[numberAdd]=iRow;
1480        addElement[numberAdd++]=1.0;
1481        if (network) {
1482          addRow[numberAdd]=numberRows;
1483          addElement[numberAdd++]=-1.0;
1484        }
1485        addCost[numberArtificials]=penalty;
1486        numberArtificials++;
1487        addStarts[numberArtificials]=numberAdd;
1488      } else if (upper[iRow]<rowSolution[iRow]-1.0e-8) {
1489        addRow[numberAdd]=iRow;
1490        addElement[numberAdd++]=-1.0;
1491        if (network) {
1492          addRow[numberAdd]=numberRows;
1493          addElement[numberAdd++]=1.0;
1494        }
1495        addCost[numberArtificials]=penalty;
1496        numberArtificials++;
1497        addStarts[numberArtificials]=numberAdd;
1498      }
1499    }
1500    if (numberArtificials) {
1501      // need copy so as not to disturb original
1502      model2 = new ClpSimplex(*model2);
1503      if (network) {
1504        // network - add a null row
1505        model2->addRow(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX);
1506        numberRows++;
1507      }
1508      model2->addColumns(numberArtificials,NULL,NULL,addCost,
1509                         addStarts,addRow,addElement);
1510    }
1511    delete [] addStarts;
1512    delete [] addRow;
1513    delete [] addElement;
1514    delete [] addCost;
1515    // look at rhs to see if to perturb
1516    double largest=0.0;
1517    double smallest = 1.0e30;
1518    for (iRow=0;iRow<numberRows;iRow++) {
1519      double value;
1520      value = fabs(model2->rowLower_[iRow]);
1521      if (value&&value<1.0e30) {
1522        largest = CoinMax(largest,value);
1523        smallest=CoinMin(smallest,value);
1524      }
1525      value = fabs(model2->rowUpper_[iRow]);
1526      if (value&&value<1.0e30) {
1527        largest = CoinMax(largest,value);
1528        smallest=CoinMin(smallest,value);
1529      }
1530    }
1531    double * saveLower = NULL;
1532    double * saveUpper = NULL;
1533    if (largest<2.01*smallest) {
1534      // perturb - so switch off standard
1535      model2->setPerturbation(100);
1536      saveLower = new double[numberRows];
1537      CoinMemcpyN(model2->rowLower_,numberRows,saveLower);
1538      saveUpper = new double[numberRows];
1539      CoinMemcpyN(model2->rowUpper_,numberRows,saveUpper);
1540      double * lower = model2->rowLower();
1541      double * upper = model2->rowUpper();
1542      for (iRow=0;iRow<numberRows;iRow++) {
1543        double lowerValue=lower[iRow], upperValue=upper[iRow];
1544        double value = randomNumberGenerator_.randomDouble();
1545        if (upperValue>lowerValue+primalTolerance_) {
1546          if (lowerValue>-1.0e20&&lowerValue)
1547            lowerValue -= value * 1.0e-4*fabs(lowerValue); 
1548          if (upperValue<1.0e20&&upperValue)
1549            upperValue += value * 1.0e-4*fabs(upperValue); 
1550        } else if (upperValue>0.0) {
1551          upperValue -= value * 1.0e-4*fabs(lowerValue); 
1552          lowerValue -= value * 1.0e-4*fabs(lowerValue); 
1553        } else if (upperValue<0.0) {
1554          upperValue += value * 1.0e-4*fabs(lowerValue); 
1555          lowerValue += value * 1.0e-4*fabs(lowerValue); 
1556        } else {
1557        }
1558        lower[iRow]=lowerValue;
1559        upper[iRow]=upperValue;
1560      }
1561    }
1562    int i;
1563    // Just do this number of passes in Sprint
1564    if (doSprint>0)
1565      maxSprintPass=options.getExtraInfo(1);
1566    // but if big use to get ratio
1567    double ratio=3;
1568    if (maxSprintPass>1000) {
1569      ratio = static_cast<double> (maxSprintPass)*0.0001;
1570      ratio = CoinMax(ratio,1.1);
1571      maxSprintPass= maxSprintPass %1000;
1572#ifdef COIN_DEVELOP
1573      printf("%d passes wanted with ratio of %g\n",maxSprintPass,ratio);
1574#endif
1575    }
1576    // Just take this number of columns in small problem
1577    int smallNumberColumns = static_cast<int> (CoinMin(ratio*numberRows,static_cast<double> (numberColumns)));
1578    smallNumberColumns = CoinMax(smallNumberColumns,3000);
1579    smallNumberColumns = CoinMin(smallNumberColumns,numberColumns);
1580    //int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
1581    //smallNumberColumns = CoinMax(smallNumberColumns,3000);
1582    //smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
1583    // redo as may have changed
1584    columnLower = model2->columnLower();
1585    columnUpper = model2->columnUpper();
1586    columnSolution = model2->primalColumnSolution();
1587    // Set up initial list
1588    numberSort=0;
1589    if (numberArtificials) {
1590      numberSort=numberArtificials;
1591      for (i=0;i<numberSort;i++)
1592        sort[i] = i+originalNumberColumns;
1593    } 
1594    // maybe a solution there already
1595    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1596      if (model2->getColumnStatus(iColumn)==basic)
1597        sort[numberSort++]=iColumn;
1598    }
1599    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
1600      if (model2->getColumnStatus(iColumn)!=basic) {
1601        if (columnSolution[iColumn]>columnLower[iColumn]&&
1602            columnSolution[iColumn]<columnUpper[iColumn]&&
1603            columnSolution[iColumn])
1604          sort[numberSort++]=iColumn;
1605      }
1606    }
1607    numberSort = CoinMin(numberSort,smallNumberColumns);
1608   
1609    int numberColumns = model2->numberColumns();
1610    double * fullSolution = model2->primalColumnSolution();
1611   
1612   
1613    int iPass;
1614    double lastObjective=1.0e31;
1615    // It will be safe to allow dense
1616    model2->setInitialDenseFactorization(true);
1617   
1618    // We will be using all rows
1619    int * whichRows = new int [numberRows];
1620    for (iRow=0;iRow<numberRows;iRow++)
1621      whichRows[iRow]=iRow;
1622    double originalOffset;
1623    model2->getDblParam(ClpObjOffset,originalOffset);
1624    int totalIterations=0;
1625    double lastSumArtificials=COIN_DBL_MAX;
1626    int originalMaxSprintPass=maxSprintPass;
1627    maxSprintPass=20; // so we do that many if infeasible
1628    for (iPass=0;iPass<maxSprintPass;iPass++) {
1629      //printf("Bug until submodel new version\n");
1630      //CoinSort_2(sort,sort+numberSort,weight);
1631      // Create small problem
1632      ClpSimplex small(model2,numberRows,whichRows,numberSort,sort);
1633      small.setPerturbation(model2->perturbation());
1634      small.setInfeasibilityCost(model2->infeasibilityCost());
1635      if (model2->factorizationFrequency()==200) {
1636        // User did not touch preset
1637        small.defaultFactorizationFrequency();
1638      }
1639      // now see what variables left out do to row solution
1640      double * rowSolution = model2->primalRowSolution();
1641      double * sumFixed = new double[numberRows];
1642      CoinZeroN (sumFixed,numberRows);
1643      int iRow,iColumn;
1644      // zero out ones in small problem
1645      for (iColumn=0;iColumn<numberSort;iColumn++) {
1646        int kColumn = sort[iColumn];
1647        fullSolution[kColumn]=0.0;
1648      }
1649      // Get objective offset
1650      double offset=0.0;
1651      const double * objective = model2->objective();
1652      for (iColumn=0;iColumn<numberColumns;iColumn++) 
1653        offset += fullSolution[iColumn]*objective[iColumn];
1654      small.setDblParam(ClpObjOffset,originalOffset-offset);
1655      model2->clpMatrix()->times(1.0,fullSolution,sumFixed);
1656     
1657      double * lower = small.rowLower();
1658      double * upper = small.rowUpper();
1659      for (iRow=0;iRow<numberRows;iRow++) {
1660        if (lower[iRow]>-1.0e50) 
1661          lower[iRow] -= sumFixed[iRow];
1662        if (upper[iRow]<1.0e50)
1663          upper[iRow] -= sumFixed[iRow];
1664        rowSolution[iRow] -= sumFixed[iRow];
1665      }
1666      delete [] sumFixed;
1667      // Solve
1668      if (interrupt)
1669        currentModel = &small;
1670      small.defaultFactorizationFrequency();
1671      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
1672        // See if original wanted vector
1673        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
1674        ClpMatrixBase * matrix = small.clpMatrix();
1675        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
1676          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
1677          clpMatrix->makeSpecialColumnCopy();
1678          small.primal(1);
1679          clpMatrix->releaseSpecialColumnCopy();
1680        } else {
1681#if 1
1682          small.primal(1);
1683#else
1684          int numberColumns = small.numberColumns();
1685          int numberRows = small.numberRows();
1686          // Use dual region
1687          double * rhs = small.dualRowSolution();
1688          int * whichRow = new int[3*numberRows];
1689          int * whichColumn = new int[2*numberColumns];
1690          int nBound;
1691          ClpSimplex * small2 = ((ClpSimplexOther *) (&small))->crunch(rhs,whichRow,whichColumn,
1692                                                                        nBound,false,false);
1693          if (small2) {
1694            small2->primal(1);
1695            if (small2->problemStatus()==0) {
1696              small.setProblemStatus(0);
1697              ((ClpSimplexOther *) (&small))->afterCrunch(*small2,whichRow,whichColumn,nBound);
1698            } else {
1699              small2->primal(1);
1700              if (small2->problemStatus())
1701                small.primal(1);
1702            }
1703            delete small2;
1704          } else {
1705            small.primal(1);
1706          }
1707          delete [] whichRow;
1708          delete [] whichColumn;
1709#endif
1710        }
1711      } else {
1712        small.primal(1);
1713      }
1714      totalIterations += small.numberIterations();
1715      // move solution back
1716      const double * solution = small.primalColumnSolution();
1717      for (iColumn=0;iColumn<numberSort;iColumn++) {
1718        int kColumn = sort[iColumn];
1719        model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
1720        fullSolution[kColumn]=solution[iColumn];
1721      }
1722      for (iRow=0;iRow<numberRows;iRow++) 
1723        model2->setRowStatus(iRow,small.getRowStatus(iRow));
1724      CoinMemcpyN(small.primalRowSolution(),
1725             numberRows,model2->primalRowSolution());
1726      double sumArtificials = 0.0;
1727      for (i=0;i<numberArtificials;i++)
1728        sumArtificials += fullSolution[i + originalNumberColumns];
1729      if (sumArtificials&&iPass>5&&sumArtificials>=lastSumArtificials) {
1730        // increase costs
1731        double * cost = model2->objective()+originalNumberColumns;
1732        double newCost = CoinMin(1.0e10,cost[0]*1.5);
1733        for (i=0;i<numberArtificials;i++)
1734          cost[i]=newCost;
1735      }
1736      lastSumArtificials = sumArtificials;
1737      // get reduced cost for large problem
1738      double * djs = model2->dualColumnSolution();
1739      CoinMemcpyN(model2->objective(),numberColumns,djs);
1740      model2->clpMatrix()->transposeTimes(-1.0,small.dualRowSolution(),djs);
1741      int numberNegative=0;
1742      double sumNegative = 0.0;
1743      // now massage weight so all basic in plus good djs
1744      // first count and do basic
1745      numberSort=0;
1746      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1747        double dj = djs[iColumn]*optimizationDirection_;
1748        double value = fullSolution[iColumn];
1749        if (model2->getColumnStatus(iColumn)==ClpSimplex::basic) {
1750          sort[numberSort++] = iColumn;
1751        } else if (dj<-dualTolerance_&&value<columnUpper[iColumn]) {
1752          numberNegative++;
1753          sumNegative -= dj;
1754        } else if (dj>dualTolerance_&&value>columnLower[iColumn]) {
1755          numberNegative++;
1756          sumNegative += dj;
1757        }
1758      }
1759      handler_->message(CLP_SPRINT,messages_)
1760        <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative
1761        <<numberNegative
1762        <<CoinMessageEol;
1763      if (sumArtificials<1.0e-8&&originalMaxSprintPass>=0) {
1764        maxSprintPass = iPass+originalMaxSprintPass;
1765        originalMaxSprintPass=-1;
1766      }
1767      if (iPass>20)
1768        sumArtificials=0.0;
1769      if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5&&sumArtificials<1.0e-8)||
1770          (!small.numberIterations()&&iPass)||
1771          iPass==maxSprintPass-1||small.status()==3) {
1772
1773        break; // finished
1774      } else {
1775        lastObjective = small.objectiveValue()*optimizationDirection_;
1776        double tolerance;
1777        double averageNegDj = sumNegative/static_cast<double> (numberNegative+1);
1778        if (numberNegative+numberSort>smallNumberColumns)
1779          tolerance = -dualTolerance_;
1780        else 
1781          tolerance = 10.0*averageNegDj;
1782        int saveN = numberSort;
1783        for (iColumn=0;iColumn<numberColumns;iColumn++) {
1784          double dj = djs[iColumn]*optimizationDirection_;
1785          double value = fullSolution[iColumn];
1786          if (model2->getColumnStatus(iColumn)!=ClpSimplex::basic) {
1787            if (dj<-dualTolerance_&&value<columnUpper[iColumn])
1788              dj = dj;
1789            else if (dj>dualTolerance_&&value>columnLower[iColumn])
1790              dj = -dj;
1791            else if (columnUpper[iColumn]>columnLower[iColumn])
1792              dj = fabs(dj);
1793            else
1794              dj = 1.0e50;
1795            if (dj<tolerance) {
1796              weight[numberSort] = dj;
1797              sort[numberSort++] = iColumn;
1798            }
1799          }
1800        }
1801        // sort
1802        CoinSort_2(weight+saveN,weight+numberSort,sort+saveN);
1803        numberSort = CoinMin(smallNumberColumns,numberSort);
1804      }
1805    }
1806    if (interrupt) 
1807      currentModel = model2;
1808    for (i=0;i<numberArtificials;i++)
1809      sort[i] = i + originalNumberColumns;
1810    model2->deleteColumns(numberArtificials,sort);
1811    if (network) {
1812      int iRow=numberRows-1;
1813      model2->deleteRows(1,&iRow);
1814    }
1815    delete [] weight;
1816    delete [] sort;
1817    delete [] whichRows;
1818    if (saveLower) {
1819      // unperturb and clean
1820      for (iRow=0;iRow<numberRows;iRow++) {
1821        double diffLower = saveLower[iRow]-model2->rowLower_[iRow];
1822        double diffUpper = saveUpper[iRow]-model2->rowUpper_[iRow];
1823        model2->rowLower_[iRow]=saveLower[iRow];
1824        model2->rowUpper_[iRow]=saveUpper[iRow];
1825        if (diffLower) 
1826          assert (!diffUpper||fabs(diffLower-diffUpper)<1.0e-5);
1827        else
1828          diffLower = diffUpper;
1829        model2->rowActivity_[iRow] += diffLower;
1830      }
1831      delete [] saveLower;
1832      delete [] saveUpper;
1833    }
1834    model2->primal(1);
1835    model2->setPerturbation(savePerturbation);
1836    if (model2!=originalModel2) {
1837      originalModel2->moveInfo(*model2);
1838      delete model2;
1839      model2 = originalModel2;
1840    }
1841    time2 = CoinCpuTime();
1842    timeCore = time2-timeX;
1843    handler_->message(CLP_INTERVAL_TIMING,messages_)
1844      <<"Sprint"<<timeCore<<time2-time1
1845      <<CoinMessageEol;
1846    timeX=time2;
1847    model2->setNumberIterations(model2->numberIterations()+totalIterations);
1848  } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
1849#ifndef SLIM_CLP
1850    //printf("***** experimental pretty crude barrier\n");
1851    //#define SAVEIT 2
1852#ifndef SAVEIT
1853#define BORROW
1854#endif
1855#ifdef BORROW
1856    ClpInterior barrier;
1857    barrier.borrowModel(*model2);
1858#else
1859    ClpInterior barrier(*model2);
1860#endif
1861    if (interrupt) 
1862      currentModel2 = &barrier;
1863    int barrierOptions = options.getSpecialOption(4);
1864    int aggressiveGamma=0;
1865    bool presolveInCrossover=false;
1866    bool scale=false;
1867    bool doKKT=false;
1868    bool forceFixing=false;
1869    int speed=0;
1870    if (barrierOptions&16) {
1871      barrierOptions &= ~16;
1872      doKKT=true;
1873    }
1874    if (barrierOptions&(32+64+128)) {
1875      aggressiveGamma=(barrierOptions&(32+64+128))>>5;
1876      barrierOptions &= ~(32+64+128);
1877    }
1878    if (barrierOptions&256) {
1879      barrierOptions &= ~256;
1880      presolveInCrossover=true;
1881    }
1882    if (barrierOptions&512) {
1883      barrierOptions &= ~512;
1884      forceFixing=true;
1885    }
1886    if (barrierOptions&1024) {
1887      barrierOptions &= ~1024;
1888      barrier.setProjectionTolerance(1.0e-9);
1889    }
1890    if (barrierOptions&(2048|4096)) {
1891      speed = (barrierOptions&(2048|4096))>>11;
1892      barrierOptions &= ~(2048|4096);
1893    }
1894    if (barrierOptions&8) {
1895      barrierOptions &= ~8;
1896      scale=true;
1897    }
1898    // If quadratic force KKT
1899    if (quadraticObj) {
1900      doKKT=true;
1901    }
1902    switch (barrierOptions) {
1903    case 0:
1904    default:
1905      if (!doKKT) {
1906        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1907        cholesky->setIntegerParameter(0,speed);
1908        barrier.setCholesky(cholesky);
1909      } else {
1910        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1911        cholesky->setKKT(true);
1912        barrier.setCholesky(cholesky);
1913      }
1914      break;
1915    case 1:
1916      if (!doKKT) {
1917        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1918        barrier.setCholesky(cholesky);
1919      } else {
1920        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1921        cholesky->setKKT(true);
1922        barrier.setCholesky(cholesky);
1923      }
1924      break;
1925#ifdef WSSMP_BARRIER
1926    case 2:
1927      {
1928        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
1929        barrier.setCholesky(cholesky);
1930        assert (!doKKT);
1931      }
1932      break;
1933    case 3:
1934      if (!doKKT) {
1935        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
1936        barrier.setCholesky(cholesky);
1937      } else {
1938        ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
1939        barrier.setCholesky(cholesky);
1940      }
1941      break;
1942#endif
1943#ifdef UFL_BARRIER
1944    case 4:
1945      if (!doKKT) {
1946        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1947        barrier.setCholesky(cholesky);
1948      } else {
1949        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1950        cholesky->setKKT(true);
1951        barrier.setCholesky(cholesky);
1952      }
1953      break;
1954#endif
1955#ifdef TAUCS_BARRIER
1956    case 5:
1957      {
1958        ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
1959        barrier.setCholesky(cholesky);
1960        assert (!doKKT);
1961      }
1962      break;
1963#endif
1964#ifdef MUMPS_BARRIER
1965    case 6:
1966      {
1967        ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
1968        barrier.setCholesky(cholesky);
1969        assert (!doKKT);
1970      }
1971      break;
1972#endif
1973    }
1974    int numberRows = model2->numberRows();
1975    int numberColumns = model2->numberColumns();
1976    int saveMaxIts = model2->maximumIterations();
1977    if (saveMaxIts<1000) {
1978      barrier.setMaximumBarrierIterations(saveMaxIts);
1979      model2->setMaximumIterations(10000000);
1980    }
1981#ifndef SAVEIT
1982    //barrier.setDiagonalPerturbation(1.0e-25);
1983    if (aggressiveGamma) {
1984      switch (aggressiveGamma) {
1985      case 1:
1986        barrier.setGamma(1.0e-5);
1987        barrier.setDelta(1.0e-5);
1988        break;
1989      case 2:
1990        barrier.setGamma(1.0e-7);
1991        break;
1992      case 3:
1993        barrier.setDelta(1.0e-5);
1994        break;
1995      case 4:
1996        barrier.setGamma(1.0e-3);
1997        barrier.setDelta(1.0e-3);
1998        break;
1999      case 5:
2000        barrier.setGamma(1.0e-3);
2001        break;
2002      case 6:
2003        barrier.setDelta(1.0e-3);
2004        break;
2005      }
2006    }
2007    if (scale)
2008      barrier.scaling(1);
2009    else
2010      barrier.scaling(0);
2011    barrier.primalDual();
2012#elif SAVEIT==1
2013    barrier.primalDual();
2014#else
2015    model2->restoreModel("xx.save");
2016    // move solutions
2017    CoinMemcpyN(model2->primalRowSolution(),
2018                numberRows,barrier.primalRowSolution());
2019    CoinMemcpyN(model2->dualRowSolution(),
2020                numberRows,barrier.dualRowSolution());
2021    CoinMemcpyN(model2->primalColumnSolution(),
2022                numberColumns,barrier.primalColumnSolution());
2023    CoinMemcpyN(model2->dualColumnSolution(),
2024                numberColumns,barrier.dualColumnSolution());
2025#endif
2026    time2 = CoinCpuTime();
2027    timeCore = time2-timeX;
2028    handler_->message(CLP_INTERVAL_TIMING,messages_)
2029      <<"Barrier"<<timeCore<<time2-time1
2030      <<CoinMessageEol;
2031    timeX=time2;
2032    int maxIts = barrier.maximumBarrierIterations();
2033    int barrierStatus=barrier.status();
2034    double gap = barrier.complementarityGap();
2035    // get which variables are fixed
2036    double * saveLower=NULL;
2037    double * saveUpper=NULL;
2038    ClpPresolve pinfo2;
2039    ClpSimplex * saveModel2=NULL;
2040    bool extraPresolve=false;
2041    int numberFixed = barrier.numberFixed();
2042    if (numberFixed) {
2043      int numberRows = barrier.numberRows();
2044      int numberColumns = barrier.numberColumns();
2045      int numberTotal = numberRows+numberColumns;
2046      saveLower = new double [numberTotal];
2047      saveUpper = new double [numberTotal];
2048      CoinMemcpyN(barrier.columnLower(),numberColumns,saveLower);
2049      CoinMemcpyN(barrier.rowLower(),numberRows,saveLower+numberColumns);
2050      CoinMemcpyN(barrier.columnUpper(),numberColumns,saveUpper);
2051      CoinMemcpyN(barrier.rowUpper(),numberRows,saveUpper+numberColumns);
2052    }
2053    if (((numberFixed*20>barrier.numberRows()&&numberFixed>5000)||forceFixing)&&
2054        presolveInCrossover) {
2055      // may as well do presolve
2056      if (!forceFixing) {
2057        barrier.fixFixed();
2058      } else {
2059        // Fix
2060        int n=barrier.numberColumns();
2061        double * lower = barrier.columnLower();
2062        double * upper = barrier.columnUpper();
2063        double * solution = barrier.primalColumnSolution();
2064        int nFix=0;
2065        for (int i=0;i<n;i++) {
2066          if (barrier.fixedOrFree(i)&&lower[i]<upper[i]) {
2067            double value = solution[i];
2068            if (value<lower[i]+1.0e-6&&value-lower[i]<upper[i]-value) {
2069              solution[i]=lower[i];
2070              upper[i]=lower[i];
2071              nFix++;
2072            } else if (value>upper[i]-1.0e-6&&value-lower[i]>upper[i]-value) {
2073              solution[i]=upper[i];
2074              lower[i]=upper[i];
2075              nFix++;
2076            }
2077          }
2078        }
2079#ifdef CLP_INVESTIGATE
2080        printf("%d columns fixed\n",nFix);
2081#endif
2082        int nr=barrier.numberRows();
2083        lower = barrier.rowLower();
2084        upper = barrier.rowUpper();
2085        solution = barrier.primalRowSolution();
2086        nFix=0;
2087        for (int i=0;i<nr;i++) {
2088          if (barrier.fixedOrFree(i+n)&&lower[i]<upper[i]) {
2089            double value = solution[i];
2090            if (value<lower[i]+1.0e-6&&value-lower[i]<upper[i]-value) {
2091              solution[i]=lower[i];
2092              upper[i]=lower[i];
2093              nFix++;
2094            } else if (value>upper[i]-1.0e-6&&value-lower[i]>upper[i]-value) {
2095              solution[i]=upper[i];
2096              lower[i]=upper[i];
2097              nFix++;
2098            }
2099          }
2100        }
2101#ifdef CLP_INVESTIGATE
2102        printf("%d row slacks fixed\n",nFix);
2103#endif
2104      }
2105      saveModel2=model2;
2106      extraPresolve=true;
2107    } else if (numberFixed) {
2108      // Set fixed to bounds (may have restored earlier solution)
2109      if (!forceFixing) {
2110        barrier.fixFixed(false);
2111      } else {
2112        // Fix
2113        int n=barrier.numberColumns();
2114        double * lower = barrier.columnLower();
2115        double * upper = barrier.columnUpper();
2116        double * solution = barrier.primalColumnSolution();
2117        int nFix=0;
2118        for (int i=0;i<n;i++) {
2119          if (barrier.fixedOrFree(i)&&lower[i]<upper[i]) {
2120            double value = solution[i];
2121            if (value<lower[i]+1.0e-8&&value-lower[i]<upper[i]-value) {
2122              solution[i]=lower[i];
2123              upper[i]=lower[i];
2124              nFix++;
2125            } else if (value>upper[i]-1.0e-8&&value-lower[i]>upper[i]-value) {
2126              solution[i]=upper[i];
2127              lower[i]=upper[i];
2128              nFix++;
2129            } else {
2130              //printf("fixcol %d %g <= %g <= %g\n",
2131              //     i,lower[i],solution[i],upper[i]);
2132            }
2133          }
2134        }
2135#ifdef CLP_INVESTIGATE
2136        printf("%d columns fixed\n",nFix);
2137#endif
2138        int nr=barrier.numberRows();
2139        lower = barrier.rowLower();
2140        upper = barrier.rowUpper();
2141        solution = barrier.primalRowSolution();
2142        nFix=0;
2143        for (int i=0;i<nr;i++) {
2144          if (barrier.fixedOrFree(i+n)&&lower[i]<upper[i]) {
2145            double value = solution[i];
2146            if (value<lower[i]+1.0e-5&&value-lower[i]<upper[i]-value) {
2147              solution[i]=lower[i];
2148              upper[i]=lower[i];
2149              nFix++;
2150            } else if (value>upper[i]-1.0e-5&&value-lower[i]>upper[i]-value) {
2151              solution[i]=upper[i];
2152              lower[i]=upper[i];
2153              nFix++;
2154            } else {
2155              //printf("fixrow %d %g <= %g <= %g\n",
2156              //     i,lower[i],solution[i],upper[i]);
2157            }
2158          }
2159        }
2160#ifdef CLP_INVESTIGATE
2161        printf("%d row slacks fixed\n",nFix);
2162#endif
2163      }
2164    }
2165#ifdef BORROW   
2166    barrier.returnModel(*model2);
2167    double * rowPrimal = new double [numberRows];
2168    double * columnPrimal = new double [numberColumns];
2169    double * rowDual = new double [numberRows];
2170    double * columnDual = new double [numberColumns];
2171    // move solutions other way
2172    CoinMemcpyN(model2->primalRowSolution(),
2173                numberRows,rowPrimal);
2174    CoinMemcpyN(model2->dualRowSolution(),
2175                numberRows,rowDual);
2176    CoinMemcpyN(model2->primalColumnSolution(),
2177                numberColumns,columnPrimal);
2178    CoinMemcpyN(model2->dualColumnSolution(),
2179                  numberColumns,columnDual);
2180#else
2181    double * rowPrimal = barrier.primalRowSolution();
2182    double * columnPrimal = barrier.primalColumnSolution();
2183    double * rowDual = barrier.dualRowSolution();
2184    double * columnDual = barrier.dualColumnSolution();
2185    // move solutions
2186    CoinMemcpyN(rowPrimal,
2187                numberRows,model2->primalRowSolution());
2188    CoinMemcpyN(rowDual,
2189                numberRows,model2->dualRowSolution());
2190    CoinMemcpyN(columnPrimal,
2191                numberColumns,model2->primalColumnSolution());
2192    CoinMemcpyN(columnDual,
2193                  numberColumns,model2->dualColumnSolution());
2194#endif
2195    if (saveModel2) {
2196      // do presolve
2197      model2 = pinfo2.presolvedModel(*model2,dblParam_[ClpPresolveTolerance],
2198                                    false,5,true);
2199      if (!model2) {
2200        model2=saveModel2;
2201        saveModel2=NULL;
2202        int numberRows = model2->numberRows();
2203        int numberColumns = model2->numberColumns();
2204        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2205        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2206        delete [] saveLower;
2207        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2208        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2209        delete [] saveUpper;
2210        saveLower=NULL;
2211        saveUpper=NULL;
2212      }
2213    }
2214    if (method==ClpSolve::useBarrier) {
2215      if (maxIts&&barrierStatus<4&&!quadraticObj) {
2216        //printf("***** crossover - needs more thought on difficult models\n");
2217#if SAVEIT==1
2218        model2->ClpSimplex::saveModel("xx.save");
2219#endif
2220        // make sure no status left
2221        model2->createStatus();
2222        // solve
2223        if (!forceFixing)
2224          model2->setPerturbation(100);
2225        if (model2->factorizationFrequency()==200) {
2226          // User did not touch preset
2227          model2->defaultFactorizationFrequency();
2228        }
2229#if 1
2230        // throw some into basis
2231        if(!forceFixing) {
2232          int numberRows = model2->numberRows();
2233          int numberColumns = model2->numberColumns();
2234          double * dsort = new double[numberColumns];
2235          int * sort = new int[numberColumns];
2236          int n=0;
2237          const double * columnLower = model2->columnLower();
2238          const double * columnUpper = model2->columnUpper();
2239          double * primalSolution = model2->primalColumnSolution();
2240          const double * dualSolution = model2->dualColumnSolution();
2241          double tolerance = 10.0*primalTolerance_;
2242          int i;
2243          for ( i=0;i<numberRows;i++) 
2244            model2->setRowStatus(i,superBasic);
2245          for ( i=0;i<numberColumns;i++) {
2246            double distance = CoinMin(columnUpper[i]-primalSolution[i],
2247                                      primalSolution[i]-columnLower[i]);
2248            if (distance>tolerance) {
2249              if (fabs(dualSolution[i])<1.0e-5)
2250                distance *= 100.0;
2251              dsort[n]=-distance;
2252              sort[n++]=i;
2253              model2->setStatus(i,superBasic);
2254            } else if (distance>primalTolerance_) {
2255              model2->setStatus(i,superBasic);
2256            } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
2257              model2->setStatus(i,atLowerBound);
2258              primalSolution[i]=columnLower[i];
2259            } else {
2260              model2->setStatus(i,atUpperBound);
2261              primalSolution[i]=columnUpper[i];
2262            }
2263          }
2264          CoinSort_2(dsort,dsort+n,sort);
2265          n = CoinMin(numberRows,n);
2266          for ( i=0;i<n;i++) {
2267            int iColumn = sort[i];
2268            model2->setStatus(iColumn,basic);
2269          }
2270          delete [] sort;
2271          delete [] dsort;
2272          // model2->allSlackBasis();
2273          if (gap<1.0e-3*static_cast<double> (numberRows+numberColumns)) {
2274            if (saveUpper) {
2275              int numberRows = model2->numberRows();
2276              int numberColumns = model2->numberColumns();
2277              CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2278              CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2279              CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2280              CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2281              delete [] saveLower;
2282              delete [] saveUpper;
2283              saveLower=NULL;
2284              saveUpper=NULL;
2285            }
2286            int numberRows = model2->numberRows();
2287            int numberColumns = model2->numberColumns();
2288            // just primal values pass
2289            double saveScale = model2->objectiveScale();
2290            model2->setObjectiveScale(1.0e-3);
2291            model2->primal(2);
2292            model2->setObjectiveScale(saveScale);
2293            // save primal solution and copy back dual
2294            CoinMemcpyN(model2->primalRowSolution(),
2295                        numberRows,rowPrimal);
2296            CoinMemcpyN(rowDual,
2297                        numberRows,model2->dualRowSolution());
2298            CoinMemcpyN(model2->primalColumnSolution(),
2299                        numberColumns,columnPrimal);
2300            CoinMemcpyN(columnDual,
2301                        numberColumns,model2->dualColumnSolution());
2302            //model2->primal(1);
2303            // clean up reduced costs and flag variables
2304            {
2305              double * dj = model2->dualColumnSolution();
2306              double * cost = model2->objective();
2307              double * saveCost = new double[numberColumns];
2308              CoinMemcpyN(cost,numberColumns,saveCost);
2309              double * saveLower = new double[numberColumns];
2310              double * lower = model2->columnLower();
2311              CoinMemcpyN(lower,numberColumns,saveLower);
2312              double * saveUpper = new double[numberColumns];
2313              double * upper = model2->columnUpper();
2314              CoinMemcpyN(upper,numberColumns,saveUpper);
2315              int i;
2316              double tolerance = 10.0*dualTolerance_;
2317              for ( i=0;i<numberColumns;i++) {
2318                if (model2->getStatus(i)==basic) {
2319                  dj[i]=0.0;
2320                } else if (model2->getStatus(i)==atLowerBound) {
2321                  if (optimizationDirection_*dj[i]<tolerance) {
2322                    if (optimizationDirection_*dj[i]<0.0) {
2323                      //if (dj[i]<-1.0e-3)
2324                      //printf("bad dj at lb %d %g\n",i,dj[i]);
2325                      cost[i] -= dj[i];
2326                      dj[i]=0.0;
2327                    }
2328                  } else {
2329                    upper[i]=lower[i];
2330                  }
2331                } else if (model2->getStatus(i)==atUpperBound) {
2332                  if (optimizationDirection_*dj[i]>tolerance) {
2333                    if (optimizationDirection_*dj[i]>0.0) {
2334                      //if (dj[i]>1.0e-3)
2335                      //printf("bad dj at ub %d %g\n",i,dj[i]);
2336                      cost[i] -= dj[i];
2337                      dj[i]=0.0;
2338                    }
2339                  } else {
2340                    lower[i]=upper[i];
2341                  }
2342                }
2343              }
2344              // just dual values pass
2345              //model2->setLogLevel(63);
2346              //model2->setFactorizationFrequency(1);
2347              model2->dual(2);
2348              CoinMemcpyN(saveCost,numberColumns,cost);
2349              delete [] saveCost;
2350              CoinMemcpyN(saveLower,numberColumns,lower);
2351              delete [] saveLower;
2352              CoinMemcpyN(saveUpper,numberColumns,upper);
2353              delete [] saveUpper;
2354            }
2355          }
2356          // and finish
2357          // move solutions
2358          CoinMemcpyN(rowPrimal,
2359                      numberRows,model2->primalRowSolution());
2360          CoinMemcpyN(columnPrimal,
2361                      numberColumns,model2->primalColumnSolution());
2362        }
2363        double saveScale = model2->objectiveScale();
2364        model2->setObjectiveScale(1.0e-3);
2365        model2->primal(2);
2366        model2->setObjectiveScale(saveScale);
2367        model2->primal(1);
2368#else
2369        // just primal
2370        model2->primal(1);
2371#endif
2372      } else if (barrierStatus==4) {
2373        // memory problems
2374        model2->setPerturbation(savePerturbation);
2375        model2->createStatus();
2376        model2->dual();
2377      } else if (maxIts&&quadraticObj) {
2378        // make sure no status left
2379        model2->createStatus();
2380        // solve
2381        model2->setPerturbation(100);
2382        model2->reducedGradient(1);
2383      }
2384    }
2385   
2386    //model2->setMaximumIterations(saveMaxIts);
2387#ifdef BORROW
2388    delete [] rowPrimal;
2389    delete [] columnPrimal;
2390    delete [] rowDual;
2391    delete [] columnDual;
2392#endif
2393    if (extraPresolve) {
2394      pinfo2.postsolve(true);
2395      delete model2;
2396      model2=saveModel2;
2397    }
2398    if (saveUpper) {
2399      if (!forceFixing) {
2400        int numberRows = model2->numberRows();
2401        int numberColumns = model2->numberColumns();
2402        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2403        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2404        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2405        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2406      }
2407      delete [] saveLower;
2408      delete [] saveUpper;
2409      saveLower=NULL;
2410      saveUpper=NULL;
2411      if (method!=ClpSolve::useBarrierNoCross) 
2412        model2->primal(1);
2413    }
2414    model2->setPerturbation(savePerturbation);
2415    time2 = CoinCpuTime();
2416    timeCore = time2-timeX;
2417    handler_->message(CLP_INTERVAL_TIMING,messages_)
2418      <<"Crossover"<<timeCore<<time2-time1
2419      <<CoinMessageEol;
2420    timeX=time2;
2421#else
2422    abort();
2423#endif
2424  } else {
2425    assert (method!=ClpSolve::automatic); // later
2426    time2=0.0;
2427  }
2428  if (saveMatrix) {
2429    if (model2==this) {
2430      // delete and replace
2431      delete model2->clpMatrix();
2432      model2->replaceMatrix(saveMatrix);
2433    } else {
2434      delete saveMatrix;
2435    }
2436  }
2437  numberIterations = model2->numberIterations();
2438  finalStatus=model2->status();
2439  int finalSecondaryStatus = model2->secondaryStatus();
2440  if (presolve==ClpSolve::presolveOn) {
2441    int saveLevel = logLevel();
2442    if ((specialOptions_&1024)==0)
2443      setLogLevel(CoinMin(1,saveLevel));
2444    else
2445      setLogLevel(CoinMin(0,saveLevel));
2446    pinfo.postsolve(true);
2447    factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
2448    time2 = CoinCpuTime();
2449    timePresolve += time2-timeX;
2450    handler_->message(CLP_INTERVAL_TIMING,messages_)
2451      <<"Postsolve"<<time2-timeX<<time2-time1
2452      <<CoinMessageEol;
2453    timeX=time2;
2454    if (!presolveToFile)
2455      delete model2;
2456    if (interrupt)
2457      currentModel = this;
2458    // checkSolution(); already done by postSolve
2459    setLogLevel(saveLevel);
2460    if (finalStatus!=3&&(finalStatus||status()==-1)) {
2461      int savePerturbation = perturbation();
2462      if (!finalStatus||(moreSpecialOptions_&2)==0) {
2463        if (finalStatus==2) {
2464          // unbounded - get feasible first
2465          double save = optimizationDirection_;
2466          optimizationDirection_=0.0;
2467          primal(1);
2468          optimizationDirection_=save;
2469          primal(1);
2470        } else if (finalStatus==1) {
2471          dual();
2472        } else {
2473          setPerturbation(100);
2474          primal(1);
2475        }
2476      } else {
2477        // just set status
2478        problemStatus_=finalStatus;
2479      }
2480      setPerturbation(savePerturbation);
2481      numberIterations += numberIterations_;
2482      numberIterations_ = numberIterations;
2483      finalStatus=status();
2484      time2 = CoinCpuTime();
2485      handler_->message(CLP_INTERVAL_TIMING,messages_)
2486      <<"Cleanup"<<time2-timeX<<time2-time1
2487      <<CoinMessageEol;
2488      timeX=time2;
2489    } else {
2490      secondaryStatus_=finalSecondaryStatus;
2491    }
2492  } else if (model2!=this) {
2493    // not presolved - but different model used (sprint probably)
2494    CoinMemcpyN(model2->primalRowSolution(),
2495                numberRows_,this->primalRowSolution());
2496    CoinMemcpyN(model2->dualRowSolution(),
2497                numberRows_,this->dualRowSolution());
2498    CoinMemcpyN(model2->primalColumnSolution(),
2499                numberColumns_,this->primalColumnSolution());
2500    CoinMemcpyN(model2->dualColumnSolution(),
2501                numberColumns_,this->dualColumnSolution());
2502    CoinMemcpyN(model2->statusArray(),
2503                numberColumns_+numberRows_,this->statusArray());
2504    objectiveValue_=model2->objectiveValue_;
2505    numberIterations_ =model2->numberIterations_;
2506    problemStatus_ =model2->problemStatus_;
2507    secondaryStatus_ =model2->secondaryStatus_;
2508    delete model2;
2509  }
2510  if (method!=ClpSolve::useBarrierNoCross&&
2511      method!=ClpSolve::useBarrier)
2512    setMaximumIterations(saveMaxIterations);
2513  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
2514                               "Errors","User stopped"};
2515  assert (finalStatus>=-1&&finalStatus<=5);
2516  handler_->message(CLP_TIMING,messages_)
2517    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
2518  handler_->printing(presolve==ClpSolve::presolveOn)
2519    <<timePresolve;
2520  handler_->printing(timeIdiot!=0.0)
2521    <<timeIdiot;
2522  handler_->message()<<CoinMessageEol;
2523  if (interrupt) 
2524    signal(SIGINT,saveSignal);
2525  perturbation_=savePerturbation;
2526  scalingFlag_=saveScaling;
2527  // If faking objective - put back correct one
2528  if (savedObjective) {
2529    delete objective_;
2530    objective_=savedObjective;
2531  }
2532  return finalStatus;
2533}
2534// General solve
2535int 
2536ClpSimplex::initialSolve()
2537{
2538  // Default so use dual
2539  ClpSolve options;
2540  return initialSolve(options);
2541}
2542// General dual solve
2543int 
2544ClpSimplex::initialDualSolve()
2545{
2546  ClpSolve options;
2547  // Use dual
2548  options.setSolveType(ClpSolve::useDual);
2549  return initialSolve(options);
2550}
2551// General dual solve
2552int 
2553ClpSimplex::initialPrimalSolve()
2554{
2555  ClpSolve options;
2556  // Use primal
2557  options.setSolveType(ClpSolve::usePrimal);
2558  return initialSolve(options);
2559}
2560// barrier solve, not to be followed by crossover
2561int 
2562ClpSimplex::initialBarrierNoCrossSolve()
2563{
2564  ClpSolve options;
2565  // Use primal
2566  options.setSolveType(ClpSolve::useBarrierNoCross);
2567  return initialSolve(options);
2568}
2569
2570// General barrier solve
2571int 
2572ClpSimplex::initialBarrierSolve()
2573{
2574  ClpSolve options;
2575  // Use primal
2576  options.setSolveType(ClpSolve::useBarrier);
2577  return initialSolve(options);
2578}
2579
2580// Default constructor
2581ClpSolve::ClpSolve (  )
2582{
2583  method_ = automatic;
2584  presolveType_=presolveOn;
2585  numberPasses_=5;
2586  int i;
2587  for (i=0;i<7;i++)
2588    options_[i]=0;
2589  // say no +-1 matrix
2590  options_[3]=1;
2591  for (i=0;i<7;i++)
2592    extraInfo_[i]=-1;
2593  independentOptions_[0]=0;
2594  // But switch off slacks
2595  independentOptions_[1]=512;
2596  // Substitute up to 3
2597  independentOptions_[2]=3;
2598 
2599}
2600// Constructor when you really know what you are doing
2601ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
2602             int numberPasses, int options[6],
2603             int extraInfo[6], int independentOptions[3])
2604{
2605  method_ = method;
2606  presolveType_=presolveType;
2607  numberPasses_=numberPasses;
2608  int i;
2609  for (i=0;i<6;i++)
2610    options_[i]=options[i];
2611  options_[6]=0;
2612  for (i=0;i<6;i++)
2613    extraInfo_[i]=extraInfo[i];
2614  extraInfo_[6]=0;
2615  for (i=0;i<3;i++)
2616    independentOptions_[i]=independentOptions[i];
2617}
2618
2619// Copy constructor.
2620ClpSolve::ClpSolve(const ClpSolve & rhs)
2621{
2622  method_ = rhs.method_;
2623  presolveType_=rhs.presolveType_;
2624  numberPasses_=rhs.numberPasses_;
2625  int i;
2626  for ( i=0;i<7;i++)
2627    options_[i]=rhs.options_[i];
2628  for ( i=0;i<7;i++)
2629    extraInfo_[i]=rhs.extraInfo_[i];
2630  for (i=0;i<3;i++)
2631    independentOptions_[i]=rhs.independentOptions_[i];
2632}
2633// Assignment operator. This copies the data
2634ClpSolve & 
2635ClpSolve::operator=(const ClpSolve & rhs)
2636{
2637  if (this != &rhs) {
2638    method_ = rhs.method_;
2639    presolveType_=rhs.presolveType_;
2640    numberPasses_=rhs.numberPasses_;
2641    int i;
2642    for (i=0;i<7;i++)
2643      options_[i]=rhs.options_[i];
2644    for (i=0;i<7;i++)
2645      extraInfo_[i]=rhs.extraInfo_[i];
2646    for (i=0;i<3;i++)
2647      independentOptions_[i]=rhs.independentOptions_[i];
2648  }
2649  return *this;
2650
2651}
2652// Destructor
2653ClpSolve::~ClpSolve (  )
2654{
2655}
2656// See header file for deatils
2657void 
2658ClpSolve::setSpecialOption(int which,int value,int extraInfo)
2659{
2660  options_[which]=value;
2661  extraInfo_[which]=extraInfo;
2662}
2663int 
2664ClpSolve::getSpecialOption(int which) const
2665{
2666  return options_[which];
2667}
2668
2669// Solve types
2670void 
2671ClpSolve::setSolveType(SolveType method, int extraInfo)
2672{
2673  method_=method;
2674}
2675
2676ClpSolve::SolveType
2677ClpSolve::getSolveType()
2678{
2679  return method_;
2680}
2681
2682// Presolve types
2683void 
2684ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
2685{
2686  presolveType_=amount;
2687  numberPasses_=extraInfo;
2688}
2689ClpSolve::PresolveType
2690ClpSolve::getPresolveType()
2691{
2692  return presolveType_;
2693}
2694// Extra info for idiot (or sprint)
2695int 
2696ClpSolve::getExtraInfo(int which) const
2697{
2698  return extraInfo_[which];
2699}
2700int 
2701ClpSolve::getPresolvePasses() const
2702{
2703  return numberPasses_;
2704}
2705/* Say to return at once if infeasible,
2706   default is to solve */
2707void 
2708ClpSolve::setInfeasibleReturn(bool trueFalse)
2709{
2710  independentOptions_[0]= trueFalse ? 1 : 0;
2711}
2712#include <string>
2713// Generates code for above constructor
2714void 
2715ClpSolve::generateCpp(FILE * fp)
2716{
2717  std::string solveType[] = {
2718    "ClpSolve::useDual",
2719    "ClpSolve::usePrimal",
2720    "ClpSolve::usePrimalorSprint",
2721    "ClpSolve::useBarrier",
2722    "ClpSolve::useBarrierNoCross",
2723    "ClpSolve::automatic",
2724    "ClpSolve::notImplemented"
2725  };
2726  std::string presolveType[] =  {
2727    "ClpSolve::presolveOn",
2728    "ClpSolve::presolveOff",
2729    "ClpSolve::presolveNumber",
2730    "ClpSolve::presolveNumberCost"
2731  };
2732  fprintf(fp,"3  ClpSolve::SolveType method = %s;\n",solveType[method_].c_str());
2733  fprintf(fp,"3  ClpSolve::PresolveType presolveType = %s;\n",
2734    presolveType[presolveType_].c_str());
2735  fprintf(fp,"3  int numberPasses = %d;\n",numberPasses_);
2736  fprintf(fp,"3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
2737    options_[0],options_[1],options_[2],
2738    options_[3],options_[4],options_[5]);
2739  fprintf(fp,"3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
2740    extraInfo_[0],extraInfo_[1],extraInfo_[2],
2741    extraInfo_[3],extraInfo_[4],extraInfo_[5]);
2742  fprintf(fp,"3  int independentOptions[] = {%d,%d,%d};\n",
2743    independentOptions_[0],independentOptions_[1],independentOptions_[2]);
2744  fprintf(fp,"3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
2745  fprintf(fp,"3                    options,extraInfo,independentOptions);\n");
2746}
2747//#############################################################################
2748#include "ClpNonLinearCost.hpp"
2749
2750ClpSimplexProgress::ClpSimplexProgress () 
2751{
2752  int i;
2753  for (i=0;i<CLP_PROGRESS;i++) {
2754    objective_[i] = COIN_DBL_MAX;
2755    infeasibility_[i] = -1.0; // set to an impossible value
2756    realInfeasibility_[i] = COIN_DBL_MAX;
2757    numberInfeasibilities_[i]=-1; 
2758    iterationNumber_[i]=-1;
2759  }
2760#ifdef CLP_PROGRESS_WEIGHT
2761  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2762    objectiveWeight_[i] = COIN_DBL_MAX;
2763    infeasibilityWeight_[i] = -1.0; // set to an impossible value
2764    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
2765    numberInfeasibilitiesWeight_[i]=-1; 
2766    iterationNumberWeight_[i]=-1;
2767  }
2768  drop_ =0.0;
2769  best_ =0.0;
2770#endif
2771  initialWeight_=0.0;
2772  for (i=0;i<CLP_CYCLE;i++) {
2773    //obj_[i]=COIN_DBL_MAX;
2774    in_[i]=-1;
2775    out_[i]=-1;
2776    way_[i]=0;
2777  }
2778  numberTimes_ = 0;
2779  numberBadTimes_ = 0;
2780  numberReallyBadTimes_ = 0;
2781  numberTimesFlagged_ = 0;
2782  model_ = NULL;
2783  oddState_=0;
2784}
2785
2786
2787//-----------------------------------------------------------------------------
2788
2789ClpSimplexProgress::~ClpSimplexProgress ()
2790{
2791}
2792// Copy constructor.
2793ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
2794{
2795  int i;
2796  for (i=0;i<CLP_PROGRESS;i++) {
2797    objective_[i] = rhs.objective_[i];
2798    infeasibility_[i] = rhs.infeasibility_[i];
2799    realInfeasibility_[i] = rhs.realInfeasibility_[i];
2800    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2801    iterationNumber_[i]=rhs.iterationNumber_[i];
2802  }
2803#ifdef CLP_PROGRESS_WEIGHT
2804  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2805    objectiveWeight_[i] = rhs.objectiveWeight_[i];
2806    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2807    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2808    numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i]; 
2809    iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
2810  }
2811  drop_ = rhs.drop_;
2812  best_ = rhs.best_;
2813#endif
2814  initialWeight_ = rhs.initialWeight_;
2815  for (i=0;i<CLP_CYCLE;i++) {
2816    //obj_[i]=rhs.obj_[i];
2817    in_[i]=rhs.in_[i];
2818    out_[i]=rhs.out_[i];
2819    way_[i]=rhs.way_[i];
2820  }
2821  numberTimes_ = rhs.numberTimes_;
2822  numberBadTimes_ = rhs.numberBadTimes_;
2823  numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
2824  numberTimesFlagged_ = rhs.numberTimesFlagged_;
2825  model_ = rhs.model_;
2826  oddState_=rhs.oddState_;
2827}
2828// Copy constructor.from model
2829ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
2830{
2831  model_ = model;
2832  reset();
2833  initialWeight_=0.0;
2834}
2835// Fill from model
2836void 
2837ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
2838{
2839  model_ = model;
2840  reset();
2841  initialWeight_=0.0;
2842}
2843// Assignment operator. This copies the data
2844ClpSimplexProgress & 
2845ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2846{
2847  if (this != &rhs) {
2848    int i;
2849    for (i=0;i<CLP_PROGRESS;i++) {
2850      objective_[i] = rhs.objective_[i];
2851      infeasibility_[i] = rhs.infeasibility_[i];
2852      realInfeasibility_[i] = rhs.realInfeasibility_[i];
2853      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2854      iterationNumber_[i]=rhs.iterationNumber_[i];
2855    }
2856#ifdef CLP_PROGRESS_WEIGHT
2857    for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2858      objectiveWeight_[i] = rhs.objectiveWeight_[i];
2859      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2860      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2861      numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i]; 
2862      iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
2863    }
2864    drop_ = rhs.drop_;
2865    best_ = rhs.best_;
2866#endif
2867    initialWeight_ = rhs.initialWeight_;
2868    for (i=0;i<CLP_CYCLE;i++) {
2869      //obj_[i]=rhs.obj_[i];
2870      in_[i]=rhs.in_[i];
2871      out_[i]=rhs.out_[i];
2872      way_[i]=rhs.way_[i];
2873    }
2874    numberTimes_ = rhs.numberTimes_;
2875    numberBadTimes_ = rhs.numberBadTimes_;
2876    numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
2877    numberTimesFlagged_ = rhs.numberTimesFlagged_;
2878    model_ = rhs.model_;
2879    oddState_=rhs.oddState_;
2880  }
2881  return *this;
2882}
2883// Seems to be something odd about exact comparison of doubles on linux
2884static bool equalDouble(double value1, double value2)
2885{
2886
2887  union { double d; int i[2]; } v1,v2;
2888  v1.d = value1;
2889  v2.d = value2;
2890  if (sizeof(int)*2==sizeof(double)) 
2891    return (v1.i[0]==v2.i[0]&&v1.i[1]==v2.i[1]);
2892  else
2893    return (v1.i[0]==v2.i[0]);
2894}
2895int
2896ClpSimplexProgress::looping()
2897{
2898  if (!model_)
2899    return -1;
2900  double objective = model_->rawObjectiveValue();
2901  if (model_->algorithm()<0)
2902    objective -= model_->bestPossibleImprovement();
2903  double infeasibility;
2904  double realInfeasibility=0.0;
2905  int numberInfeasibilities;
2906  int iterationNumber = model_->numberIterations();
2907  numberTimesFlagged_=0;
2908  if (model_->algorithm()<0) {
2909    // dual
2910    infeasibility = model_->sumPrimalInfeasibilities();
2911    numberInfeasibilities = model_->numberPrimalInfeasibilities();
2912  } else {
2913    //primal
2914    infeasibility = model_->sumDualInfeasibilities();
2915    realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
2916    numberInfeasibilities = model_->numberDualInfeasibilities();
2917  }
2918  int i;
2919  int numberMatched=0;
2920  int matched=0;
2921  int nsame=0;
2922  for (i=0;i<CLP_PROGRESS;i++) {
2923    bool matchedOnObjective = equalDouble(objective,objective_[i]);
2924    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
2925    bool matchedOnInfeasibilities = 
2926      (numberInfeasibilities==numberInfeasibilities_[i]);
2927   
2928    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
2929      matched |= (1<<i);
2930      // Check not same iteration
2931      if (iterationNumber!=iterationNumber_[i]) {
2932        numberMatched++;
2933        // here mainly to get over compiler bug?
2934        if (model_->messageHandler()->logLevel()>10)
2935          printf("%d %d %d %d %d loop check\n",i,numberMatched,
2936                 matchedOnObjective, matchedOnInfeasibility, 
2937                 matchedOnInfeasibilities);
2938      } else {
2939        // stuck but code should notice
2940        nsame++;
2941      }
2942    }
2943    if (i) {
2944      objective_[i-1] = objective_[i];
2945      infeasibility_[i-1] = infeasibility_[i];
2946      realInfeasibility_[i-1] = realInfeasibility_[i];
2947      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
2948      iterationNumber_[i-1]=iterationNumber_[i];
2949    }
2950  }
2951  objective_[CLP_PROGRESS-1] = objective;
2952  infeasibility_[CLP_PROGRESS-1] = infeasibility;
2953  realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
2954  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
2955  iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
2956  if (nsame==CLP_PROGRESS)
2957    numberMatched=CLP_PROGRESS; // really stuck
2958  if (model_->progressFlag())
2959    numberMatched=0;
2960  numberTimes_++;
2961  if (numberTimes_<10)
2962    numberMatched=0;
2963  // skip if just last time as may be checking something
2964  if (matched==(1<<(CLP_PROGRESS-1)))
2965    numberMatched=0;
2966  if (numberMatched) {
2967    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
2968      <<numberMatched
2969      <<matched
2970      <<numberTimes_
2971      <<CoinMessageEol;
2972    numberBadTimes_++;
2973    if (numberBadTimes_<10) {
2974      // make factorize every iteration
2975      model_->forceFactorization(1);
2976      if (numberBadTimes_<2) {
2977        startCheck(); // clear other loop check
2978        if (model_->algorithm()<0) {
2979          // dual - change tolerance
2980          model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
2981          // if infeasible increase dual bound
2982          if (model_->dualBound()<1.0e17) {
2983            model_->setDualBound(model_->dualBound()*1.1);
2984            static_cast<ClpSimplexDual *>(model_)->resetFakeBounds(0);
2985          }
2986        } else {
2987          // primal - change tolerance
2988          if (numberBadTimes_>3)
2989            model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
2990          // if infeasible increase infeasibility cost
2991          if (model_->nonLinearCost()->numberInfeasibilities()&&
2992              model_->infeasibilityCost()<1.0e17) {
2993            model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
2994          }
2995        }
2996      } else {
2997        // flag
2998        int iSequence;
2999        if (model_->algorithm()<0) {
3000          // dual
3001          if (model_->dualBound()>1.0e14) 
3002            model_->setDualBound(1.0e14);
3003          iSequence=in_[CLP_CYCLE-1];
3004        } else {
3005          // primal
3006          if (model_->infeasibilityCost()>1.0e14) 
3007            model_->setInfeasibilityCost(1.0e14);
3008          iSequence=out_[CLP_CYCLE-1];
3009        }
3010        if (iSequence>=0) {
3011          char x = model_->isColumn(iSequence) ? 'C' :'R';
3012          if (model_->messageHandler()->logLevel()>=63)
3013            model_->messageHandler()->message(CLP_SIMPLEX_FLAG,model_->messages())
3014              <<x<<model_->sequenceWithin(iSequence)
3015              <<CoinMessageEol;
3016          // if Gub then needs to be sequenceIn_
3017          int save=model_->sequenceIn();
3018          model_->setSequenceIn(iSequence);
3019          model_->setFlagged(iSequence);
3020          model_->setSequenceIn(save);
3021          //printf("flagging %d from loop\n",iSequence);
3022          startCheck();
3023        } else {
3024          // Give up
3025          if (model_->messageHandler()->logLevel()>=63)
3026            printf("***** All flagged?\n");
3027          return 4;
3028        }
3029        // reset
3030        numberBadTimes_=2;
3031      }
3032      return -2;
3033    } else {
3034      // look at solution and maybe declare victory
3035      if (infeasibility<1.0e-4) {
3036        return 0;
3037      } else {
3038        model_->messageHandler()->message(CLP_LOOP,model_->messages())
3039          <<CoinMessageEol;
3040#ifndef NDEBUG
3041        printf("debug loop ClpSimplex A\n");
3042        abort();
3043#endif
3044        return 3;
3045      }
3046    }
3047  }
3048  return -1;
3049}
3050// Resets as much as possible
3051void 
3052ClpSimplexProgress::reset()
3053{
3054  int i;
3055  for (i=0;i<CLP_PROGRESS;i++) {
3056    if (model_->algorithm()>=0)
3057      objective_[i] = COIN_DBL_MAX;
3058    else
3059      objective_[i] = -COIN_DBL_MAX;
3060    infeasibility_[i] = -1.0; // set to an impossible value
3061    realInfeasibility_[i] = COIN_DBL_MAX;
3062    numberInfeasibilities_[i]=-1; 
3063    iterationNumber_[i]=-1;
3064  }
3065#ifdef CLP_PROGRESS_WEIGHT
3066  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
3067    objectiveWeight_[i] = COIN_DBL_MAX;
3068    infeasibilityWeight_[i] = -1.0; // set to an impossible value
3069    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3070    numberInfeasibilitiesWeight_[i]=-1; 
3071    iterationNumberWeight_[i]=-1;
3072  }
3073  drop_ =0.0;
3074  best_ =0.0;
3075#endif
3076  for (i=0;i<CLP_CYCLE;i++) {
3077    //obj_[i]=COIN_DBL_MAX;
3078    in_[i]=-1;
3079    out_[i]=-1;
3080    way_[i]=0;
3081  }
3082  numberTimes_ = 0;
3083  numberBadTimes_ = 0;
3084  numberReallyBadTimes_ = 0;
3085  numberTimesFlagged_ = 0;
3086  oddState_=0;
3087}
3088// Returns previous objective (if -1) - current if (0)
3089double 
3090ClpSimplexProgress::lastObjective(int back) const
3091{
3092  return objective_[CLP_PROGRESS-1-back];
3093}
3094// Returns previous infeasibility (if -1) - current if (0)
3095double 
3096ClpSimplexProgress::lastInfeasibility(int back) const
3097{
3098  return realInfeasibility_[CLP_PROGRESS-1-back];
3099}
3100// Sets real primal infeasibility
3101void
3102ClpSimplexProgress::setInfeasibility(double value)
3103{
3104  for (int i=1;i<CLP_PROGRESS;i++) 
3105    realInfeasibility_[i-1] = realInfeasibility_[i];
3106  realInfeasibility_[CLP_PROGRESS-1]=value;
3107}
3108// Modify objective e.g. if dual infeasible in dual
3109void 
3110ClpSimplexProgress::modifyObjective(double value)
3111{
3112  objective_[CLP_PROGRESS-1]=value;
3113}
3114// Returns previous iteration number (if -1) - current if (0)
3115int 
3116ClpSimplexProgress::lastIterationNumber(int back) const
3117{
3118  return iterationNumber_[CLP_PROGRESS-1-back];
3119}
3120// clears iteration numbers (to switch off panic)
3121void 
3122ClpSimplexProgress::clearIterationNumbers()
3123{
3124  for (int i=0;i<CLP_PROGRESS;i++) 
3125    iterationNumber_[i]=-1;
3126}
3127// Start check at beginning of whileIterating
3128void 
3129ClpSimplexProgress::startCheck()
3130{
3131  int i;
3132  for (i=0;i<CLP_CYCLE;i++) {
3133    //obj_[i]=COIN_DBL_MAX;
3134    in_[i]=-1;
3135    out_[i]=-1;
3136    way_[i]=0;
3137  }
3138}
3139// Returns cycle length in whileIterating
3140int 
3141ClpSimplexProgress::cycle(int in, int out,int wayIn,int wayOut)
3142{
3143  int i;
3144#if 0
3145  if (model_->numberIterations()>206571) {
3146    printf("in %d out %d\n",in,out);
3147    for (i=0;i<CLP_CYCLE;i++)
3148      printf("cy %d in %d out %d\n",i,in_[i],out_[i]);
3149  }
3150#endif
3151  int matched=0;
3152  // first see if in matches any out
3153  for (i=1;i<CLP_CYCLE;i++) {
3154    if (in==out_[i]) {
3155      // even if flip then suspicious
3156      matched=-1;
3157      break;
3158    }
3159  }
3160#if 0
3161  if (!matched||in_[0]<0) {
3162    // can't be cycle
3163    for (i=0;i<CLP_CYCLE-1;i++) {
3164      //obj_[i]=obj_[i+1];
3165      in_[i]=in_[i+1];
3166      out_[i]=out_[i+1];
3167      way_[i]=way_[i+1];
3168    }
3169  } else {
3170    // possible cycle
3171    matched=0;
3172    for (i=0;i<CLP_CYCLE-1;i++) {
3173      int k;
3174      char wayThis = way_[i];
3175      int inThis = in_[i];
3176      int outThis = out_[i];
3177      //double objThis = obj_[i];
3178      for(k=i+1;k<CLP_CYCLE;k++) {
3179        if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
3180          int distance = k-i;
3181          if (k+distance<CLP_CYCLE) {
3182            // See if repeats
3183            int j=k+distance;
3184            if (inThis==in_[j]&&outThis==out_[j]&&wayThis==way_[j]) {
3185              matched=distance;
3186              break;
3187            }
3188          } else {
3189            matched=distance;
3190            break;
3191          }
3192        }
3193      }
3194      //obj_[i]=obj_[i+1];
3195      in_[i]=in_[i+1];
3196      out_[i]=out_[i+1];
3197      way_[i]=way_[i+1];
3198    }
3199  }
3200#else
3201  if (matched&&in_[0]>=0) {
3202    // possible cycle - only check [0] against all
3203    matched=0;
3204    int nMatched=0;
3205    char way0 = way_[0];
3206    int in0 = in_[0];
3207    int out0 = out_[0];
3208    //double obj0 = obj_[i];
3209    for(int k=1;k<CLP_CYCLE-4;k++) {
3210      if (in0==in_[k]&&out0==out_[k]&&way0==way_[k]) {
3211        nMatched++;
3212        // See if repeats
3213        int end = CLP_CYCLE-k;
3214        int j;
3215        for ( j=1;j<end;j++) {
3216          if (in_[j+k]!=in_[j]||out_[j+k]!=out_[j]||way_[j+k]!=way_[j])
3217            break;
3218        }
3219        if (j==end) {
3220          matched=k;
3221          break;
3222        }
3223      }
3224    }
3225    // If three times then that is too much even if not regular
3226    if (matched<=0&&nMatched>1)
3227      matched=100;
3228  }
3229  for (i=0;i<CLP_CYCLE-1;i++) {
3230    //obj_[i]=obj_[i+1];
3231    in_[i]=in_[i+1];
3232    out_[i]=out_[i+1];
3233    way_[i]=way_[i+1];
3234  }
3235#endif
3236  int way = 1-wayIn+4*(1-wayOut);
3237  //obj_[i]=model_->objectiveValue();
3238  in_[CLP_CYCLE-1]=in;
3239  out_[CLP_CYCLE-1]=out;
3240  way_[CLP_CYCLE-1]=static_cast<char>(way);
3241  return matched;
3242}
3243#include "CoinStructuredModel.hpp"
3244// Solve using structure of model and maybe in parallel
3245int 
3246ClpSimplex::solve(CoinStructuredModel * model)
3247{
3248  // analyze structure
3249  int numberRowBlocks=model->numberRowBlocks();
3250  int numberColumnBlocks = model->numberColumnBlocks();
3251  int numberElementBlocks = model->numberElementBlocks();
3252  if (numberElementBlocks==1) {
3253    loadProblem(*model,false);
3254    return dual();
3255  }
3256  // For now just get top level structure
3257  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3258  for (int i=0;i<numberElementBlocks;i++) {
3259    CoinStructuredModel * subModel = 
3260      dynamic_cast<CoinStructuredModel *>(model->block(i));
3261    CoinModel * thisBlock;
3262    if (subModel) {
3263      thisBlock=subModel->coinModelBlock(blockInfo[i]);
3264      model->setCoinModel(thisBlock,i);
3265    } else {
3266      thisBlock = dynamic_cast<CoinModel *>(model->block(i));
3267      assert (thisBlock);
3268      // just fill in info
3269      CoinModelBlockInfo info = CoinModelBlockInfo();
3270      int whatsSet = thisBlock->whatIsSet();
3271      info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
3272      info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
3273      info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
3274      info.integer = ((whatsSet&32)!=0) ? 1 : 0;
3275      info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
3276      info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
3277      // Which block
3278      int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
3279      info.rowBlock=iRowBlock;
3280      int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
3281      info.columnBlock=iColumnBlock;
3282      blockInfo[i] = info;
3283    }
3284  }
3285  int * rowCounts = new int [numberRowBlocks];
3286  CoinZeroN(rowCounts,numberRowBlocks);
3287  int * columnCounts = new int [numberColumnBlocks+1];
3288  CoinZeroN(columnCounts,numberColumnBlocks);
3289  int decomposeType=0;
3290  for (int i=0;i<numberElementBlocks;i++) {
3291    int iRowBlock = blockInfo[i].rowBlock;
3292    int iColumnBlock = blockInfo[i].columnBlock;
3293    rowCounts[iRowBlock]++;
3294    columnCounts[iColumnBlock]++;
3295  }
3296  if (numberRowBlocks==numberColumnBlocks||
3297      numberRowBlocks==numberColumnBlocks+1) {
3298    // could be Dantzig-Wolfe
3299    int numberG1=0;
3300    for (int i=0;i<numberRowBlocks;i++) {
3301      if (rowCounts[i]>1)
3302        numberG1++;
3303    }
3304    bool masterColumns = (numberColumnBlocks==numberRowBlocks);
3305    if ((masterColumns&&numberElementBlocks==2*numberRowBlocks-1)
3306        ||(!masterColumns&&numberElementBlocks==2*numberRowBlocks)) {
3307      if (numberG1<2)
3308        decomposeType=1;
3309    }
3310  }
3311  if (!decomposeType&&(numberRowBlocks==numberColumnBlocks||
3312                       numberRowBlocks==numberColumnBlocks-1)) {
3313    // could be Benders
3314    int numberG1=0;
3315    for (int i=0;i<numberColumnBlocks;i++) {
3316      if (columnCounts[i]>1)
3317        numberG1++;
3318    }
3319    bool masterRows = (numberColumnBlocks==numberRowBlocks);
3320    if ((masterRows&&numberElementBlocks==2*numberColumnBlocks-1)
3321        ||(!masterRows&&numberElementBlocks==2*numberColumnBlocks)) {
3322      if (numberG1<2)
3323        decomposeType=2;
3324    }
3325  }
3326  delete [] rowCounts;
3327  delete [] columnCounts;
3328  delete [] blockInfo;
3329  // decide what to do
3330  switch (decomposeType) {
3331    // No good
3332  case 0:
3333    loadProblem(*model,false);
3334    return dual();
3335    // DW
3336  case 1:
3337    return solveDW(model);
3338    // Benders
3339  case 2:
3340    return solveBenders(model);
3341  }
3342  return 0; // to stop compiler warning
3343}
3344/* This loads a model from a CoinStructuredModel object - returns number of errors.
3345   If originalOrder then keep to order stored in blocks,
3346   otherwise first column/rows correspond to first block - etc.
3347   If keepSolution true and size is same as current then
3348   keeps current status and solution
3349*/
3350int 
3351ClpSimplex::loadProblem (  CoinStructuredModel & coinModel,
3352                           bool originalOrder,
3353                           bool keepSolution)
3354{
3355  unsigned char * status = NULL;
3356  double * psol = NULL;
3357  double * dsol = NULL;
3358  int numberRows=coinModel.numberRows();
3359  int numberColumns = coinModel.numberColumns();
3360  int numberRowBlocks=coinModel.numberRowBlocks();
3361  int numberColumnBlocks = coinModel.numberColumnBlocks();
3362  int numberElementBlocks = coinModel.numberElementBlocks();
3363  if (status_&&numberRows_&&numberRows_==numberRows&&
3364      numberColumns_==numberColumns&&keepSolution) {
3365    status = new unsigned char [numberRows_+numberColumns_];
3366    CoinMemcpyN(status_,numberRows_+numberColumns_,status);
3367    psol = new double [numberRows_+numberColumns_];
3368    CoinMemcpyN(columnActivity_,numberColumns_,psol);
3369    CoinMemcpyN(rowActivity_,numberRows_,psol+numberColumns_);
3370    dsol = new double [numberRows_+numberColumns_];
3371    CoinMemcpyN(reducedCost_,numberColumns_,dsol);
3372    CoinMemcpyN(dual_,numberRows_,dsol+numberColumns_);
3373  }
3374  int returnCode=0;
3375  double * rowLower = new double [numberRows];
3376  double * rowUpper = new double [numberRows];
3377  double * columnLower = new double [numberColumns];
3378  double * columnUpper = new double [numberColumns];
3379  double * objective = new double [numberColumns];
3380  int * integerType = new int [numberColumns];
3381  CoinBigIndex numberElements=0;
3382  // Bases for blocks
3383  int * rowBase = new int[numberRowBlocks];
3384  CoinFillN(rowBase,numberRowBlocks,-1);
3385  // And row to put it
3386  int * whichRow = new int [numberRows];
3387  int * columnBase = new int[numberColumnBlocks];
3388  CoinFillN(columnBase,numberColumnBlocks,-1);
3389  // And column to put it
3390  int * whichColumn = new int [numberColumns];
3391  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
3392    CoinModel * block = coinModel.coinBlock(iBlock);
3393    numberElements += block->numberElements();
3394    //and set up elements etc
3395    double * associated = block->associatedArray();
3396    // If strings then do copies
3397    if (block->stringsExist()) 
3398      returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
3399                                          objective, integerType,associated);
3400    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3401    int iRowBlock = info.rowBlock;
3402    int iColumnBlock = info.columnBlock;
3403    if (rowBase[iRowBlock]<0) {
3404      rowBase[iRowBlock]=block->numberRows();
3405      // Save block number
3406      whichRow[numberRows-numberRowBlocks+iRowBlock]= iBlock;
3407    } else {
3408      assert(rowBase[iRowBlock]==block->numberRows());
3409    }
3410    if (columnBase[iColumnBlock]<0) {
3411      columnBase[iColumnBlock]=block->numberColumns();
3412      // Save block number
3413      whichColumn[numberColumns-numberColumnBlocks+iColumnBlock]= iBlock;
3414    } else {
3415      assert(columnBase[iColumnBlock]==block->numberColumns());
3416    }
3417  }
3418  // Fill arrays with defaults
3419  CoinFillN(rowLower,numberRows,-COIN_DBL_MAX);
3420  CoinFillN(rowUpper,numberRows,COIN_DBL_MAX);
3421  CoinFillN(columnLower,numberColumns,0.0);
3422  CoinFillN(columnUpper,numberColumns,COIN_DBL_MAX);
3423  CoinFillN(objective,numberColumns,0.0);
3424  CoinFillN(integerType,numberColumns,0);
3425  int n=0;
3426  for (int iBlock=0;iBlock<numberRowBlocks;iBlock++) {
3427    int k = rowBase[iBlock];
3428    rowBase[iBlock]=n;
3429    assert (k>=0);
3430    // block number
3431    int jBlock = whichRow[numberRows-numberRowBlocks+iBlock];
3432    if (originalOrder) {
3433      memcpy(whichRow+n,coinModel.coinBlock(jBlock)->originalRows(),k*sizeof(int));
3434    } else {
3435      CoinIotaN(whichRow+n,k,n);
3436    }
3437    n+=k;
3438  }
3439  assert (n==numberRows);
3440  n=0;
3441  for (int iBlock=0;iBlock<numberColumnBlocks;iBlock++) {
3442    int k = columnBase[iBlock];
3443    columnBase[iBlock]=n;
3444    assert (k>=0);
3445    if (k) {
3446      // block number
3447      int jBlock = whichColumn[numberColumns-numberColumnBlocks+iBlock];
3448      if (originalOrder) {
3449        memcpy(whichColumn+n,coinModel.coinBlock(jBlock)->originalColumns(),
3450               k*sizeof(int));
3451      } else {
3452        CoinIotaN(whichColumn+n,k,n);
3453      }
3454      n+=k;
3455    }
3456  }
3457  assert (n==numberColumns);
3458  bool gotIntegers=false;
3459  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
3460    CoinModel * block = coinModel.coinBlock(iBlock);
3461    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3462    int iRowBlock = info.rowBlock;
3463    int iRowBase = rowBase[iRowBlock];
3464    int iColumnBlock = info.columnBlock;
3465    int iColumnBase = columnBase[iColumnBlock];
3466    if (info.rhs) {
3467      int nRows = block->numberRows();
3468      const double * lower = block->rowLowerArray();
3469      const double * upper = block->rowUpperArray();
3470      for (int i=0;i<nRows;i++) {
3471        int put = whichRow[i+iRowBase];
3472        rowLower[put] = lower[i];
3473        rowUpper[put] = upper[i];
3474      }
3475    }
3476    if (info.bounds) {
3477      int nColumns = block->numberColumns();
3478      const double * lower = block->columnLowerArray();
3479      const double * upper = block->columnUpperArray();
3480      const double * obj = block->objectiveArray();
3481      for (int i=0;i<nColumns;i++) {
3482        int put = whichColumn[i+iColumnBase];
3483        columnLower[put] = lower[i];
3484        columnUpper[put] = upper[i];
3485        objective[put] = obj[i];
3486      }
3487    }
3488    if (info.integer) {
3489      gotIntegers=true;
3490      int nColumns = block->numberColumns();
3491      const int * type = block->integerTypeArray();
3492      for (int i=0;i<nColumns;i++) {
3493        int put = whichColumn[i+iColumnBase];
3494        integerType[put] = type[i];
3495      }
3496    }
3497  }
3498  gutsOfLoadModel(numberRows, numberColumns,
3499                  columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
3500  delete [] rowLower;
3501  delete [] rowUpper;
3502  delete [] columnLower;
3503  delete [] columnUpper;
3504  delete [] objective;
3505  // Do integers if wanted
3506  if (gotIntegers) {
3507    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3508      if (integerType[iColumn])
3509        setInteger(iColumn);
3510    }
3511  }
3512  delete [] integerType;
3513  setObjectiveOffset(coinModel.objectiveOffset());
3514  // Space for elements
3515  int * row = new int[numberElements];
3516  int * column = new int[numberElements];
3517  double * element = new double[numberElements];
3518  numberElements=0;
3519  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
3520    CoinModel * block = coinModel.coinBlock(iBlock);
3521    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3522    int iRowBlock = info.rowBlock;
3523    int iRowBase = rowBase[iRowBlock];
3524    int iColumnBlock = info.columnBlock;
3525    int iColumnBase = columnBase[iColumnBlock];
3526    if (info.rowName) {
3527      int numberItems = block->rowNames()->numberItems();
3528      assert( block->numberRows()>=numberItems);
3529      if (numberItems) {
3530        const char *const * rowNames=block->rowNames()->names();
3531        for (int i=0;i<numberItems;i++) {
3532          int put = whichRow[i+iRowBase];
3533          std::string name = rowNames[i];
3534          setRowName(put,name);
3535        }
3536      }
3537    }
3538    if (info.columnName) {
3539      int numberItems = block->columnNames()->numberItems();
3540      assert( block->numberColumns()>=numberItems);
3541      if (numberItems) {
3542        const char *const * columnNames=block->columnNames()->names();
3543        for (int i=0;i<numberItems;i++) {
3544          int put = whichColumn[i+iColumnBase];
3545          std::string name = columnNames[i];
3546          setColumnName(put,name);
3547        }
3548      }
3549    }
3550    if (info.matrix) {
3551      CoinPackedMatrix matrix2;
3552      const CoinPackedMatrix * matrix = block->packedMatrix();
3553      if (!matrix) {
3554        double * associated = block->associatedArray();
3555        block->createPackedMatrix(matrix2,associated);
3556        matrix = &matrix2;
3557      }
3558      // get matrix data pointers
3559      const int * row2 = matrix->getIndices();
3560      const CoinBigIndex * columnStart = matrix->getVectorStarts();
3561      const double * elementByColumn = matrix->getElements();
3562      const int * columnLength = matrix->getVectorLengths(); 
3563      int n = matrix->getNumCols();
3564      assert (matrix->isColOrdered());
3565      for (int iColumn=0;iColumn<n;iColumn++) {
3566        CoinBigIndex j;
3567        int jColumn = whichColumn[iColumn+iColumnBase];
3568        for (j=columnStart[iColumn];
3569             j<columnStart[iColumn]+columnLength[iColumn];j++) {
3570          row[numberElements]=whichRow[row2[j]+iRowBase];
3571          column[numberElements]=jColumn;
3572          element[numberElements++]=elementByColumn[j];
3573        }
3574      }
3575    }
3576  }
3577  delete [] whichRow;
3578  delete [] whichColumn;
3579  delete [] rowBase;
3580  delete [] columnBase;
3581  CoinPackedMatrix * matrix =
3582    new CoinPackedMatrix (true,row,column,element,numberElements);
3583  matrix_ = new ClpPackedMatrix(matrix);
3584  matrix_->setDimensions(numberRows,numberColumns);
3585  delete [] row;
3586  delete [] column;
3587  delete [] element;
3588  createStatus();
3589  if (status) {
3590    // copy back
3591    CoinMemcpyN(status,numberRows_+numberColumns_,status_);
3592    CoinMemcpyN(psol,numberColumns_,columnActivity_);
3593    CoinMemcpyN(psol+numberColumns_,numberRows_,rowActivity_);
3594    CoinMemcpyN(dsol,numberColumns_,reducedCost_);
3595    CoinMemcpyN(dsol+numberColumns_,numberRows_,dual_);
3596    delete [] status;
3597    delete [] psol;
3598    delete [] dsol;
3599  }
3600  optimizationDirection_ = coinModel.optimizationDirection(); 
3601  return returnCode;
3602}
3603/*  If input negative scales objective so maximum <= -value
3604    and returns scale factor used.  If positive unscales and also
3605    redoes dual stuff
3606*/
3607double 
3608ClpSimplex::scaleObjective(double value)
3609{
3610  double * obj = objective(); 
3611  double largest=0.0;
3612  if (value<0.0) {
3613    value = - value;
3614    for (int i=0;i<numberColumns_;i++) {
3615      largest = CoinMax(largest,fabs(obj[i]));
3616    }
3617    if (largest>value) {
3618      double scaleFactor = value/largest;
3619      for (int i=0;i<numberColumns_;i++) {
3620        obj[i] *= scaleFactor;
3621        reducedCost_[i] *= scaleFactor;
3622      }
3623      for (int i=0;i<numberRows_;i++) {
3624        dual_[i] *= scaleFactor;
3625      }
3626      largest /= value;
3627    } else {
3628      // no need
3629      largest=1.0;
3630    }
3631  } else {
3632    // at end
3633    if (value!=1.0) {
3634      for (int i=0;i<numberColumns_;i++) {
3635        obj[i] *= value;
3636        reducedCost_[i] *= value;
3637      }
3638      for (int i=0;i<numberRows_;i++) {
3639        dual_[i] *= value;
3640      }
3641      computeObjectiveValue();
3642    }
3643  }
3644  return largest;
3645}
3646// Solve using Dantzig-Wolfe decomposition and maybe in parallel
3647int 
3648ClpSimplex::solveDW(CoinStructuredModel * model)
3649{
3650  double time1 = CoinCpuTime();
3651  int numberColumns = model->numberColumns();
3652  int numberRowBlocks=model->numberRowBlocks();
3653  int numberColumnBlocks = model->numberColumnBlocks();
3654  int numberElementBlocks = model->numberElementBlocks();
3655  // We already have top level structure
3656  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3657  for (int i=0;i<numberElementBlocks;i++) {
3658    CoinModel * thisBlock = model->coinBlock(i);
3659    assert (thisBlock);
3660    // just fill in info
3661    CoinModelBlockInfo info = CoinModelBlockInfo();
3662    int whatsSet = thisBlock->whatIsSet();
3663    info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
3664    info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
3665    info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
3666    info.integer = ((whatsSet&32)!=0) ? 1 : 0;
3667    info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
3668    info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
3669    // Which block
3670    int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
3671    info.rowBlock=iRowBlock;
3672    int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
3673    info.columnBlock=iColumnBlock;
3674    blockInfo[i] = info;
3675  }
3676  // make up problems
3677  int numberBlocks=numberRowBlocks-1;
3678  // Find master rows and columns
3679  int * rowCounts = new int [numberRowBlocks];
3680  CoinZeroN(rowCounts,numberRowBlocks);
3681  int * columnCounts = new int [numberColumnBlocks+1];
3682  CoinZeroN(columnCounts,numberColumnBlocks);
3683  int iBlock;
3684  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
3685    int iRowBlock = blockInfo[iBlock].rowBlock;
3686    rowCounts[iRowBlock]++;
3687    int iColumnBlock =blockInfo[iBlock].columnBlock;
3688    columnCounts[iColumnBlock]++;
3689  }
3690  int * whichBlock = new int [numberElementBlocks];
3691  int masterRowBlock=-1;
3692  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
3693    if (rowCounts[iBlock]>1) {
3694      if (masterRowBlock==-1) {
3695        masterRowBlock=iBlock;
3696      } else {
3697        // Can't decode
3698        masterRowBlock=-2;
3699        break;
3700      }
3701    }
3702  }
3703  int masterColumnBlock=-1;
3704  int kBlock=0;
3705  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
3706    int count=columnCounts[iBlock];
3707    columnCounts[iBlock]=kBlock;
3708    kBlock += count;
3709  }
3710  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
3711    int iColumnBlock = blockInfo[iBlock].columnBlock;
3712    whichBlock[columnCounts[iColumnBlock]]=iBlock;
3713    columnCounts[iColumnBlock]++;
3714  }
3715  for (iBlock = numberColumnBlocks-1;iBlock>=0;iBlock--) 
3716    columnCounts[iBlock+1]=columnCounts[iBlock];
3717  columnCounts[0]=0;
3718  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
3719    int count=columnCounts[iBlock+1]-columnCounts[iBlock];
3720    if (count==1) {
3721      int kBlock = whichBlock[columnCounts[iBlock]];
3722      int iRowBlock = blockInfo[kBlock].rowBlock;
3723      if (iRowBlock==masterRowBlock) {
3724        if (masterColumnBlock==-1) {
3725          masterColumnBlock=iBlock;
3726        } else {
3727          // Can't decode
3728          masterColumnBlock=-2;
3729          break;
3730        }
3731      }
3732    }
3733  }
3734  if (masterRowBlock<0||masterColumnBlock==-2) {
3735    // What now
3736    abort();
3737  }
3738  delete [] rowCounts;
3739  // create all data
3740  const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
3741  ClpSimplex * sub = new ClpSimplex [numberBlocks];
3742  ClpSimplex master;
3743  // Set offset
3744  master.setObjectiveOffset(model->objectiveOffset());
3745  kBlock=0;
3746  int masterBlock=-1;
3747  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
3748    top[kBlock]=NULL;
3749    int start=columnCounts[iBlock];
3750    int end=columnCounts[iBlock+1];
3751    assert (end-start<=2);
3752    for (int j=start;j<end;j++) {
3753      int jBlock = whichBlock[j];
3754      int iRowBlock = blockInfo[jBlock].rowBlock;
3755      int iColumnBlock =blockInfo[jBlock].columnBlock;
3756      assert (iColumnBlock==iBlock);
3757      if (iColumnBlock!=masterColumnBlock&&iRowBlock==masterRowBlock) {
3758        // top matrix
3759        top[kBlock]=model->coinBlock(jBlock)->packedMatrix();
3760      } else {
3761        const CoinPackedMatrix * matrix
3762          = model->coinBlock(jBlock)->packedMatrix();
3763        // Get pointers to arrays
3764        const double * rowLower;
3765        const double * rowUpper;
3766        const double * columnLower;
3767        const double * columnUpper;
3768        const double * objective;
3769        model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
3770                     columnLower,columnUpper,objective);
3771        if (iColumnBlock!=masterColumnBlock) {
3772          // diagonal block
3773          sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
3774                                  objective,rowLower,rowUpper);
3775          if (true) {
3776            double * lower = sub[kBlock].columnLower();
3777            double * upper = sub[kBlock].columnUpper();
3778            int n = sub[kBlock].numberColumns();
3779            for (int i=0;i<n;i++) {
3780              lower[i]=CoinMax(-1.0e8,lower[i]);
3781              upper[i]=CoinMin(1.0e8,upper[i]);
3782            }
3783          }
3784          if (optimizationDirection_<0.0) {
3785            double * obj = sub[kBlock].objective();
3786            int n = sub[kBlock].numberColumns();
3787            for (int i=0;i<n;i++)
3788              obj[i] = - obj[i];
3789          }
3790          if (this->factorizationFrequency()==200) {
3791            // User did not touch preset
3792            sub[kBlock].defaultFactorizationFrequency();
3793          } else {
3794            // make sure model has correct value
3795            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
3796          }
3797          sub[kBlock].setPerturbation(50);
3798          // Set columnCounts to be diagonal block index for cleanup
3799          columnCounts[kBlock]=jBlock;
3800        } else {
3801          // master
3802          masterBlock=jBlock;
3803          master.loadProblem(*matrix,columnLower,columnUpper,
3804                             objective,rowLower,rowUpper);
3805          if (optimizationDirection_<0.0) {
3806            double * obj = master.objective();
3807            int n = master.numberColumns();
3808            for (int i=0;i<n;i++)
3809              obj[i] = - obj[i];
3810          }
3811        }
3812      }
3813    }
3814    if (iBlock!=masterColumnBlock) 
3815      kBlock++;
3816  }
3817  delete [] whichBlock;
3818  delete [] blockInfo;
3819  // For now master must have been defined (does not have to have columns)
3820  assert (master.numberRows());
3821  assert (masterBlock>=0);
3822  int numberMasterRows = master.numberRows();
3823  // Overkill in terms of space
3824  int spaceNeeded = CoinMax(numberBlocks*(numberMasterRows+1),
3825                            2*numberMasterRows);
3826  int * rowAdd = new int[spaceNeeded];
3827  double * elementAdd = new double[spaceNeeded];
3828  spaceNeeded = numberBlocks;
3829  int * columnAdd = new int[spaceNeeded+1];
3830  double * objective = new double[spaceNeeded];
3831  // Add in costed slacks
3832  int firstArtificial = master.numberColumns();
3833  int lastArtificial = firstArtificial;
3834  if (true) {
3835    const double * lower = master.rowLower();
3836    const double * upper = master.rowUpper();
3837    int kCol=0;
3838    for (int iRow=0;iRow<numberMasterRows;iRow++) {
3839      if (lower[iRow]>-1.0e10) {
3840        rowAdd[kCol]=iRow;
3841        elementAdd[kCol++]=1.0;
3842      }
3843      if (upper[iRow]<1.0e10) {
3844        rowAdd[kCol]=iRow;
3845        elementAdd[kCol++]=-1.0;
3846      }
3847    }
3848    if (kCol>spaceNeeded) {
3849      spaceNeeded = kCol;
3850      delete [] columnAdd;
3851      delete [] objective;
3852      columnAdd = new int[spaceNeeded+1];
3853      objective = new double[spaceNeeded];
3854    }
3855    for (int i=0;i<kCol;i++) {
3856      columnAdd[i]=i;
3857      objective[i]=1.0e13;
3858    }
3859    columnAdd[kCol]=kCol;
3860    master.addColumns(kCol,NULL,NULL,objective,
3861                      columnAdd,rowAdd,elementAdd);
3862    lastArtificial = master.numberColumns();
3863  }
3864  int maxPass=500;
3865  int iPass;
3866  double lastObjective=1.0e31;
3867  // Create convexity rows for proposals
3868  int numberMasterColumns = master.numberColumns();
3869  master.resize(numberMasterRows+numberBlocks,numberMasterColumns);
3870  if (this->factorizationFrequency()==200) {
3871    // User did not touch preset
3872    master.defaultFactorizationFrequency();
3873  } else {
3874    // make sure model has correct value
3875    master.setFactorizationFrequency(this->factorizationFrequency());
3876  }
3877  master.setPerturbation(50);
3878  // Arrays to say which block and when created
3879  int maximumColumns = 2*numberMasterRows+10*numberBlocks;
3880  whichBlock = new int[maximumColumns];
3881  int * when = new int[maximumColumns];
3882  int numberColumnsGenerated=numberBlocks;
3883  // fill in rhs and add in artificials
3884  {
3885    double * rowLower = master.rowLower();
3886    double * rowUpper = master.rowUpper();
3887    int iBlock;
3888    columnAdd[0] = 0;
3889    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
3890      int iRow = iBlock + numberMasterRows;;
3891      rowLower[iRow]=1.0;
3892      rowUpper[iRow]=1.0;
3893      rowAdd[iBlock] = iRow;
3894      elementAdd[iBlock] = 1.0;
3895      objective[iBlock] = 1.0e13;
3896      columnAdd[iBlock+1] = iBlock+1;
3897      when[iBlock]=-1;
3898      whichBlock[iBlock] = iBlock;
3899    }
3900    master.addColumns(numberBlocks,NULL,NULL,objective,
3901                      columnAdd, rowAdd, elementAdd);
3902  }
3903  // and resize matrix to double check clp will be happy
3904  //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
3905  //                     numberMasterColumns+numberBlocks);
3906  std::cout<<"Time to decompose "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
3907  for (iPass=0;iPass<maxPass;iPass++) {
3908    printf("Start of pass %d\n",iPass);
3909    // Solve master - may be infeasible
3910    //master.scaling(0);
3911    if (0) {
3912      master.writeMps("yy.mps");
3913    }
3914    // Correct artificials
3915    double sumArtificials=0.0;
3916    if (iPass) {
3917      double * upper = master.columnUpper();
3918      double * solution = master.primalColumnSolution();
3919      double * obj = master.objective();
3920      sumArtificials=0.0;
3921      for (int i=firstArtificial;i<lastArtificial;i++) {
3922        sumArtificials += solution[i];
3923        //assert (solution[i]>-1.0e-2);
3924        if (solution[i]<1.0e-6) {
3925#if 0
3926          // Could take out
3927          obj[i]=0.0;
3928          upper[i]=0.0;
3929#else
3930          obj[i]=1.0e7;
3931          upper[i]=1.0e-1;
3932#endif
3933          solution[i]=0.0;
3934          master.setColumnStatus(i,isFixed);
3935        } else {
3936          upper[i]=solution[i]+1.0e-5*(1.0+solution[i]);
3937        }
3938      }
3939      printf("Sum of artificials before solve is %g\n",sumArtificials);
3940    }
3941    // scale objective to be reasonable
3942    double scaleFactor = master.scaleObjective(-1.0e9);
3943    {
3944      double * dual = master.dualRowSolution();
3945      int n=master.numberRows();
3946      memset(dual,0,n*sizeof(double));
3947      double * solution = master.primalColumnSolution();
3948      master.clpMatrix()->times(1.0,solution,dual);
3949      double sum=0.0;
3950      double * lower = master.rowLower();
3951      double * upper = master.rowUpper();
3952      for (int iRow=0;iRow<n;iRow++) {
3953        double value = dual[iRow];
3954        if (value>upper[iRow])
3955          sum += value-upper[iRow];
3956        else if (value<lower[iRow])
3957          sum -= value-lower[iRow];
3958      }
3959      printf("suminf %g\n",sum);
3960      lower = master.columnLower();
3961      upper = master.columnUpper();
3962      n=master.numberColumns();
3963      for (int iColumn=0;iColumn<n;iColumn++) {
3964        double value = solution[iColumn];
3965        if (value>upper[iColumn]+1.0e-5)
3966          sum += value-upper[iColumn];
3967        else if (value<lower[iColumn]-1.0e-5)
3968          sum -= value-lower[iColumn];
3969      }
3970      printf("suminf %g\n",sum);
3971    }
3972    master.primal(1);
3973    // Correct artificials
3974    sumArtificials=0.0;
3975    {
3976      double * solution = master.primalColumnSolution();
3977      for (int i=firstArtificial;i<lastArtificial;i++) {
3978        sumArtificials += solution[i];
3979      }
3980      printf("Sum of artificials after solve is %g\n",sumArtificials);
3981    }
3982    master.scaleObjective(scaleFactor);
3983    int problemStatus = master.status(); // do here as can change (delcols)
3984    if (master.numberIterations()==0&&iPass)
3985      break; // finished
3986    if (master.objectiveValue()>lastObjective-1.0e-7&&iPass>555)
3987      break; // finished
3988    lastObjective = master.objectiveValue();
3989    // mark basic ones and delete if necessary
3990    int iColumn;
3991    numberColumnsGenerated=master.numberColumns()-numberMasterColumns;
3992    for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
3993      if (master.getStatus(iColumn+numberMasterColumns)==ClpSimplex::basic)
3994        when[iColumn]=iPass;
3995    }
3996    if (numberColumnsGenerated+numberBlocks>maximumColumns) {
3997      // delete
3998      int numberKeep=0;
3999      int numberDelete=0;
4000      int * whichDelete = new int[numberColumnsGenerated];
4001      for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
4002        if (when[iColumn]>iPass-7) {
4003          // keep
4004          when[numberKeep]=when[iColumn];
4005          whichBlock[numberKeep++]=whichBlock[iColumn];
4006        } else {
4007          // delete
4008          whichDelete[numberDelete++]=iColumn+numberMasterColumns;
4009        }
4010      }
4011      numberColumnsGenerated -= numberDelete;
4012      master.deleteColumns(numberDelete,whichDelete);
4013      delete [] whichDelete;
4014    }
4015    const double * dual=NULL;
4016    bool deleteDual=false;
4017    if (problemStatus==0) {
4018      dual = master.dualRowSolution();
4019    } else if (problemStatus==1) {
4020      // could do composite objective
4021      dual = master.infeasibilityRay();
4022      deleteDual = true;
4023      printf("The sum of infeasibilities is %g\n",
4024             master.sumPrimalInfeasibilities());
4025    } else if (!master.numberColumns()) {
4026      assert(!iPass);
4027      dual = master.dualRowSolution();
4028      memset(master.dualRowSolution(),
4029             0,(numberMasterRows+numberBlocks)*sizeof(double));
4030    } else {
4031      abort();
4032    }
4033    // Scale back on first time
4034    if (!iPass) {
4035      double * dual2 = master.dualRowSolution();
4036      for (int i=0;i<numberMasterRows+numberBlocks;i++) {
4037        dual2[i] *= 1.0e-7;
4038      }
4039      dual = master.dualRowSolution();
4040    }
4041    // Create objective for sub problems and solve
4042    columnAdd[0]=0;
4043    int numberProposals=0;
4044    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4045      int numberColumns2 = sub[iBlock].numberColumns();
4046      double * saveObj = new double [numberColumns2];
4047      double * objective2 = sub[iBlock].objective();
4048      memcpy(saveObj,objective2,numberColumns2*sizeof(double));
4049      // new objective
4050      top[iBlock]->transposeTimes(dual,objective2);
4051      int i;
4052      if (problemStatus==0) {
4053        for (i=0;i<numberColumns2;i++)
4054          objective2[i] = saveObj[i]-objective2[i];
4055      } else {
4056        for (i=0;i<numberColumns2;i++)
4057          objective2[i] = -objective2[i];
4058      }
4059      // scale objective to be reasonable
4060      double scaleFactor = 
4061        sub[iBlock].scaleObjective((sumArtificials>1.0e-5) ? -1.0e-4 :-1.0e9);
4062      if (iPass) {
4063        sub[iBlock].primal();
4064      } else {
4065        sub[iBlock].dual();
4066      }
4067      sub[iBlock].scaleObjective(scaleFactor);
4068      if (!sub[iBlock].isProvenOptimal()&& 
4069          !sub[iBlock].isProvenDualInfeasible()) {
4070        memset(objective2,0,numberColumns2*sizeof(double));
4071        sub[iBlock].primal();
4072        if (problemStatus==0) {
4073          for (i=0;i<numberColumns2;i++)
4074            objective2[i] = saveObj[i]-objective2[i];
4075        } else {
4076          for (i=0;i<numberColumns2;i++)
4077            objective2[i] = -objective2[i];
4078        }
4079        double scaleFactor = sub[iBlock].scaleObjective(-1.0e9);
4080        sub[iBlock].primal(1);
4081        sub[iBlock].scaleObjective(scaleFactor);
4082      }
4083      memcpy(objective2,saveObj,numberColumns2*sizeof(double));
4084      // get proposal
4085      if (sub[iBlock].numberIterations()||!iPass) {
4086        double objValue =0.0;
4087        int start = columnAdd[numberProposals];
4088        // proposal
4089        if (sub[iBlock].isProvenOptimal()) {
4090          const double * solution = sub[iBlock].primalColumnSolution();
4091          top[iBlock]->times(solution,elementAdd+start);
4092          for (i=0;i<numberColumns2;i++)
4093            objValue += solution[i]*saveObj[i];
4094          // See if good dj and pack down
4095          int number=start;
4096          double dj = objValue;
4097          if (problemStatus) 
4098            dj=0.0;
4099          double smallest=1.0e100;
4100          double largest=0.0;
4101          for (i=0;i<numberMasterRows;i++) {
4102            double value = elementAdd[start+i];
4103            if (fabs(value)>1.0e-15) {
4104              dj -= dual[i]*value;
4105              smallest = CoinMin(smallest,fabs(value));
4106              largest = CoinMax(largest,fabs(value));
4107              rowAdd[number]=i;
4108              elementAdd[number++]=value;
4109            }
4110          }
4111          // and convexity
4112          dj -= dual[numberMasterRows+iBlock];
4113          rowAdd[number]=numberMasterRows+iBlock;
4114          elementAdd[number++]=1.0;
4115          // if elements large then scale?
4116          //if (largest>1.0e8||smallest<1.0e-8)
4117          printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
4118                 iBlock,smallest,largest,dj);
4119          if (dj<-1.0e-6||!iPass) {
4120            // take
4121            objective[numberProposals]=objValue;
4122            columnAdd[++numberProposals]=number;
4123            when[numberColumnsGenerated]=iPass;
4124            whichBlock[numberColumnsGenerated++]=iBlock;
4125          }
4126        } else if (sub[iBlock].isProvenDualInfeasible()) {
4127          // use ray
4128          const double * solution = sub[iBlock].unboundedRay();
4129          top[iBlock]->times(solution,elementAdd+start);
4130          for (i=0;i<numberColumns2;i++)
4131            objValue += solution[i]*saveObj[i];
4132          // See if good dj and pack down
4133          int number=start;
4134          double dj = objValue;
4135          double smallest=1.0e100;
4136          double largest=0.0;
4137          for (i=0;i<numberMasterRows;i++) {
4138            double value = elementAdd[start+i];
4139            if (fabs(value)>1.0e-15) {
4140              dj -= dual[i]*value;
4141              smallest = CoinMin(smallest,fabs(value));
4142              largest = CoinMax(largest,fabs(value));
4143              rowAdd[number]=i;
4144              elementAdd[number++]=value;
4145            }
4146          }
4147          // if elements large or small then scale?
4148          //if (largest>1.0e8||smallest<1.0e-8)
4149          printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
4150                 iBlock,smallest,largest,dj);
4151          if (dj<-1.0e-6) {
4152            // take
4153            objective[numberProposals]=objValue;
4154            columnAdd[++numberProposals]=number;
4155            when[numberColumnsGenerated]=iPass;
4156            whichBlock[numberColumnsGenerated++]=iBlock;
4157          }
4158        } else {
4159          abort();
4160        }
4161      }
4162      delete [] saveObj;
4163    }
4164    if (deleteDual)
4165      delete [] dual;
4166    if (numberProposals) 
4167      master.addColumns(numberProposals,NULL,NULL,objective,
4168                        columnAdd,rowAdd,elementAdd);
4169  }
4170  std::cout<<"Time at end of D-W "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4171  //master.scaling(0);
4172  //master.primal(1);
4173  loadProblem(*model);
4174  // now put back a good solution
4175  double * lower = new double[numberMasterRows+numberBlocks];
4176  double * upper = new double[numberMasterRows+numberBlocks];
4177  numberColumnsGenerated  += numberMasterColumns;
4178  double * sol = new double[numberColumnsGenerated];
4179  const double * solution = master.primalColumnSolution();
4180  const double * masterLower = master.rowLower();
4181  const double * masterUpper = master.rowUpper();
4182  double * fullSolution = primalColumnSolution();
4183  const double * fullLower = columnLower();
4184  const double * fullUpper = columnUpper();
4185  const double * rowSolution = master.primalRowSolution();
4186  double * fullRowSolution = primalRowSolution();
4187  const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4188  int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4189  const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4190  int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4191  for (int iRow=0;iRow<numberRows2;iRow++) {
4192    int kRow = rowBack[iRow];
4193    setRowStatus(kRow,master.getRowStatus(iRow));
4194    fullRowSolution[kRow]=rowSolution[iRow];
4195  }
4196  for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4197    int kColumn = columnBack[iColumn];
4198    setStatus(kColumn,master.getStatus(iColumn));
4199    fullSolution[kColumn]=solution[iColumn];
4200  }
4201  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4202    // move basis
4203    int kBlock = columnCounts[iBlock];
4204    const int * rowBack = model->coinBlock(kBlock)->originalRows();
4205    int numberRows2 = model->coinBlock(kBlock)->numberRows();
4206    const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4207    int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4208    for (int iRow=0;iRow<numberRows2;iRow++) {
4209      int kRow = rowBack[iRow];
4210      setRowStatus(kRow,sub[iBlock].getRowStatus(iRow));
4211    }
4212    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4213      int kColumn = columnBack[iColumn];
4214      setStatus(kColumn,sub[iBlock].getStatus(iColumn));
4215    }
4216    // convert top bit to by rows
4217    CoinPackedMatrix topMatrix = *top[iBlock];
4218    topMatrix.reverseOrdering();
4219    // zero solution
4220    memset(sol,0,numberColumnsGenerated*sizeof(double));
4221   
4222    for (int i=numberMasterColumns;i<numberColumnsGenerated;i++) {
4223      if (whichBlock[i-numberMasterColumns]==iBlock)
4224        sol[i] = solution[i];
4225    }
4226    memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double));
4227    master.clpMatrix()->times(1.0,sol,lower);
4228    for (int iRow=0;iRow<numberMasterRows;iRow++) {
4229      double value = lower[iRow];
4230      if (masterUpper[iRow]<1.0e20) 
4231        upper[iRow] = value;
4232      else
4233        upper[iRow]=COIN_DBL_MAX;
4234      if (masterLower[iRow]>-1.0e20) 
4235        lower[iRow] = value;
4236      else
4237        lower[iRow]=-COIN_DBL_MAX;
4238    }
4239    sub[iBlock].addRows(numberMasterRows,lower,upper,
4240                        topMatrix.getVectorStarts(),
4241                        topMatrix.getVectorLengths(),
4242                        topMatrix.getIndices(),
4243                        topMatrix.getElements());
4244    sub[iBlock].primal(1);
4245    const double * subSolution = sub[iBlock].primalColumnSolution();
4246    const double * subRowSolution = sub[iBlock].primalRowSolution();
4247    // move solution
4248    for (int iRow=0;iRow<numberRows2;iRow++) {
4249      int kRow = rowBack[iRow];
4250      fullRowSolution[kRow]=subRowSolution[iRow];
4251    }
4252    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4253      int kColumn = columnBack[iColumn];
4254      fullSolution[kColumn]=subSolution[iColumn];
4255    }
4256  }
4257  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
4258    if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&&
4259        fullSolution[iColumn]>fullLower[iColumn]+1.0e-8) {
4260      if (getStatus(iColumn)!=ClpSimplex::basic) {
4261        if (columnLower_[iColumn]>-1.0e30||
4262            columnUpper_[iColumn]<1.0e30)
4263          setStatus(iColumn,ClpSimplex::superBasic);
4264        else
4265          setStatus(iColumn,ClpSimplex::isFree);
4266      }
4267    } else if (fullSolution[iColumn]>=fullUpper[iColumn]-1.0e-8) {
4268      // may help to make rest non basic
4269      if (getStatus(iColumn)!=ClpSimplex::basic) 
4270        setStatus(iColumn,ClpSimplex::atUpperBound);
4271    } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e-8) {
4272      // may help to make rest non basic
4273      if (getStatus(iColumn)!=ClpSimplex::basic) 
4274        setStatus(iColumn,ClpSimplex::atLowerBound);
4275    }
4276  }
4277  //int numberRows=model->numberRows();
4278  //for (int iRow=0;iRow<numberRows;iRow++)
4279  //setRowStatus(iRow,ClpSimplex::superBasic);
4280  std::cout<<"Time before cleanup of full model "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4281  primal(1);
4282  std::cout<<"Total time "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4283  delete [] columnCounts;
4284  delete [] sol;
4285  delete [] lower;
4286  delete [] upper;
4287  delete [] whichBlock;
4288  delete [] when;
4289  delete [] columnAdd;
4290  delete [] rowAdd;
4291  delete [] elementAdd;
4292  delete [] objective;
4293  delete [] top;
4294  delete [] sub;
4295  return 0;
4296}
4297// Solve using Benders decomposition and maybe in parallel
4298int 
4299ClpSimplex::solveBenders(CoinStructuredModel * model)
4300{
4301  double time1 = CoinCpuTime();
4302  //int numberColumns = model->numberColumns();
4303  int numberRowBlocks=model->numberRowBlocks();
4304  int numberColumnBlocks = model->numberColumnBlocks();
4305  int numberElementBlocks = model->numberElementBlocks();
4306  // We already have top level structure
4307  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4308  for (int i=0;i<numberElementBlocks;i++) {
4309    CoinModel * thisBlock = model->coinBlock(i);
4310    assert (thisBlock);
4311    // just fill in info
4312    CoinModelBlockInfo info = CoinModelBlockInfo();
4313    int whatsSet = thisBlock->whatIsSet();
4314    info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
4315    info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
4316    info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
4317    info.integer = ((whatsSet&32)!=0) ? 1 : 0;
4318    info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
4319    info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
4320    // Which block
4321    int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
4322    info.rowBlock=iRowBlock;
4323    int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
4324    info.columnBlock=iColumnBlock;
4325    blockInfo[i] = info;
4326  }
4327  // make up problems
4328  int numberBlocks=numberColumnBlocks-1;
4329  // Find master columns and rows
4330  int * columnCounts = new int [numberColumnBlocks];
4331  CoinZeroN(columnCounts,numberColumnBlocks);
4332  int * rowCounts = new int [numberRowBlocks+1];
4333  CoinZeroN(rowCounts,numberRowBlocks);
4334  int iBlock;
4335  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
4336    int iColumnBlock = blockInfo[iBlock].columnBlock;
4337    columnCounts[iColumnBlock]++;
4338    int iRowBlock =blockInfo[iBlock].rowBlock;
4339    rowCounts[iRowBlock]++;
4340  }
4341  int * whichBlock = new int [numberElementBlocks];
4342  int masterColumnBlock=-1;
4343  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
4344    if (columnCounts[iBlock]>1) {
4345      if (masterColumnBlock==-1) {
4346        masterColumnBlock=iBlock;
4347      } else {
4348        // Can't decode
4349        masterColumnBlock=-2;
4350        break;
4351      }
4352    }
4353  }
4354  int masterRowBlock=-1;
4355  int kBlock=0;
4356  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
4357    int count=rowCounts[iBlock];
4358    rowCounts[iBlock]=kBlock;
4359    kBlock += count;
4360  }
4361  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
4362    int iRowBlock = blockInfo[iBlock].rowBlock;
4363    whichBlock[rowCounts[iRowBlock]]=iBlock;
4364    rowCounts[iRowBlock]++;
4365  }
4366  for (iBlock = numberRowBlocks-1;iBlock>=0;iBlock--) 
4367    rowCounts[iBlock+1]=rowCounts[iBlock];
4368  rowCounts[0]=0;
4369  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
4370    int count=rowCounts[iBlock+1]-rowCounts[iBlock];
4371    if (count==1) {
4372      int kBlock = whichBlock[rowCounts[iBlock]];
4373      int iColumnBlock = blockInfo[kBlock].columnBlock;
4374      if (iColumnBlock==masterColumnBlock) {
4375        if (masterRowBlock==-1) {
4376          masterRowBlock=iBlock;
4377        } else {
4378          // Can't decode
4379          masterRowBlock=-2;
4380          break;
4381        }
4382      }
4383    }
4384  }
4385  if (masterColumnBlock<0||masterRowBlock==-2) {
4386    // What now
4387    abort();
4388  }
4389  delete [] columnCounts;
4390  // create all data
4391  const CoinPackedMatrix ** first = new const CoinPackedMatrix * [numberRowBlocks];
4392  ClpSimplex * sub = new ClpSimplex [numberBlocks];
4393  ClpSimplex master;
4394  // Set offset
4395  master.setObjectiveOffset(model->objectiveOffset());
4396  kBlock=0;
4397  int masterBlock=-1;
4398  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
4399    first[kBlock]=NULL;
4400    int start=rowCounts[iBlock];
4401    int end=rowCounts[iBlock+1];
4402    assert (end-start<=2);
4403    for (int j=start;j<end;j++) {
4404      int jBlock = whichBlock[j];
4405      int iColumnBlock = blockInfo[jBlock].columnBlock;
4406      int iRowBlock =blockInfo[jBlock].rowBlock;
4407      assert (iRowBlock==iBlock);
4408      if (iRowBlock!=masterRowBlock&&iColumnBlock==masterColumnBlock) {
4409        // first matrix
4410        first[kBlock]=model->coinBlock(jBlock)->packedMatrix();
4411      } else {
4412        const CoinPackedMatrix * matrix
4413          = model->coinBlock(jBlock)->packedMatrix();
4414        // Get pointers to arrays
4415        const double * columnLower;
4416        const double * columnUpper;
4417        const double * rowLower;
4418        const double * rowUpper;
4419        const double * objective;
4420        model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
4421                     columnLower,columnUpper,objective);
4422        if (iRowBlock!=masterRowBlock) {
4423          // diagonal block
4424          sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
4425                                  objective,rowLower,rowUpper);
4426          if (optimizationDirection_<0.0) {
4427            double * obj = sub[kBlock].objective();
4428            int n = sub[kBlock].numberColumns();
4429            for (int i=0;i<n;i++)
4430              obj[i] = - obj[i];
4431          }
4432          if (this->factorizationFrequency()==200) {
4433            // User did not touch preset
4434            sub[kBlock].defaultFactorizationFrequency();
4435          } else {
4436            // make sure model has correct value
4437            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
4438          }
4439          sub[kBlock].setPerturbation(50);
4440          // Set rowCounts to be diagonal block index for cleanup
4441          rowCounts[kBlock]=jBlock;
4442        } else {
4443          // master
4444          masterBlock=jBlock;
4445          master.loadProblem(*matrix,columnLower,columnUpper,
4446                             objective,rowLower,rowUpper);
4447          if (optimizationDirection_<0.0) {
4448            double * obj = master.objective();
4449            int n = master.numberColumns();
4450            for (int i=0;i<n;i++)
4451              obj[i] = - obj[i];
4452          }
4453        }
4454      }
4455    }
4456    if (iBlock!=masterRowBlock) 
4457      kBlock++;
4458  }
4459  delete [] whichBlock;
4460  delete [] blockInfo;
4461  // For now master must have been defined (does not have to have rows)
4462  assert (master.numberColumns());
4463  assert (masterBlock>=0);
4464  int numberMasterColumns = master.numberColumns();
4465  // Overkill in terms of space
4466  int spaceNeeded = CoinMax(numberBlocks*(numberMasterColumns+1),
4467                            2*numberMasterColumns);
4468  int * columnAdd = new int[spaceNeeded];
4469  double * elementAdd = new double[spaceNeeded];
4470  spaceNeeded = numberBlocks;
4471  int * rowAdd = new int[spaceNeeded+1];
4472  double * objective = new double[spaceNeeded];
4473  int maxPass=500;
4474  int iPass;
4475  double lastObjective=-1.0e31;
4476  // Create columns for proposals
4477  int numberMasterRows = master.numberRows();
4478  master.resize(numberMasterColumns+numberBlocks,numberMasterRows);
4479  if (this->factorizationFrequency()==200) {
4480    // User did not touch preset
4481    master.defaultFactorizationFrequency();
4482  } else {
4483    // make sure model has correct value
4484    master.setFactorizationFrequency(this->factorizationFrequency());
4485  }
4486  master.setPerturbation(50);
4487  // Arrays to say which block and when created
4488  int maximumRows = 2*numberMasterColumns+10*numberBlocks;
4489  whichBlock = new int[maximumRows];
4490  int * when = new int[maximumRows];
4491  int numberRowsGenerated=numberBlocks;
4492  // Add extra variables
4493  {
4494    int iBlock;
4495    columnAdd[0] = 0;
4496    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4497      objective[iBlock] = 1.0;
4498      columnAdd[iBlock+1] = 0;
4499      when[iBlock]=-1;
4500      whichBlock[iBlock] = iBlock;
4501    }
4502    master.addColumns(numberBlocks,NULL,NULL,objective,
4503                      columnAdd, rowAdd, elementAdd);
4504  }
4505  std::cout<<"Time to decompose "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4506  for (iPass=0;iPass<maxPass;iPass++) {
4507    printf("Start of pass %d\n",iPass);
4508    // Solve master - may be unbounded
4509    //master.scaling(0);
4510    if (1) {
4511      master.writeMps("yy.mps");
4512    }
4513    master.dual();
4514    int problemStatus = master.status(); // do here as can change (delcols)
4515    if (master.numberIterations()==0&&iPass)
4516      break; // finished
4517    if (master.objectiveValue()<lastObjective+1.0e-7&&iPass>555)
4518      break; // finished
4519    lastObjective = master.objectiveValue();
4520    // mark non-basic rows and delete if necessary
4521    int iRow;
4522    numberRowsGenerated=master.numberRows()-numberMasterRows;
4523    for (iRow=0;iRow<numberRowsGenerated;iRow++) {
4524      if (master.getStatus(iRow+numberMasterRows)!=ClpSimplex::basic)
4525        when[iRow]=iPass;
4526    }
4527    if (numberRowsGenerated>maximumRows) {
4528      // delete
4529      int numberKeep=0;
4530      int numberDelete=0;
4531      int * whichDelete = new int[numberRowsGenerated];
4532      for (iRow=0;iRow<numberRowsGenerated;iRow++) {
4533        if (when[iRow]>iPass-7) {
4534          // keep
4535          when[numberKeep]=when[iRow];
4536          whichBlock[numberKeep++]=whichBlock[iRow];
4537        } else {
4538          // delete
4539          whichDelete[numberDelete++]=iRow+numberMasterRows;
4540        }
4541      }
4542      numberRowsGenerated -= numberDelete;
4543      master.deleteRows(numberDelete,whichDelete);
4544      delete [] whichDelete;
4545    }
4546    const double * primal=NULL;
4547    bool deletePrimal=false;
4548    if (problemStatus==0) {
4549      primal = master.primalColumnSolution();
4550    } else if (problemStatus==2) {
4551      // could do composite objective
4552      primal = master.unboundedRay();
4553      deletePrimal = true;
4554      printf("The sum of infeasibilities is %g\n",
4555             master.sumPrimalInfeasibilities());
4556    } else if (!master.numberRows()) {
4557      assert(!iPass);
4558      primal = master.primalColumnSolution();
4559      memset(master.primalColumnSolution(),
4560             0,numberMasterColumns*sizeof(double));
4561    } else {
4562      abort();
4563    }
4564    // Create rhs for sub problems and solve
4565    rowAdd[0]=0;
4566    int numberProposals=0;
4567    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4568      int numberRows2 = sub[iBlock].numberRows();
4569      double * saveLower = new double [numberRows2];
4570      double * lower2 = sub[iBlock].rowLower();
4571      double * saveUpper = new double [numberRows2];
4572      double * upper2 = sub[iBlock].rowUpper();
4573      // new rhs
4574      CoinZeroN(saveUpper,numberRows2);
4575      first[iBlock]->times(primal,saveUpper);
4576      for (int i=0;i<numberRows2;i++) {
4577        double value = saveUpper[i];
4578        saveLower[i]=lower2[i];
4579        saveUpper[i]=upper2[i];
4580        if (saveLower[i]>-1.0e30)
4581          lower2[i] -= value;
4582        if (saveUpper[i]<1.0e30)
4583          upper2[i] -= value;
4584      }
4585      sub[iBlock].dual();
4586      memcpy(lower2,saveLower,numberRows2*sizeof(double));
4587      memcpy(upper2,saveUpper,numberRows2*sizeof(double));
4588      // get proposal
4589      if (sub[iBlock].numberIterations()||!iPass) {
4590        double objValue =0.0;
4591        int start = rowAdd[numberProposals];
4592        // proposal
4593        if (sub[iBlock].isProvenOptimal()) {
4594          const double * solution = sub[iBlock].dualRowSolution();
4595          first[iBlock]->transposeTimes(solution,elementAdd+start);
4596          for (int i=0;i<numberRows2;i++) {
4597            if (solution[i]<-dualTolerance_) {
4598              // at upper
4599              assert (saveUpper[i]<1.0e30);
4600              objValue += solution[i]*saveUpper[i];
4601            } else if (solution[i]>dualTolerance_) {
4602              // at lower
4603              assert (saveLower[i]>-1.0e30);
4604              objValue += solution[i]*saveLower[i];
4605            }
4606          }
4607         
4608          // See if cuts off and pack down
4609          int number=start;
4610          double infeas = objValue;
4611          double smallest=1.0e100;
4612          double largest=0.0;
4613          for (int i=0;i<numberMasterColumns;i++) {
4614            double value = elementAdd[start+i];
4615            if (fabs(value)>1.0e-15) {
4616              infeas -= primal[i]*value;
4617              smallest = CoinMin(smallest,fabs(value));
4618              largest = CoinMax(largest,fabs(value));
4619              columnAdd[number]=i;
4620              elementAdd[number++]=-value;
4621            }
4622          }
4623          columnAdd[number]=numberMasterColumns+iBlock;
4624          elementAdd[number++]=-1.0;
4625          // if elements large then scale?
4626          //if (largest>1.0e8||smallest<1.0e-8)
4627          printf("For subproblem %d smallest - %g, largest %g - infeas %g\n",
4628                 iBlock,smallest,largest,infeas);
4629          if (infeas<-1.0e-6||!iPass) {
4630            // take
4631            objective[numberProposals]=objValue;
4632            rowAdd[++numberProposals]=number;
4633            when[numberRowsGenerated]=iPass;
4634            whichBlock[numberRowsGenerated++]=iBlock;
4635          }
4636        } else if (sub[iBlock].isProvenPrimalInfeasible()) {
4637          // use ray
4638          const double * solution = sub[iBlock].infeasibilityRay();
4639          first[iBlock]->transposeTimes(solution,elementAdd+start);
4640          for (int i=0;i<numberRows2;i++) {
4641            if (solution[i]<-dualTolerance_) {
4642              // at upper
4643              assert (saveUpper[i]<1.0e30);
4644              objValue += solution[i]*saveUpper[i];
4645            } else if (solution[i]>dualTolerance_) {
4646              // at lower
4647              assert (saveLower[i]>-1.0e30);
4648              objValue += solution[i]*saveLower[i];
4649            }
4650          }
4651          // See if good infeas and pack down
4652          int number=start;
4653          double infeas = objValue;
4654          double smallest=1.0e100;
4655          double largest=0.0;
4656          for (int i=0;i<numberMasterColumns;i++) {
4657            double value = elementAdd[start+i];
4658            if (fabs(value)>1.0e-15) {
4659              infeas -= primal[i]*value;
4660              smallest = CoinMin(smallest,fabs(value));
4661              largest = CoinMax(largest,fabs(value));
4662              columnAdd[number]=i;
4663              elementAdd[number++]=-value;
4664            }
4665          }
4666          // if elements large or small then scale?
4667          //if (largest>1.0e8||smallest<1.0e-8)
4668          printf("For subproblem ray %d smallest - %g, largest %g - infeas %g\n",
4669                 iBlock,smallest,largest,infeas);
4670          if (infeas<-1.0e-6) {
4671            // take
4672            objective[numberProposals]=objValue;
4673            rowAdd[++numberProposals]=number;
4674            when[numberRowsGenerated]=iPass;
4675            whichBlock[numberRowsGenerated++]=iBlock;
4676          }
4677        } else {
4678          abort();
4679        }
4680      }
4681      delete [] saveLower;
4682      delete [] saveUpper;
4683    }
4684    if (deletePrimal)
4685      delete [] primal;
4686    if (numberProposals) {
4687      master.addRows(numberProposals,NULL,objective,
4688                        rowAdd,columnAdd,elementAdd);
4689    }
4690  }
4691  std::cout<<"Time at end of Benders "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4692  //master.scaling(0);
4693  //master.primal(1);
4694  loadProblem(*model);
4695  // now put back a good solution
4696  const double * columnSolution = master.primalColumnSolution();
4697  double * fullColumnSolution = primalColumnSolution();
4698  const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4699  int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4700  const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4701  int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4702  for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4703    int kColumn = columnBack[iColumn];
4704    setColumnStatus(kColumn,master.getColumnStatus(iColumn));
4705    fullColumnSolution[kColumn]=columnSolution[iColumn];
4706  }
4707  for (int iRow=0;iRow<numberRows2;iRow++) {
4708    int kRow = rowBack[iRow];
4709    setStatus(kRow,master.getStatus(iRow));
4710    //fullSolution[kRow]=solution[iRow];
4711  }
4712  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4713    // move basis
4714    int kBlock = rowCounts[iBlock];
4715    const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4716    int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4717    const int * rowBack = model->coinBlock(kBlock)->originalRows();
4718    int numberRows2 = model->coinBlock(kBlock)->numberRows();
4719    const double * subColumnSolution = sub[iBlock].primalColumnSolution();
4720    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4721      int kColumn = columnBack[iColumn];
4722      setColumnStatus(kColumn,sub[iBlock].getColumnStatus(iColumn));
4723      fullColumnSolution[kColumn]=subColumnSolution[iColumn];
4724    }
4725    for (int iRow=0;iRow<numberRows2;iRow++) {
4726      int kRow = rowBack[iRow];
4727      setStatus(kRow,sub[iBlock].getStatus(iRow));
4728      setStatus(kRow,atLowerBound);
4729    }
4730  }
4731  double * fullSolution = primalRowSolution();
4732  CoinZeroN(fullSolution,numberRows_);
4733  times(1.0,fullColumnSolution,fullSolution);
4734  //int numberColumns=model->numberColumns();
4735  //for (int iColumn=0;iColumn<numberColumns;iColumn++)
4736  //setColumnStatus(iColumn,ClpSimplex::superBasic);
4737  std::cout<<"Time before cleanup of full model "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4738  this->primal(1);
4739  std::cout<<"Total time "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4740  delete [] rowCounts;
4741  //delete [] sol;
4742  //delete [] lower;
4743  //delete [] upper;
4744  delete [] whichBlock;
4745  delete [] when;
4746  delete [] rowAdd;
4747  delete [] columnAdd;
4748  delete [] elementAdd;
4749  delete [] objective;
4750  delete [] first;
4751  delete [] sub;
4752  return 0;
4753}
Note: See TracBrowser for help on using the repository browser.