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

Last change on this file since 2044 was 2044, checked in by stefan, 5 years ago

use Sleep on windows and usleep on non-windows

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