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

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

try and improve quadratic and barrier

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 146.3 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    bool forceFixing=false;
1878    if (barrierOptions&16) {
1879      barrierOptions &= ~16;
1880      doKKT=true;
1881    }
1882    if (barrierOptions&(32+64+128)) {
1883      aggressiveGamma=(barrierOptions&(32+64+128))>>5;
1884      barrierOptions &= ~(32+64+128);
1885    }
1886    if (barrierOptions&256) {
1887      barrierOptions &= ~256;
1888      presolveInCrossover=true;
1889    }
1890    if (barrierOptions&512) {
1891      barrierOptions &= ~512;
1892      forceFixing=true;
1893    }
1894    if (barrierOptions&8) {
1895      barrierOptions &= ~8;
1896      scale=true;
1897    }
1898#ifdef COIN_DEVELOP
1899#ifndef FAST_BARRIER
1900    if (!numberBarrier)
1901      std::cout<<"Warning - the default ordering is just on row counts! "
1902               <<"The factorization is being improved"<<std::endl;
1903    numberBarrier++;
1904#endif
1905#endif
1906    // If quadratic force KKT
1907    if (quadraticObj) {
1908      doKKT=true;
1909    }
1910    switch (barrierOptions) {
1911    case 0:
1912    default:
1913      if (!doKKT) {
1914        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1915        barrier.setCholesky(cholesky);
1916      } else {
1917        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
1918        cholesky->setKKT(true);
1919        barrier.setCholesky(cholesky);
1920      }
1921      break;
1922    case 1:
1923      if (!doKKT) {
1924        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1925        barrier.setCholesky(cholesky);
1926      } else {
1927        ClpCholeskyDense * cholesky = new ClpCholeskyDense();
1928        cholesky->setKKT(true);
1929        barrier.setCholesky(cholesky);
1930      }
1931      break;
1932#ifdef WSSMP_BARRIER
1933    case 2:
1934      {
1935        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
1936        barrier.setCholesky(cholesky);
1937        assert (!doKKT);
1938      }
1939      break;
1940    case 3:
1941      if (!doKKT) {
1942        ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
1943        barrier.setCholesky(cholesky);
1944      } else {
1945        ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
1946        barrier.setCholesky(cholesky);
1947      }
1948      break;
1949#endif
1950#ifdef UFL_BARRIER
1951    case 4:
1952      if (!doKKT) {
1953        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1954        barrier.setCholesky(cholesky);
1955      } else {
1956        ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
1957        cholesky->setKKT(true);
1958        barrier.setCholesky(cholesky);
1959      }
1960      break;
1961#endif
1962#ifdef TAUCS_BARRIER
1963    case 5:
1964      {
1965        ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
1966        barrier.setCholesky(cholesky);
1967        assert (!doKKT);
1968      }
1969      break;
1970#endif
1971#ifdef MUMPS_BARRIER
1972    case 6:
1973      {
1974        ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
1975        barrier.setCholesky(cholesky);
1976        assert (!doKKT);
1977      }
1978      break;
1979#endif
1980    }
1981    int numberRows = model2->numberRows();
1982    int numberColumns = model2->numberColumns();
1983    int saveMaxIts = model2->maximumIterations();
1984    if (saveMaxIts<1000) {
1985      barrier.setMaximumBarrierIterations(saveMaxIts);
1986      model2->setMaximumIterations(1000000);
1987    }
1988#ifndef SAVEIT
1989    //barrier.setDiagonalPerturbation(1.0e-25);
1990    if (aggressiveGamma) {
1991      switch (aggressiveGamma) {
1992      case 1:
1993        barrier.setGamma(1.0e-5);
1994        barrier.setDelta(1.0e-5);
1995        break;
1996      case 2:
1997        barrier.setGamma(1.0e-5);
1998        break;
1999      case 3:
2000        barrier.setDelta(1.0e-5);
2001        break;
2002      case 4:
2003        barrier.setGamma(1.0e-3);
2004        barrier.setDelta(1.0e-3);
2005        break;
2006      case 5:
2007        barrier.setGamma(1.0e-3);
2008        break;
2009      case 6:
2010        barrier.setDelta(1.0e-3);
2011        break;
2012      }
2013    }
2014    if (scale)
2015      barrier.scaling(1);
2016    else
2017      barrier.scaling(0);
2018    barrier.primalDual();
2019#elif SAVEIT==1
2020    barrier.primalDual();
2021#else
2022    model2->restoreModel("xx.save");
2023    // move solutions
2024    CoinMemcpyN(model2->primalRowSolution(),
2025                numberRows,barrier.primalRowSolution());
2026    CoinMemcpyN(model2->dualRowSolution(),
2027                numberRows,barrier.dualRowSolution());
2028    CoinMemcpyN(model2->primalColumnSolution(),
2029                numberColumns,barrier.primalColumnSolution());
2030    CoinMemcpyN(model2->dualColumnSolution(),
2031                numberColumns,barrier.dualColumnSolution());
2032#endif
2033    time2 = CoinCpuTime();
2034    timeCore = time2-timeX;
2035    handler_->message(CLP_INTERVAL_TIMING,messages_)
2036      <<"Barrier"<<timeCore<<time2-time1
2037      <<CoinMessageEol;
2038    timeX=time2;
2039    int maxIts = barrier.maximumBarrierIterations();
2040    int barrierStatus=barrier.status();
2041    double gap = barrier.complementarityGap();
2042    // get which variables are fixed
2043    double * saveLower=NULL;
2044    double * saveUpper=NULL;
2045    ClpPresolve pinfo2;
2046    ClpSimplex * saveModel2=NULL;
2047    bool extraPresolve=false;
2048    int numberFixed = barrier.numberFixed();
2049    if (numberFixed) {
2050      int numberRows = barrier.numberRows();
2051      int numberColumns = barrier.numberColumns();
2052      int numberTotal = numberRows+numberColumns;
2053      saveLower = new double [numberTotal];
2054      saveUpper = new double [numberTotal];
2055      CoinMemcpyN(barrier.columnLower(),numberColumns,saveLower);
2056      CoinMemcpyN(barrier.rowLower(),numberRows,saveLower+numberColumns);
2057      CoinMemcpyN(barrier.columnUpper(),numberColumns,saveUpper);
2058      CoinMemcpyN(barrier.rowUpper(),numberRows,saveUpper+numberColumns);
2059    }
2060    if (((numberFixed*20>barrier.numberRows()&&numberFixed>5000)||forceFixing)&&
2061        presolveInCrossover) {
2062      // may as well do presolve
2063      if (!forceFixing) {
2064        barrier.fixFixed();
2065      } else {
2066        // Fix
2067        int n=barrier.numberColumns();
2068        double * lower = barrier.columnLower();
2069        double * upper = barrier.columnUpper();
2070        double * solution = barrier.primalColumnSolution();
2071        int nFix=0;
2072        for (int i=0;i<n;i++) {
2073          if (barrier.fixedOrFree(i)&&lower[i]<upper[i]) {
2074            double value = solution[i];
2075            if (value<lower[i]+1.0e-6&&value-lower[i]<upper[i]-value) {
2076              solution[i]=lower[i];
2077              upper[i]=lower[i];
2078              nFix++;
2079            } else if (value>upper[i]-1.0e-6&&value-lower[i]>upper[i]-value) {
2080              solution[i]=upper[i];
2081              lower[i]=upper[i];
2082              nFix++;
2083            }
2084          }
2085        }
2086        printf("%d columns fixed\n",nFix);
2087        int nr=barrier.numberRows();
2088        lower = barrier.rowLower();
2089        upper = barrier.rowUpper();
2090        solution = barrier.primalRowSolution();
2091        nFix=0;
2092        for (int i=0;i<nr;i++) {
2093          if (barrier.fixedOrFree(i+n)&&lower[i]<upper[i]) {
2094            double value = solution[i];
2095            if (value<lower[i]+1.0e-6&&value-lower[i]<upper[i]-value) {
2096              solution[i]=lower[i];
2097              upper[i]=lower[i];
2098              nFix++;
2099            } else if (value>upper[i]-1.0e-6&&value-lower[i]>upper[i]-value) {
2100              solution[i]=upper[i];
2101              lower[i]=upper[i];
2102              nFix++;
2103            }
2104          }
2105        }
2106        printf("%d row slacks fixed\n",nFix);
2107      }
2108      saveModel2=model2;
2109      extraPresolve=true;
2110    } else if (numberFixed) {
2111      // Set fixed to bounds (may have restored earlier solution)
2112      if (!forceFixing) {
2113        barrier.fixFixed(false);
2114      } else {
2115        // Fix
2116        int n=barrier.numberColumns();
2117        double * lower = barrier.columnLower();
2118        double * upper = barrier.columnUpper();
2119        double * solution = barrier.primalColumnSolution();
2120        int nFix=0;
2121        for (int i=0;i<n;i++) {
2122          if (barrier.fixedOrFree(i)&&lower[i]<upper[i]) {
2123            double value = solution[i];
2124            if (value<lower[i]+1.0e-8&&value-lower[i]<upper[i]-value) {
2125              solution[i]=lower[i];
2126              upper[i]=lower[i];
2127              nFix++;
2128            } else if (value>upper[i]-1.0e-8&&value-lower[i]>upper[i]-value) {
2129              solution[i]=upper[i];
2130              lower[i]=upper[i];
2131              nFix++;
2132            } else {
2133              //printf("fixcol %d %g <= %g <= %g\n",
2134              //     i,lower[i],solution[i],upper[i]);
2135            }
2136          }
2137        }
2138        printf("%d columns fixed\n",nFix);
2139        int nr=barrier.numberRows();
2140        lower = barrier.rowLower();
2141        upper = barrier.rowUpper();
2142        solution = barrier.primalRowSolution();
2143        nFix=0;
2144        for (int i=0;i<nr;i++) {
2145          if (barrier.fixedOrFree(i+n)&&lower[i]<upper[i]) {
2146            double value = solution[i];
2147            if (value<lower[i]+1.0e-5&&value-lower[i]<upper[i]-value) {
2148              solution[i]=lower[i];
2149              upper[i]=lower[i];
2150              nFix++;
2151            } else if (value>upper[i]-1.0e-5&&value-lower[i]>upper[i]-value) {
2152              solution[i]=upper[i];
2153              lower[i]=upper[i];
2154              nFix++;
2155            } else {
2156              printf("fixrow %d %g <= %g <= %g\n",
2157                     i,lower[i],solution[i],upper[i]);
2158            }
2159          }
2160        }
2161        printf("%d row slacks fixed\n",nFix);
2162      }
2163    }
2164#ifdef BORROW   
2165    barrier.returnModel(*model2);
2166    double * rowPrimal = new double [numberRows];
2167    double * columnPrimal = new double [numberColumns];
2168    double * rowDual = new double [numberRows];
2169    double * columnDual = new double [numberColumns];
2170    // move solutions other way
2171    CoinMemcpyN(model2->primalRowSolution(),
2172                numberRows,rowPrimal);
2173    CoinMemcpyN(model2->dualRowSolution(),
2174                numberRows,rowDual);
2175    CoinMemcpyN(model2->primalColumnSolution(),
2176                numberColumns,columnPrimal);
2177    CoinMemcpyN(model2->dualColumnSolution(),
2178                  numberColumns,columnDual);
2179#else
2180    double * rowPrimal = barrier.primalRowSolution();
2181    double * columnPrimal = barrier.primalColumnSolution();
2182    double * rowDual = barrier.dualRowSolution();
2183    double * columnDual = barrier.dualColumnSolution();
2184    // move solutions
2185    CoinMemcpyN(rowPrimal,
2186                numberRows,model2->primalRowSolution());
2187    CoinMemcpyN(rowDual,
2188                numberRows,model2->dualRowSolution());
2189    CoinMemcpyN(columnPrimal,
2190                numberColumns,model2->primalColumnSolution());
2191    CoinMemcpyN(columnDual,
2192                  numberColumns,model2->dualColumnSolution());
2193#endif
2194    if (saveModel2) {
2195      // do presolve
2196      model2 = pinfo2.presolvedModel(*model2,dblParam_[ClpPresolveTolerance],
2197                                    false,5,true);
2198      if (!model2) {
2199        model2=saveModel2;
2200        saveModel2=NULL;
2201        int numberRows = model2->numberRows();
2202        int numberColumns = model2->numberColumns();
2203        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2204        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2205        delete [] saveLower;
2206        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2207        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2208        delete [] saveUpper;
2209        saveLower=NULL;
2210        saveUpper=NULL;
2211      }
2212    }
2213    if (method==ClpSolve::useBarrier) {
2214      if (maxIts&&barrierStatus<4&&!quadraticObj) {
2215        //printf("***** crossover - needs more thought on difficult models\n");
2216#if SAVEIT==1
2217        model2->ClpSimplex::saveModel("xx.save");
2218#endif
2219        // make sure no status left
2220        model2->createStatus();
2221        // solve
2222        if (!forceFixing)
2223          model2->setPerturbation(100);
2224        if (model2->factorizationFrequency()==200) {
2225          // User did not touch preset
2226          model2->defaultFactorizationFrequency();
2227        }
2228#if 1
2229        // throw some into basis
2230        if(!forceFixing) {
2231          int numberRows = model2->numberRows();
2232          int numberColumns = model2->numberColumns();
2233          double * dsort = new double[numberColumns];
2234          int * sort = new int[numberColumns];
2235          int n=0;
2236          const double * columnLower = model2->columnLower();
2237          const double * columnUpper = model2->columnUpper();
2238          double * primalSolution = model2->primalColumnSolution();
2239          const double * dualSolution = model2->dualColumnSolution();
2240          double tolerance = 10.0*primalTolerance_;
2241          int i;
2242          for ( i=0;i<numberRows;i++) 
2243            model2->setRowStatus(i,superBasic);
2244          for ( i=0;i<numberColumns;i++) {
2245            double distance = CoinMin(columnUpper[i]-primalSolution[i],
2246                                      primalSolution[i]-columnLower[i]);
2247            if (distance>tolerance) {
2248              if (fabs(dualSolution[i])<1.0e-5)
2249                distance *= 100.0;
2250              dsort[n]=-distance;
2251              sort[n++]=i;
2252              model2->setStatus(i,superBasic);
2253            } else if (distance>primalTolerance_) {
2254              model2->setStatus(i,superBasic);
2255            } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
2256              model2->setStatus(i,atLowerBound);
2257              primalSolution[i]=columnLower[i];
2258            } else {
2259              model2->setStatus(i,atUpperBound);
2260              primalSolution[i]=columnUpper[i];
2261            }
2262          }
2263          CoinSort_2(dsort,dsort+n,sort);
2264          n = CoinMin(numberRows,n);
2265          for ( i=0;i<n;i++) {
2266            int iColumn = sort[i];
2267            model2->setStatus(iColumn,basic);
2268          }
2269          delete [] sort;
2270          delete [] dsort;
2271          // model2->allSlackBasis();
2272          if (gap<1.0e-3*static_cast<double> (numberRows+numberColumns)) {
2273            if (saveUpper) {
2274              int numberRows = model2->numberRows();
2275              int numberColumns = model2->numberColumns();
2276              CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2277              CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2278              CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2279              CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2280              delete [] saveLower;
2281              delete [] saveUpper;
2282              saveLower=NULL;
2283              saveUpper=NULL;
2284            }
2285            int numberRows = model2->numberRows();
2286            int numberColumns = model2->numberColumns();
2287            // just primal values pass
2288            double saveScale = model2->objectiveScale();
2289            model2->setObjectiveScale(1.0e-3);
2290            model2->primal(2);
2291            model2->setObjectiveScale(saveScale);
2292            // save primal solution and copy back dual
2293            CoinMemcpyN(model2->primalRowSolution(),
2294                        numberRows,rowPrimal);
2295            CoinMemcpyN(rowDual,
2296                        numberRows,model2->dualRowSolution());
2297            CoinMemcpyN(model2->primalColumnSolution(),
2298                        numberColumns,columnPrimal);
2299            CoinMemcpyN(columnDual,
2300                        numberColumns,model2->dualColumnSolution());
2301            //model2->primal(1);
2302            // clean up reduced costs and flag variables
2303            {
2304              double * dj = model2->dualColumnSolution();
2305              double * cost = model2->objective();
2306              double * saveCost = new double[numberColumns];
2307              CoinMemcpyN(cost,numberColumns,saveCost);
2308              double * saveLower = new double[numberColumns];
2309              double * lower = model2->columnLower();
2310              CoinMemcpyN(lower,numberColumns,saveLower);
2311              double * saveUpper = new double[numberColumns];
2312              double * upper = model2->columnUpper();
2313              CoinMemcpyN(upper,numberColumns,saveUpper);
2314              int i;
2315              double tolerance = 10.0*dualTolerance_;
2316              for ( i=0;i<numberColumns;i++) {
2317                if (model2->getStatus(i)==basic) {
2318                  dj[i]=0.0;
2319                } else if (model2->getStatus(i)==atLowerBound) {
2320                  if (optimizationDirection_*dj[i]<tolerance) {
2321                    if (optimizationDirection_*dj[i]<0.0) {
2322                      //if (dj[i]<-1.0e-3)
2323                      //printf("bad dj at lb %d %g\n",i,dj[i]);
2324                      cost[i] -= dj[i];
2325                      dj[i]=0.0;
2326                    }
2327                  } else {
2328                    upper[i]=lower[i];
2329                  }
2330                } else if (model2->getStatus(i)==atUpperBound) {
2331                  if (optimizationDirection_*dj[i]>tolerance) {
2332                    if (optimizationDirection_*dj[i]>0.0) {
2333                      //if (dj[i]>1.0e-3)
2334                      //printf("bad dj at ub %d %g\n",i,dj[i]);
2335                      cost[i] -= dj[i];
2336                      dj[i]=0.0;
2337                    }
2338                  } else {
2339                    lower[i]=upper[i];
2340                  }
2341                }
2342              }
2343              // just dual values pass
2344              //model2->setLogLevel(63);
2345              //model2->setFactorizationFrequency(1);
2346              model2->dual(2);
2347              CoinMemcpyN(saveCost,numberColumns,cost);
2348              delete [] saveCost;
2349              CoinMemcpyN(saveLower,numberColumns,lower);
2350              delete [] saveLower;
2351              CoinMemcpyN(saveUpper,numberColumns,upper);
2352              delete [] saveUpper;
2353            }
2354          }
2355          // and finish
2356          // move solutions
2357          CoinMemcpyN(rowPrimal,
2358                      numberRows,model2->primalRowSolution());
2359          CoinMemcpyN(columnPrimal,
2360                      numberColumns,model2->primalColumnSolution());
2361        }
2362        double saveScale = model2->objectiveScale();
2363        model2->setObjectiveScale(1.0e-3);
2364        model2->primal(2);
2365        model2->setObjectiveScale(saveScale);
2366        model2->primal(1);
2367#else
2368        // just primal
2369        model2->primal(1);
2370#endif
2371      } else if (barrierStatus==4) {
2372        // memory problems
2373        model2->setPerturbation(savePerturbation);
2374        model2->createStatus();
2375        model2->dual();
2376      } else if (maxIts&&quadraticObj) {
2377        // make sure no status left
2378        model2->createStatus();
2379        // solve
2380        model2->setPerturbation(100);
2381        model2->reducedGradient(1);
2382      }
2383    }
2384    model2->setMaximumIterations(saveMaxIts);
2385#ifdef BORROW
2386    delete [] rowPrimal;
2387    delete [] columnPrimal;
2388    delete [] rowDual;
2389    delete [] columnDual;
2390#endif
2391    if (extraPresolve) {
2392      pinfo2.postsolve(true);
2393      delete model2;
2394      model2=saveModel2;
2395    }
2396    if (saveUpper) {
2397      if (!forceFixing) {
2398        int numberRows = model2->numberRows();
2399        int numberColumns = model2->numberColumns();
2400        CoinMemcpyN(saveLower,numberColumns,model2->columnLower());
2401        CoinMemcpyN(saveLower+numberColumns,numberRows,model2->rowLower());
2402        CoinMemcpyN(saveUpper,numberColumns,model2->columnUpper());
2403        CoinMemcpyN(saveUpper+numberColumns,numberRows,model2->rowUpper());
2404      }
2405      delete [] saveLower;
2406      delete [] saveUpper;
2407      saveLower=NULL;
2408      saveUpper=NULL;
2409      if (method!=ClpSolve::useBarrierNoCross) 
2410        model2->primal(1);
2411    }
2412    model2->setPerturbation(savePerturbation);
2413    time2 = CoinCpuTime();
2414    timeCore = time2-timeX;
2415    handler_->message(CLP_INTERVAL_TIMING,messages_)
2416      <<"Crossover"<<timeCore<<time2-time1
2417      <<CoinMessageEol;
2418    timeX=time2;
2419#else
2420    abort();
2421#endif
2422  } else {
2423    assert (method!=ClpSolve::automatic); // later
2424    time2=0.0;
2425  }
2426  if (saveMatrix) {
2427    if (model2==this) {
2428      // delete and replace
2429      delete model2->clpMatrix();
2430      model2->replaceMatrix(saveMatrix);
2431    } else {
2432      delete saveMatrix;
2433    }
2434  }
2435  numberIterations = model2->numberIterations();
2436  finalStatus=model2->status();
2437  int finalSecondaryStatus = model2->secondaryStatus();
2438  if (presolve==ClpSolve::presolveOn) {
2439    int saveLevel = logLevel();
2440    if ((specialOptions_&1024)==0)
2441      setLogLevel(CoinMin(1,saveLevel));
2442    else
2443      setLogLevel(CoinMin(0,saveLevel));
2444    pinfo.postsolve(true);
2445    factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
2446    time2 = CoinCpuTime();
2447    timePresolve += time2-timeX;
2448    handler_->message(CLP_INTERVAL_TIMING,messages_)
2449      <<"Postsolve"<<time2-timeX<<time2-time1
2450      <<CoinMessageEol;
2451    timeX=time2;
2452    if (!presolveToFile)
2453      delete model2;
2454    if (interrupt)
2455      currentModel = this;
2456    // checkSolution(); already done by postSolve
2457    setLogLevel(saveLevel);
2458    if (finalStatus!=3&&(finalStatus||status()==-1)) {
2459      int savePerturbation = perturbation();
2460      if (!finalStatus||(moreSpecialOptions_&2)==0) {
2461        if (finalStatus==2) {
2462          // unbounded - get feasible first
2463          double save = optimizationDirection_;
2464          optimizationDirection_=0.0;
2465          primal(1);
2466          optimizationDirection_=save;
2467          primal(1);
2468        } else if (finalStatus==1) {
2469          dual();
2470        } else {
2471          setPerturbation(100);
2472          primal(1);
2473        }
2474      } else {
2475        // just set status
2476        problemStatus_=finalStatus;
2477      }
2478      setPerturbation(savePerturbation);
2479      numberIterations += numberIterations_;
2480      numberIterations_ = numberIterations;
2481      finalStatus=status();
2482      time2 = CoinCpuTime();
2483      handler_->message(CLP_INTERVAL_TIMING,messages_)
2484      <<"Cleanup"<<time2-timeX<<time2-time1
2485      <<CoinMessageEol;
2486      timeX=time2;
2487    } else {
2488      secondaryStatus_=finalSecondaryStatus;
2489    }
2490  } else if (model2!=this) {
2491    // not presolved - but different model used (sprint probably)
2492    CoinMemcpyN(model2->primalRowSolution(),
2493                numberRows_,this->primalRowSolution());
2494    CoinMemcpyN(model2->dualRowSolution(),
2495                numberRows_,this->dualRowSolution());
2496    CoinMemcpyN(model2->primalColumnSolution(),
2497                numberColumns_,this->primalColumnSolution());
2498    CoinMemcpyN(model2->dualColumnSolution(),
2499                numberColumns_,this->dualColumnSolution());
2500    CoinMemcpyN(model2->statusArray(),
2501                numberColumns_+numberRows_,this->statusArray());
2502    objectiveValue_=model2->objectiveValue_;
2503    numberIterations_ =model2->numberIterations_;
2504    problemStatus_ =model2->problemStatus_;
2505    secondaryStatus_ =model2->secondaryStatus_;
2506    delete model2;
2507  }
2508  setMaximumIterations(saveMaxIterations);
2509  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
2510                               "Errors","User stopped"};
2511  assert (finalStatus>=-1&&finalStatus<=5);
2512  handler_->message(CLP_TIMING,messages_)
2513    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
2514  handler_->printing(presolve==ClpSolve::presolveOn)
2515    <<timePresolve;
2516  handler_->printing(timeIdiot!=0.0)
2517    <<timeIdiot;
2518  handler_->message()<<CoinMessageEol;
2519  if (interrupt) 
2520    signal(SIGINT,saveSignal);
2521  perturbation_=savePerturbation;
2522  scalingFlag_=saveScaling;
2523  // If faking objective - put back correct one
2524  if (savedObjective) {
2525    delete objective_;
2526    objective_=savedObjective;
2527  }
2528  return finalStatus;
2529}
2530// General solve
2531int 
2532ClpSimplex::initialSolve()
2533{
2534  // Default so use dual
2535  ClpSolve options;
2536  return initialSolve(options);
2537}
2538// General dual solve
2539int 
2540ClpSimplex::initialDualSolve()
2541{
2542  ClpSolve options;
2543  // Use dual
2544  options.setSolveType(ClpSolve::useDual);
2545  return initialSolve(options);
2546}
2547// General dual solve
2548int 
2549ClpSimplex::initialPrimalSolve()
2550{
2551  ClpSolve options;
2552  // Use primal
2553  options.setSolveType(ClpSolve::usePrimal);
2554  return initialSolve(options);
2555}
2556// barrier solve, not to be followed by crossover
2557int 
2558ClpSimplex::initialBarrierNoCrossSolve()
2559{
2560  ClpSolve options;
2561  // Use primal
2562  options.setSolveType(ClpSolve::useBarrierNoCross);
2563  return initialSolve(options);
2564}
2565
2566// General barrier solve
2567int 
2568ClpSimplex::initialBarrierSolve()
2569{
2570  ClpSolve options;
2571  // Use primal
2572  options.setSolveType(ClpSolve::useBarrier);
2573  return initialSolve(options);
2574}
2575
2576// Default constructor
2577ClpSolve::ClpSolve (  )
2578{
2579  method_ = automatic;
2580  presolveType_=presolveOn;
2581  numberPasses_=5;
2582  int i;
2583  for (i=0;i<7;i++)
2584    options_[i]=0;
2585  // say no +-1 matrix
2586  options_[3]=1;
2587  for (i=0;i<7;i++)
2588    extraInfo_[i]=-1;
2589  independentOptions_[0]=0;
2590  // But switch off slacks
2591  independentOptions_[1]=512;
2592  // Substitute up to 3
2593  independentOptions_[2]=3;
2594 
2595}
2596// Constructor when you really know what you are doing
2597ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
2598             int numberPasses, int options[6],
2599             int extraInfo[6], int independentOptions[3])
2600{
2601  method_ = method;
2602  presolveType_=presolveType;
2603  numberPasses_=numberPasses;
2604  int i;
2605  for (i=0;i<6;i++)
2606    options_[i]=options[i];
2607  options_[6]=0;
2608  for (i=0;i<6;i++)
2609    extraInfo_[i]=extraInfo[i];
2610  extraInfo_[6]=0;
2611  for (i=0;i<3;i++)
2612    independentOptions_[i]=independentOptions[i];
2613}
2614
2615// Copy constructor.
2616ClpSolve::ClpSolve(const ClpSolve & rhs)
2617{
2618  method_ = rhs.method_;
2619  presolveType_=rhs.presolveType_;
2620  numberPasses_=rhs.numberPasses_;
2621  int i;
2622  for ( i=0;i<7;i++)
2623    options_[i]=rhs.options_[i];
2624  for ( i=0;i<7;i++)
2625    extraInfo_[i]=rhs.extraInfo_[i];
2626  for (i=0;i<3;i++)
2627    independentOptions_[i]=rhs.independentOptions_[i];
2628}
2629// Assignment operator. This copies the data
2630ClpSolve & 
2631ClpSolve::operator=(const ClpSolve & rhs)
2632{
2633  if (this != &rhs) {
2634    method_ = rhs.method_;
2635    presolveType_=rhs.presolveType_;
2636    numberPasses_=rhs.numberPasses_;
2637    int i;
2638    for (i=0;i<7;i++)
2639      options_[i]=rhs.options_[i];
2640    for (i=0;i<7;i++)
2641      extraInfo_[i]=rhs.extraInfo_[i];
2642    for (i=0;i<3;i++)
2643      independentOptions_[i]=rhs.independentOptions_[i];
2644  }
2645  return *this;
2646
2647}
2648// Destructor
2649ClpSolve::~ClpSolve (  )
2650{
2651}
2652// See header file for deatils
2653void 
2654ClpSolve::setSpecialOption(int which,int value,int extraInfo)
2655{
2656  options_[which]=value;
2657  extraInfo_[which]=extraInfo;
2658}
2659int 
2660ClpSolve::getSpecialOption(int which) const
2661{
2662  return options_[which];
2663}
2664
2665// Solve types
2666void 
2667ClpSolve::setSolveType(SolveType method, int extraInfo)
2668{
2669  method_=method;
2670}
2671
2672ClpSolve::SolveType
2673ClpSolve::getSolveType()
2674{
2675  return method_;
2676}
2677
2678// Presolve types
2679void 
2680ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
2681{
2682  presolveType_=amount;
2683  numberPasses_=extraInfo;
2684}
2685ClpSolve::PresolveType
2686ClpSolve::getPresolveType()
2687{
2688  return presolveType_;
2689}
2690// Extra info for idiot (or sprint)
2691int 
2692ClpSolve::getExtraInfo(int which) const
2693{
2694  return extraInfo_[which];
2695}
2696int 
2697ClpSolve::getPresolvePasses() const
2698{
2699  return numberPasses_;
2700}
2701/* Say to return at once if infeasible,
2702   default is to solve */
2703void 
2704ClpSolve::setInfeasibleReturn(bool trueFalse)
2705{
2706  independentOptions_[0]= trueFalse ? 1 : 0;
2707}
2708#include <string>
2709// Generates code for above constructor
2710void 
2711ClpSolve::generateCpp(FILE * fp)
2712{
2713  std::string solveType[] = {
2714    "ClpSolve::useDual",
2715    "ClpSolve::usePrimal",
2716    "ClpSolve::usePrimalorSprint",
2717    "ClpSolve::useBarrier",
2718    "ClpSolve::useBarrierNoCross",
2719    "ClpSolve::automatic",
2720    "ClpSolve::notImplemented"
2721  };
2722  std::string presolveType[] =  {
2723    "ClpSolve::presolveOn",
2724    "ClpSolve::presolveOff",
2725    "ClpSolve::presolveNumber",
2726    "ClpSolve::presolveNumberCost"
2727  };
2728  fprintf(fp,"3  ClpSolve::SolveType method = %s;\n",solveType[method_].c_str());
2729  fprintf(fp,"3  ClpSolve::PresolveType presolveType = %s;\n",
2730    presolveType[presolveType_].c_str());
2731  fprintf(fp,"3  int numberPasses = %d;\n",numberPasses_);
2732  fprintf(fp,"3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
2733    options_[0],options_[1],options_[2],
2734    options_[3],options_[4],options_[5]);
2735  fprintf(fp,"3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
2736    extraInfo_[0],extraInfo_[1],extraInfo_[2],
2737    extraInfo_[3],extraInfo_[4],extraInfo_[5]);
2738  fprintf(fp,"3  int independentOptions[] = {%d,%d,%d};\n",
2739    independentOptions_[0],independentOptions_[1],independentOptions_[2]);
2740  fprintf(fp,"3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
2741  fprintf(fp,"3                    options,extraInfo,independentOptions);\n");
2742}
2743//#############################################################################
2744#include "ClpNonLinearCost.hpp"
2745
2746ClpSimplexProgress::ClpSimplexProgress () 
2747{
2748  int i;
2749  for (i=0;i<CLP_PROGRESS;i++) {
2750    objective_[i] = COIN_DBL_MAX;
2751    infeasibility_[i] = -1.0; // set to an impossible value
2752    realInfeasibility_[i] = COIN_DBL_MAX;
2753    numberInfeasibilities_[i]=-1; 
2754    iterationNumber_[i]=-1;
2755  }
2756#ifdef CLP_PROGRESS_WEIGHT
2757  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2758    objectiveWeight_[i] = COIN_DBL_MAX;
2759    infeasibilityWeight_[i] = -1.0; // set to an impossible value
2760    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
2761    numberInfeasibilitiesWeight_[i]=-1; 
2762    iterationNumberWeight_[i]=-1;
2763  }
2764  drop_ =0.0;
2765  best_ =0.0;
2766#endif
2767  initialWeight_=0.0;
2768  for (i=0;i<CLP_CYCLE;i++) {
2769    //obj_[i]=COIN_DBL_MAX;
2770    in_[i]=-1;
2771    out_[i]=-1;
2772    way_[i]=0;
2773  }
2774  numberTimes_ = 0;
2775  numberBadTimes_ = 0;
2776  numberReallyBadTimes_ = 0;
2777  numberTimesFlagged_ = 0;
2778  model_ = NULL;
2779  oddState_=0;
2780}
2781
2782
2783//-----------------------------------------------------------------------------
2784
2785ClpSimplexProgress::~ClpSimplexProgress ()
2786{
2787}
2788// Copy constructor.
2789ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
2790{
2791  int i;
2792  for (i=0;i<CLP_PROGRESS;i++) {
2793    objective_[i] = rhs.objective_[i];
2794    infeasibility_[i] = rhs.infeasibility_[i];
2795    realInfeasibility_[i] = rhs.realInfeasibility_[i];
2796    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2797    iterationNumber_[i]=rhs.iterationNumber_[i];
2798  }
2799#ifdef CLP_PROGRESS_WEIGHT
2800  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2801    objectiveWeight_[i] = rhs.objectiveWeight_[i];
2802    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2803    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2804    numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i]; 
2805    iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
2806  }
2807  drop_ = rhs.drop_;
2808  best_ = rhs.best_;
2809#endif
2810  initialWeight_ = rhs.initialWeight_;
2811  for (i=0;i<CLP_CYCLE;i++) {
2812    //obj_[i]=rhs.obj_[i];
2813    in_[i]=rhs.in_[i];
2814    out_[i]=rhs.out_[i];
2815    way_[i]=rhs.way_[i];
2816  }
2817  numberTimes_ = rhs.numberTimes_;
2818  numberBadTimes_ = rhs.numberBadTimes_;
2819  numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
2820  numberTimesFlagged_ = rhs.numberTimesFlagged_;
2821  model_ = rhs.model_;
2822  oddState_=rhs.oddState_;
2823}
2824// Copy constructor.from model
2825ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
2826{
2827  model_ = model;
2828  reset();
2829  initialWeight_=0.0;
2830}
2831// Fill from model
2832void 
2833ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
2834{
2835  model_ = model;
2836  reset();
2837  initialWeight_=0.0;
2838}
2839// Assignment operator. This copies the data
2840ClpSimplexProgress & 
2841ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2842{
2843  if (this != &rhs) {
2844    int i;
2845    for (i=0;i<CLP_PROGRESS;i++) {
2846      objective_[i] = rhs.objective_[i];
2847      infeasibility_[i] = rhs.infeasibility_[i];
2848      realInfeasibility_[i] = rhs.realInfeasibility_[i];
2849      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2850      iterationNumber_[i]=rhs.iterationNumber_[i];
2851    }
2852#ifdef CLP_PROGRESS_WEIGHT
2853    for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
2854      objectiveWeight_[i] = rhs.objectiveWeight_[i];
2855      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2856      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2857      numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i]; 
2858      iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
2859    }
2860    drop_ = rhs.drop_;
2861    best_ = rhs.best_;
2862#endif
2863    initialWeight_ = rhs.initialWeight_;
2864    for (i=0;i<CLP_CYCLE;i++) {
2865      //obj_[i]=rhs.obj_[i];
2866      in_[i]=rhs.in_[i];
2867      out_[i]=rhs.out_[i];
2868      way_[i]=rhs.way_[i];
2869    }
2870    numberTimes_ = rhs.numberTimes_;
2871    numberBadTimes_ = rhs.numberBadTimes_;
2872    numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
2873    numberTimesFlagged_ = rhs.numberTimesFlagged_;
2874    model_ = rhs.model_;
2875    oddState_=rhs.oddState_;
2876  }
2877  return *this;
2878}
2879// Seems to be something odd about exact comparison of doubles on linux
2880static bool equalDouble(double value1, double value2)
2881{
2882
2883  union { double d; int i[2]; } v1,v2;
2884  v1.d = value1;
2885  v2.d = value2;
2886  if (sizeof(int)*2==sizeof(double)) 
2887    return (v1.i[0]==v2.i[0]&&v1.i[1]==v2.i[1]);
2888  else
2889    return (v1.i[0]==v2.i[0]);
2890}
2891int
2892ClpSimplexProgress::looping()
2893{
2894  if (!model_)
2895    return -1;
2896  double objective = model_->rawObjectiveValue();
2897  if (model_->algorithm()<0)
2898    objective -= model_->bestPossibleImprovement();
2899  double infeasibility;
2900  double realInfeasibility=0.0;
2901  int numberInfeasibilities;
2902  int iterationNumber = model_->numberIterations();
2903  numberTimesFlagged_=0;
2904  if (model_->algorithm()<0) {
2905    // dual
2906    infeasibility = model_->sumPrimalInfeasibilities();
2907    numberInfeasibilities = model_->numberPrimalInfeasibilities();
2908  } else {
2909    //primal
2910    infeasibility = model_->sumDualInfeasibilities();
2911    realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
2912    numberInfeasibilities = model_->numberDualInfeasibilities();
2913  }
2914  int i;
2915  int numberMatched=0;
2916  int matched=0;
2917  int nsame=0;
2918  for (i=0;i<CLP_PROGRESS;i++) {
2919    bool matchedOnObjective = equalDouble(objective,objective_[i]);
2920    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
2921    bool matchedOnInfeasibilities = 
2922      (numberInfeasibilities==numberInfeasibilities_[i]);
2923   
2924    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
2925      matched |= (1<<i);
2926      // Check not same iteration
2927      if (iterationNumber!=iterationNumber_[i]) {
2928        numberMatched++;
2929        // here mainly to get over compiler bug?
2930        if (model_->messageHandler()->logLevel()>10)
2931          printf("%d %d %d %d %d loop check\n",i,numberMatched,
2932                 matchedOnObjective, matchedOnInfeasibility, 
2933                 matchedOnInfeasibilities);
2934      } else {
2935        // stuck but code should notice
2936        nsame++;
2937      }
2938    }
2939    if (i) {
2940      objective_[i-1] = objective_[i];
2941      infeasibility_[i-1] = infeasibility_[i];
2942      realInfeasibility_[i-1] = realInfeasibility_[i];
2943      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
2944      iterationNumber_[i-1]=iterationNumber_[i];
2945    }
2946  }
2947  objective_[CLP_PROGRESS-1] = objective;
2948  infeasibility_[CLP_PROGRESS-1] = infeasibility;
2949  realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
2950  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
2951  iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
2952  if (nsame==CLP_PROGRESS)
2953    numberMatched=CLP_PROGRESS; // really stuck
2954  if (model_->progressFlag())
2955    numberMatched=0;
2956  numberTimes_++;
2957  if (numberTimes_<10)
2958    numberMatched=0;
2959  // skip if just last time as may be checking something
2960  if (matched==(1<<(CLP_PROGRESS-1)))
2961    numberMatched=0;
2962  if (numberMatched) {
2963    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
2964      <<numberMatched
2965      <<matched
2966      <<numberTimes_
2967      <<CoinMessageEol;
2968    numberBadTimes_++;
2969    if (numberBadTimes_<10) {
2970      // make factorize every iteration
2971      model_->forceFactorization(1);
2972      if (numberBadTimes_<2) {
2973        startCheck(); // clear other loop check
2974        if (model_->algorithm()<0) {
2975          // dual - change tolerance
2976          model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
2977          // if infeasible increase dual bound
2978          if (model_->dualBound()<1.0e17) {
2979            model_->setDualBound(model_->dualBound()*1.1);
2980            static_cast<ClpSimplexDual *>(model_)->resetFakeBounds(0);
2981          }
2982        } else {
2983          // primal - change tolerance
2984          if (numberBadTimes_>3)
2985            model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
2986          // if infeasible increase infeasibility cost
2987          if (model_->nonLinearCost()->numberInfeasibilities()&&
2988              model_->infeasibilityCost()<1.0e17) {
2989            model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
2990          }
2991        }
2992      } else {
2993        // flag
2994        int iSequence;
2995        if (model_->algorithm()<0) {
2996          // dual
2997          if (model_->dualBound()>1.0e14) 
2998            model_->setDualBound(1.0e14);
2999          iSequence=in_[CLP_CYCLE-1];
3000        } else {
3001          // primal
3002          if (model_->infeasibilityCost()>1.0e14) 
3003            model_->setInfeasibilityCost(1.0e14);
3004          iSequence=out_[CLP_CYCLE-1];
3005        }
3006        if (iSequence>=0) {
3007          char x = model_->isColumn(iSequence) ? 'C' :'R';
3008          if (model_->messageHandler()->logLevel()>=63)
3009            model_->messageHandler()->message(CLP_SIMPLEX_FLAG,model_->messages())
3010              <<x<<model_->sequenceWithin(iSequence)
3011              <<CoinMessageEol;
3012          // if Gub then needs to be sequenceIn_
3013          int save=model_->sequenceIn();
3014          model_->setSequenceIn(iSequence);
3015          model_->setFlagged(iSequence);
3016          model_->setSequenceIn(save);
3017          //printf("flagging %d from loop\n",iSequence);
3018          startCheck();
3019        } else {
3020          // Give up
3021          if (model_->messageHandler()->logLevel()>=63)
3022            printf("***** All flagged?\n");
3023          return 4;
3024        }
3025        // reset
3026        numberBadTimes_=2;
3027      }
3028      return -2;
3029    } else {
3030      // look at solution and maybe declare victory
3031      if (infeasibility<1.0e-4) {
3032        return 0;
3033      } else {
3034        model_->messageHandler()->message(CLP_LOOP,model_->messages())
3035          <<CoinMessageEol;
3036#ifndef NDEBUG
3037        printf("debug loop ClpSimplex A\n");
3038        abort();
3039#endif
3040        return 3;
3041      }
3042    }
3043  }
3044  return -1;
3045}
3046// Resets as much as possible
3047void 
3048ClpSimplexProgress::reset()
3049{
3050  int i;
3051  for (i=0;i<CLP_PROGRESS;i++) {
3052    if (model_->algorithm()>=0)
3053      objective_[i] = COIN_DBL_MAX;
3054    else
3055      objective_[i] = -COIN_DBL_MAX;
3056    infeasibility_[i] = -1.0; // set to an impossible value
3057    realInfeasibility_[i] = COIN_DBL_MAX;
3058    numberInfeasibilities_[i]=-1; 
3059    iterationNumber_[i]=-1;
3060  }
3061#ifdef CLP_PROGRESS_WEIGHT
3062  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
3063    objectiveWeight_[i] = COIN_DBL_MAX;
3064    infeasibilityWeight_[i] = -1.0; // set to an impossible value
3065    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3066    numberInfeasibilitiesWeight_[i]=-1; 
3067    iterationNumberWeight_[i]=-1;
3068  }
3069  drop_ =0.0;
3070  best_ =0.0;
3071#endif
3072  for (i=0;i<CLP_CYCLE;i++) {
3073    //obj_[i]=COIN_DBL_MAX;
3074    in_[i]=-1;
3075    out_[i]=-1;
3076    way_[i]=0;
3077  }
3078  numberTimes_ = 0;
3079  numberBadTimes_ = 0;
3080  numberReallyBadTimes_ = 0;
3081  numberTimesFlagged_ = 0;
3082  oddState_=0;
3083}
3084// Returns previous objective (if -1) - current if (0)
3085double 
3086ClpSimplexProgress::lastObjective(int back) const
3087{
3088  return objective_[CLP_PROGRESS-1-back];
3089}
3090// Returns previous infeasibility (if -1) - current if (0)
3091double 
3092ClpSimplexProgress::lastInfeasibility(int back) const
3093{
3094  return realInfeasibility_[CLP_PROGRESS-1-back];
3095}
3096// Sets real primal infeasibility
3097void
3098ClpSimplexProgress::setInfeasibility(double value)
3099{
3100  for (int i=1;i<CLP_PROGRESS;i++) 
3101    realInfeasibility_[i-1] = realInfeasibility_[i];
3102  realInfeasibility_[CLP_PROGRESS-1]=value;
3103}
3104// Modify objective e.g. if dual infeasible in dual
3105void 
3106ClpSimplexProgress::modifyObjective(double value)
3107{
3108  objective_[CLP_PROGRESS-1]=value;
3109}
3110// Returns previous iteration number (if -1) - current if (0)
3111int 
3112ClpSimplexProgress::lastIterationNumber(int back) const
3113{
3114  return iterationNumber_[CLP_PROGRESS-1-back];
3115}
3116// clears iteration numbers (to switch off panic)
3117void 
3118ClpSimplexProgress::clearIterationNumbers()
3119{
3120  for (int i=0;i<CLP_PROGRESS;i++) 
3121    iterationNumber_[i]=-1;
3122}
3123// Start check at beginning of whileIterating
3124void 
3125ClpSimplexProgress::startCheck()
3126{
3127  int i;
3128  for (i=0;i<CLP_CYCLE;i++) {
3129    //obj_[i]=COIN_DBL_MAX;
3130    in_[i]=-1;
3131    out_[i]=-1;
3132    way_[i]=0;
3133  }
3134}
3135// Returns cycle length in whileIterating
3136int 
3137ClpSimplexProgress::cycle(int in, int out,int wayIn,int wayOut)
3138{
3139  int i;
3140#if 0
3141  if (model_->numberIterations()>206571) {
3142    printf("in %d out %d\n",in,out);
3143    for (i=0;i<CLP_CYCLE;i++)
3144      printf("cy %d in %d out %d\n",i,in_[i],out_[i]);
3145  }
3146#endif
3147  int matched=0;
3148  // first see if in matches any out
3149  for (i=1;i<CLP_CYCLE;i++) {
3150    if (in==out_[i]) {
3151      // even if flip then suspicious
3152      matched=-1;
3153      break;
3154    }
3155  }
3156#if 0
3157  if (!matched||in_[0]<0) {
3158    // can't be cycle
3159    for (i=0;i<CLP_CYCLE-1;i++) {
3160      //obj_[i]=obj_[i+1];
3161      in_[i]=in_[i+1];
3162      out_[i]=out_[i+1];
3163      way_[i]=way_[i+1];
3164    }
3165  } else {
3166    // possible cycle
3167    matched=0;
3168    for (i=0;i<CLP_CYCLE-1;i++) {
3169      int k;
3170      char wayThis = way_[i];
3171      int inThis = in_[i];
3172      int outThis = out_[i];
3173      //double objThis = obj_[i];
3174      for(k=i+1;k<CLP_CYCLE;k++) {
3175        if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
3176          int distance = k-i;
3177          if (k+distance<CLP_CYCLE) {
3178            // See if repeats
3179            int j=k+distance;
3180            if (inThis==in_[j]&&outThis==out_[j]&&wayThis==way_[j]) {
3181              matched=distance;
3182              break;
3183            }
3184          } else {
3185            matched=distance;
3186            break;
3187          }
3188        }
3189      }
3190      //obj_[i]=obj_[i+1];
3191      in_[i]=in_[i+1];
3192      out_[i]=out_[i+1];
3193      way_[i]=way_[i+1];
3194    }
3195  }
3196#else
3197  if (matched&&in_[0]>=0) {
3198    // possible cycle - only check [0] against all
3199    matched=0;
3200    int nMatched=0;
3201    char way0 = way_[0];
3202    int in0 = in_[0];
3203    int out0 = out_[0];
3204    //double obj0 = obj_[i];
3205    for(int k=1;k<CLP_CYCLE-4;k++) {
3206      if (in0==in_[k]&&out0==out_[k]&&way0==way_[k]) {
3207        nMatched++;
3208        // See if repeats
3209        int end = CLP_CYCLE-k;
3210        int j;
3211        for ( j=1;j<end;j++) {
3212          if (in_[j+k]!=in_[j]||out_[j+k]!=out_[j]||way_[j+k]!=way_[j])
3213            break;
3214        }
3215        if (j==end) {
3216          matched=k;
3217          break;
3218        }
3219      }
3220    }
3221    // If three times then that is too much even if not regular
3222    if (matched<=0&&nMatched>1)
3223      matched=100;
3224  }
3225  for (i=0;i<CLP_CYCLE-1;i++) {
3226    //obj_[i]=obj_[i+1];
3227    in_[i]=in_[i+1];
3228    out_[i]=out_[i+1];
3229    way_[i]=way_[i+1];
3230  }
3231#endif
3232  int way = 1-wayIn+4*(1-wayOut);
3233  //obj_[i]=model_->objectiveValue();
3234  in_[CLP_CYCLE-1]=in;
3235  out_[CLP_CYCLE-1]=out;
3236  way_[CLP_CYCLE-1]=static_cast<char>(way);
3237  return matched;
3238}
3239#include "CoinStructuredModel.hpp"
3240// Solve using structure of model and maybe in parallel
3241int 
3242ClpSimplex::solve(CoinStructuredModel * model)
3243{
3244  // analyze structure
3245  int numberRowBlocks=model->numberRowBlocks();
3246  int numberColumnBlocks = model->numberColumnBlocks();
3247  int numberElementBlocks = model->numberElementBlocks();
3248  if (numberElementBlocks==1) {
3249    loadProblem(*model,false);
3250    return dual();
3251  }
3252  // For now just get top level structure
3253  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3254  for (int i=0;i<numberElementBlocks;i++) {
3255    CoinStructuredModel * subModel = 
3256      dynamic_cast<CoinStructuredModel *>(model->block(i));
3257    CoinModel * thisBlock;
3258    if (subModel) {
3259      thisBlock=subModel->coinModelBlock(blockInfo[i]);
3260      model->setCoinModel(thisBlock,i);
3261    } else {
3262      thisBlock = dynamic_cast<CoinModel *>(model->block(i));
3263      assert (thisBlock);
3264      // just fill in info
3265      CoinModelBlockInfo info = CoinModelBlockInfo();
3266      int whatsSet = thisBlock->whatIsSet();
3267      info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
3268      info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
3269      info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
3270      info.integer = ((whatsSet&32)!=0) ? 1 : 0;
3271      info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
3272      info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
3273      // Which block
3274      int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
3275      info.rowBlock=iRowBlock;
3276      int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
3277      info.columnBlock=iColumnBlock;
3278      blockInfo[i] = info;
3279    }
3280  }
3281  int * rowCounts = new int [numberRowBlocks];
3282  CoinZeroN(rowCounts,numberRowBlocks);
3283  int * columnCounts = new int [numberColumnBlocks+1];
3284  CoinZeroN(columnCounts,numberColumnBlocks);
3285  int decomposeType=0;
3286  for (int i=0;i<numberElementBlocks;i++) {
3287    int iRowBlock = blockInfo[i].rowBlock;
3288    int iColumnBlock = blockInfo[i].columnBlock;
3289    rowCounts[iRowBlock]++;
3290    columnCounts[iColumnBlock]++;
3291  }
3292  if (numberRowBlocks==numberColumnBlocks||
3293      numberRowBlocks==numberColumnBlocks+1) {
3294    // could be Dantzig-Wolfe
3295    int numberG1=0;
3296    for (int i=0;i<numberRowBlocks;i++) {
3297      if (rowCounts[i]>1)
3298        numberG1++;
3299    }
3300    bool masterColumns = (numberColumnBlocks==numberRowBlocks);
3301    if ((masterColumns&&numberElementBlocks==2*numberRowBlocks-1)
3302        ||(!masterColumns&&numberElementBlocks==2*numberRowBlocks)) {
3303      if (numberG1<2)
3304        decomposeType=1;
3305    }
3306  }
3307  if (!decomposeType&&(numberRowBlocks==numberColumnBlocks||
3308                       numberRowBlocks==numberColumnBlocks-1)) {
3309    // could be Benders
3310    int numberG1=0;
3311    for (int i=0;i<numberColumnBlocks;i++) {
3312      if (columnCounts[i]>1)
3313        numberG1++;
3314    }
3315    bool masterRows = (numberColumnBlocks==numberRowBlocks);
3316    if ((masterRows&&numberElementBlocks==2*numberColumnBlocks-1)
3317        ||(!masterRows&&numberElementBlocks==2*numberColumnBlocks)) {
3318      if (numberG1<2)
3319        decomposeType=2;
3320    }
3321  }
3322  delete [] rowCounts;
3323  delete [] columnCounts;
3324  delete [] blockInfo;
3325  // decide what to do
3326  switch (decomposeType) {
3327    // No good
3328  case 0:
3329    loadProblem(*model,false);
3330    return dual();
3331    // DW
3332  case 1:
3333    return solveDW(model);
3334    // Benders
3335  case 2:
3336    return solveBenders(model);
3337  }
3338  return 0; // to stop compiler warning
3339}
3340/* This loads a model from a CoinStructuredModel object - returns number of errors.
3341   If originalOrder then keep to order stored in blocks,
3342   otherwise first column/rows correspond to first block - etc.
3343   If keepSolution true and size is same as current then
3344   keeps current status and solution
3345*/
3346int 
3347ClpSimplex::loadProblem (  CoinStructuredModel & coinModel,
3348                           bool originalOrder,
3349                           bool keepSolution)
3350{
3351  unsigned char * status = NULL;
3352  double * psol = NULL;
3353  double * dsol = NULL;
3354  int numberRows=coinModel.numberRows();
3355  int numberColumns = coinModel.numberColumns();
3356  int numberRowBlocks=coinModel.numberRowBlocks();
3357  int numberColumnBlocks = coinModel.numberColumnBlocks();
3358  int numberElementBlocks = coinModel.numberElementBlocks();
3359  if (status_&&numberRows_&&numberRows_==numberRows&&
3360      numberColumns_==numberColumns&&keepSolution) {
3361    status = new unsigned char [numberRows_+numberColumns_];
3362    CoinMemcpyN(status_,numberRows_+numberColumns_,status);
3363    psol = new double [numberRows_+numberColumns_];
3364    CoinMemcpyN(columnActivity_,numberColumns_,psol);
3365    CoinMemcpyN(rowActivity_,numberRows_,psol+numberColumns_);
3366    dsol = new double [numberRows_+numberColumns_];
3367    CoinMemcpyN(reducedCost_,numberColumns_,dsol);
3368    CoinMemcpyN(dual_,numberRows_,dsol+numberColumns_);
3369  }
3370  int returnCode=0;
3371  double * rowLower = new double [numberRows];
3372  double * rowUpper = new double [numberRows];
3373  double * columnLower = new double [numberColumns];
3374  double * columnUpper = new double [numberColumns];
3375  double * objective = new double [numberColumns];
3376  int * integerType = new int [numberColumns];
3377  CoinBigIndex numberElements=0;
3378  // Bases for blocks
3379  int * rowBase = new int[numberRowBlocks];
3380  CoinFillN(rowBase,numberRowBlocks,-1);
3381  // And row to put it
3382  int * whichRow = new int [numberRows];
3383  int * columnBase = new int[numberColumnBlocks];
3384  CoinFillN(columnBase,numberColumnBlocks,-1);
3385  // And column to put it
3386  int * whichColumn = new int [numberColumns];
3387  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
3388    CoinModel * block = coinModel.coinBlock(iBlock);
3389    numberElements += block->numberElements();
3390    //and set up elements etc
3391    double * associated = block->associatedArray();
3392    // If strings then do copies
3393    if (block->stringsExist()) 
3394      returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
3395                                          objective, integerType,associated);
3396    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3397    int iRowBlock = info.rowBlock;
3398    int iColumnBlock = info.columnBlock;
3399    if (rowBase[iRowBlock]<0) {
3400      rowBase[iRowBlock]=block->numberRows();
3401      // Save block number
3402      whichRow[numberRows-numberRowBlocks+iRowBlock]= iBlock;
3403    } else {
3404      assert(rowBase[iRowBlock]==block->numberRows());
3405    }
3406    if (columnBase[iColumnBlock]<0) {
3407      columnBase[iColumnBlock]=block->numberColumns();
3408      // Save block number
3409      whichColumn[numberColumns-numberColumnBlocks+iColumnBlock]= iBlock;
3410    } else {
3411      assert(columnBase[iColumnBlock]==block->numberColumns());
3412    }
3413  }
3414  // Fill arrays with defaults
3415  CoinFillN(rowLower,numberRows,-COIN_DBL_MAX);
3416  CoinFillN(rowUpper,numberRows,COIN_DBL_MAX);
3417  CoinFillN(columnLower,numberColumns,0.0);
3418  CoinFillN(columnUpper,numberColumns,COIN_DBL_MAX);
3419  CoinFillN(objective,numberColumns,0.0);
3420  CoinFillN(integerType,numberColumns,0);
3421  int n=0;
3422  for (int iBlock=0;iBlock<numberRowBlocks;iBlock++) {
3423    int k = rowBase[iBlock];
3424    rowBase[iBlock]=n;
3425    assert (k>=0);
3426    // block number
3427    int jBlock = whichRow[numberRows-numberRowBlocks+iBlock];
3428    if (originalOrder) {
3429      memcpy(whichRow+n,coinModel.coinBlock(jBlock)->originalRows(),k*sizeof(int));
3430    } else {
3431      CoinIotaN(whichRow+n,k,n);
3432    }
3433    n+=k;
3434  }
3435  assert (n==numberRows);
3436  n=0;
3437  for (int iBlock=0;iBlock<numberColumnBlocks;iBlock++) {
3438    int k = columnBase[iBlock];
3439    columnBase[iBlock]=n;
3440    assert (k>=0);
3441    if (k) {
3442      // block number
3443      int jBlock = whichColumn[numberColumns-numberColumnBlocks+iBlock];
3444      if (originalOrder) {
3445        memcpy(whichColumn+n,coinModel.coinBlock(jBlock)->originalColumns(),
3446               k*sizeof(int));
3447      } else {
3448        CoinIotaN(whichColumn+n,k,n);
3449      }
3450      n+=k;
3451    }
3452  }
3453  assert (n==numberColumns);
3454  bool gotIntegers=false;
3455  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
3456    CoinModel * block = coinModel.coinBlock(iBlock);
3457    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3458    int iRowBlock = info.rowBlock;
3459    int iRowBase = rowBase[iRowBlock];
3460    int iColumnBlock = info.columnBlock;
3461    int iColumnBase = columnBase[iColumnBlock];
3462    if (info.rhs) {
3463      int nRows = block->numberRows();
3464      const double * lower = block->rowLowerArray();
3465      const double * upper = block->rowUpperArray();
3466      for (int i=0;i<nRows;i++) {
3467        int put = whichRow[i+iRowBase];
3468        rowLower[put] = lower[i];
3469        rowUpper[put] = upper[i];
3470      }
3471    }
3472    if (info.bounds) {
3473      int nColumns = block->numberColumns();
3474      const double * lower = block->columnLowerArray();
3475      const double * upper = block->columnUpperArray();
3476      const double * obj = block->objectiveArray();
3477      for (int i=0;i<nColumns;i++) {
3478        int put = whichColumn[i+iColumnBase];
3479        columnLower[put] = lower[i];
3480        columnUpper[put] = upper[i];
3481        objective[put] = obj[i];
3482      }
3483    }
3484    if (info.integer) {
3485      gotIntegers=true;
3486      int nColumns = block->numberColumns();
3487      const int * type = block->integerTypeArray();
3488      for (int i=0;i<nColumns;i++) {
3489        int put = whichColumn[i+iColumnBase];
3490        integerType[put] = type[i];
3491      }
3492    }
3493  }
3494  gutsOfLoadModel(numberRows, numberColumns,
3495                  columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
3496  delete [] rowLower;
3497  delete [] rowUpper;
3498  delete [] columnLower;
3499  delete [] columnUpper;
3500  delete [] objective;
3501  // Do integers if wanted
3502  if (gotIntegers) {
3503    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
3504      if (integerType[iColumn])
3505        setInteger(iColumn);
3506    }
3507  }
3508  delete [] integerType;
3509  setObjectiveOffset(coinModel.objectiveOffset());
3510  // Space for elements
3511  int * row = new int[numberElements];
3512  int * column = new int[numberElements];
3513  double * element = new double[numberElements];
3514  numberElements=0;
3515  for (int iBlock=0;iBlock<numberElementBlocks;iBlock++) {
3516    CoinModel * block = coinModel.coinBlock(iBlock);
3517    const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3518    int iRowBlock = info.rowBlock;
3519    int iRowBase = rowBase[iRowBlock];
3520    int iColumnBlock = info.columnBlock;
3521    int iColumnBase = columnBase[iColumnBlock];
3522    if (info.rowName) {
3523      int numberItems = block->rowNames()->numberItems();
3524      assert( block->numberRows()>=numberItems);
3525      if (numberItems) {
3526        const char *const * rowNames=block->rowNames()->names();
3527        for (int i=0;i<numberItems;i++) {
3528          int put = whichRow[i+iRowBase];
3529          std::string name = rowNames[i];
3530          setRowName(put,name);
3531        }
3532      }
3533    }
3534    if (info.columnName) {
3535      int numberItems = block->columnNames()->numberItems();
3536      assert( block->numberColumns()>=numberItems);
3537      if (numberItems) {
3538        const char *const * columnNames=block->columnNames()->names();
3539        for (int i=0;i<numberItems;i++) {
3540          int put = whichColumn[i+iColumnBase];
3541          std::string name = columnNames[i];
3542          setColumnName(put,name);
3543        }
3544      }
3545    }
3546    if (info.matrix) {
3547      CoinPackedMatrix matrix2;
3548      const CoinPackedMatrix * matrix = block->packedMatrix();
3549      if (!matrix) {
3550        double * associated = block->associatedArray();
3551        block->createPackedMatrix(matrix2,associated);
3552        matrix = &matrix2;
3553      }
3554      // get matrix data pointers
3555      const int * row2 = matrix->getIndices();
3556      const CoinBigIndex * columnStart = matrix->getVectorStarts();
3557      const double * elementByColumn = matrix->getElements();
3558      const int * columnLength = matrix->getVectorLengths(); 
3559      int n = matrix->getNumCols();
3560      assert (matrix->isColOrdered());
3561      for (int iColumn=0;iColumn<n;iColumn++) {
3562        CoinBigIndex j;
3563        int jColumn = whichColumn[iColumn+iColumnBase];
3564        for (j=columnStart[iColumn];
3565             j<columnStart[iColumn]+columnLength[iColumn];j++) {
3566          row[numberElements]=whichRow[row2[j]+iRowBase];
3567          column[numberElements]=jColumn;
3568          element[numberElements++]=elementByColumn[j];
3569        }
3570      }
3571    }
3572  }
3573  delete [] whichRow;
3574  delete [] whichColumn;
3575  delete [] rowBase;
3576  delete [] columnBase;
3577  CoinPackedMatrix * matrix =
3578    new CoinPackedMatrix (true,row,column,element,numberElements);
3579  matrix_ = new ClpPackedMatrix(matrix);
3580  matrix_->setDimensions(numberRows,numberColumns);
3581  delete [] row;
3582  delete [] column;
3583  delete [] element;
3584  createStatus();
3585  if (status) {
3586    // copy back
3587    CoinMemcpyN(status,numberRows_+numberColumns_,status_);
3588    CoinMemcpyN(psol,numberColumns_,columnActivity_);
3589    CoinMemcpyN(psol+numberColumns_,numberRows_,rowActivity_);
3590    CoinMemcpyN(dsol,numberColumns_,reducedCost_);
3591    CoinMemcpyN(dsol+numberColumns_,numberRows_,dual_);
3592    delete [] status;
3593    delete [] psol;
3594    delete [] dsol;
3595  }
3596  optimizationDirection_ = coinModel.optimizationDirection(); 
3597  return returnCode;
3598}
3599/*  If input negative scales objective so maximum <= -value
3600    and returns scale factor used.  If positive unscales and also
3601    redoes dual stuff
3602*/
3603double 
3604ClpSimplex::scaleObjective(double value)
3605{
3606  double * obj = objective(); 
3607  double largest=0.0;
3608  if (value<0.0) {
3609    value = - value;
3610    for (int i=0;i<numberColumns_;i++) {
3611      largest = CoinMax(largest,fabs(obj[i]));
3612    }
3613    if (largest>value) {
3614      double scaleFactor = value/largest;
3615      for (int i=0;i<numberColumns_;i++) {
3616        obj[i] *= scaleFactor;
3617        reducedCost_[i] *= scaleFactor;
3618      }
3619      for (int i=0;i<numberRows_;i++) {
3620        dual_[i] *= scaleFactor;
3621      }
3622      largest /= value;
3623    } else {
3624      // no need
3625      largest=1.0;
3626    }
3627  } else {
3628    // at end
3629    if (value!=1.0) {
3630      for (int i=0;i<numberColumns_;i++) {
3631        obj[i] *= value;
3632        reducedCost_[i] *= value;
3633      }
3634      for (int i=0;i<numberRows_;i++) {
3635        dual_[i] *= value;
3636      }
3637      computeObjectiveValue();
3638    }
3639  }
3640  return largest;
3641}
3642// Solve using Dantzig-Wolfe decomposition and maybe in parallel
3643int 
3644ClpSimplex::solveDW(CoinStructuredModel * model)
3645{
3646  double time1 = CoinCpuTime();
3647  int numberColumns = model->numberColumns();
3648  int numberRowBlocks=model->numberRowBlocks();
3649  int numberColumnBlocks = model->numberColumnBlocks();
3650  int numberElementBlocks = model->numberElementBlocks();
3651  // We already have top level structure
3652  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3653  for (int i=0;i<numberElementBlocks;i++) {
3654    CoinModel * thisBlock = model->coinBlock(i);
3655    assert (thisBlock);
3656    // just fill in info
3657    CoinModelBlockInfo info = CoinModelBlockInfo();
3658    int whatsSet = thisBlock->whatIsSet();
3659    info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
3660    info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
3661    info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
3662    info.integer = ((whatsSet&32)!=0) ? 1 : 0;
3663    info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
3664    info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
3665    // Which block
3666    int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
3667    info.rowBlock=iRowBlock;
3668    int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
3669    info.columnBlock=iColumnBlock;
3670    blockInfo[i] = info;
3671  }
3672  // make up problems
3673  int numberBlocks=numberRowBlocks-1;
3674  // Find master rows and columns
3675  int * rowCounts = new int [numberRowBlocks];
3676  CoinZeroN(rowCounts,numberRowBlocks);
3677  int * columnCounts = new int [numberColumnBlocks+1];
3678  CoinZeroN(columnCounts,numberColumnBlocks);
3679  int iBlock;
3680  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
3681    int iRowBlock = blockInfo[iBlock].rowBlock;
3682    rowCounts[iRowBlock]++;
3683    int iColumnBlock =blockInfo[iBlock].columnBlock;
3684    columnCounts[iColumnBlock]++;
3685  }
3686  int * whichBlock = new int [numberElementBlocks];
3687  int masterRowBlock=-1;
3688  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
3689    if (rowCounts[iBlock]>1) {
3690      if (masterRowBlock==-1) {
3691        masterRowBlock=iBlock;
3692      } else {
3693        // Can't decode
3694        masterRowBlock=-2;
3695        break;
3696      }
3697    }
3698  }
3699  int masterColumnBlock=-1;
3700  int kBlock=0;
3701  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
3702    int count=columnCounts[iBlock];
3703    columnCounts[iBlock]=kBlock;
3704    kBlock += count;
3705  }
3706  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
3707    int iColumnBlock = blockInfo[iBlock].columnBlock;
3708    whichBlock[columnCounts[iColumnBlock]]=iBlock;
3709    columnCounts[iColumnBlock]++;
3710  }
3711  for (iBlock = numberColumnBlocks-1;iBlock>=0;iBlock--) 
3712    columnCounts[iBlock+1]=columnCounts[iBlock];
3713  columnCounts[0]=0;
3714  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
3715    int count=columnCounts[iBlock+1]-columnCounts[iBlock];
3716    if (count==1) {
3717      int kBlock = whichBlock[columnCounts[iBlock]];
3718      int iRowBlock = blockInfo[kBlock].rowBlock;
3719      if (iRowBlock==masterRowBlock) {
3720        if (masterColumnBlock==-1) {
3721          masterColumnBlock=iBlock;
3722        } else {
3723          // Can't decode
3724          masterColumnBlock=-2;
3725          break;
3726        }
3727      }
3728    }
3729  }
3730  if (masterRowBlock<0||masterColumnBlock==-2) {
3731    // What now
3732    abort();
3733  }
3734  delete [] rowCounts;
3735  // create all data
3736  const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
3737  ClpSimplex * sub = new ClpSimplex [numberBlocks];
3738  ClpSimplex master;
3739  // Set offset
3740  master.setObjectiveOffset(model->objectiveOffset());
3741  kBlock=0;
3742  int masterBlock=-1;
3743  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
3744    top[kBlock]=NULL;
3745    int start=columnCounts[iBlock];
3746    int end=columnCounts[iBlock+1];
3747    assert (end-start<=2);
3748    for (int j=start;j<end;j++) {
3749      int jBlock = whichBlock[j];
3750      int iRowBlock = blockInfo[jBlock].rowBlock;
3751      int iColumnBlock =blockInfo[jBlock].columnBlock;
3752      assert (iColumnBlock==iBlock);
3753      if (iColumnBlock!=masterColumnBlock&&iRowBlock==masterRowBlock) {
3754        // top matrix
3755        top[kBlock]=model->coinBlock(jBlock)->packedMatrix();
3756      } else {
3757        const CoinPackedMatrix * matrix
3758          = model->coinBlock(jBlock)->packedMatrix();
3759        // Get pointers to arrays
3760        const double * rowLower;
3761        const double * rowUpper;
3762        const double * columnLower;
3763        const double * columnUpper;
3764        const double * objective;
3765        model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
3766                     columnLower,columnUpper,objective);
3767        if (iColumnBlock!=masterColumnBlock) {
3768          // diagonal block
3769          sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
3770                                  objective,rowLower,rowUpper);
3771          if (true) {
3772            double * lower = sub[kBlock].columnLower();
3773            double * upper = sub[kBlock].columnUpper();
3774            int n = sub[kBlock].numberColumns();
3775            for (int i=0;i<n;i++) {
3776              lower[i]=CoinMax(-1.0e8,lower[i]);
3777              upper[i]=CoinMin(1.0e8,upper[i]);
3778            }
3779          }
3780          if (optimizationDirection_<0.0) {
3781            double * obj = sub[kBlock].objective();
3782            int n = sub[kBlock].numberColumns();
3783            for (int i=0;i<n;i++)
3784              obj[i] = - obj[i];
3785          }
3786          if (this->factorizationFrequency()==200) {
3787            // User did not touch preset
3788            sub[kBlock].defaultFactorizationFrequency();
3789          } else {
3790            // make sure model has correct value
3791            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
3792          }
3793          sub[kBlock].setPerturbation(50);
3794          // Set columnCounts to be diagonal block index for cleanup
3795          columnCounts[kBlock]=jBlock;
3796        } else {
3797          // master
3798          masterBlock=jBlock;
3799          master.loadProblem(*matrix,columnLower,columnUpper,
3800                             objective,rowLower,rowUpper);
3801          if (optimizationDirection_<0.0) {
3802            double * obj = master.objective();
3803            int n = master.numberColumns();
3804            for (int i=0;i<n;i++)
3805              obj[i] = - obj[i];
3806          }
3807        }
3808      }
3809    }
3810    if (iBlock!=masterColumnBlock) 
3811      kBlock++;
3812  }
3813  delete [] whichBlock;
3814  delete [] blockInfo;
3815  // For now master must have been defined (does not have to have columns)
3816  assert (master.numberRows());
3817  assert (masterBlock>=0);
3818  int numberMasterRows = master.numberRows();
3819  // Overkill in terms of space
3820  int spaceNeeded = CoinMax(numberBlocks*(numberMasterRows+1),
3821                            2*numberMasterRows);
3822  int * rowAdd = new int[spaceNeeded];
3823  double * elementAdd = new double[spaceNeeded];
3824  spaceNeeded = numberBlocks;
3825  int * columnAdd = new int[spaceNeeded+1];
3826  double * objective = new double[spaceNeeded];
3827  // Add in costed slacks
3828  int firstArtificial = master.numberColumns();
3829  int lastArtificial = firstArtificial;
3830  if (true) {
3831    const double * lower = master.rowLower();
3832    const double * upper = master.rowUpper();
3833    int kCol=0;
3834    for (int iRow=0;iRow<numberMasterRows;iRow++) {
3835      if (lower[iRow]>-1.0e10) {
3836        rowAdd[kCol]=iRow;
3837        elementAdd[kCol++]=1.0;
3838      }
3839      if (upper[iRow]<1.0e10) {
3840        rowAdd[kCol]=iRow;
3841        elementAdd[kCol++]=-1.0;
3842      }
3843    }
3844    if (kCol>spaceNeeded) {
3845      spaceNeeded = kCol;
3846      delete [] columnAdd;
3847      delete [] objective;
3848      columnAdd = new int[spaceNeeded+1];
3849      objective = new double[spaceNeeded];
3850    }
3851    for (int i=0;i<kCol;i++) {
3852      columnAdd[i]=i;
3853      objective[i]=1.0e13;
3854    }
3855    columnAdd[kCol]=kCol;
3856    master.addColumns(kCol,NULL,NULL,objective,
3857                      columnAdd,rowAdd,elementAdd);
3858    lastArtificial = master.numberColumns();
3859  }
3860  int maxPass=500;
3861  int iPass;
3862  double lastObjective=1.0e31;
3863  // Create convexity rows for proposals
3864  int numberMasterColumns = master.numberColumns();
3865  master.resize(numberMasterRows+numberBlocks,numberMasterColumns);
3866  if (this->factorizationFrequency()==200) {
3867    // User did not touch preset
3868    master.defaultFactorizationFrequency();
3869  } else {
3870    // make sure model has correct value
3871    master.setFactorizationFrequency(this->factorizationFrequency());
3872  }
3873  master.setPerturbation(50);
3874  // Arrays to say which block and when created
3875  int maximumColumns = 2*numberMasterRows+10*numberBlocks;
3876  whichBlock = new int[maximumColumns];
3877  int * when = new int[maximumColumns];
3878  int numberColumnsGenerated=numberBlocks;
3879  // fill in rhs and add in artificials
3880  {
3881    double * rowLower = master.rowLower();
3882    double * rowUpper = master.rowUpper();
3883    int iBlock;
3884    columnAdd[0] = 0;
3885    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
3886      int iRow = iBlock + numberMasterRows;;
3887      rowLower[iRow]=1.0;
3888      rowUpper[iRow]=1.0;
3889      rowAdd[iBlock] = iRow;
3890      elementAdd[iBlock] = 1.0;
3891      objective[iBlock] = 1.0e13;
3892      columnAdd[iBlock+1] = iBlock+1;
3893      when[iBlock]=-1;
3894      whichBlock[iBlock] = iBlock;
3895    }
3896    master.addColumns(numberBlocks,NULL,NULL,objective,
3897                      columnAdd, rowAdd, elementAdd);
3898  }
3899  // and resize matrix to double check clp will be happy
3900  //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
3901  //                     numberMasterColumns+numberBlocks);
3902  std::cout<<"Time to decompose "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
3903  for (iPass=0;iPass<maxPass;iPass++) {
3904    printf("Start of pass %d\n",iPass);
3905    // Solve master - may be infeasible
3906    //master.scaling(0);
3907    if (0) {
3908      master.writeMps("yy.mps");
3909    }
3910    // Correct artificials
3911    double sumArtificials=0.0;
3912    if (iPass) {
3913      double * upper = master.columnUpper();
3914      double * solution = master.primalColumnSolution();
3915      double * obj = master.objective();
3916      sumArtificials=0.0;
3917      for (int i=firstArtificial;i<lastArtificial;i++) {
3918        sumArtificials += solution[i];
3919        //assert (solution[i]>-1.0e-2);
3920        if (solution[i]<1.0e-6) {
3921#if 0
3922          // Could take out
3923          obj[i]=0.0;
3924          upper[i]=0.0;
3925#else
3926          obj[i]=1.0e7;
3927          upper[i]=1.0e-1;
3928#endif
3929          solution[i]=0.0;
3930          master.setColumnStatus(i,isFixed);
3931        } else {
3932          upper[i]=solution[i]+1.0e-5*(1.0+solution[i]);
3933        }
3934      }
3935      printf("Sum of artificials before solve is %g\n",sumArtificials);
3936    }
3937    // scale objective to be reasonable
3938    double scaleFactor = master.scaleObjective(-1.0e9);
3939    {
3940      double * dual = master.dualRowSolution();
3941      int n=master.numberRows();
3942      memset(dual,0,n*sizeof(double));
3943      double * solution = master.primalColumnSolution();
3944      master.clpMatrix()->times(1.0,solution,dual);
3945      double sum=0.0;
3946      double * lower = master.rowLower();
3947      double * upper = master.rowUpper();
3948      for (int iRow=0;iRow<n;iRow++) {
3949        double value = dual[iRow];
3950        if (value>upper[iRow])
3951          sum += value-upper[iRow];
3952        else if (value<lower[iRow])
3953          sum -= value-lower[iRow];
3954      }
3955      printf("suminf %g\n",sum);
3956      lower = master.columnLower();
3957      upper = master.columnUpper();
3958      n=master.numberColumns();
3959      for (int iColumn=0;iColumn<n;iColumn++) {
3960        double value = solution[iColumn];
3961        if (value>upper[iColumn]+1.0e-5)
3962          sum += value-upper[iColumn];
3963        else if (value<lower[iColumn]-1.0e-5)
3964          sum -= value-lower[iColumn];
3965      }
3966      printf("suminf %g\n",sum);
3967    }
3968    master.primal(1);
3969    // Correct artificials
3970    sumArtificials=0.0;
3971    {
3972      double * solution = master.primalColumnSolution();
3973      for (int i=firstArtificial;i<lastArtificial;i++) {
3974        sumArtificials += solution[i];
3975      }
3976      printf("Sum of artificials after solve is %g\n",sumArtificials);
3977    }
3978    master.scaleObjective(scaleFactor);
3979    int problemStatus = master.status(); // do here as can change (delcols)
3980    if (master.numberIterations()==0&&iPass)
3981      break; // finished
3982    if (master.objectiveValue()>lastObjective-1.0e-7&&iPass>555)
3983      break; // finished
3984    lastObjective = master.objectiveValue();
3985    // mark basic ones and delete if necessary
3986    int iColumn;
3987    numberColumnsGenerated=master.numberColumns()-numberMasterColumns;
3988    for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
3989      if (master.getStatus(iColumn+numberMasterColumns)==ClpSimplex::basic)
3990        when[iColumn]=iPass;
3991    }
3992    if (numberColumnsGenerated+numberBlocks>maximumColumns) {
3993      // delete
3994      int numberKeep=0;
3995      int numberDelete=0;
3996      int * whichDelete = new int[numberColumnsGenerated];
3997      for (iColumn=0;iColumn<numberColumnsGenerated;iColumn++) {
3998        if (when[iColumn]>iPass-7) {
3999          // keep
4000          when[numberKeep]=when[iColumn];
4001          whichBlock[numberKeep++]=whichBlock[iColumn];
4002        } else {
4003          // delete
4004          whichDelete[numberDelete++]=iColumn+numberMasterColumns;
4005        }
4006      }
4007      numberColumnsGenerated -= numberDelete;
4008      master.deleteColumns(numberDelete,whichDelete);
4009      delete [] whichDelete;
4010    }
4011    const double * dual=NULL;
4012    bool deleteDual=false;
4013    if (problemStatus==0) {
4014      dual = master.dualRowSolution();
4015    } else if (problemStatus==1) {
4016      // could do composite objective
4017      dual = master.infeasibilityRay();
4018      deleteDual = true;
4019      printf("The sum of infeasibilities is %g\n",
4020             master.sumPrimalInfeasibilities());
4021    } else if (!master.numberColumns()) {
4022      assert(!iPass);
4023      dual = master.dualRowSolution();
4024      memset(master.dualRowSolution(),
4025             0,(numberMasterRows+numberBlocks)*sizeof(double));
4026    } else {
4027      abort();
4028    }
4029    // Scale back on first time
4030    if (!iPass) {
4031      double * dual2 = master.dualRowSolution();
4032      for (int i=0;i<numberMasterRows+numberBlocks;i++) {
4033        dual2[i] *= 1.0e-7;
4034      }
4035      dual = master.dualRowSolution();
4036    }
4037    // Create objective for sub problems and solve
4038    columnAdd[0]=0;
4039    int numberProposals=0;
4040    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4041      int numberColumns2 = sub[iBlock].numberColumns();
4042      double * saveObj = new double [numberColumns2];
4043      double * objective2 = sub[iBlock].objective();
4044      memcpy(saveObj,objective2,numberColumns2*sizeof(double));
4045      // new objective
4046      top[iBlock]->transposeTimes(dual,objective2);
4047      int i;
4048      if (problemStatus==0) {
4049        for (i=0;i<numberColumns2;i++)
4050          objective2[i] = saveObj[i]-objective2[i];
4051      } else {
4052        for (i=0;i<numberColumns2;i++)
4053          objective2[i] = -objective2[i];
4054      }
4055      // scale objective to be reasonable
4056      double scaleFactor = 
4057        sub[iBlock].scaleObjective((sumArtificials>1.0e-5) ? -1.0e-4 :-1.0e9);
4058      if (iPass) {
4059        sub[iBlock].primal();
4060      } else {
4061        sub[iBlock].dual();
4062      }
4063      sub[iBlock].scaleObjective(scaleFactor);
4064      if (!sub[iBlock].isProvenOptimal()&& 
4065          !sub[iBlock].isProvenDualInfeasible()) {
4066        memset(objective2,0,numberColumns2*sizeof(double));
4067        sub[iBlock].primal();
4068        if (problemStatus==0) {
4069          for (i=0;i<numberColumns2;i++)
4070            objective2[i] = saveObj[i]-objective2[i];
4071        } else {
4072          for (i=0;i<numberColumns2;i++)
4073            objective2[i] = -objective2[i];
4074        }
4075        double scaleFactor = sub[iBlock].scaleObjective(-1.0e9);
4076        sub[iBlock].primal(1);
4077        sub[iBlock].scaleObjective(scaleFactor);
4078      }
4079      memcpy(objective2,saveObj,numberColumns2*sizeof(double));
4080      // get proposal
4081      if (sub[iBlock].numberIterations()||!iPass) {
4082        double objValue =0.0;
4083        int start = columnAdd[numberProposals];
4084        // proposal
4085        if (sub[iBlock].isProvenOptimal()) {
4086          const double * solution = sub[iBlock].primalColumnSolution();
4087          top[iBlock]->times(solution,elementAdd+start);
4088          for (i=0;i<numberColumns2;i++)
4089            objValue += solution[i]*saveObj[i];
4090          // See if good dj and pack down
4091          int number=start;
4092          double dj = objValue;
4093          if (problemStatus) 
4094            dj=0.0;
4095          double smallest=1.0e100;
4096          double largest=0.0;
4097          for (i=0;i<numberMasterRows;i++) {
4098            double value = elementAdd[start+i];
4099            if (fabs(value)>1.0e-15) {
4100              dj -= dual[i]*value;
4101              smallest = CoinMin(smallest,fabs(value));
4102              largest = CoinMax(largest,fabs(value));
4103              rowAdd[number]=i;
4104              elementAdd[number++]=value;
4105            }
4106          }
4107          // and convexity
4108          dj -= dual[numberMasterRows+iBlock];
4109          rowAdd[number]=numberMasterRows+iBlock;
4110          elementAdd[number++]=1.0;
4111          // if elements large then scale?
4112          //if (largest>1.0e8||smallest<1.0e-8)
4113          printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
4114                 iBlock,smallest,largest,dj);
4115          if (dj<-1.0e-6||!iPass) {
4116            // take
4117            objective[numberProposals]=objValue;
4118            columnAdd[++numberProposals]=number;
4119            when[numberColumnsGenerated]=iPass;
4120            whichBlock[numberColumnsGenerated++]=iBlock;
4121          }
4122        } else if (sub[iBlock].isProvenDualInfeasible()) {
4123          // use ray
4124          const double * solution = sub[iBlock].unboundedRay();
4125          top[iBlock]->times(solution,elementAdd+start);
4126          for (i=0;i<numberColumns2;i++)
4127            objValue += solution[i]*saveObj[i];
4128          // See if good dj and pack down
4129          int number=start;
4130          double dj = objValue;
4131          double smallest=1.0e100;
4132          double largest=0.0;
4133          for (i=0;i<numberMasterRows;i++) {
4134            double value = elementAdd[start+i];
4135            if (fabs(value)>1.0e-15) {
4136              dj -= dual[i]*value;
4137              smallest = CoinMin(smallest,fabs(value));
4138              largest = CoinMax(largest,fabs(value));
4139              rowAdd[number]=i;
4140              elementAdd[number++]=value;
4141            }
4142          }
4143          // if elements large or small then scale?
4144          //if (largest>1.0e8||smallest<1.0e-8)
4145          printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
4146                 iBlock,smallest,largest,dj);
4147          if (dj<-1.0e-6) {
4148            // take
4149            objective[numberProposals]=objValue;
4150            columnAdd[++numberProposals]=number;
4151            when[numberColumnsGenerated]=iPass;
4152            whichBlock[numberColumnsGenerated++]=iBlock;
4153          }
4154        } else {
4155          abort();
4156        }
4157      }
4158      delete [] saveObj;
4159    }
4160    if (deleteDual)
4161      delete [] dual;
4162    if (numberProposals) 
4163      master.addColumns(numberProposals,NULL,NULL,objective,
4164                        columnAdd,rowAdd,elementAdd);
4165  }
4166  std::cout<<"Time at end of D-W "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4167  //master.scaling(0);
4168  //master.primal(1);
4169  loadProblem(*model);
4170  // now put back a good solution
4171  double * lower = new double[numberMasterRows+numberBlocks];
4172  double * upper = new double[numberMasterRows+numberBlocks];
4173  numberColumnsGenerated  += numberMasterColumns;
4174  double * sol = new double[numberColumnsGenerated];
4175  const double * solution = master.primalColumnSolution();
4176  const double * masterLower = master.rowLower();
4177  const double * masterUpper = master.rowUpper();
4178  double * fullSolution = primalColumnSolution();
4179  const double * fullLower = columnLower();
4180  const double * fullUpper = columnUpper();
4181  const double * rowSolution = master.primalRowSolution();
4182  double * fullRowSolution = primalRowSolution();
4183  const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4184  int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4185  const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4186  int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4187  for (int iRow=0;iRow<numberRows2;iRow++) {
4188    int kRow = rowBack[iRow];
4189    setRowStatus(kRow,master.getRowStatus(iRow));
4190    fullRowSolution[kRow]=rowSolution[iRow];
4191  }
4192  for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4193    int kColumn = columnBack[iColumn];
4194    setStatus(kColumn,master.getStatus(iColumn));
4195    fullSolution[kColumn]=solution[iColumn];
4196  }
4197  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4198    // move basis
4199    int kBlock = columnCounts[iBlock];
4200    const int * rowBack = model->coinBlock(kBlock)->originalRows();
4201    int numberRows2 = model->coinBlock(kBlock)->numberRows();
4202    const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4203    int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4204    for (int iRow=0;iRow<numberRows2;iRow++) {
4205      int kRow = rowBack[iRow];
4206      setRowStatus(kRow,sub[iBlock].getRowStatus(iRow));
4207    }
4208    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4209      int kColumn = columnBack[iColumn];
4210      setStatus(kColumn,sub[iBlock].getStatus(iColumn));
4211    }
4212    // convert top bit to by rows
4213    CoinPackedMatrix topMatrix = *top[iBlock];
4214    topMatrix.reverseOrdering();
4215    // zero solution
4216    memset(sol,0,numberColumnsGenerated*sizeof(double));
4217   
4218    for (int i=numberMasterColumns;i<numberColumnsGenerated;i++) {
4219      if (whichBlock[i-numberMasterColumns]==iBlock)
4220        sol[i] = solution[i];
4221    }
4222    memset(lower,0,(numberMasterRows+numberBlocks)*sizeof(double));
4223    master.clpMatrix()->times(1.0,sol,lower);
4224    for (int iRow=0;iRow<numberMasterRows;iRow++) {
4225      double value = lower[iRow];
4226      if (masterUpper[iRow]<1.0e20) 
4227        upper[iRow] = value;
4228      else
4229        upper[iRow]=COIN_DBL_MAX;
4230      if (masterLower[iRow]>-1.0e20) 
4231        lower[iRow] = value;
4232      else
4233        lower[iRow]=-COIN_DBL_MAX;
4234    }
4235    sub[iBlock].addRows(numberMasterRows,lower,upper,
4236                        topMatrix.getVectorStarts(),
4237                        topMatrix.getVectorLengths(),
4238                        topMatrix.getIndices(),
4239                        topMatrix.getElements());
4240    sub[iBlock].primal(1);
4241    const double * subSolution = sub[iBlock].primalColumnSolution();
4242    const double * subRowSolution = sub[iBlock].primalRowSolution();
4243    // move solution
4244    for (int iRow=0;iRow<numberRows2;iRow++) {
4245      int kRow = rowBack[iRow];
4246      fullRowSolution[kRow]=subRowSolution[iRow];
4247    }
4248    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4249      int kColumn = columnBack[iColumn];
4250      fullSolution[kColumn]=subSolution[iColumn];
4251    }
4252  }
4253  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
4254    if (fullSolution[iColumn]<fullUpper[iColumn]-1.0e-8&&
4255        fullSolution[iColumn]>fullLower[iColumn]+1.0e-8) {
4256      if (getStatus(iColumn)!=ClpSimplex::basic) {
4257        if (columnLower_[iColumn]>-1.0e30||
4258            columnUpper_[iColumn]<1.0e30)
4259          setStatus(iColumn,ClpSimplex::superBasic);
4260        else
4261          setStatus(iColumn,ClpSimplex::isFree);
4262      }
4263    } else if (fullSolution[iColumn]>=fullUpper[iColumn]-1.0e-8) {
4264      // may help to make rest non basic
4265      if (getStatus(iColumn)!=ClpSimplex::basic) 
4266        setStatus(iColumn,ClpSimplex::atUpperBound);
4267    } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e-8) {
4268      // may help to make rest non basic
4269      if (getStatus(iColumn)!=ClpSimplex::basic) 
4270        setStatus(iColumn,ClpSimplex::atLowerBound);
4271    }
4272  }
4273  //int numberRows=model->numberRows();
4274  //for (int iRow=0;iRow<numberRows;iRow++)
4275  //setRowStatus(iRow,ClpSimplex::superBasic);
4276  std::cout<<"Time before cleanup of full model "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4277  primal(1);
4278  std::cout<<"Total time "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4279  delete [] columnCounts;
4280  delete [] sol;
4281  delete [] lower;
4282  delete [] upper;
4283  delete [] whichBlock;
4284  delete [] when;
4285  delete [] columnAdd;
4286  delete [] rowAdd;
4287  delete [] elementAdd;
4288  delete [] objective;
4289  delete [] top;
4290  delete [] sub;
4291  return 0;
4292}
4293// Solve using Benders decomposition and maybe in parallel
4294int 
4295ClpSimplex::solveBenders(CoinStructuredModel * model)
4296{
4297  double time1 = CoinCpuTime();
4298  //int numberColumns = model->numberColumns();
4299  int numberRowBlocks=model->numberRowBlocks();
4300  int numberColumnBlocks = model->numberColumnBlocks();
4301  int numberElementBlocks = model->numberElementBlocks();
4302  // We already have top level structure
4303  CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4304  for (int i=0;i<numberElementBlocks;i++) {
4305    CoinModel * thisBlock = model->coinBlock(i);
4306    assert (thisBlock);
4307    // just fill in info
4308    CoinModelBlockInfo info = CoinModelBlockInfo();
4309    int whatsSet = thisBlock->whatIsSet();
4310    info.matrix = ((whatsSet&1)!=0) ? 1 : 0;
4311    info.rhs = ((whatsSet&2)!=0) ? 1 : 0;
4312    info.rowName = ((whatsSet&4)!=0) ? 1 : 0;
4313    info.integer = ((whatsSet&32)!=0) ? 1 : 0;
4314    info.bounds = ((whatsSet&8)!=0) ? 1 : 0;
4315    info.columnName = ((whatsSet&16)!=0) ? 1 : 0;
4316    // Which block
4317    int iRowBlock=model->rowBlock(thisBlock->getRowBlock());
4318    info.rowBlock=iRowBlock;
4319    int iColumnBlock=model->columnBlock(thisBlock->getColumnBlock());
4320    info.columnBlock=iColumnBlock;
4321    blockInfo[i] = info;
4322  }
4323  // make up problems
4324  int numberBlocks=numberColumnBlocks-1;
4325  // Find master columns and rows
4326  int * columnCounts = new int [numberColumnBlocks];
4327  CoinZeroN(columnCounts,numberColumnBlocks);
4328  int * rowCounts = new int [numberRowBlocks+1];
4329  CoinZeroN(rowCounts,numberRowBlocks);
4330  int iBlock;
4331  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
4332    int iColumnBlock = blockInfo[iBlock].columnBlock;
4333    columnCounts[iColumnBlock]++;
4334    int iRowBlock =blockInfo[iBlock].rowBlock;
4335    rowCounts[iRowBlock]++;
4336  }
4337  int * whichBlock = new int [numberElementBlocks];
4338  int masterColumnBlock=-1;
4339  for (iBlock = 0;iBlock<numberColumnBlocks;iBlock++) {
4340    if (columnCounts[iBlock]>1) {
4341      if (masterColumnBlock==-1) {
4342        masterColumnBlock=iBlock;
4343      } else {
4344        // Can't decode
4345        masterColumnBlock=-2;
4346        break;
4347      }
4348    }
4349  }
4350  int masterRowBlock=-1;
4351  int kBlock=0;
4352  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
4353    int count=rowCounts[iBlock];
4354    rowCounts[iBlock]=kBlock;
4355    kBlock += count;
4356  }
4357  for (iBlock = 0;iBlock<numberElementBlocks;iBlock++) {
4358    int iRowBlock = blockInfo[iBlock].rowBlock;
4359    whichBlock[rowCounts[iRowBlock]]=iBlock;
4360    rowCounts[iRowBlock]++;
4361  }
4362  for (iBlock = numberRowBlocks-1;iBlock>=0;iBlock--) 
4363    rowCounts[iBlock+1]=rowCounts[iBlock];
4364  rowCounts[0]=0;
4365  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
4366    int count=rowCounts[iBlock+1]-rowCounts[iBlock];
4367    if (count==1) {
4368      int kBlock = whichBlock[rowCounts[iBlock]];
4369      int iColumnBlock = blockInfo[kBlock].columnBlock;
4370      if (iColumnBlock==masterColumnBlock) {
4371        if (masterRowBlock==-1) {
4372          masterRowBlock=iBlock;
4373        } else {
4374          // Can't decode
4375          masterRowBlock=-2;
4376          break;
4377        }
4378      }
4379    }
4380  }
4381  if (masterColumnBlock<0||masterRowBlock==-2) {
4382    // What now
4383    abort();
4384  }
4385  delete [] columnCounts;
4386  // create all data
4387  const CoinPackedMatrix ** first = new const CoinPackedMatrix * [numberRowBlocks];
4388  ClpSimplex * sub = new ClpSimplex [numberBlocks];
4389  ClpSimplex master;
4390  // Set offset
4391  master.setObjectiveOffset(model->objectiveOffset());
4392  kBlock=0;
4393  int masterBlock=-1;
4394  for (iBlock = 0;iBlock<numberRowBlocks;iBlock++) {
4395    first[kBlock]=NULL;
4396    int start=rowCounts[iBlock];
4397    int end=rowCounts[iBlock+1];
4398    assert (end-start<=2);
4399    for (int j=start;j<end;j++) {
4400      int jBlock = whichBlock[j];
4401      int iColumnBlock = blockInfo[jBlock].columnBlock;
4402      int iRowBlock =blockInfo[jBlock].rowBlock;
4403      assert (iRowBlock==iBlock);
4404      if (iRowBlock!=masterRowBlock&&iColumnBlock==masterColumnBlock) {
4405        // first matrix
4406        first[kBlock]=model->coinBlock(jBlock)->packedMatrix();
4407      } else {
4408        const CoinPackedMatrix * matrix
4409          = model->coinBlock(jBlock)->packedMatrix();
4410        // Get pointers to arrays
4411        const double * columnLower;
4412        const double * columnUpper;
4413        const double * rowLower;
4414        const double * rowUpper;
4415        const double * objective;
4416        model->block(iRowBlock,iColumnBlock,rowLower,rowUpper,
4417                     columnLower,columnUpper,objective);
4418        if (iRowBlock!=masterRowBlock) {
4419          // diagonal block
4420          sub[kBlock].loadProblem(*matrix,columnLower,columnUpper,
4421                                  objective,rowLower,rowUpper);
4422          if (optimizationDirection_<0.0) {
4423            double * obj = sub[kBlock].objective();
4424            int n = sub[kBlock].numberColumns();
4425            for (int i=0;i<n;i++)
4426              obj[i] = - obj[i];
4427          }
4428          if (this->factorizationFrequency()==200) {
4429            // User did not touch preset
4430            sub[kBlock].defaultFactorizationFrequency();
4431          } else {
4432            // make sure model has correct value
4433            sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
4434          }
4435          sub[kBlock].setPerturbation(50);
4436          // Set rowCounts to be diagonal block index for cleanup
4437          rowCounts[kBlock]=jBlock;
4438        } else {
4439          // master
4440          masterBlock=jBlock;
4441          master.loadProblem(*matrix,columnLower,columnUpper,
4442                             objective,rowLower,rowUpper);
4443          if (optimizationDirection_<0.0) {
4444            double * obj = master.objective();
4445            int n = master.numberColumns();
4446            for (int i=0;i<n;i++)
4447              obj[i] = - obj[i];
4448          }
4449        }
4450      }
4451    }
4452    if (iBlock!=masterRowBlock) 
4453      kBlock++;
4454  }
4455  delete [] whichBlock;
4456  delete [] blockInfo;
4457  // For now master must have been defined (does not have to have rows)
4458  assert (master.numberColumns());
4459  assert (masterBlock>=0);
4460  int numberMasterColumns = master.numberColumns();
4461  // Overkill in terms of space
4462  int spaceNeeded = CoinMax(numberBlocks*(numberMasterColumns+1),
4463                            2*numberMasterColumns);
4464  int * columnAdd = new int[spaceNeeded];
4465  double * elementAdd = new double[spaceNeeded];
4466  spaceNeeded = numberBlocks;
4467  int * rowAdd = new int[spaceNeeded+1];
4468  double * objective = new double[spaceNeeded];
4469  int maxPass=500;
4470  int iPass;
4471  double lastObjective=-1.0e31;
4472  // Create columns for proposals
4473  int numberMasterRows = master.numberRows();
4474  master.resize(numberMasterColumns+numberBlocks,numberMasterRows);
4475  if (this->factorizationFrequency()==200) {
4476    // User did not touch preset
4477    master.defaultFactorizationFrequency();
4478  } else {
4479    // make sure model has correct value
4480    master.setFactorizationFrequency(this->factorizationFrequency());
4481  }
4482  master.setPerturbation(50);
4483  // Arrays to say which block and when created
4484  int maximumRows = 2*numberMasterColumns+10*numberBlocks;
4485  whichBlock = new int[maximumRows];
4486  int * when = new int[maximumRows];
4487  int numberRowsGenerated=numberBlocks;
4488  // Add extra variables
4489  {
4490    int iBlock;
4491    columnAdd[0] = 0;
4492    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4493      objective[iBlock] = 1.0;
4494      columnAdd[iBlock+1] = 0;
4495      when[iBlock]=-1;
4496      whichBlock[iBlock] = iBlock;
4497    }
4498    master.addColumns(numberBlocks,NULL,NULL,objective,
4499                      columnAdd, rowAdd, elementAdd);
4500  }
4501  std::cout<<"Time to decompose "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4502  for (iPass=0;iPass<maxPass;iPass++) {
4503    printf("Start of pass %d\n",iPass);
4504    // Solve master - may be unbounded
4505    //master.scaling(0);
4506    if (1) {
4507      master.writeMps("yy.mps");
4508    }
4509    master.dual();
4510    int problemStatus = master.status(); // do here as can change (delcols)
4511    if (master.numberIterations()==0&&iPass)
4512      break; // finished
4513    if (master.objectiveValue()<lastObjective+1.0e-7&&iPass>555)
4514      break; // finished
4515    lastObjective = master.objectiveValue();
4516    // mark non-basic rows and delete if necessary
4517    int iRow;
4518    numberRowsGenerated=master.numberRows()-numberMasterRows;
4519    for (iRow=0;iRow<numberRowsGenerated;iRow++) {
4520      if (master.getStatus(iRow+numberMasterRows)!=ClpSimplex::basic)
4521        when[iRow]=iPass;
4522    }
4523    if (numberRowsGenerated>maximumRows) {
4524      // delete
4525      int numberKeep=0;
4526      int numberDelete=0;
4527      int * whichDelete = new int[numberRowsGenerated];
4528      for (iRow=0;iRow<numberRowsGenerated;iRow++) {
4529        if (when[iRow]>iPass-7) {
4530          // keep
4531          when[numberKeep]=when[iRow];
4532          whichBlock[numberKeep++]=whichBlock[iRow];
4533        } else {
4534          // delete
4535          whichDelete[numberDelete++]=iRow+numberMasterRows;
4536        }
4537      }
4538      numberRowsGenerated -= numberDelete;
4539      master.deleteRows(numberDelete,whichDelete);
4540      delete [] whichDelete;
4541    }
4542    const double * primal=NULL;
4543    bool deletePrimal=false;
4544    if (problemStatus==0) {
4545      primal = master.primalColumnSolution();
4546    } else if (problemStatus==2) {
4547      // could do composite objective
4548      primal = master.unboundedRay();
4549      deletePrimal = true;
4550      printf("The sum of infeasibilities is %g\n",
4551             master.sumPrimalInfeasibilities());
4552    } else if (!master.numberRows()) {
4553      assert(!iPass);
4554      primal = master.primalColumnSolution();
4555      memset(master.primalColumnSolution(),
4556             0,numberMasterColumns*sizeof(double));
4557    } else {
4558      abort();
4559    }
4560    // Create rhs for sub problems and solve
4561    rowAdd[0]=0;
4562    int numberProposals=0;
4563    for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4564      int numberRows2 = sub[iBlock].numberRows();
4565      double * saveLower = new double [numberRows2];
4566      double * lower2 = sub[iBlock].rowLower();
4567      double * saveUpper = new double [numberRows2];
4568      double * upper2 = sub[iBlock].rowUpper();
4569      // new rhs
4570      CoinZeroN(saveUpper,numberRows2);
4571      first[iBlock]->times(primal,saveUpper);
4572      for (int i=0;i<numberRows2;i++) {
4573        double value = saveUpper[i];
4574        saveLower[i]=lower2[i];
4575        saveUpper[i]=upper2[i];
4576        if (saveLower[i]>-1.0e30)
4577          lower2[i] -= value;
4578        if (saveUpper[i]<1.0e30)
4579          upper2[i] -= value;
4580      }
4581      sub[iBlock].dual();
4582      memcpy(lower2,saveLower,numberRows2*sizeof(double));
4583      memcpy(upper2,saveUpper,numberRows2*sizeof(double));
4584      // get proposal
4585      if (sub[iBlock].numberIterations()||!iPass) {
4586        double objValue =0.0;
4587        int start = rowAdd[numberProposals];
4588        // proposal
4589        if (sub[iBlock].isProvenOptimal()) {
4590          const double * solution = sub[iBlock].dualRowSolution();
4591          first[iBlock]->transposeTimes(solution,elementAdd+start);
4592          for (int i=0;i<numberRows2;i++) {
4593            if (solution[i]<-dualTolerance_) {
4594              // at upper
4595              assert (saveUpper[i]<1.0e30);
4596              objValue += solution[i]*saveUpper[i];
4597            } else if (solution[i]>dualTolerance_) {
4598              // at lower
4599              assert (saveLower[i]>-1.0e30);
4600              objValue += solution[i]*saveLower[i];
4601            }
4602          }
4603         
4604          // See if cuts off and pack down
4605          int number=start;
4606          double infeas = objValue;
4607          double smallest=1.0e100;
4608          double largest=0.0;
4609          for (int i=0;i<numberMasterColumns;i++) {
4610            double value = elementAdd[start+i];
4611            if (fabs(value)>1.0e-15) {
4612              infeas -= primal[i]*value;
4613              smallest = CoinMin(smallest,fabs(value));
4614              largest = CoinMax(largest,fabs(value));
4615              columnAdd[number]=i;
4616              elementAdd[number++]=-value;
4617            }
4618          }
4619          columnAdd[number]=numberMasterColumns+iBlock;
4620          elementAdd[number++]=-1.0;
4621          // if elements large then scale?
4622          //if (largest>1.0e8||smallest<1.0e-8)
4623          printf("For subproblem %d smallest - %g, largest %g - infeas %g\n",
4624                 iBlock,smallest,largest,infeas);
4625          if (infeas<-1.0e-6||!iPass) {
4626            // take
4627            objective[numberProposals]=objValue;
4628            rowAdd[++numberProposals]=number;
4629            when[numberRowsGenerated]=iPass;
4630            whichBlock[numberRowsGenerated++]=iBlock;
4631          }
4632        } else if (sub[iBlock].isProvenPrimalInfeasible()) {
4633          // use ray
4634          const double * solution = sub[iBlock].infeasibilityRay();
4635          first[iBlock]->transposeTimes(solution,elementAdd+start);
4636          for (int i=0;i<numberRows2;i++) {
4637            if (solution[i]<-dualTolerance_) {
4638              // at upper
4639              assert (saveUpper[i]<1.0e30);
4640              objValue += solution[i]*saveUpper[i];
4641            } else if (solution[i]>dualTolerance_) {
4642              // at lower
4643              assert (saveLower[i]>-1.0e30);
4644              objValue += solution[i]*saveLower[i];
4645            }
4646          }
4647          // See if good infeas and pack down
4648          int number=start;
4649          double infeas = objValue;
4650          double smallest=1.0e100;
4651          double largest=0.0;
4652          for (int i=0;i<numberMasterColumns;i++) {
4653            double value = elementAdd[start+i];
4654            if (fabs(value)>1.0e-15) {
4655              infeas -= primal[i]*value;
4656              smallest = CoinMin(smallest,fabs(value));
4657              largest = CoinMax(largest,fabs(value));
4658              columnAdd[number]=i;
4659              elementAdd[number++]=-value;
4660            }
4661          }
4662          // if elements large or small then scale?
4663          //if (largest>1.0e8||smallest<1.0e-8)
4664          printf("For subproblem ray %d smallest - %g, largest %g - infeas %g\n",
4665                 iBlock,smallest,largest,infeas);
4666          if (infeas<-1.0e-6) {
4667            // take
4668            objective[numberProposals]=objValue;
4669            rowAdd[++numberProposals]=number;
4670            when[numberRowsGenerated]=iPass;
4671            whichBlock[numberRowsGenerated++]=iBlock;
4672          }
4673        } else {
4674          abort();
4675        }
4676      }
4677      delete [] saveLower;
4678      delete [] saveUpper;
4679    }
4680    if (deletePrimal)
4681      delete [] primal;
4682    if (numberProposals) {
4683      master.addRows(numberProposals,NULL,objective,
4684                        rowAdd,columnAdd,elementAdd);
4685    }
4686  }
4687  std::cout<<"Time at end of Benders "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4688  //master.scaling(0);
4689  //master.primal(1);
4690  loadProblem(*model);
4691  // now put back a good solution
4692  const double * columnSolution = master.primalColumnSolution();
4693  double * fullColumnSolution = primalColumnSolution();
4694  const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4695  int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4696  const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4697  int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4698  for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4699    int kColumn = columnBack[iColumn];
4700    setColumnStatus(kColumn,master.getColumnStatus(iColumn));
4701    fullColumnSolution[kColumn]=columnSolution[iColumn];
4702  }
4703  for (int iRow=0;iRow<numberRows2;iRow++) {
4704    int kRow = rowBack[iRow];
4705    setStatus(kRow,master.getStatus(iRow));
4706    //fullSolution[kRow]=solution[iRow];
4707  }
4708  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
4709    // move basis
4710    int kBlock = rowCounts[iBlock];
4711    const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4712    int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4713    const int * rowBack = model->coinBlock(kBlock)->originalRows();
4714    int numberRows2 = model->coinBlock(kBlock)->numberRows();
4715    const double * subColumnSolution = sub[iBlock].primalColumnSolution();
4716    for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
4717      int kColumn = columnBack[iColumn];
4718      setColumnStatus(kColumn,sub[iBlock].getColumnStatus(iColumn));
4719      fullColumnSolution[kColumn]=subColumnSolution[iColumn];
4720    }
4721    for (int iRow=0;iRow<numberRows2;iRow++) {
4722      int kRow = rowBack[iRow];
4723      setStatus(kRow,sub[iBlock].getStatus(iRow));
4724      setStatus(kRow,atLowerBound);
4725    }
4726  }
4727  double * fullSolution = primalRowSolution();
4728  CoinZeroN(fullSolution,numberRows_);
4729  times(1.0,fullColumnSolution,fullSolution);
4730  //int numberColumns=model->numberColumns();
4731  //for (int iColumn=0;iColumn<numberColumns;iColumn++)
4732  //setColumnStatus(iColumn,ClpSimplex::superBasic);
4733  std::cout<<"Time before cleanup of full model "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4734  this->primal(1);
4735  std::cout<<"Total time "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
4736  delete [] rowCounts;
4737  //delete [] sol;
4738  //delete [] lower;
4739  //delete [] upper;
4740  delete [] whichBlock;
4741  delete [] when;
4742  delete [] rowAdd;
4743  delete [] columnAdd;
4744  delete [] elementAdd;
4745  delete [] objective;
4746  delete [] first;
4747  delete [] sub;
4748  return 0;
4749}
Note: See TracBrowser for help on using the repository browser.