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

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

get rid of compiler warnings

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