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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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