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

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

fix bug in decomposition

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