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