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

Last change on this file since 2321 was 2321, checked in by forrest, 20 months ago

primal if any superbasic

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