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

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

for build pthreads for use in Cbc analyze

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