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

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

changes for idiot code

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