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

Last change on this file since 1453 was 1453, checked in by forrest, 10 years ago

second try at changes for seg fault and idiot

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