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

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

minor mistake

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 201.6 KB
Line 
1/* $Id: ClpSolve.cpp 1612 2010-09-30 16:17:02Z 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          int oldStatus=problemStatus_;
2486          setProblemStatus(finalStatus);
2487          setSecondaryStatus(finalSecondaryStatus);
2488          int rcode=eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
2489          if (finalStatus != 3 && rcode < 0 && (finalStatus || oldStatus == -1)) {
2490               int savePerturbation = perturbation();
2491               if (!finalStatus || (moreSpecialOptions_ & 2) == 0) {
2492                    if (finalStatus == 2) {
2493                         // unbounded - get feasible first
2494                         double save = optimizationDirection_;
2495                         optimizationDirection_ = 0.0;
2496                         primal(1);
2497                         optimizationDirection_ = save;
2498                         primal(1);
2499                    } else if (finalStatus == 1) {
2500                         dual();
2501                    } else {
2502                         setPerturbation(100);
2503                         primal(1);
2504                    }
2505               } else {
2506                    // just set status
2507                    problemStatus_ = finalStatus;
2508               }
2509               setPerturbation(savePerturbation);
2510               numberIterations += numberIterations_;
2511               numberIterations_ = numberIterations;
2512               finalStatus = status();
2513               time2 = CoinCpuTime();
2514               handler_->message(CLP_INTERVAL_TIMING, messages_)
2515                         << "Cleanup" << time2 - timeX << time2 - time1
2516                         << CoinMessageEol;
2517               timeX = time2;
2518          } else if (rcode >= 0) {
2519              primal(1);
2520          } else {
2521               secondaryStatus_ = finalSecondaryStatus;
2522          }
2523     } else if (model2 != this) {
2524          // not presolved - but different model used (sprint probably)
2525          CoinMemcpyN(model2->primalRowSolution(),
2526                      numberRows_, this->primalRowSolution());
2527          CoinMemcpyN(model2->dualRowSolution(),
2528                      numberRows_, this->dualRowSolution());
2529          CoinMemcpyN(model2->primalColumnSolution(),
2530                      numberColumns_, this->primalColumnSolution());
2531          CoinMemcpyN(model2->dualColumnSolution(),
2532                      numberColumns_, this->dualColumnSolution());
2533          CoinMemcpyN(model2->statusArray(),
2534                      numberColumns_ + numberRows_, this->statusArray());
2535          objectiveValue_ = model2->objectiveValue_;
2536          numberIterations_ = model2->numberIterations_;
2537          problemStatus_ = model2->problemStatus_;
2538          secondaryStatus_ = model2->secondaryStatus_;
2539          delete model2;
2540     }
2541     if (method != ClpSolve::useBarrierNoCross &&
2542               method != ClpSolve::useBarrier)
2543          setMaximumIterations(saveMaxIterations);
2544     std::string statusMessage[] = {"Unknown", "Optimal", "PrimalInfeasible", "DualInfeasible", "Stopped",
2545                                    "Errors", "User stopped"
2546                                   };
2547     assert (finalStatus >= -1 && finalStatus <= 5);
2548     handler_->message(CLP_TIMING, messages_)
2549               << statusMessage[finalStatus+1] << objectiveValue() << numberIterations << time2 - time1;
2550     handler_->printing(presolve == ClpSolve::presolveOn)
2551               << timePresolve;
2552     handler_->printing(timeIdiot != 0.0)
2553               << timeIdiot;
2554     handler_->message() << CoinMessageEol;
2555     if (interrupt)
2556          signal(SIGINT, saveSignal);
2557     perturbation_ = savePerturbation;
2558     scalingFlag_ = saveScaling;
2559     // If faking objective - put back correct one
2560     if (savedObjective) {
2561          delete objective_;
2562          objective_ = savedObjective;
2563     }
2564     if (options.getSpecialOption(1) == 2 &&
2565               options.getExtraInfo(1) > 1000000) {
2566          ClpObjective * savedObjective = objective_;
2567          // make up zero objective
2568          double * obj = new double[numberColumns_];
2569          for (int i = 0; i < numberColumns_; i++)
2570               obj[i] = 0.0;
2571          objective_ = new ClpLinearObjective(obj, numberColumns_);
2572          delete [] obj;
2573          primal(1);
2574          delete objective_;
2575          objective_ = savedObjective;
2576          finalStatus = status();
2577     }
2578     eventHandler()->event(ClpEventHandler::presolveEnd);
2579     return finalStatus;
2580}
2581// General solve
2582int
2583ClpSimplex::initialSolve()
2584{
2585     // Default so use dual
2586     ClpSolve options;
2587     return initialSolve(options);
2588}
2589// General dual solve
2590int
2591ClpSimplex::initialDualSolve()
2592{
2593     ClpSolve options;
2594     // Use dual
2595     options.setSolveType(ClpSolve::useDual);
2596     return initialSolve(options);
2597}
2598// General dual solve
2599int
2600ClpSimplex::initialPrimalSolve()
2601{
2602     ClpSolve options;
2603     // Use primal
2604     options.setSolveType(ClpSolve::usePrimal);
2605     return initialSolve(options);
2606}
2607// barrier solve, not to be followed by crossover
2608int
2609ClpSimplex::initialBarrierNoCrossSolve()
2610{
2611     ClpSolve options;
2612     // Use primal
2613     options.setSolveType(ClpSolve::useBarrierNoCross);
2614     return initialSolve(options);
2615}
2616
2617// General barrier solve
2618int
2619ClpSimplex::initialBarrierSolve()
2620{
2621     ClpSolve options;
2622     // Use primal
2623     options.setSolveType(ClpSolve::useBarrier);
2624     return initialSolve(options);
2625}
2626
2627// Default constructor
2628ClpSolve::ClpSolve (  )
2629{
2630     method_ = automatic;
2631     presolveType_ = presolveOn;
2632     numberPasses_ = 5;
2633     int i;
2634     for (i = 0; i < 7; i++)
2635          options_[i] = 0;
2636     // say no +-1 matrix
2637     options_[3] = 1;
2638     for (i = 0; i < 7; i++)
2639          extraInfo_[i] = -1;
2640     independentOptions_[0] = 0;
2641     // But switch off slacks
2642     independentOptions_[1] = 512;
2643     // Substitute up to 3
2644     independentOptions_[2] = 3;
2645
2646}
2647// Constructor when you really know what you are doing
2648ClpSolve::ClpSolve ( SolveType method, PresolveType presolveType,
2649                     int numberPasses, int options[6],
2650                     int extraInfo[6], int independentOptions[3])
2651{
2652     method_ = method;
2653     presolveType_ = presolveType;
2654     numberPasses_ = numberPasses;
2655     int i;
2656     for (i = 0; i < 6; i++)
2657          options_[i] = options[i];
2658     options_[6] = 0;
2659     for (i = 0; i < 6; i++)
2660          extraInfo_[i] = extraInfo[i];
2661     extraInfo_[6] = 0;
2662     for (i = 0; i < 3; i++)
2663          independentOptions_[i] = independentOptions[i];
2664}
2665
2666// Copy constructor.
2667ClpSolve::ClpSolve(const ClpSolve & rhs)
2668{
2669     method_ = rhs.method_;
2670     presolveType_ = rhs.presolveType_;
2671     numberPasses_ = rhs.numberPasses_;
2672     int i;
2673     for ( i = 0; i < 7; i++)
2674          options_[i] = rhs.options_[i];
2675     for ( i = 0; i < 7; i++)
2676          extraInfo_[i] = rhs.extraInfo_[i];
2677     for (i = 0; i < 3; i++)
2678          independentOptions_[i] = rhs.independentOptions_[i];
2679}
2680// Assignment operator. This copies the data
2681ClpSolve &
2682ClpSolve::operator=(const ClpSolve & rhs)
2683{
2684     if (this != &rhs) {
2685          method_ = rhs.method_;
2686          presolveType_ = rhs.presolveType_;
2687          numberPasses_ = rhs.numberPasses_;
2688          int i;
2689          for (i = 0; i < 7; i++)
2690               options_[i] = rhs.options_[i];
2691          for (i = 0; i < 7; i++)
2692               extraInfo_[i] = rhs.extraInfo_[i];
2693          for (i = 0; i < 3; i++)
2694               independentOptions_[i] = rhs.independentOptions_[i];
2695     }
2696     return *this;
2697
2698}
2699// Destructor
2700ClpSolve::~ClpSolve (  )
2701{
2702}
2703// See header file for deatils
2704void
2705ClpSolve::setSpecialOption(int which, int value, int extraInfo)
2706{
2707     options_[which] = value;
2708     extraInfo_[which] = extraInfo;
2709}
2710int
2711ClpSolve::getSpecialOption(int which) const
2712{
2713     return options_[which];
2714}
2715
2716// Solve types
2717void
2718ClpSolve::setSolveType(SolveType method, int /*extraInfo*/)
2719{
2720     method_ = method;
2721}
2722
2723ClpSolve::SolveType
2724ClpSolve::getSolveType()
2725{
2726     return method_;
2727}
2728
2729// Presolve types
2730void
2731ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
2732{
2733     presolveType_ = amount;
2734     numberPasses_ = extraInfo;
2735}
2736ClpSolve::PresolveType
2737ClpSolve::getPresolveType()
2738{
2739     return presolveType_;
2740}
2741// Extra info for idiot (or sprint)
2742int
2743ClpSolve::getExtraInfo(int which) const
2744{
2745     return extraInfo_[which];
2746}
2747int
2748ClpSolve::getPresolvePasses() const
2749{
2750     return numberPasses_;
2751}
2752/* Say to return at once if infeasible,
2753   default is to solve */
2754void
2755ClpSolve::setInfeasibleReturn(bool trueFalse)
2756{
2757     independentOptions_[0] = trueFalse ? 1 : 0;
2758}
2759#include <string>
2760// Generates code for above constructor
2761void
2762ClpSolve::generateCpp(FILE * fp)
2763{
2764     std::string solveType[] = {
2765          "ClpSolve::useDual",
2766          "ClpSolve::usePrimal",
2767          "ClpSolve::usePrimalorSprint",
2768          "ClpSolve::useBarrier",
2769          "ClpSolve::useBarrierNoCross",
2770          "ClpSolve::automatic",
2771          "ClpSolve::notImplemented"
2772     };
2773     std::string presolveType[] =  {
2774          "ClpSolve::presolveOn",
2775          "ClpSolve::presolveOff",
2776          "ClpSolve::presolveNumber",
2777          "ClpSolve::presolveNumberCost"
2778     };
2779     fprintf(fp, "3  ClpSolve::SolveType method = %s;\n", solveType[method_].c_str());
2780     fprintf(fp, "3  ClpSolve::PresolveType presolveType = %s;\n",
2781             presolveType[presolveType_].c_str());
2782     fprintf(fp, "3  int numberPasses = %d;\n", numberPasses_);
2783     fprintf(fp, "3  int options[] = {%d,%d,%d,%d,%d,%d};\n",
2784             options_[0], options_[1], options_[2],
2785             options_[3], options_[4], options_[5]);
2786     fprintf(fp, "3  int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
2787             extraInfo_[0], extraInfo_[1], extraInfo_[2],
2788             extraInfo_[3], extraInfo_[4], extraInfo_[5]);
2789     fprintf(fp, "3  int independentOptions[] = {%d,%d,%d};\n",
2790             independentOptions_[0], independentOptions_[1], independentOptions_[2]);
2791     fprintf(fp, "3  ClpSolve clpSolve(method,presolveType,numberPasses,\n");
2792     fprintf(fp, "3                    options,extraInfo,independentOptions);\n");
2793}
2794//#############################################################################
2795#include "ClpNonLinearCost.hpp"
2796
2797ClpSimplexProgress::ClpSimplexProgress ()
2798{
2799     int i;
2800     for (i = 0; i < CLP_PROGRESS; i++) {
2801          objective_[i] = COIN_DBL_MAX;
2802          infeasibility_[i] = -1.0; // set to an impossible value
2803          realInfeasibility_[i] = COIN_DBL_MAX;
2804          numberInfeasibilities_[i] = -1;
2805          iterationNumber_[i] = -1;
2806     }
2807#ifdef CLP_PROGRESS_WEIGHT
2808     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
2809          objectiveWeight_[i] = COIN_DBL_MAX;
2810          infeasibilityWeight_[i] = -1.0; // set to an impossible value
2811          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
2812          numberInfeasibilitiesWeight_[i] = -1;
2813          iterationNumberWeight_[i] = -1;
2814     }
2815     drop_ = 0.0;
2816     best_ = 0.0;
2817#endif
2818     initialWeight_ = 0.0;
2819     for (i = 0; i < CLP_CYCLE; i++) {
2820          //obj_[i]=COIN_DBL_MAX;
2821          in_[i] = -1;
2822          out_[i] = -1;
2823          way_[i] = 0;
2824     }
2825     numberTimes_ = 0;
2826     numberBadTimes_ = 0;
2827     numberReallyBadTimes_ = 0;
2828     numberTimesFlagged_ = 0;
2829     model_ = NULL;
2830     oddState_ = 0;
2831}
2832
2833
2834//-----------------------------------------------------------------------------
2835
2836ClpSimplexProgress::~ClpSimplexProgress ()
2837{
2838}
2839// Copy constructor.
2840ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
2841{
2842     int i;
2843     for (i = 0; i < CLP_PROGRESS; i++) {
2844          objective_[i] = rhs.objective_[i];
2845          infeasibility_[i] = rhs.infeasibility_[i];
2846          realInfeasibility_[i] = rhs.realInfeasibility_[i];
2847          numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
2848          iterationNumber_[i] = rhs.iterationNumber_[i];
2849     }
2850#ifdef CLP_PROGRESS_WEIGHT
2851     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
2852          objectiveWeight_[i] = rhs.objectiveWeight_[i];
2853          infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2854          realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2855          numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
2856          iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
2857     }
2858     drop_ = rhs.drop_;
2859     best_ = rhs.best_;
2860#endif
2861     initialWeight_ = rhs.initialWeight_;
2862     for (i = 0; i < CLP_CYCLE; i++) {
2863          //obj_[i]=rhs.obj_[i];
2864          in_[i] = rhs.in_[i];
2865          out_[i] = rhs.out_[i];
2866          way_[i] = rhs.way_[i];
2867     }
2868     numberTimes_ = rhs.numberTimes_;
2869     numberBadTimes_ = rhs.numberBadTimes_;
2870     numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
2871     numberTimesFlagged_ = rhs.numberTimesFlagged_;
2872     model_ = rhs.model_;
2873     oddState_ = rhs.oddState_;
2874}
2875// Copy constructor.from model
2876ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model)
2877{
2878     model_ = model;
2879     reset();
2880     initialWeight_ = 0.0;
2881}
2882// Fill from model
2883void
2884ClpSimplexProgress::fillFromModel ( ClpSimplex * model )
2885{
2886     model_ = model;
2887     reset();
2888     initialWeight_ = 0.0;
2889}
2890// Assignment operator. This copies the data
2891ClpSimplexProgress &
2892ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2893{
2894     if (this != &rhs) {
2895          int i;
2896          for (i = 0; i < CLP_PROGRESS; i++) {
2897               objective_[i] = rhs.objective_[i];
2898               infeasibility_[i] = rhs.infeasibility_[i];
2899               realInfeasibility_[i] = rhs.realInfeasibility_[i];
2900               numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
2901               iterationNumber_[i] = rhs.iterationNumber_[i];
2902          }
2903#ifdef CLP_PROGRESS_WEIGHT
2904          for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
2905               objectiveWeight_[i] = rhs.objectiveWeight_[i];
2906               infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
2907               realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
2908               numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
2909               iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
2910          }
2911          drop_ = rhs.drop_;
2912          best_ = rhs.best_;
2913#endif
2914          initialWeight_ = rhs.initialWeight_;
2915          for (i = 0; i < CLP_CYCLE; i++) {
2916               //obj_[i]=rhs.obj_[i];
2917               in_[i] = rhs.in_[i];
2918               out_[i] = rhs.out_[i];
2919               way_[i] = rhs.way_[i];
2920          }
2921          numberTimes_ = rhs.numberTimes_;
2922          numberBadTimes_ = rhs.numberBadTimes_;
2923          numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
2924          numberTimesFlagged_ = rhs.numberTimesFlagged_;
2925          model_ = rhs.model_;
2926          oddState_ = rhs.oddState_;
2927     }
2928     return *this;
2929}
2930// Seems to be something odd about exact comparison of doubles on linux
2931static bool equalDouble(double value1, double value2)
2932{
2933
2934     union {
2935          double d;
2936          int i[2];
2937     } v1, v2;
2938     v1.d = value1;
2939     v2.d = value2;
2940     if (sizeof(int) * 2 == sizeof(double))
2941          return (v1.i[0] == v2.i[0] && v1.i[1] == v2.i[1]);
2942     else
2943          return (v1.i[0] == v2.i[0]);
2944}
2945int
2946ClpSimplexProgress::looping()
2947{
2948     if (!model_)
2949          return -1;
2950     double objective = model_->rawObjectiveValue();
2951     if (model_->algorithm() < 0)
2952          objective -= model_->bestPossibleImprovement();
2953     double infeasibility;
2954     double realInfeasibility = 0.0;
2955     int numberInfeasibilities;
2956     int iterationNumber = model_->numberIterations();
2957     numberTimesFlagged_ = 0;
2958     if (model_->algorithm() < 0) {
2959          // dual
2960          infeasibility = model_->sumPrimalInfeasibilities();
2961          numberInfeasibilities = model_->numberPrimalInfeasibilities();
2962     } else {
2963          //primal
2964          infeasibility = model_->sumDualInfeasibilities();
2965          realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
2966          numberInfeasibilities = model_->numberDualInfeasibilities();
2967     }
2968     int i;
2969     int numberMatched = 0;
2970     int matched = 0;
2971     int nsame = 0;
2972     for (i = 0; i < CLP_PROGRESS; i++) {
2973          bool matchedOnObjective = equalDouble(objective, objective_[i]);
2974          bool matchedOnInfeasibility = equalDouble(infeasibility, infeasibility_[i]);
2975          bool matchedOnInfeasibilities =
2976               (numberInfeasibilities == numberInfeasibilities_[i]);
2977
2978          if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
2979               matched |= (1 << i);
2980               // Check not same iteration
2981               if (iterationNumber != iterationNumber_[i]) {
2982                    numberMatched++;
2983                    // here mainly to get over compiler bug?
2984                    if (model_->messageHandler()->logLevel() > 10)
2985                         printf("%d %d %d %d %d loop check\n", i, numberMatched,
2986                                matchedOnObjective, matchedOnInfeasibility,
2987                                matchedOnInfeasibilities);
2988               } else {
2989                    // stuck but code should notice
2990                    nsame++;
2991               }
2992          }
2993          if (i) {
2994               objective_[i-1] = objective_[i];
2995               infeasibility_[i-1] = infeasibility_[i];
2996               realInfeasibility_[i-1] = realInfeasibility_[i];
2997               numberInfeasibilities_[i-1] = numberInfeasibilities_[i];
2998               iterationNumber_[i-1] = iterationNumber_[i];
2999          }
3000     }
3001     objective_[CLP_PROGRESS-1] = objective;
3002     infeasibility_[CLP_PROGRESS-1] = infeasibility;
3003     realInfeasibility_[CLP_PROGRESS-1] = realInfeasibility;
3004     numberInfeasibilities_[CLP_PROGRESS-1] = numberInfeasibilities;
3005     iterationNumber_[CLP_PROGRESS-1] = iterationNumber;
3006     if (nsame == CLP_PROGRESS)
3007          numberMatched = CLP_PROGRESS; // really stuck
3008     if (model_->progressFlag())
3009          numberMatched = 0;
3010     numberTimes_++;
3011     if (numberTimes_ < 10)
3012          numberMatched = 0;
3013     // skip if just last time as may be checking something
3014     if (matched == (1 << (CLP_PROGRESS - 1)))
3015          numberMatched = 0;
3016     if (numberMatched) {
3017          model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
3018                    << numberMatched
3019                    << matched
3020                    << numberTimes_
3021                    << CoinMessageEol;
3022          numberBadTimes_++;
3023          if (numberBadTimes_ < 10) {
3024               // make factorize every iteration
3025               model_->forceFactorization(1);
3026               if (numberBadTimes_ < 2) {
3027                    startCheck(); // clear other loop check
3028                    if (model_->algorithm() < 0) {
3029                         // dual - change tolerance
3030                         model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
3031                         // if infeasible increase dual bound
3032                         if (model_->dualBound() < 1.0e17) {
3033                              model_->setDualBound(model_->dualBound() * 1.1);
3034                              static_cast<ClpSimplexDual *>(model_)->resetFakeBounds(0);
3035                         }
3036                    } else {
3037                         // primal - change tolerance
3038                         if (numberBadTimes_ > 3)
3039                              model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
3040                         // if infeasible increase infeasibility cost
3041                         if (model_->nonLinearCost()->numberInfeasibilities() &&
3042                                   model_->infeasibilityCost() < 1.0e17) {
3043                              model_->setInfeasibilityCost(model_->infeasibilityCost() * 1.1);
3044                         }
3045                    }
3046               } else {
3047                    // flag
3048                    int iSequence;
3049                    if (model_->algorithm() < 0) {
3050                         // dual
3051                         if (model_->dualBound() > 1.0e14)
3052                              model_->setDualBound(1.0e14);
3053                         iSequence = in_[CLP_CYCLE-1];
3054                    } else {
3055                         // primal
3056                         if (model_->infeasibilityCost() > 1.0e14)
3057                              model_->setInfeasibilityCost(1.0e14);
3058                         iSequence = out_[CLP_CYCLE-1];
3059                    }
3060                    if (iSequence >= 0) {
3061                         char x = model_->isColumn(iSequence) ? 'C' : 'R';
3062                         if (model_->messageHandler()->logLevel() >= 63)
3063                              model_->messageHandler()->message(CLP_SIMPLEX_FLAG, model_->messages())
3064                                        << x << model_->sequenceWithin(iSequence)
3065                                        << CoinMessageEol;
3066                         // if Gub then needs to be sequenceIn_
3067                         int save = model_->sequenceIn();
3068                         model_->setSequenceIn(iSequence);
3069                         model_->setFlagged(iSequence);
3070                         model_->setSequenceIn(save);
3071                         //printf("flagging %d from loop\n",iSequence);
3072                         startCheck();
3073                    } else {
3074                         // Give up
3075                         if (model_->messageHandler()->logLevel() >= 63)
3076                              printf("***** All flagged?\n");
3077                         return 4;
3078                    }
3079                    // reset
3080                    numberBadTimes_ = 2;
3081               }
3082               return -2;
3083          } else {
3084               // look at solution and maybe declare victory
3085               if (infeasibility < 1.0e-4) {
3086                    return 0;
3087               } else {
3088                    model_->messageHandler()->message(CLP_LOOP, model_->messages())
3089                              << CoinMessageEol;
3090#ifndef NDEBUG
3091                    printf("debug loop ClpSimplex A\n");
3092                    abort();
3093#endif
3094                    return 3;
3095               }
3096          }
3097     }
3098     return -1;
3099}
3100// Resets as much as possible
3101void
3102ClpSimplexProgress::reset()
3103{
3104     int i;
3105     for (i = 0; i < CLP_PROGRESS; i++) {
3106          if (model_->algorithm() >= 0)
3107               objective_[i] = COIN_DBL_MAX;
3108          else
3109               objective_[i] = -COIN_DBL_MAX;
3110          infeasibility_[i] = -1.0; // set to an impossible value
3111          realInfeasibility_[i] = COIN_DBL_MAX;
3112          numberInfeasibilities_[i] = -1;
3113          iterationNumber_[i] = -1;
3114     }
3115#ifdef CLP_PROGRESS_WEIGHT
3116     for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
3117          objectiveWeight_[i] = COIN_DBL_MAX;
3118          infeasibilityWeight_[i] = -1.0; // set to an impossible value
3119          realInfeasibilityWeight_[i] = COIN_DBL_MAX;
3120          numberInfeasibilitiesWeight_[i] = -1;
3121          iterationNumberWeight_[i] = -1;
3122     }
3123     drop_ = 0.0;
3124     best_ = 0.0;
3125#endif
3126     for (i = 0; i < CLP_CYCLE; i++) {
3127          //obj_[i]=COIN_DBL_MAX;
3128          in_[i] = -1;
3129          out_[i] = -1;
3130          way_[i] = 0;
3131     }
3132     numberTimes_ = 0;
3133     numberBadTimes_ = 0;
3134     numberReallyBadTimes_ = 0;
3135     numberTimesFlagged_ = 0;
3136     oddState_ = 0;
3137}
3138// Returns previous objective (if -1) - current if (0)
3139double
3140ClpSimplexProgress::lastObjective(int back) const
3141{
3142     return objective_[CLP_PROGRESS-1-back];
3143}
3144// Returns previous infeasibility (if -1) - current if (0)
3145double
3146ClpSimplexProgress::lastInfeasibility(int back) const
3147{
3148     return realInfeasibility_[CLP_PROGRESS-1-back];
3149}
3150// Sets real primal infeasibility
3151void
3152ClpSimplexProgress::setInfeasibility(double value)
3153{
3154     for (int i = 1; i < CLP_PROGRESS; i++)
3155          realInfeasibility_[i-1] = realInfeasibility_[i];
3156     realInfeasibility_[CLP_PROGRESS-1] = value;
3157}
3158// Modify objective e.g. if dual infeasible in dual
3159void
3160ClpSimplexProgress::modifyObjective(double value)
3161{
3162     objective_[CLP_PROGRESS-1] = value;
3163}
3164// Returns previous iteration number (if -1) - current if (0)
3165int
3166ClpSimplexProgress::lastIterationNumber(int back) const
3167{
3168     return iterationNumber_[CLP_PROGRESS-1-back];
3169}
3170// clears iteration numbers (to switch off panic)
3171void
3172ClpSimplexProgress::clearIterationNumbers()
3173{
3174     for (int i = 0; i < CLP_PROGRESS; i++)
3175          iterationNumber_[i] = -1;
3176}
3177// Start check at beginning of whileIterating
3178void
3179ClpSimplexProgress::startCheck()
3180{
3181     int i;
3182     for (i = 0; i < CLP_CYCLE; i++) {
3183          //obj_[i]=COIN_DBL_MAX;
3184          in_[i] = -1;
3185          out_[i] = -1;
3186          way_[i] = 0;
3187     }
3188}
3189// Returns cycle length in whileIterating
3190int
3191ClpSimplexProgress::cycle(int in, int out, int wayIn, int wayOut)
3192{
3193     int i;
3194#if 0
3195     if (model_->numberIterations() > 206571) {
3196          printf("in %d out %d\n", in, out);
3197          for (i = 0; i < CLP_CYCLE; i++)
3198               printf("cy %d in %d out %d\n", i, in_[i], out_[i]);
3199     }
3200#endif
3201     int matched = 0;
3202     // first see if in matches any out
3203     for (i = 1; i < CLP_CYCLE; i++) {
3204          if (in == out_[i]) {
3205               // even if flip then suspicious
3206               matched = -1;
3207               break;
3208          }
3209     }
3210#if 0
3211     if (!matched || in_[0] < 0) {
3212          // can't be cycle
3213          for (i = 0; i < CLP_CYCLE - 1; i++) {
3214               //obj_[i]=obj_[i+1];
3215               in_[i] = in_[i+1];
3216               out_[i] = out_[i+1];
3217               way_[i] = way_[i+1];
3218          }
3219     } else {
3220          // possible cycle
3221          matched = 0;
3222          for (i = 0; i < CLP_CYCLE - 1; i++) {
3223               int k;
3224               char wayThis = way_[i];
3225               int inThis = in_[i];
3226               int outThis = out_[i];
3227               //double objThis = obj_[i];
3228               for(k = i + 1; k < CLP_CYCLE; k++) {
3229                    if (inThis == in_[k] && outThis == out_[k] && wayThis == way_[k]) {
3230                         int distance = k - i;
3231                         if (k + distance < CLP_CYCLE) {
3232                              // See if repeats
3233                              int j = k + distance;
3234                              if (inThis == in_[j] && outThis == out_[j] && wayThis == way_[j]) {
3235                                   matched = distance;
3236                                   break;
3237                              }
3238                         } else {
3239                              matched = distance;
3240                              break;
3241                         }
3242                    }
3243               }
3244               //obj_[i]=obj_[i+1];
3245               in_[i] = in_[i+1];
3246               out_[i] = out_[i+1];
3247               way_[i] = way_[i+1];
3248          }
3249     }
3250#else
3251     if (matched && in_[0] >= 0) {
3252          // possible cycle - only check [0] against all
3253          matched = 0;
3254          int nMatched = 0;
3255          char way0 = way_[0];
3256          int in0 = in_[0];
3257          int out0 = out_[0];
3258          //double obj0 = obj_[i];
3259          for(int k = 1; k < CLP_CYCLE - 4; k++) {
3260               if (in0 == in_[k] && out0 == out_[k] && way0 == way_[k]) {
3261                    nMatched++;
3262                    // See if repeats
3263                    int end = CLP_CYCLE - k;
3264                    int j;
3265                    for ( j = 1; j < end; j++) {
3266                         if (in_[j+k] != in_[j] || out_[j+k] != out_[j] || way_[j+k] != way_[j])
3267                              break;
3268                    }
3269                    if (j == end) {
3270                         matched = k;
3271                         break;
3272                    }
3273               }
3274          }
3275          // If three times then that is too much even if not regular
3276          if (matched <= 0 && nMatched > 1)
3277               matched = 100;
3278     }
3279     for (i = 0; i < CLP_CYCLE - 1; i++) {
3280          //obj_[i]=obj_[i+1];
3281          in_[i] = in_[i+1];
3282          out_[i] = out_[i+1];
3283          way_[i] = way_[i+1];
3284     }
3285#endif
3286     int way = 1 - wayIn + 4 * (1 - wayOut);
3287     //obj_[i]=model_->objectiveValue();
3288     in_[CLP_CYCLE-1] = in;
3289     out_[CLP_CYCLE-1] = out;
3290     way_[CLP_CYCLE-1] = static_cast<char>(way);
3291     return matched;
3292}
3293#include "CoinStructuredModel.hpp"
3294// Solve using structure of model and maybe in parallel
3295int
3296ClpSimplex::solve(CoinStructuredModel * model)
3297{
3298     // analyze structure
3299     int numberRowBlocks = model->numberRowBlocks();
3300     int numberColumnBlocks = model->numberColumnBlocks();
3301     int numberElementBlocks = model->numberElementBlocks();
3302     if (numberElementBlocks == 1) {
3303          loadProblem(*model, false);
3304          return dual();
3305     }
3306     // For now just get top level structure
3307     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3308     for (int i = 0; i < numberElementBlocks; i++) {
3309          CoinStructuredModel * subModel =
3310               dynamic_cast<CoinStructuredModel *>(model->block(i));
3311          CoinModel * thisBlock;
3312          if (subModel) {
3313               thisBlock = subModel->coinModelBlock(blockInfo[i]);
3314               model->setCoinModel(thisBlock, i);
3315          } else {
3316               thisBlock = dynamic_cast<CoinModel *>(model->block(i));
3317               assert (thisBlock);
3318               // just fill in info
3319               CoinModelBlockInfo info = CoinModelBlockInfo();
3320               int whatsSet = thisBlock->whatIsSet();
3321               info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
3322               info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
3323               info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
3324               info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
3325               info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
3326               info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
3327               // Which block
3328               int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
3329               info.rowBlock = iRowBlock;
3330               int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
3331               info.columnBlock = iColumnBlock;
3332               blockInfo[i] = info;
3333          }
3334     }
3335     int * rowCounts = new int [numberRowBlocks];
3336     CoinZeroN(rowCounts, numberRowBlocks);
3337     int * columnCounts = new int [numberColumnBlocks+1];
3338     CoinZeroN(columnCounts, numberColumnBlocks);
3339     int decomposeType = 0;
3340     for (int i = 0; i < numberElementBlocks; i++) {
3341          int iRowBlock = blockInfo[i].rowBlock;
3342          int iColumnBlock = blockInfo[i].columnBlock;
3343          rowCounts[iRowBlock]++;
3344          columnCounts[iColumnBlock]++;
3345     }
3346     if (numberRowBlocks == numberColumnBlocks ||
3347               numberRowBlocks == numberColumnBlocks + 1) {
3348          // could be Dantzig-Wolfe
3349          int numberG1 = 0;
3350          for (int i = 0; i < numberRowBlocks; i++) {
3351               if (rowCounts[i] > 1)
3352                    numberG1++;
3353          }
3354          bool masterColumns = (numberColumnBlocks == numberRowBlocks);
3355          if ((masterColumns && numberElementBlocks == 2 * numberRowBlocks - 1)
3356                    || (!masterColumns && numberElementBlocks == 2 * numberRowBlocks)) {
3357               if (numberG1 < 2)
3358                    decomposeType = 1;
3359          }
3360     }
3361     if (!decomposeType && (numberRowBlocks == numberColumnBlocks ||
3362                            numberRowBlocks == numberColumnBlocks - 1)) {
3363          // could be Benders
3364          int numberG1 = 0;
3365          for (int i = 0; i < numberColumnBlocks; i++) {
3366               if (columnCounts[i] > 1)
3367                    numberG1++;
3368          }
3369          bool masterRows = (numberColumnBlocks == numberRowBlocks);
3370          if ((masterRows && numberElementBlocks == 2 * numberColumnBlocks - 1)
3371                    || (!masterRows && numberElementBlocks == 2 * numberColumnBlocks)) {
3372               if (numberG1 < 2)
3373                    decomposeType = 2;
3374          }
3375     }
3376     delete [] rowCounts;
3377     delete [] columnCounts;
3378     delete [] blockInfo;
3379     // decide what to do
3380     switch (decomposeType) {
3381          // No good
3382     case 0:
3383          loadProblem(*model, false);
3384          return dual();
3385          // DW
3386     case 1:
3387          return solveDW(model);
3388          // Benders
3389     case 2:
3390          return solveBenders(model);
3391     }
3392     return 0; // to stop compiler warning
3393}
3394/* This loads a model from a CoinStructuredModel object - returns number of errors.
3395   If originalOrder then keep to order stored in blocks,
3396   otherwise first column/rows correspond to first block - etc.
3397   If keepSolution true and size is same as current then
3398   keeps current status and solution
3399*/
3400int
3401ClpSimplex::loadProblem (  CoinStructuredModel & coinModel,
3402                           bool originalOrder,
3403                           bool keepSolution)
3404{
3405     unsigned char * status = NULL;
3406     double * psol = NULL;
3407     double * dsol = NULL;
3408     int numberRows = coinModel.numberRows();
3409     int numberColumns = coinModel.numberColumns();
3410     int numberRowBlocks = coinModel.numberRowBlocks();
3411     int numberColumnBlocks = coinModel.numberColumnBlocks();
3412     int numberElementBlocks = coinModel.numberElementBlocks();
3413     if (status_ && numberRows_ && numberRows_ == numberRows &&
3414               numberColumns_ == numberColumns && keepSolution) {
3415          status = new unsigned char [numberRows_+numberColumns_];
3416          CoinMemcpyN(status_, numberRows_ + numberColumns_, status);
3417          psol = new double [numberRows_+numberColumns_];
3418          CoinMemcpyN(columnActivity_, numberColumns_, psol);
3419          CoinMemcpyN(rowActivity_, numberRows_, psol + numberColumns_);
3420          dsol = new double [numberRows_+numberColumns_];
3421          CoinMemcpyN(reducedCost_, numberColumns_, dsol);
3422          CoinMemcpyN(dual_, numberRows_, dsol + numberColumns_);
3423     }
3424     int returnCode = 0;
3425     double * rowLower = new double [numberRows];
3426     double * rowUpper = new double [numberRows];
3427     double * columnLower = new double [numberColumns];
3428     double * columnUpper = new double [numberColumns];
3429     double * objective = new double [numberColumns];
3430     int * integerType = new int [numberColumns];
3431     CoinBigIndex numberElements = 0;
3432     // Bases for blocks
3433     int * rowBase = new int[numberRowBlocks];
3434     CoinFillN(rowBase, numberRowBlocks, -1);
3435     // And row to put it
3436     int * whichRow = new int [numberRows];
3437     int * columnBase = new int[numberColumnBlocks];
3438     CoinFillN(columnBase, numberColumnBlocks, -1);
3439     // And column to put it
3440     int * whichColumn = new int [numberColumns];
3441     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
3442          CoinModel * block = coinModel.coinBlock(iBlock);
3443          numberElements += block->numberElements();
3444          //and set up elements etc
3445          double * associated = block->associatedArray();
3446          // If strings then do copies
3447          if (block->stringsExist())
3448               returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
3449                                                 objective, integerType, associated);
3450          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3451          int iRowBlock = info.rowBlock;
3452          int iColumnBlock = info.columnBlock;
3453          if (rowBase[iRowBlock] < 0) {
3454               rowBase[iRowBlock] = block->numberRows();
3455               // Save block number
3456               whichRow[numberRows-numberRowBlocks+iRowBlock] = iBlock;
3457          } else {
3458               assert(rowBase[iRowBlock] == block->numberRows());
3459          }
3460          if (columnBase[iColumnBlock] < 0) {
3461               columnBase[iColumnBlock] = block->numberColumns();
3462               // Save block number
3463               whichColumn[numberColumns-numberColumnBlocks+iColumnBlock] = iBlock;
3464          } else {
3465               assert(columnBase[iColumnBlock] == block->numberColumns());
3466          }
3467     }
3468     // Fill arrays with defaults
3469     CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
3470     CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
3471     CoinFillN(columnLower, numberColumns, 0.0);
3472     CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
3473     CoinFillN(objective, numberColumns, 0.0);
3474     CoinFillN(integerType, numberColumns, 0);
3475     int n = 0;
3476     for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
3477          int k = rowBase[iBlock];
3478          rowBase[iBlock] = n;
3479          assert (k >= 0);
3480          // block number
3481          int jBlock = whichRow[numberRows-numberRowBlocks+iBlock];
3482          if (originalOrder) {
3483               memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
3484          } else {
3485               CoinIotaN(whichRow + n, k, n);
3486          }
3487          n += k;
3488     }
3489     assert (n == numberRows);
3490     n = 0;
3491     for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
3492          int k = columnBase[iBlock];
3493          columnBase[iBlock] = n;
3494          assert (k >= 0);
3495          if (k) {
3496               // block number
3497               int jBlock = whichColumn[numberColumns-numberColumnBlocks+iBlock];
3498               if (originalOrder) {
3499                    memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
3500                           k * sizeof(int));
3501               } else {
3502                    CoinIotaN(whichColumn + n, k, n);
3503               }
3504               n += k;
3505          }
3506     }
3507     assert (n == numberColumns);
3508     bool gotIntegers = false;
3509     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
3510          CoinModel * block = coinModel.coinBlock(iBlock);
3511          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3512          int iRowBlock = info.rowBlock;
3513          int iRowBase = rowBase[iRowBlock];
3514          int iColumnBlock = info.columnBlock;
3515          int iColumnBase = columnBase[iColumnBlock];
3516          if (info.rhs) {
3517               int nRows = block->numberRows();
3518               const double * lower = block->rowLowerArray();
3519               const double * upper = block->rowUpperArray();
3520               for (int i = 0; i < nRows; i++) {
3521                    int put = whichRow[i+iRowBase];
3522                    rowLower[put] = lower[i];
3523                    rowUpper[put] = upper[i];
3524               }
3525          }
3526          if (info.bounds) {
3527               int nColumns = block->numberColumns();
3528               const double * lower = block->columnLowerArray();
3529               const double * upper = block->columnUpperArray();
3530               const double * obj = block->objectiveArray();
3531               for (int i = 0; i < nColumns; i++) {
3532                    int put = whichColumn[i+iColumnBase];
3533                    columnLower[put] = lower[i];
3534                    columnUpper[put] = upper[i];
3535                    objective[put] = obj[i];
3536               }
3537          }
3538          if (info.integer) {
3539               gotIntegers = true;
3540               int nColumns = block->numberColumns();
3541               const int * type = block->integerTypeArray();
3542               for (int i = 0; i < nColumns; i++) {
3543                    int put = whichColumn[i+iColumnBase];
3544                    integerType[put] = type[i];
3545               }
3546          }
3547     }
3548     gutsOfLoadModel(numberRows, numberColumns,
3549                     columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
3550     delete [] rowLower;
3551     delete [] rowUpper;
3552     delete [] columnLower;
3553     delete [] columnUpper;
3554     delete [] objective;
3555     // Do integers if wanted
3556     if (gotIntegers) {
3557          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3558               if (integerType[iColumn])
3559                    setInteger(iColumn);
3560          }
3561     }
3562     delete [] integerType;
3563     setObjectiveOffset(coinModel.objectiveOffset());
3564     // Space for elements
3565     int * row = new int[numberElements];
3566     int * column = new int[numberElements];
3567     double * element = new double[numberElements];
3568     numberElements = 0;
3569     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
3570          CoinModel * block = coinModel.coinBlock(iBlock);
3571          const CoinModelBlockInfo & info = coinModel.blockType(iBlock);
3572          int iRowBlock = info.rowBlock;
3573          int iRowBase = rowBase[iRowBlock];
3574          int iColumnBlock = info.columnBlock;
3575          int iColumnBase = columnBase[iColumnBlock];
3576          if (info.rowName) {
3577               int numberItems = block->rowNames()->numberItems();
3578               assert( block->numberRows() >= numberItems);
3579               if (numberItems) {
3580                    const char *const * rowNames = block->rowNames()->names();
3581                    for (int i = 0; i < numberItems; i++) {
3582                         int put = whichRow[i+iRowBase];
3583                         std::string name = rowNames[i];
3584                         setRowName(put, name);
3585                    }
3586               }
3587          }
3588          if (info.columnName) {
3589               int numberItems = block->columnNames()->numberItems();
3590               assert( block->numberColumns() >= numberItems);
3591               if (numberItems) {
3592                    const char *const * columnNames = block->columnNames()->names();
3593                    for (int i = 0; i < numberItems; i++) {
3594                         int put = whichColumn[i+iColumnBase];
3595                         std::string name = columnNames[i];
3596                         setColumnName(put, name);
3597                    }
3598               }
3599          }
3600          if (info.matrix) {
3601               CoinPackedMatrix matrix2;
3602               const CoinPackedMatrix * matrix = block->packedMatrix();
3603               if (!matrix) {
3604                    double * associated = block->associatedArray();
3605                    block->createPackedMatrix(matrix2, associated);
3606                    matrix = &matrix2;
3607               }
3608               // get matrix data pointers
3609               const int * row2 = matrix->getIndices();
3610               const CoinBigIndex * columnStart = matrix->getVectorStarts();
3611               const double * elementByColumn = matrix->getElements();
3612               const int * columnLength = matrix->getVectorLengths();
3613               int n = matrix->getNumCols();
3614               assert (matrix->isColOrdered());
3615               for (int iColumn = 0; iColumn < n; iColumn++) {
3616                    CoinBigIndex j;
3617                    int jColumn = whichColumn[iColumn+iColumnBase];
3618                    for (j = columnStart[iColumn];
3619                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3620                         row[numberElements] = whichRow[row2[j] + iRowBase];
3621                         column[numberElements] = jColumn;
3622                         element[numberElements++] = elementByColumn[j];
3623                    }
3624               }
3625          }
3626     }
3627     delete [] whichRow;
3628     delete [] whichColumn;
3629     delete [] rowBase;
3630     delete [] columnBase;
3631     CoinPackedMatrix * matrix =
3632          new CoinPackedMatrix (true, row, column, element, numberElements);
3633     matrix_ = new ClpPackedMatrix(matrix);
3634     matrix_->setDimensions(numberRows, numberColumns);
3635     delete [] row;
3636     delete [] column;
3637     delete [] element;
3638     createStatus();
3639     if (status) {
3640          // copy back
3641          CoinMemcpyN(status, numberRows_ + numberColumns_, status_);
3642          CoinMemcpyN(psol, numberColumns_, columnActivity_);
3643          CoinMemcpyN(psol + numberColumns_, numberRows_, rowActivity_);
3644          CoinMemcpyN(dsol, numberColumns_, reducedCost_);
3645          CoinMemcpyN(dsol + numberColumns_, numberRows_, dual_);
3646          delete [] status;
3647          delete [] psol;
3648          delete [] dsol;
3649     }
3650     optimizationDirection_ = coinModel.optimizationDirection();
3651     return returnCode;
3652}
3653/*  If input negative scales objective so maximum <= -value
3654    and returns scale factor used.  If positive unscales and also
3655    redoes dual stuff
3656*/
3657double
3658ClpSimplex::scaleObjective(double value)
3659{
3660     double * obj = objective();
3661     double largest = 0.0;
3662     if (value < 0.0) {
3663          value = - value;
3664          for (int i = 0; i < numberColumns_; i++) {
3665               largest = CoinMax(largest, fabs(obj[i]));
3666          }
3667          if (largest > value) {
3668               double scaleFactor = value / largest;
3669               for (int i = 0; i < numberColumns_; i++) {
3670                    obj[i] *= scaleFactor;
3671                    reducedCost_[i] *= scaleFactor;
3672               }
3673               for (int i = 0; i < numberRows_; i++) {
3674                    dual_[i] *= scaleFactor;
3675               }
3676               largest /= value;
3677          } else {
3678               // no need
3679               largest = 1.0;
3680          }
3681     } else {
3682          // at end
3683          if (value != 1.0) {
3684               for (int i = 0; i < numberColumns_; i++) {
3685                    obj[i] *= value;
3686                    reducedCost_[i] *= value;
3687               }
3688               for (int i = 0; i < numberRows_; i++) {
3689                    dual_[i] *= value;
3690               }
3691               computeObjectiveValue();
3692          }
3693     }
3694     return largest;
3695}
3696// Solve using Dantzig-Wolfe decomposition and maybe in parallel
3697int
3698ClpSimplex::solveDW(CoinStructuredModel * model)
3699{
3700     double time1 = CoinCpuTime();
3701     int numberColumns = model->numberColumns();
3702     int numberRowBlocks = model->numberRowBlocks();
3703     int numberColumnBlocks = model->numberColumnBlocks();
3704     int numberElementBlocks = model->numberElementBlocks();
3705     // We already have top level structure
3706     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
3707     for (int i = 0; i < numberElementBlocks; i++) {
3708          CoinModel * thisBlock = model->coinBlock(i);
3709          assert (thisBlock);
3710          // just fill in info
3711          CoinModelBlockInfo info = CoinModelBlockInfo();
3712          int whatsSet = thisBlock->whatIsSet();
3713          info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
3714          info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
3715          info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
3716          info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
3717          info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
3718          info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
3719          // Which block
3720          int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
3721          info.rowBlock = iRowBlock;
3722          int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
3723          info.columnBlock = iColumnBlock;
3724          blockInfo[i] = info;
3725     }
3726     // make up problems
3727     int numberBlocks = numberRowBlocks - 1;
3728     // Find master rows and columns
3729     int * rowCounts = new int [numberRowBlocks];
3730     CoinZeroN(rowCounts, numberRowBlocks);
3731     int * columnCounts = new int [numberColumnBlocks+1];
3732     CoinZeroN(columnCounts, numberColumnBlocks);
3733     int iBlock;
3734     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
3735          int iRowBlock = blockInfo[iBlock].rowBlock;
3736          rowCounts[iRowBlock]++;
3737          int iColumnBlock = blockInfo[iBlock].columnBlock;
3738          columnCounts[iColumnBlock]++;
3739     }
3740     int * whichBlock = new int [numberElementBlocks];
3741     int masterRowBlock = -1;
3742     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
3743          if (rowCounts[iBlock] > 1) {
3744               if (masterRowBlock == -1) {
3745                    masterRowBlock = iBlock;
3746               } else {
3747                    // Can't decode
3748                    masterRowBlock = -2;
3749                    break;
3750               }
3751          }
3752     }
3753     int masterColumnBlock = -1;
3754     int kBlock = 0;
3755     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
3756          int count = columnCounts[iBlock];
3757          columnCounts[iBlock] = kBlock;
3758          kBlock += count;
3759     }
3760     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
3761          int iColumnBlock = blockInfo[iBlock].columnBlock;
3762          whichBlock[columnCounts[iColumnBlock]] = iBlock;
3763          columnCounts[iColumnBlock]++;
3764     }
3765     for (iBlock = numberColumnBlocks - 1; iBlock >= 0; iBlock--)
3766          columnCounts[iBlock+1] = columnCounts[iBlock];
3767     columnCounts[0] = 0;
3768     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
3769          int count = columnCounts[iBlock+1] - columnCounts[iBlock];
3770          if (count == 1) {
3771               int kBlock = whichBlock[columnCounts[iBlock]];
3772               int iRowBlock = blockInfo[kBlock].rowBlock;
3773               if (iRowBlock == masterRowBlock) {
3774                    if (masterColumnBlock == -1) {
3775                         masterColumnBlock = iBlock;
3776                    } else {
3777                         // Can't decode
3778                         masterColumnBlock = -2;
3779                         break;
3780                    }
3781               }
3782          }
3783     }
3784     if (masterRowBlock < 0 || masterColumnBlock == -2) {
3785          // What now
3786          abort();
3787     }
3788     delete [] rowCounts;
3789     // create all data
3790     const CoinPackedMatrix ** top = new const CoinPackedMatrix * [numberColumnBlocks];
3791     ClpSimplex * sub = new ClpSimplex [numberBlocks];
3792     ClpSimplex master;
3793     // Set offset
3794     master.setObjectiveOffset(model->objectiveOffset());
3795     kBlock = 0;
3796     int masterBlock = -1;
3797     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
3798          top[kBlock] = NULL;
3799          int start = columnCounts[iBlock];
3800          int end = columnCounts[iBlock+1];
3801          assert (end - start <= 2);
3802          for (int j = start; j < end; j++) {
3803               int jBlock = whichBlock[j];
3804               int iRowBlock = blockInfo[jBlock].rowBlock;
3805               int iColumnBlock = blockInfo[jBlock].columnBlock;
3806               assert (iColumnBlock == iBlock);
3807               if (iColumnBlock != masterColumnBlock && iRowBlock == masterRowBlock) {
3808                    // top matrix
3809                    top[kBlock] = model->coinBlock(jBlock)->packedMatrix();
3810               } else {
3811                    const CoinPackedMatrix * matrix
3812                    = model->coinBlock(jBlock)->packedMatrix();
3813                    // Get pointers to arrays
3814                    const double * rowLower;
3815                    const double * rowUpper;
3816                    const double * columnLower;
3817                    const double * columnUpper;
3818                    const double * objective;
3819                    model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
3820                                 columnLower, columnUpper, objective);
3821                    if (iColumnBlock != masterColumnBlock) {
3822                         // diagonal block
3823                         sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
3824                                                 objective, rowLower, rowUpper);
3825                         if (true) {
3826                              double * lower = sub[kBlock].columnLower();
3827                              double * upper = sub[kBlock].columnUpper();
3828                              int n = sub[kBlock].numberColumns();
3829                              for (int i = 0; i < n; i++) {
3830                                   lower[i] = CoinMax(-1.0e8, lower[i]);
3831                                   upper[i] = CoinMin(1.0e8, upper[i]);
3832                              }
3833                         }
3834                         if (optimizationDirection_ < 0.0) {
3835                              double * obj = sub[kBlock].objective();
3836                              int n = sub[kBlock].numberColumns();
3837                              for (int i = 0; i < n; i++)
3838                                   obj[i] = - obj[i];
3839                         }
3840                         if (this->factorizationFrequency() == 200) {
3841                              // User did not touch preset
3842                              sub[kBlock].defaultFactorizationFrequency();
3843                         } else {
3844                              // make sure model has correct value
3845                              sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
3846                         }
3847                         sub[kBlock].setPerturbation(50);
3848                         // Set columnCounts to be diagonal block index for cleanup
3849                         columnCounts[kBlock] = jBlock;
3850                    } else {
3851                         // master
3852                         masterBlock = jBlock;
3853                         master.loadProblem(*matrix, columnLower, columnUpper,
3854                                            objective, rowLower, rowUpper);
3855                         if (optimizationDirection_ < 0.0) {
3856                              double * obj = master.objective();
3857                              int n = master.numberColumns();
3858                              for (int i = 0; i < n; i++)
3859                                   obj[i] = - obj[i];
3860                         }
3861                    }
3862               }
3863          }
3864          if (iBlock != masterColumnBlock)
3865               kBlock++;
3866     }
3867     delete [] whichBlock;
3868     delete [] blockInfo;
3869     // For now master must have been defined (does not have to have columns)
3870     assert (master.numberRows());
3871     assert (masterBlock >= 0);
3872     int numberMasterRows = master.numberRows();
3873     // Overkill in terms of space
3874     int spaceNeeded = CoinMax(numberBlocks * (numberMasterRows + 1),
3875                               2 * numberMasterRows);
3876     int * rowAdd = new int[spaceNeeded];
3877     double * elementAdd = new double[spaceNeeded];
3878     spaceNeeded = numberBlocks;
3879     int * columnAdd = new int[spaceNeeded+1];
3880     double * objective = new double[spaceNeeded];
3881     // Add in costed slacks
3882     int firstArtificial = master.numberColumns();
3883     int lastArtificial = firstArtificial;
3884     if (true) {
3885          const double * lower = master.rowLower();
3886          const double * upper = master.rowUpper();
3887          int kCol = 0;
3888          for (int iRow = 0; iRow < numberMasterRows; iRow++) {
3889               if (lower[iRow] > -1.0e10) {
3890                    rowAdd[kCol] = iRow;
3891                    elementAdd[kCol++] = 1.0;
3892               }
3893               if (upper[iRow] < 1.0e10) {
3894                    rowAdd[kCol] = iRow;
3895                    elementAdd[kCol++] = -1.0;
3896               }
3897          }
3898          if (kCol > spaceNeeded) {
3899               spaceNeeded = kCol;
3900               delete [] columnAdd;
3901               delete [] objective;
3902               columnAdd = new int[spaceNeeded+1];
3903               objective = new double[spaceNeeded];
3904          }
3905          for (int i = 0; i < kCol; i++) {
3906               columnAdd[i] = i;
3907               objective[i] = 1.0e13;
3908          }
3909          columnAdd[kCol] = kCol;
3910          master.addColumns(kCol, NULL, NULL, objective,
3911                            columnAdd, rowAdd, elementAdd);
3912          lastArtificial = master.numberColumns();
3913     }
3914     int maxPass = 500;
3915     int iPass;
3916     double lastObjective = 1.0e31;
3917     // Create convexity rows for proposals
3918     int numberMasterColumns = master.numberColumns();
3919     master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
3920     if (this->factorizationFrequency() == 200) {
3921          // User did not touch preset
3922          master.defaultFactorizationFrequency();
3923     } else {
3924          // make sure model has correct value
3925          master.setFactorizationFrequency(this->factorizationFrequency());
3926     }
3927     master.setPerturbation(50);
3928     // Arrays to say which block and when created
3929     int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
3930     whichBlock = new int[maximumColumns];
3931     int * when = new int[maximumColumns];
3932     int numberColumnsGenerated = numberBlocks;
3933     // fill in rhs and add in artificials
3934     {
3935          double * rowLower = master.rowLower();
3936          double * rowUpper = master.rowUpper();
3937          int iBlock;
3938          columnAdd[0] = 0;
3939          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
3940               int iRow = iBlock + numberMasterRows;;
3941               rowLower[iRow] = 1.0;
3942               rowUpper[iRow] = 1.0;
3943               rowAdd[iBlock] = iRow;
3944               elementAdd[iBlock] = 1.0;
3945               objective[iBlock] = 1.0e13;
3946               columnAdd[iBlock+1] = iBlock + 1;
3947               when[iBlock] = -1;
3948               whichBlock[iBlock] = iBlock;
3949          }
3950          master.addColumns(numberBlocks, NULL, NULL, objective,
3951                            columnAdd, rowAdd, elementAdd);
3952     }
3953     // and resize matrix to double check clp will be happy
3954     //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
3955     //                  numberMasterColumns+numberBlocks);
3956     std::cout << "Time to decompose " << CoinCpuTime() - time1 << " seconds" << std::endl;
3957     for (iPass = 0; iPass < maxPass; iPass++) {
3958          printf("Start of pass %d\n", iPass);
3959          // Solve master - may be infeasible
3960          //master.scaling(0);
3961          if (0) {
3962               master.writeMps("yy.mps");
3963          }
3964          // Correct artificials
3965          double sumArtificials = 0.0;
3966          if (iPass) {
3967               double * upper = master.columnUpper();
3968               double * solution = master.primalColumnSolution();
3969               double * obj = master.objective();
3970               sumArtificials = 0.0;
3971               for (int i = firstArtificial; i < lastArtificial; i++) {
3972                    sumArtificials += solution[i];
3973                    //assert (solution[i]>-1.0e-2);
3974                    if (solution[i] < 1.0e-6) {
3975#if 0
3976                         // Could take out
3977                         obj[i] = 0.0;
3978                         upper[i] = 0.0;
3979#else
3980                         obj[i] = 1.0e7;
3981                         upper[i] = 1.0e-1;
3982#endif
3983                         solution[i] = 0.0;
3984                         master.setColumnStatus(i, isFixed);
3985                    } else {
3986                         upper[i] = solution[i] + 1.0e-5 * (1.0 + solution[i]);
3987                    }
3988               }
3989               printf("Sum of artificials before solve is %g\n", sumArtificials);
3990          }
3991          // scale objective to be reasonable
3992          double scaleFactor = master.scaleObjective(-1.0e9);
3993          {
3994               double * dual = master.dualRowSolution();
3995               int n = master.numberRows();
3996               memset(dual, 0, n * sizeof(double));
3997               double * solution = master.primalColumnSolution();
3998               master.clpMatrix()->times(1.0, solution, dual);
3999               double sum = 0.0;
4000               double * lower = master.rowLower();
4001               double * upper = master.rowUpper();
4002               for (int iRow = 0; iRow < n; iRow++) {
4003                    double value = dual[iRow];
4004                    if (value > upper[iRow])
4005                         sum += value - upper[iRow];
4006                    else if (value < lower[iRow])
4007                         sum -= value - lower[iRow];
4008               }
4009               printf("suminf %g\n", sum);
4010               lower = master.columnLower();
4011               upper = master.columnUpper();
4012               n = master.numberColumns();
4013               for (int iColumn = 0; iColumn < n; iColumn++) {
4014                    double value = solution[iColumn];
4015                    if (value > upper[iColumn] + 1.0e-5)
4016                         sum += value - upper[iColumn];
4017                    else if (value < lower[iColumn] - 1.0e-5)
4018                         sum -= value - lower[iColumn];
4019               }
4020               printf("suminf %g\n", sum);
4021          }
4022          master.primal(1);
4023          // Correct artificials
4024          sumArtificials = 0.0;
4025          {
4026               double * solution = master.primalColumnSolution();
4027               for (int i = firstArtificial; i < lastArtificial; i++) {
4028                    sumArtificials += solution[i];
4029               }
4030               printf("Sum of artificials after solve is %g\n", sumArtificials);
4031          }
4032          master.scaleObjective(scaleFactor);
4033          int problemStatus = master.status(); // do here as can change (delcols)
4034          if (master.numberIterations() == 0 && iPass)
4035               break; // finished
4036          if (master.objectiveValue() > lastObjective - 1.0e-7 && iPass > 555)
4037               break; // finished
4038          lastObjective = master.objectiveValue();
4039          // mark basic ones and delete if necessary
4040          int iColumn;
4041          numberColumnsGenerated = master.numberColumns() - numberMasterColumns;
4042          for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
4043               if (master.getStatus(iColumn + numberMasterColumns) == ClpSimplex::basic)
4044                    when[iColumn] = iPass;
4045          }
4046          if (numberColumnsGenerated + numberBlocks > maximumColumns) {
4047               // delete
4048               int numberKeep = 0;
4049               int numberDelete = 0;
4050               int * whichDelete = new int[numberColumnsGenerated];
4051               for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
4052                    if (when[iColumn] > iPass - 7) {
4053                         // keep
4054                         when[numberKeep] = when[iColumn];
4055                         whichBlock[numberKeep++] = whichBlock[iColumn];
4056                    } else {
4057                         // delete
4058                         whichDelete[numberDelete++] = iColumn + numberMasterColumns;
4059                    }
4060               }
4061               numberColumnsGenerated -= numberDelete;
4062               master.deleteColumns(numberDelete, whichDelete);
4063               delete [] whichDelete;
4064          }
4065          const double * dual = NULL;
4066          bool deleteDual = false;
4067          if (problemStatus == 0) {
4068               dual = master.dualRowSolution();
4069          } else if (problemStatus == 1) {
4070               // could do composite objective
4071               dual = master.infeasibilityRay();
4072               deleteDual = true;
4073               printf("The sum of infeasibilities is %g\n",
4074                      master.sumPrimalInfeasibilities());
4075          } else if (!master.numberColumns()) {
4076               assert(!iPass);
4077               dual = master.dualRowSolution();
4078               memset(master.dualRowSolution(),
4079                      0, (numberMasterRows + numberBlocks)*sizeof(double));
4080          } else {
4081               abort();
4082          }
4083          // Scale back on first time
4084          if (!iPass) {
4085               double * dual2 = master.dualRowSolution();
4086               for (int i = 0; i < numberMasterRows + numberBlocks; i++) {
4087                    dual2[i] *= 1.0e-7;
4088               }
4089               dual = master.dualRowSolution();
4090          }
4091          // Create objective for sub problems and solve
4092          columnAdd[0] = 0;
4093          int numberProposals = 0;
4094          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4095               int numberColumns2 = sub[iBlock].numberColumns();
4096               double * saveObj = new double [numberColumns2];
4097               double * objective2 = sub[iBlock].objective();
4098               memcpy(saveObj, objective2, numberColumns2 * sizeof(double));
4099               // new objective
4100               top[iBlock]->transposeTimes(dual, objective2);
4101               int i;
4102               if (problemStatus == 0) {
4103                    for (i = 0; i < numberColumns2; i++)
4104                         objective2[i] = saveObj[i] - objective2[i];
4105               } else {
4106                    for (i = 0; i < numberColumns2; i++)
4107                         objective2[i] = -objective2[i];
4108               }
4109               // scale objective to be reasonable
4110               double scaleFactor =
4111                    sub[iBlock].scaleObjective((sumArtificials > 1.0e-5) ? -1.0e-4 : -1.0e9);
4112               if (iPass) {
4113                    sub[iBlock].primal();
4114               } else {
4115                    sub[iBlock].dual();
4116               }
4117               sub[iBlock].scaleObjective(scaleFactor);
4118               if (!sub[iBlock].isProvenOptimal() &&
4119                         !sub[iBlock].isProvenDualInfeasible()) {
4120                    memset(objective2, 0, numberColumns2 * sizeof(double));
4121                    sub[iBlock].primal();
4122                    if (problemStatus == 0) {
4123                         for (i = 0; i < numberColumns2; i++)
4124                              objective2[i] = saveObj[i] - objective2[i];
4125                    } else {
4126                         for (i = 0; i < numberColumns2; i++)
4127                              objective2[i] = -objective2[i];
4128                    }
4129                    double scaleFactor = sub[iBlock].scaleObjective(-1.0e9);
4130                    sub[iBlock].primal(1);
4131                    sub[iBlock].scaleObjective(scaleFactor);
4132               }
4133               memcpy(objective2, saveObj, numberColumns2 * sizeof(double));
4134               // get proposal
4135               if (sub[iBlock].numberIterations() || !iPass) {
4136                    double objValue = 0.0;
4137                    int start = columnAdd[numberProposals];
4138                    // proposal
4139                    if (sub[iBlock].isProvenOptimal()) {
4140                         const double * solution = sub[iBlock].primalColumnSolution();
4141                         top[iBlock]->times(solution, elementAdd + start);
4142                         for (i = 0; i < numberColumns2; i++)
4143                              objValue += solution[i] * saveObj[i];
4144                         // See if good dj and pack down
4145                         int number = start;
4146                         double dj = objValue;
4147                         if (problemStatus)
4148                              dj = 0.0;
4149                         double smallest = 1.0e100;
4150                         double largest = 0.0;
4151                         for (i = 0; i < numberMasterRows; i++) {
4152                              double value = elementAdd[start+i];
4153                              if (fabs(value) > 1.0e-15) {
4154                                   dj -= dual[i] * value;
4155                                   smallest = CoinMin(smallest, fabs(value));
4156                                   largest = CoinMax(largest, fabs(value));
4157                                   rowAdd[number] = i;
4158                                   elementAdd[number++] = value;
4159                              }
4160                         }
4161                         // and convexity
4162                         dj -= dual[numberMasterRows+iBlock];
4163                         rowAdd[number] = numberMasterRows + iBlock;
4164                         elementAdd[number++] = 1.0;
4165                         // if elements large then scale?
4166                         //if (largest>1.0e8||smallest<1.0e-8)
4167                         printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
4168                                iBlock, smallest, largest, dj);
4169                         if (dj < -1.0e-6 || !iPass) {
4170                              // take
4171                              objective[numberProposals] = objValue;
4172                              columnAdd[++numberProposals] = number;
4173                              when[numberColumnsGenerated] = iPass;
4174                              whichBlock[numberColumnsGenerated++] = iBlock;
4175                         }
4176                    } else if (sub[iBlock].isProvenDualInfeasible()) {
4177                         // use ray
4178                         const double * solution = sub[iBlock].unboundedRay();
4179                         top[iBlock]->times(solution, elementAdd + start);
4180                         for (i = 0; i < numberColumns2; i++)
4181                              objValue += solution[i] * saveObj[i];
4182                         // See if good dj and pack down
4183                         int number = start;
4184                         double dj = objValue;
4185                         double smallest = 1.0e100;
4186                         double largest = 0.0;
4187                         for (i = 0; i < numberMasterRows; i++) {
4188                              double value = elementAdd[start+i];
4189                              if (fabs(value) > 1.0e-15) {
4190                                   dj -= dual[i] * value;
4191                                   smallest = CoinMin(smallest, fabs(value));
4192                                   largest = CoinMax(largest, fabs(value));
4193                                   rowAdd[number] = i;
4194                                   elementAdd[number++] = value;
4195                              }
4196                         }
4197                         // if elements large or small then scale?
4198                         //if (largest>1.0e8||smallest<1.0e-8)
4199                         printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
4200                                iBlock, smallest, largest, dj);
4201                         if (dj < -1.0e-6) {
4202                              // take
4203                              objective[numberProposals] = objValue;
4204                              columnAdd[++numberProposals] = number;
4205                              when[numberColumnsGenerated] = iPass;
4206                              whichBlock[numberColumnsGenerated++] = iBlock;
4207                         }
4208                    } else {
4209                         abort();
4210                    }
4211               }
4212               delete [] saveObj;
4213          }
4214          if (deleteDual)
4215               delete [] dual;
4216          if (numberProposals)
4217               master.addColumns(numberProposals, NULL, NULL, objective,
4218                                 columnAdd, rowAdd, elementAdd);
4219     }
4220     std::cout << "Time at end of D-W " << CoinCpuTime() - time1 << " seconds" << std::endl;
4221     //master.scaling(0);
4222     //master.primal(1);
4223     loadProblem(*model);
4224     // now put back a good solution
4225     double * lower = new double[numberMasterRows+numberBlocks];
4226     double * upper = new double[numberMasterRows+numberBlocks];
4227     numberColumnsGenerated  += numberMasterColumns;
4228     double * sol = new double[numberColumnsGenerated];
4229     const double * solution = master.primalColumnSolution();
4230     const double * masterLower = master.rowLower();
4231     const double * masterUpper = master.rowUpper();
4232     double * fullSolution = primalColumnSolution();
4233     const double * fullLower = columnLower();
4234     const double * fullUpper = columnUpper();
4235     const double * rowSolution = master.primalRowSolution();
4236     double * fullRowSolution = primalRowSolution();
4237     const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4238     int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4239     const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4240     int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4241     for (int iRow = 0; iRow < numberRows2; iRow++) {
4242          int kRow = rowBack[iRow];
4243          setRowStatus(kRow, master.getRowStatus(iRow));
4244          fullRowSolution[kRow] = rowSolution[iRow];
4245     }
4246     for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4247          int kColumn = columnBack[iColumn];
4248          setStatus(kColumn, master.getStatus(iColumn));
4249          fullSolution[kColumn] = solution[iColumn];
4250     }
4251     for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4252          // move basis
4253          int kBlock = columnCounts[iBlock];
4254          const int * rowBack = model->coinBlock(kBlock)->originalRows();
4255          int numberRows2 = model->coinBlock(kBlock)->numberRows();
4256          const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4257          int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4258          for (int iRow = 0; iRow < numberRows2; iRow++) {
4259               int kRow = rowBack[iRow];
4260               setRowStatus(kRow, sub[iBlock].getRowStatus(iRow));
4261          }
4262          for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4263               int kColumn = columnBack[iColumn];
4264               setStatus(kColumn, sub[iBlock].getStatus(iColumn));
4265          }
4266          // convert top bit to by rows
4267          CoinPackedMatrix topMatrix = *top[iBlock];
4268          topMatrix.reverseOrdering();
4269          // zero solution
4270          memset(sol, 0, numberColumnsGenerated * sizeof(double));
4271
4272          for (int i = numberMasterColumns; i < numberColumnsGenerated; i++) {
4273               if (whichBlock[i-numberMasterColumns] == iBlock)
4274                    sol[i] = solution[i];
4275          }
4276          memset(lower, 0, (numberMasterRows + numberBlocks)*sizeof(double));
4277          master.clpMatrix()->times(1.0, sol, lower);
4278          for (int iRow = 0; iRow < numberMasterRows; iRow++) {
4279               double value = lower[iRow];
4280               if (masterUpper[iRow] < 1.0e20)
4281                    upper[iRow] = value;
4282               else
4283                    upper[iRow] = COIN_DBL_MAX;
4284               if (masterLower[iRow] > -1.0e20)
4285                    lower[iRow] = value;
4286               else
4287                    lower[iRow] = -COIN_DBL_MAX;
4288          }
4289          sub[iBlock].addRows(numberMasterRows, lower, upper,
4290                              topMatrix.getVectorStarts(),
4291                              topMatrix.getVectorLengths(),
4292                              topMatrix.getIndices(),
4293                              topMatrix.getElements());
4294          sub[iBlock].primal(1);
4295          const double * subSolution = sub[iBlock].primalColumnSolution();
4296          const double * subRowSolution = sub[iBlock].primalRowSolution();
4297          // move solution
4298          for (int iRow = 0; iRow < numberRows2; iRow++) {
4299               int kRow = rowBack[iRow];
4300               fullRowSolution[kRow] = subRowSolution[iRow];
4301          }
4302          for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4303               int kColumn = columnBack[iColumn];
4304               fullSolution[kColumn] = subSolution[iColumn];
4305          }
4306     }
4307     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4308          if (fullSolution[iColumn] < fullUpper[iColumn] - 1.0e-8 &&
4309                    fullSolution[iColumn] > fullLower[iColumn] + 1.0e-8) {
4310               if (getStatus(iColumn) != ClpSimplex::basic) {
4311                    if (columnLower_[iColumn] > -1.0e30 ||
4312                              columnUpper_[iColumn] < 1.0e30)
4313                         setStatus(iColumn, ClpSimplex::superBasic);
4314                    else
4315                         setStatus(iColumn, ClpSimplex::isFree);
4316               }
4317          } else if (fullSolution[iColumn] >= fullUpper[iColumn] - 1.0e-8) {
4318               // may help to make rest non basic
4319               if (getStatus(iColumn) != ClpSimplex::basic)
4320                    setStatus(iColumn, ClpSimplex::atUpperBound);
4321          } else if (fullSolution[iColumn] <= fullLower[iColumn] + 1.0e-8) {
4322               // may help to make rest non basic
4323               if (getStatus(iColumn) != ClpSimplex::basic)
4324                    setStatus(iColumn, ClpSimplex::atLowerBound);
4325          }
4326     }
4327     //int numberRows=model->numberRows();
4328     //for (int iRow=0;iRow<numberRows;iRow++)
4329     //setRowStatus(iRow,ClpSimplex::superBasic);
4330     std::cout << "Time before cleanup of full model " << CoinCpuTime() - time1 << " seconds" << std::endl;
4331     primal(1);
4332     std::cout << "Total time " << CoinCpuTime() - time1 << " seconds" << std::endl;
4333     delete [] columnCounts;
4334     delete [] sol;
4335     delete [] lower;
4336     delete [] upper;
4337     delete [] whichBlock;
4338     delete [] when;
4339     delete [] columnAdd;
4340     delete [] rowAdd;
4341     delete [] elementAdd;
4342     delete [] objective;
4343     delete [] top;
4344     delete [] sub;
4345     return 0;
4346}
4347// Solve using Benders decomposition and maybe in parallel
4348int
4349ClpSimplex::solveBenders(CoinStructuredModel * model)
4350{
4351     double time1 = CoinCpuTime();
4352     //int numberColumns = model->numberColumns();
4353     int numberRowBlocks = model->numberRowBlocks();
4354     int numberColumnBlocks = model->numberColumnBlocks();
4355     int numberElementBlocks = model->numberElementBlocks();
4356     // We already have top level structure
4357     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
4358     for (int i = 0; i < numberElementBlocks; i++) {
4359          CoinModel * thisBlock = model->coinBlock(i);
4360          assert (thisBlock);
4361          // just fill in info
4362          CoinModelBlockInfo info = CoinModelBlockInfo();
4363          int whatsSet = thisBlock->whatIsSet();
4364          info.matrix = static_cast<char>(((whatsSet & 1) != 0) ? 1 : 0);
4365          info.rhs = static_cast<char>(((whatsSet & 2) != 0) ? 1 : 0);
4366          info.rowName = static_cast<char>(((whatsSet & 4) != 0) ? 1 : 0);
4367          info.integer = static_cast<char>(((whatsSet & 32) != 0) ? 1 : 0);
4368          info.bounds = static_cast<char>(((whatsSet & 8) != 0) ? 1 : 0);
4369          info.columnName = static_cast<char>(((whatsSet & 16) != 0) ? 1 : 0);
4370          // Which block
4371          int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4372          info.rowBlock = iRowBlock;
4373          int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4374          info.columnBlock = iColumnBlock;
4375          blockInfo[i] = info;
4376     }
4377     // make up problems
4378     int numberBlocks = numberColumnBlocks - 1;
4379     // Find master columns and rows
4380     int * columnCounts = new int [numberColumnBlocks];
4381     CoinZeroN(columnCounts, numberColumnBlocks);
4382     int * rowCounts = new int [numberRowBlocks+1];
4383     CoinZeroN(rowCounts, numberRowBlocks);
4384     int iBlock;
4385     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4386          int iColumnBlock = blockInfo[iBlock].columnBlock;
4387          columnCounts[iColumnBlock]++;
4388          int iRowBlock = blockInfo[iBlock].rowBlock;
4389          rowCounts[iRowBlock]++;
4390     }
4391     int * whichBlock = new int [numberElementBlocks];
4392     int masterColumnBlock = -1;
4393     for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4394          if (columnCounts[iBlock] > 1) {
4395               if (masterColumnBlock == -1) {
4396                    masterColumnBlock = iBlock;
4397               } else {
4398                    // Can't decode
4399                    masterColumnBlock = -2;
4400                    break;
4401               }
4402          }
4403     }
4404     int masterRowBlock = -1;
4405     int kBlock = 0;
4406     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4407          int count = rowCounts[iBlock];
4408          rowCounts[iBlock] = kBlock;
4409          kBlock += count;
4410     }
4411     for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4412          int iRowBlock = blockInfo[iBlock].rowBlock;
4413          whichBlock[rowCounts[iRowBlock]] = iBlock;
4414          rowCounts[iRowBlock]++;
4415     }
4416     for (iBlock = numberRowBlocks - 1; iBlock >= 0; iBlock--)
4417          rowCounts[iBlock+1] = rowCounts[iBlock];
4418     rowCounts[0] = 0;
4419     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4420          int count = rowCounts[iBlock+1] - rowCounts[iBlock];
4421          if (count == 1) {
4422               int kBlock = whichBlock[rowCounts[iBlock]];
4423               int iColumnBlock = blockInfo[kBlock].columnBlock;
4424               if (iColumnBlock == masterColumnBlock) {
4425                    if (masterRowBlock == -1) {
4426                         masterRowBlock = iBlock;
4427                    } else {
4428                         // Can't decode
4429                         masterRowBlock = -2;
4430                         break;
4431                    }
4432               }
4433          }
4434     }
4435     if (masterColumnBlock < 0 || masterRowBlock == -2) {
4436          // What now
4437          abort();
4438     }
4439     delete [] columnCounts;
4440     // create all data
4441     const CoinPackedMatrix ** first = new const CoinPackedMatrix * [numberRowBlocks];
4442     ClpSimplex * sub = new ClpSimplex [numberBlocks];
4443     ClpSimplex master;
4444     // Set offset
4445     master.setObjectiveOffset(model->objectiveOffset());
4446     kBlock = 0;
4447     int masterBlock = -1;
4448     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4449          first[kBlock] = NULL;
4450          int start = rowCounts[iBlock];
4451          int end = rowCounts[iBlock+1];
4452          assert (end - start <= 2);
4453          for (int j = start; j < end; j++) {
4454               int jBlock = whichBlock[j];
4455               int iColumnBlock = blockInfo[jBlock].columnBlock;
4456               int iRowBlock = blockInfo[jBlock].rowBlock;
4457               assert (iRowBlock == iBlock);
4458               if (iRowBlock != masterRowBlock && iColumnBlock == masterColumnBlock) {
4459                    // first matrix
4460                    first[kBlock] = model->coinBlock(jBlock)->packedMatrix();
4461               } else {
4462                    const CoinPackedMatrix * matrix
4463                    = model->coinBlock(jBlock)->packedMatrix();
4464                    // Get pointers to arrays
4465                    const double * columnLower;
4466                    const double * columnUpper;
4467                    const double * rowLower;
4468                    const double * rowUpper;
4469                    const double * objective;
4470                    model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
4471                                 columnLower, columnUpper, objective);
4472                    if (iRowBlock != masterRowBlock) {
4473                         // diagonal block
4474                         sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
4475                                                 objective, rowLower, rowUpper);
4476                         if (optimizationDirection_ < 0.0) {
4477                              double * obj = sub[kBlock].objective();
4478                              int n = sub[kBlock].numberColumns();
4479                              for (int i = 0; i < n; i++)
4480                                   obj[i] = - obj[i];
4481                         }
4482                         if (this->factorizationFrequency() == 200) {
4483                              // User did not touch preset
4484                              sub[kBlock].defaultFactorizationFrequency();
4485                         } else {
4486                              // make sure model has correct value
4487                              sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
4488                         }
4489                         sub[kBlock].setPerturbation(50);
4490                         // Set rowCounts to be diagonal block index for cleanup
4491                         rowCounts[kBlock] = jBlock;
4492                    } else {
4493                         // master
4494                         masterBlock = jBlock;
4495                         master.loadProblem(*matrix, columnLower, columnUpper,
4496                                            objective, rowLower, rowUpper);
4497                         if (optimizationDirection_ < 0.0) {
4498                              double * obj = master.objective();
4499                              int n = master.numberColumns();
4500                              for (int i = 0; i < n; i++)
4501                                   obj[i] = - obj[i];
4502                         }
4503                    }
4504               }
4505          }
4506          if (iBlock != masterRowBlock)
4507               kBlock++;
4508     }
4509     delete [] whichBlock;
4510     delete [] blockInfo;
4511     // For now master must have been defined (does not have to have rows)
4512     assert (master.numberColumns());
4513     assert (masterBlock >= 0);
4514     int numberMasterColumns = master.numberColumns();
4515     // Overkill in terms of space
4516     int spaceNeeded = CoinMax(numberBlocks * (numberMasterColumns + 1),
4517                               2 * numberMasterColumns);
4518     int * columnAdd = new int[spaceNeeded];
4519     double * elementAdd = new double[spaceNeeded];
4520     spaceNeeded = numberBlocks;
4521     int * rowAdd = new int[spaceNeeded+1];
4522     double * objective = new double[spaceNeeded];
4523     int maxPass = 500;
4524     int iPass;
4525     double lastObjective = -1.0e31;
4526     // Create columns for proposals
4527     int numberMasterRows = master.numberRows();
4528     master.resize(numberMasterColumns + numberBlocks, numberMasterRows);
4529     if (this->factorizationFrequency() == 200) {
4530          // User did not touch preset
4531          master.defaultFactorizationFrequency();
4532     } else {
4533          // make sure model has correct value
4534          master.setFactorizationFrequency(this->factorizationFrequency());
4535     }
4536     master.setPerturbation(50);
4537     // Arrays to say which block and when created
4538     int maximumRows = 2 * numberMasterColumns + 10 * numberBlocks;
4539     whichBlock = new int[maximumRows];
4540     int * when = new int[maximumRows];
4541     int numberRowsGenerated = numberBlocks;
4542     // Add extra variables
4543     {
4544          int iBlock;
4545          columnAdd[0] = 0;
4546          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4547               objective[iBlock] = 1.0;
4548               columnAdd[iBlock+1] = 0;
4549               when[iBlock] = -1;
4550               whichBlock[iBlock] = iBlock;
4551          }
4552          master.addColumns(numberBlocks, NULL, NULL, objective,
4553                            columnAdd, rowAdd, elementAdd);
4554     }
4555     std::cout << "Time to decompose " << CoinCpuTime() - time1 << " seconds" << std::endl;
4556     for (iPass = 0; iPass < maxPass; iPass++) {
4557          printf("Start of pass %d\n", iPass);
4558          // Solve master - may be unbounded
4559          //master.scaling(0);
4560          if (1) {
4561               master.writeMps("yy.mps");
4562          }
4563          master.dual();
4564          int problemStatus = master.status(); // do here as can change (delcols)
4565          if (master.numberIterations() == 0 && iPass)
4566               break; // finished
4567          if (master.objectiveValue() < lastObjective + 1.0e-7 && iPass > 555)
4568               break; // finished
4569          lastObjective = master.objectiveValue();
4570          // mark non-basic rows and delete if necessary
4571          int iRow;
4572          numberRowsGenerated = master.numberRows() - numberMasterRows;
4573          for (iRow = 0; iRow < numberRowsGenerated; iRow++) {
4574               if (master.getStatus(iRow + numberMasterRows) != ClpSimplex::basic)
4575                    when[iRow] = iPass;
4576          }
4577          if (numberRowsGenerated > maximumRows) {
4578               // delete
4579               int numberKeep = 0;
4580               int numberDelete = 0;
4581               int * whichDelete = new int[numberRowsGenerated];
4582               for (iRow = 0; iRow < numberRowsGenerated; iRow++) {
4583                    if (when[iRow] > iPass - 7) {
4584                         // keep
4585                         when[numberKeep] = when[iRow];
4586                         whichBlock[numberKeep++] = whichBlock[iRow];
4587                    } else {
4588                         // delete
4589                         whichDelete[numberDelete++] = iRow + numberMasterRows;
4590                    }
4591               }
4592               numberRowsGenerated -= numberDelete;
4593               master.deleteRows(numberDelete, whichDelete);
4594               delete [] whichDelete;
4595          }
4596          const double * primal = NULL;
4597          bool deletePrimal = false;
4598          if (problemStatus == 0) {
4599               primal = master.primalColumnSolution();
4600          } else if (problemStatus == 2) {
4601               // could do composite objective
4602               primal = master.unboundedRay();
4603               deletePrimal = true;
4604               printf("The sum of infeasibilities is %g\n",
4605                      master.sumPrimalInfeasibilities());
4606          } else if (!master.numberRows()) {
4607               assert(!iPass);
4608               primal = master.primalColumnSolution();
4609               memset(master.primalColumnSolution(),
4610                      0, numberMasterColumns * sizeof(double));
4611          } else {
4612               abort();
4613          }
4614          // Create rhs for sub problems and solve
4615          rowAdd[0] = 0;
4616          int numberProposals = 0;
4617          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4618               int numberRows2 = sub[iBlock].numberRows();
4619               double * saveLower = new double [numberRows2];
4620               double * lower2 = sub[iBlock].rowLower();
4621               double * saveUpper = new double [numberRows2];
4622               double * upper2 = sub[iBlock].rowUpper();
4623               // new rhs
4624               CoinZeroN(saveUpper, numberRows2);
4625               first[iBlock]->times(primal, saveUpper);
4626               for (int i = 0; i < numberRows2; i++) {
4627                    double value = saveUpper[i];
4628                    saveLower[i] = lower2[i];
4629                    saveUpper[i] = upper2[i];
4630                    if (saveLower[i] > -1.0e30)
4631                         lower2[i] -= value;
4632                    if (saveUpper[i] < 1.0e30)
4633                         upper2[i] -= value;
4634               }
4635               sub[iBlock].dual();
4636               memcpy(lower2, saveLower, numberRows2 * sizeof(double));
4637               memcpy(upper2, saveUpper, numberRows2 * sizeof(double));
4638               // get proposal
4639               if (sub[iBlock].numberIterations() || !iPass) {
4640                    double objValue = 0.0;
4641                    int start = rowAdd[numberProposals];
4642                    // proposal
4643                    if (sub[iBlock].isProvenOptimal()) {
4644                         const double * solution = sub[iBlock].dualRowSolution();
4645                         first[iBlock]->transposeTimes(solution, elementAdd + start);
4646                         for (int i = 0; i < numberRows2; i++) {
4647                              if (solution[i] < -dualTolerance_) {
4648                                   // at upper
4649                                   assert (saveUpper[i] < 1.0e30);
4650                                   objValue += solution[i] * saveUpper[i];
4651                              } else if (solution[i] > dualTolerance_) {
4652                                   // at lower
4653                                   assert (saveLower[i] > -1.0e30);
4654                                   objValue += solution[i] * saveLower[i];
4655                              }
4656                         }
4657
4658                         // See if cuts off and pack down
4659                         int number = start;
4660                         double infeas = objValue;
4661                         double smallest = 1.0e100;
4662                         double largest = 0.0;
4663                         for (int i = 0; i < numberMasterColumns; i++) {
4664                              double value = elementAdd[start+i];
4665                              if (fabs(value) > 1.0e-15) {
4666                                   infeas -= primal[i] * value;
4667                                   smallest = CoinMin(smallest, fabs(value));
4668                                   largest = CoinMax(largest, fabs(value));
4669                                   columnAdd[number] = i;
4670                                   elementAdd[number++] = -value;
4671                              }
4672                         }
4673                         columnAdd[number] = numberMasterColumns + iBlock;
4674                         elementAdd[number++] = -1.0;
4675                         // if elements large then scale?
4676                         //if (largest>1.0e8||smallest<1.0e-8)
4677                         printf("For subproblem %d smallest - %g, largest %g - infeas %g\n",
4678                                iBlock, smallest, largest, infeas);
4679                         if (infeas < -1.0e-6 || !iPass) {
4680                              // take
4681                              objective[numberProposals] = objValue;
4682                              rowAdd[++numberProposals] = number;
4683                              when[numberRowsGenerated] = iPass;
4684                              whichBlock[numberRowsGenerated++] = iBlock;
4685                         }
4686                    } else if (sub[iBlock].isProvenPrimalInfeasible()) {
4687                         // use ray
4688                         const double * solution = sub[iBlock].infeasibilityRay();
4689                         first[iBlock]->transposeTimes(solution, elementAdd + start);
4690                         for (int i = 0; i < numberRows2; i++) {
4691                              if (solution[i] < -dualTolerance_) {
4692                                   // at upper
4693                                   assert (saveUpper[i] < 1.0e30);
4694                                   objValue += solution[i] * saveUpper[i];
4695                              } else if (solution[i] > dualTolerance_) {
4696                                   // at lower
4697                                   assert (saveLower[i] > -1.0e30);
4698                                   objValue += solution[i] * saveLower[i];
4699                              }
4700                         }
4701                         // See if good infeas and pack down
4702                         int number = start;
4703                         double infeas = objValue;
4704                         double smallest = 1.0e100;
4705                         double largest = 0.0;
4706                         for (int i = 0; i < numberMasterColumns; i++) {
4707                              double value = elementAdd[start+i];
4708                              if (fabs(value) > 1.0e-15) {
4709                                   infeas -= primal[i] * value;
4710                                   smallest = CoinMin(smallest, fabs(value));
4711                                   largest = CoinMax(largest, fabs(value));
4712                                   columnAdd[number] = i;
4713                                   elementAdd[number++] = -value;
4714                              }
4715                         }
4716                         // if elements large or small then scale?
4717                         //if (largest>1.0e8||smallest<1.0e-8)
4718                         printf("For subproblem ray %d smallest - %g, largest %g - infeas %g\n",
4719                                iBlock, smallest, largest, infeas);
4720                         if (infeas < -1.0e-6) {
4721                              // take
4722                              objective[numberProposals] = objValue;
4723                              rowAdd[++numberProposals] = number;
4724                              when[numberRowsGenerated] = iPass;
4725                              whichBlock[numberRowsGenerated++] = iBlock;
4726                         }
4727                    } else {
4728                         abort();
4729                    }
4730               }
4731               delete [] saveLower;
4732               delete [] saveUpper;
4733          }
4734          if (deletePrimal)
4735               delete [] primal;
4736          if (numberProposals) {
4737               master.addRows(numberProposals, NULL, objective,
4738                              rowAdd, columnAdd, elementAdd);
4739          }
4740     }
4741     std::cout << "Time at end of Benders " << CoinCpuTime() - time1 << " seconds" << std::endl;
4742     //master.scaling(0);
4743     //master.primal(1);
4744     loadProblem(*model);
4745     // now put back a good solution
4746     const double * columnSolution = master.primalColumnSolution();
4747     double * fullColumnSolution = primalColumnSolution();
4748     const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
4749     int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
4750     const int * rowBack = model->coinBlock(masterBlock)->originalRows();
4751     int numberRows2 = model->coinBlock(masterBlock)->numberRows();
4752     for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4753          int kColumn = columnBack[iColumn];
4754          setColumnStatus(kColumn, master.getColumnStatus(iColumn));
4755          fullColumnSolution[kColumn] = columnSolution[iColumn];
4756     }
4757     for (int iRow = 0; iRow < numberRows2; iRow++) {
4758          int kRow = rowBack[iRow];
4759          setStatus(kRow, master.getStatus(iRow));
4760          //fullSolution[kRow]=solution[iRow];
4761     }
4762     for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
4763          // move basis
4764          int kBlock = rowCounts[iBlock];
4765          const int * columnBack = model->coinBlock(kBlock)->originalColumns();
4766          int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
4767          const int * rowBack = model->coinBlock(kBlock)->originalRows();
4768          int numberRows2 = model->coinBlock(kBlock)->numberRows();
4769          const double * subColumnSolution = sub[iBlock].primalColumnSolution();
4770          for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
4771               int kColumn = columnBack[iColumn];
4772               setColumnStatus(kColumn, sub[iBlock].getColumnStatus(iColumn));
4773               fullColumnSolution[kColumn] = subColumnSolution[iColumn];
4774          }
4775          for (int iRow = 0; iRow < numberRows2; iRow++) {
4776               int kRow = rowBack[iRow];
4777               setStatus(kRow, sub[iBlock].getStatus(iRow));
4778               setStatus(kRow, atLowerBound);
4779          }
4780     }
4781     double * fullSolution = primalRowSolution();
4782     CoinZeroN(fullSolution, numberRows_);
4783     times(1.0, fullColumnSolution, fullSolution);
4784     //int numberColumns=model->numberColumns();
4785     //for (int iColumn=0;iColumn<numberColumns;iColumn++)
4786     //setColumnStatus(iColumn,ClpSimplex::superBasic);
4787     std::cout << "Time before cleanup of full model " << CoinCpuTime() - time1 << " seconds" << std::endl;
4788     this->primal(1);
4789     std::cout << "Total time " << CoinCpuTime() - time1 << " seconds" << std::endl;
4790     delete [] rowCounts;
4791     //delete [] sol;
4792     //delete [] lower;
4793     //delete [] upper;
4794     delete [] whichBlock;
4795     delete [] when;
4796     delete [] rowAdd;
4797     delete [] columnAdd;
4798     delete [] elementAdd;
4799     delete [] objective;
4800     delete [] first;
4801     delete [] sub;
4802     return 0;
4803}
Note: See TracBrowser for help on using the repository browser.