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

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

save some memory

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