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

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

make WSMP recoginition and use similar to Ipopt's WSMP interface

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