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

Last change on this file since 1655 was 1655, checked in by mjs, 8 years ago

Relicense under the Eclipse Public License.

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