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

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

fix memory leak

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