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

Last change on this file since 2030 was 2026, checked in by forrest, 6 years ago

fix compile error for abc

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