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

Last change on this file since 1321 was 1321, checked in by forrest, 13 years ago

out compiler warnings and stability improvements

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