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

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

changes to simplex and lots of stuff and start Mumps cholesky

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