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

Last change on this file since 1326 was 1326, checked in by forrest, 12 years ago

First attempt at Benders decomposition

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