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

Last change on this file since 1611 was 1611, checked in by forrest, 9 years ago

fix typo

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