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

Last change on this file since 1691 was 1691, checked in by stefan, 9 years ago

check for AMD and CHOLMOD and GLPK's AMD in configure; compile ClpCholeskyUfl? only if AMD, CHOLMOD, or GLPK available

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