1/* $Id: ClpDynamicMatrix.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2004, 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
7#include <cstdio>
8
9#include "CoinPragma.hpp"
10#include "CoinIndexedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12
13#include "ClpSimplex.hpp"
14#include "ClpFactorization.hpp"
15#include "ClpQuadraticObjective.hpp"
16#include "ClpNonLinearCost.hpp"
17// at end to get min/max!
18#include "ClpDynamicMatrix.hpp"
19#include "ClpMessage.hpp"
20//#define CLP_DEBUG
21//#define CLP_DEBUG_PRINT
22//#############################################################################
23// Constructors / Destructor / Assignment
24//#############################################################################
25
26//-------------------------------------------------------------------
27// Default Constructor
28//-------------------------------------------------------------------
29ClpDynamicMatrix::ClpDynamicMatrix ()
30 : ClpPackedMatrix(),
31 sumDualInfeasibilities_(0.0),
32 sumPrimalInfeasibilities_(0.0),
33 sumOfRelaxedDualInfeasibilities_(0.0),
34 sumOfRelaxedPrimalInfeasibilities_(0.0),
35 savedBestGubDual_(0.0),
36 savedBestSet_(0),
37 backToPivotRow_(NULL),
38 keyVariable_(NULL),
39 toIndex_(NULL),
40 fromIndex_(NULL),
41 numberSets_(0),
42 numberActiveSets_(0),
43 objectiveOffset_(0.0),
44 lowerSet_(NULL),
45 upperSet_(NULL),
46 status_(NULL),
47 model_(NULL),
48 firstAvailable_(0),
49 firstAvailableBefore_(0),
50 firstDynamic_(0),
51 lastDynamic_(0),
52 numberStaticRows_(0),
53 numberElements_(0),
54 numberDualInfeasibilities_(0),
55 numberPrimalInfeasibilities_(0),
56 noCheck_(-1),
57 infeasibilityWeight_(0.0),
58 numberGubColumns_(0),
59 maximumGubColumns_(0),
60 maximumElements_(0),
61 startSet_(NULL),
62 next_(NULL),
63 startColumn_(NULL),
64 row_(NULL),
65 element_(NULL),
66 cost_(NULL),
67 id_(NULL),
68 dynamicStatus_(NULL),
69 columnLower_(NULL),
70 columnUpper_(NULL)
71{
72 setType(15);
73}
74
75//-------------------------------------------------------------------
76// Copy constructor
77//-------------------------------------------------------------------
78ClpDynamicMatrix::ClpDynamicMatrix (const ClpDynamicMatrix & rhs)
79 : ClpPackedMatrix(rhs)
80{
81 objectiveOffset_ = rhs.objectiveOffset_;
82 numberSets_ = rhs.numberSets_;
83 numberActiveSets_ = rhs.numberActiveSets_;
84 firstAvailable_ = rhs.firstAvailable_;
85 firstAvailableBefore_ = rhs.firstAvailableBefore_;
86 firstDynamic_ = rhs.firstDynamic_;
87 lastDynamic_ = rhs.lastDynamic_;
88 numberStaticRows_ = rhs.numberStaticRows_;
89 numberElements_ = rhs.numberElements_;
90 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_);
91 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
92 toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
93 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_);
94 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
95 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
96 status_ = ClpCopyOfArray(rhs.status_, static_cast<int>(2*numberSets_+4*sizeof(int)));
97 model_ = rhs.model_;
98 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
99 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
100 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
101 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
102 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
103 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
104 savedBestGubDual_ = rhs.savedBestGubDual_;
105 savedBestSet_ = rhs.savedBestSet_;
106 noCheck_ = rhs.noCheck_;
107 infeasibilityWeight_ = rhs.infeasibilityWeight_;
108 // Now secondary data
109 numberGubColumns_ = rhs.numberGubColumns_;
110 maximumGubColumns_ = rhs.maximumGubColumns_;
111 maximumElements_ = rhs.maximumElements_;
112 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
113 next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
114 startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
115 row_ = ClpCopyOfArray(rhs.row_, maximumElements_);
116 element_ = ClpCopyOfArray(rhs.element_, maximumElements_);
117 cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_);
118 id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
119 columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
120 columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
121 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
122}
123
124/* This is the real constructor*/
125ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex * model, int numberSets,
126 int numberGubColumns, const int * starts,
127 const double * lower, const double * upper,
128 const CoinBigIndex * startColumn, const int * row,
129 const double * element, const double * cost,
130 const double * columnLower, const double * columnUpper,
131 const unsigned char * status,
132 const unsigned char * dynamicStatus)
133 : ClpPackedMatrix()
134{
135 setType(15);
136 objectiveOffset_ = model->objectiveOffset();
137 model_ = model;
138 numberSets_ = numberSets;
139 numberGubColumns_ = numberGubColumns;
140 maximumGubColumns_ = numberGubColumns_;
141 if (numberGubColumns_)
142 maximumElements_ = startColumn[numberGubColumns_];
143 else
144 maximumElements_ = 0;
145 startSet_ = new int [numberSets_+1];
146 next_ = new int [maximumGubColumns_];
147 // fill in startSet and next
148 int iSet;
149 if (numberGubColumns_) {
150 for (iSet = 0; iSet < numberSets_; iSet++) {
151 int first = starts[iSet];
152 int last = starts[iSet+1] - 1;
153 startSet_[iSet] = first;
154 for (int i = first; i < last; i++)
155 next_[i] = i + 1;
156 next_[last] = -iSet - 1;
157 }
158 }
159 startSet_[numberSets_] = starts[numberSets_];
160 int numberColumns = model->numberColumns();
161 int numberRows = model->numberRows();
162 numberStaticRows_ = numberRows;
163 savedBestGubDual_ = 0.0;
164 savedBestSet_ = 0;
165 // Number of columns needed
166 int frequency = model->factorizationFrequency();
167 int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4;
168 // But we may have two per row + one for incoming (make it two)
169 numberGubInSmall = CoinMax(2*numberRows+2,numberGubInSmall);
170 // for small problems this could be too big
171 //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
172 int numberNeeded = numberGubInSmall + numberColumns;
173 firstAvailable_ = numberColumns;
174 firstAvailableBefore_ = firstAvailable_;
175 firstDynamic_ = numberColumns;
176 lastDynamic_ = numberNeeded;
177 startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
178 if (!numberGubColumns_) {
179 //startColumn_ = new CoinBigIndex [1];
180 startColumn_[0] = 0;
181 }
182 CoinBigIndex numberElements = startColumn_[numberGubColumns_];
183 row_ = ClpCopyOfArray(row, numberElements);
184 element_ = new double[numberElements];
185 CoinBigIndex i;
186 for (i = 0; i < numberElements; i++)
187 element_[i] = element[i];
188 cost_ = new double[numberGubColumns_];
189 for (i = 0; i < numberGubColumns_; i++) {
190 cost_[i] = cost[i];
191 // I don't think I need sorted but ...
192 CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i+1], element_ + startColumn_[i]);
193 }
194 if (columnLower) {
195 columnLower_ = new double[numberGubColumns_];
196 for (i = 0; i < numberGubColumns_; i++)
197 columnLower_[i] = columnLower[i];
198 } else {
199 columnLower_ = NULL;
200 }
201 if (columnUpper) {
202 columnUpper_ = new double[numberGubColumns_];
203 for (i = 0; i < numberGubColumns_; i++)
204 columnUpper_[i] = columnUpper[i];
205 } else {
206 columnUpper_ = NULL;
207 }
208 lowerSet_ = new double[numberSets_];
209 for (i = 0; i < numberSets_; i++) {
210 if (lower[i] > -1.0e20)
211 lowerSet_[i] = lower[i];
212 else
213 lowerSet_[i] = -1.0e30;
214 }
215 upperSet_ = new double[numberSets_];
216 for (i = 0; i < numberSets_; i++) {
217 if (upper[i] < 1.0e20)
218 upperSet_[i] = upper[i];
219 else
220 upperSet_[i] = 1.0e30;
221 }
222 id_ = new int[numberGubInSmall];
223 for (i = 0; i < numberGubInSmall; i++)
224 id_[i] = -1;
225 ClpPackedMatrix* originalMatrixA =
226 dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
227 assert (originalMatrixA);
228 CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
229 originalMatrixA->setMatrixNull(); // so can be deleted safely
230 // guess how much space needed
231 double guess = numberElements;
232 guess /= static_cast<double> (numberColumns);
233 guess *= 2 * numberGubInSmall;
234 numberElements_ = static_cast<int> (guess);
235 numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
236 matrix_ = originalMatrix;
237 //delete originalMatrixA;
238 flags_ &= ~1;
239 // resize model (matrix stays same)
240 // modify frequency
241 if (frequency>=50)
242 frequency = 50+(frequency-50)/2;
243 int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1;
244 model->resize(newRowSize, numberNeeded);
245 for (i = numberRows; i < newRowSize; i++)
246 model->setRowStatus(i, ClpSimplex::basic);
247 if (columnUpper_) {
248 // set all upper bounds so we have enough space
249 double * columnUpper = model->columnUpper();
250 for(i = firstDynamic_; i < lastDynamic_; i++)
251 columnUpper[i] = 1.0e10;
252 }
253 // resize matrix
254 // extra 1 is so can keep number of elements handy
255 originalMatrix->reserve(numberNeeded, numberElements_, true);
256 originalMatrix->reserve(numberNeeded + 1, numberElements_, false);
257 originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements();
258 originalMatrix->setDimensions(newRowSize, -1);
259 numberActiveColumns_ = firstDynamic_;
260 // redo number of columns
261 numberColumns = matrix_->getNumCols();
262 backToPivotRow_ = new int[numberNeeded];
263 keyVariable_ = new int[numberSets_];
264 if (status) {
265 status_ = ClpCopyOfArray(status, static_cast<int>(2*numberSets_+4*sizeof(int)));
266 assert (dynamicStatus);
267 dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_);
268 } else {
269 status_ = new unsigned char [2*numberSets_+4*sizeof(int)];
270 memset(status_, 0, numberSets_);
271 int i;
272 for (i = 0; i < numberSets_; i++) {
273 // make slack key
274 setStatus(i, ClpSimplex::basic);
275 }
276 dynamicStatus_ = new unsigned char [2*numberGubColumns_];
277 memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
278 for (i = 0; i < numberGubColumns_; i++)
279 setDynamicStatus(i, atLowerBound);
280 }
281 toIndex_ = new int[numberSets_];
282 for (iSet = 0; iSet < numberSets_; iSet++)
283 toIndex_[iSet] = -1;
284 fromIndex_ = new int [newRowSize-numberStaticRows_+1];
285 numberActiveSets_ = 0;
286 rhsOffset_ = NULL;
287 if (numberGubColumns_) {
288 if (!status) {
289 gubCrash();
290 } else {
291 initialProblem();
292 }
293 }
294 noCheck_ = -1;
295 infeasibilityWeight_ = 0.0;
296}
297
298//-------------------------------------------------------------------
299// Destructor
300//-------------------------------------------------------------------
301ClpDynamicMatrix::~ClpDynamicMatrix ()
302{
303 delete [] backToPivotRow_;
304 delete [] keyVariable_;
305 delete [] toIndex_;
306 delete [] fromIndex_;
307 delete [] lowerSet_;
308 delete [] upperSet_;
309 delete [] status_;
310 delete [] startSet_;
311 delete [] next_;
312 delete [] startColumn_;
313 delete [] row_;
314 delete [] element_;
315 delete [] cost_;
316 delete [] id_;
317 delete [] dynamicStatus_;
318 delete [] columnLower_;
319 delete [] columnUpper_;
320}
321
322//----------------------------------------------------------------
323// Assignment operator
324//-------------------------------------------------------------------
325ClpDynamicMatrix &
326ClpDynamicMatrix::operator=(const ClpDynamicMatrix& rhs)
327{
328 if (this != &rhs) {
329 ClpPackedMatrix::operator=(rhs);
330 delete [] backToPivotRow_;
331 delete [] keyVariable_;
332 delete [] toIndex_;
333 delete [] fromIndex_;
334 delete [] lowerSet_;
335 delete [] upperSet_;
336 delete [] status_;
337 delete [] startSet_;
338 delete [] next_;
339 delete [] startColumn_;
340 delete [] row_;
341 delete [] element_;
342 delete [] cost_;
343 delete [] id_;
344 delete [] dynamicStatus_;
345 delete [] columnLower_;
346 delete [] columnUpper_;
347 objectiveOffset_ = rhs.objectiveOffset_;
348 numberSets_ = rhs.numberSets_;
349 numberActiveSets_ = rhs.numberActiveSets_;
350 firstAvailable_ = rhs.firstAvailable_;
351 firstAvailableBefore_ = rhs.firstAvailableBefore_;
352 firstDynamic_ = rhs.firstDynamic_;
353 lastDynamic_ = rhs.lastDynamic_;
354 numberStaticRows_ = rhs.numberStaticRows_;
355 numberElements_ = rhs.numberElements_;
356 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_);
357 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
358 toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
359 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_);
360 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
361 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
362 status_ = ClpCopyOfArray(rhs.status_, static_cast<int>(2*numberSets_+4*sizeof(int)));
363 model_ = rhs.model_;
364 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
365 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
366 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
367 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
368 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
369 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
370 savedBestGubDual_ = rhs.savedBestGubDual_;
371 savedBestSet_ = rhs.savedBestSet_;
372 noCheck_ = rhs.noCheck_;
373 infeasibilityWeight_ = rhs.infeasibilityWeight_;
374 // Now secondary data
375 numberGubColumns_ = rhs.numberGubColumns_;
376 maximumGubColumns_ = rhs.maximumGubColumns_;
377 maximumElements_ = rhs.maximumElements_;
378 startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
379 next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
380 startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
381 row_ = ClpCopyOfArray(rhs.row_, maximumElements_);
382 element_ = ClpCopyOfArray(rhs.element_, maximumElements_);
383 cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_);
384 id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
385 columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
386 columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
387 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
388 }
389 return *this;
390}
391//-------------------------------------------------------------------
392// Clone
393//-------------------------------------------------------------------
394ClpMatrixBase * ClpDynamicMatrix::clone() const
395{
396 return new ClpDynamicMatrix(*this);
397}
398// Partial pricing
399void
400ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
401 int & bestSequence, int & numberWanted)
402{
403 numberWanted = currentWanted_;
404 assert(!model->rowScale());
405 if (numberSets_) {
406 // Do packed part before gub
407 // always???
408 //printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
409 ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
410 } else {
411 // no gub
412 ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
413 return;
414 }
415 if (numberWanted > 0) {
416 // and do some proportion of full set
417 int startG2 = static_cast<int> (startFraction * numberSets_);
418 int endG2 = static_cast<int> (endFraction * numberSets_ + 0.1);
419 endG2 = CoinMin(endG2, numberSets_);
420 //printf("gub price - set start %d end %d\n",
421 // startG2,endG2);
422 double tolerance = model->currentDualTolerance();
423 double * reducedCost = model->djRegion();
424 const double * duals = model->dualRowSolution();
425 double bestDj;
426 int numberRows = model->numberRows();
427 int slackOffset = lastDynamic_ + numberRows;
428 int structuralOffset = slackOffset + numberSets_;
429 // If nothing found yet can go all the way to end
430 int endAll = endG2;
431 if (bestSequence < 0 && !startG2)
432 endAll = numberSets_;
433 if (bestSequence >= 0) {
434 if (bestSequence != savedBestSequence_)
435 bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent
436 else
437 bestDj = savedBestDj_;
438 } else {
439 bestDj = tolerance;
440 }
441 int saveSequence = bestSequence;
442 double djMod = 0.0;
443 double bestDjMod = 0.0;
444 //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
445 // startG2,endG2,numberWanted);
446 int bestSet = -1;
447#if 0
448 // make sure first available is clean (in case last iteration rejected)
449 cost[firstAvailable_] = 0.0;
450 length[firstAvailable_] = 0;
451 model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
452 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
453 {
454 for (int i = firstAvailable_; i < lastDynamic_; i++)
455 assert(!cost[i]);
456 }
457#endif
458 int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
459 int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_;
460 for (int iSet = startG2; iSet < endAll; iSet++) {
461 if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) {
462 // give up
463 numberWanted = 0;
464 break;
465 } else if (iSet == endG2 && bestSequence >= 0) {
466 break;
467 }
468 int gubRow = toIndex_[iSet];
469 if (gubRow >= 0) {
470 djMod = duals[gubRow+numberStaticRows_]; // have I got sign right?
471 } else {
472 int iBasic = keyVariable_[iSet];
473 if (iBasic >= maximumGubColumns_) {
474 djMod = 0.0; // set not in
475 } else {
476 // get dj without
477 djMod = 0.0;
478 for (CoinBigIndex j = startColumn_[iBasic];
479 j < startColumn_[iBasic+1]; j++) {
480 int jRow = row_[j];
481 djMod -= duals[jRow] * element_[j];
482 }
483 djMod += cost_[iBasic];
484 // See if gub slack possible - dj is djMod
485 if (getStatus(iSet) == ClpSimplex::atLowerBound) {
486 double value = -djMod;
487 if (value > tolerance) {
488 numberWanted--;
489 if (value > bestDj) {
490 // check flagged variable and correct dj
491 if (!flagged(iSet)) {
492 bestDj = value;
493 bestSequence = slackOffset + iSet;
494 bestDjMod = djMod;
495 bestSet = iSet;
496 } else {
497 // just to make sure we don't exit before got something
498 numberWanted++;
499 abort();
500 }
501 }
502 }
503 } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
504 double value = djMod;
505 if (value > tolerance) {
506 numberWanted--;
507 if (value > bestDj) {
508 // check flagged variable and correct dj
509 if (!flagged(iSet)) {
510 bestDj = value;
511 bestSequence = slackOffset + iSet;
512 bestDjMod = djMod;
513 bestSet = iSet;
514 } else {
515 // just to make sure we don't exit before got something
516 numberWanted++;
517 abort();
518 }
519 }
520 }
521 }
522 }
523 }
524 int iSequence = startSet_[iSet];
525 while (iSequence >= 0) {
526 DynamicStatus status = getDynamicStatus(iSequence);
527 if (status == atLowerBound || status == atUpperBound) {
528 double value = cost_[iSequence] - djMod;
529 for (CoinBigIndex j = startColumn_[iSequence];
530 j < startColumn_[iSequence+1]; j++) {
531 int jRow = row_[j];
532 value -= duals[jRow] * element_[j];
533 }
534 // change sign if at lower bound
535 if (status == atLowerBound)
536 value = -value;
537 if (value > tolerance) {
538 numberWanted--;
539 if (value > bestDj) {
540 // check flagged variable and correct dj
541 if (!flagged(iSequence)) {
542 if (false/*status == atLowerBound
543 &&keyValue(iSet)<1.0e-7*/) {
544 // can't come in because
545 // of ones at ub
546 numberWanted++;
547 } else {
548
549 bestDj = value;
550 bestSequence = structuralOffset + iSequence;
551 bestDjMod = djMod;
552 bestSet = iSet;
553 }
554 } else {
555 // just to make sure we don't exit before got something
556 numberWanted++;
557 }
558 }
559 }
560 }
561 iSequence = next_[iSequence]; //onto next in set
562 }
563 if (numberWanted <= 0) {
564 numberWanted = 0;
565 break;
566 }
567 }
568 if (bestSequence != saveSequence) {
569 savedBestGubDual_ = bestDjMod;
570 savedBestDj_ = bestDj;
571 savedBestSequence_ = bestSequence;
572 savedBestSet_ = bestSet;
573 }
574 // See if may be finished
575 if (!startG2 && bestSequence < 0)
576 infeasibilityWeight_ = model_->infeasibilityCost();
577 else if (bestSequence >= 0)
578 infeasibilityWeight_ = -1.0;
579 }
580 currentWanted_ = numberWanted;
581}
582/* Returns effective RHS if it is being used. This is used for long problems
583 or big gub or anywhere where going through full columns is
584 expensive. This may re-compute */
585double *
586ClpDynamicMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh,
587 bool
588#ifdef CLP_DEBUG
589 check
590#endif
591 )
592{
593 // forceRefresh=true;printf("take out forceRefresh\n");
594 if (!model_->numberIterations())
595 forceRefresh = true;
596 //check=false;
597#ifdef CLP_DEBUG
598 double * saveE = NULL;
599 if (rhsOffset_ && check) {
600 int numberRows = model->numberRows();
601 saveE = new double[numberRows];
602 }
603#endif
604 if (rhsOffset_) {
605#ifdef CLP_DEBUG
606 if (check) {
607 // no need - but check anyway
608 int numberRows = model->numberRows();
609 double * rhs = new double[numberRows];
610 int iRow;
611 int iSet;
612 CoinZeroN(rhs, numberRows);
613 // do ones at bounds before gub
614 const double * smallSolution = model->solutionRegion();
615 const double * element = matrix_->getElements();
616 const int * row = matrix_->getIndices();
617 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
618 const int * length = matrix_->getVectorLengths();
619 int iColumn;
620 double objectiveOffset = 0.0;
621 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
622 if (model->getStatus(iColumn) != ClpSimplex::basic) {
623 double value = smallSolution[iColumn];
624 for (CoinBigIndex j = startColumn[iColumn];
625 j < startColumn[iColumn] + length[iColumn]; j++) {
626 int jRow = row[j];
627 rhs[jRow] -= value * element[j];
628 }
629 }
630 }
631 if (columnLower_ || columnUpper_) {
632 double * solution = new double [numberGubColumns_];
633 for (iSet = 0; iSet < numberSets_; iSet++) {
634 int j = startSet_[iSet];
635 while (j >= 0) {
636 double value = 0.0;
637 if (getDynamicStatus(j) != inSmall) {
638 if (getDynamicStatus(j) == atLowerBound) {
639 if (columnLower_)
640 value = columnLower_[j];
641 } else if (getDynamicStatus(j) == atUpperBound) {
642 value = columnUpper_[j];
643 } else if (getDynamicStatus(j) == soloKey) {
644 value = keyValue(iSet);
645 }
646 objectiveOffset += value * cost_[j];
647 }
648 solution[j] = value;
649 j = next_[j]; //onto next in set
650 }
651 }
652 // ones in gub and in small problem
653 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
654 if (model_->getStatus(iColumn) != ClpSimplex::basic) {
655 int jFull = id_[iColumn-firstDynamic_];
656 solution[jFull] = smallSolution[iColumn];
657 }
658 }
659 for (iSet = 0; iSet < numberSets_; iSet++) {
660 int kRow = toIndex_[iSet];
661 if (kRow >= 0)
662 kRow += numberStaticRows_;
663 int j = startSet_[iSet];
664 while (j >= 0) {
665 double value = solution[j];
666 if (value) {
667 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
668 int iRow = row_[k];
669 rhs[iRow] -= element_[k] * value;
670 }
671 if (kRow >= 0)
672 rhs[kRow] -= value;
673 }
674 j = next_[j]; //onto next in set
675 }
676 }
677 delete [] solution;
678 } else {
679 // bounds
680 ClpSimplex::Status iStatus;
681 for (int iSet = 0; iSet < numberSets_; iSet++) {
682 int kRow = toIndex_[iSet];
683 if (kRow < 0) {
684 int iColumn = keyVariable_[iSet];
685 if (iColumn < maximumGubColumns_) {
686 // key is not treated as basic
687 double b = 0.0;
688 // key is structural - where is slack
689 iStatus = getStatus(iSet);
690 assert (iStatus != ClpSimplex::basic);
691 if (iStatus == ClpSimplex::atLowerBound)
692 b = lowerSet_[iSet];
693 else
694 b = upperSet_[iSet];
695 if (b) {
696 objectiveOffset += b * cost_[iColumn];
697 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
698 int iRow = row_[j];
699 rhs[iRow] -= element_[j] * b;
700 }
701 }
702 }
703 }
704 }
705 }
706 if (fabs(model->objectiveOffset() - objectiveOffset_ - objectiveOffset) > 1.0e-1)
707 printf("old offset %g, true %g\n", model->objectiveOffset(),
708 objectiveOffset_ - objectiveOffset);
709 for (iRow = 0; iRow < numberRows; iRow++) {
710 if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
711 printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
712 }
713 CoinMemcpyN(rhs, numberRows, saveE);
714 delete [] rhs;
715 }
716#endif
717 if (forceRefresh || (refreshFrequency_ && model->numberIterations() >=
718 lastRefresh_ + refreshFrequency_)) {
719 int numberRows = model->numberRows();
720 int iSet;
721 CoinZeroN(rhsOffset_, numberRows);
722 // do ones at bounds before gub
723 const double * smallSolution = model->solutionRegion();
724 const double * element = matrix_->getElements();
725 const int * row = matrix_->getIndices();
726 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
727 const int * length = matrix_->getVectorLengths();
728 int iColumn;
729 double objectiveOffset = 0.0;
730 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
731 if (model->getStatus(iColumn) != ClpSimplex::basic) {
732 double value = smallSolution[iColumn];
733 for (CoinBigIndex j = startColumn[iColumn];
734 j < startColumn[iColumn] + length[iColumn]; j++) {
735 int jRow = row[j];
736 rhsOffset_[jRow] -= value * element[j];
737 }
738 }
739 }
740 if (columnLower_ || columnUpper_) {
741 double * solution = new double [numberGubColumns_];
742 for (iSet = 0; iSet < numberSets_; iSet++) {
743 int j = startSet_[iSet];
744 while (j >= 0) {
745 double value = 0.0;
746 if (getDynamicStatus(j) != inSmall) {
747 if (getDynamicStatus(j) == atLowerBound) {
748 if (columnLower_)
749 value = columnLower_[j];
750 } else if (getDynamicStatus(j) == atUpperBound) {
751 value = columnUpper_[j];
752 assert (value<1.0e30);
753 } else if (getDynamicStatus(j) == soloKey) {
754 value = keyValue(iSet);
755 }
756 objectiveOffset += value * cost_[j];
757 }
758 solution[j] = value;
759 j = next_[j]; //onto next in set
760 }
761 }
762 // ones in gub and in small problem
763 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
764 if (model_->getStatus(iColumn) != ClpSimplex::basic) {
765 int jFull = id_[iColumn-firstDynamic_];
766 solution[jFull] = smallSolution[iColumn];
767 }
768 }
769 for (iSet = 0; iSet < numberSets_; iSet++) {
770 int kRow = toIndex_[iSet];
771 if (kRow >= 0)
772 kRow += numberStaticRows_;
773 int j = startSet_[iSet];
774 while (j >= 0) {
775 //? use DynamicStatus status = getDynamicStatus(j);
776 double value = solution[j];
777 if (value) {
778 for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
779 int iRow = row_[k];
780 rhsOffset_[iRow] -= element_[k] * value;
781 }
782 if (kRow >= 0)
783 rhsOffset_[kRow] -= value;
784 }
785 j = next_[j]; //onto next in set
786 }
787 }
788 delete [] solution;
789 } else {
790 // bounds
791 ClpSimplex::Status iStatus;
792 for (int iSet = 0; iSet < numberSets_; iSet++) {
793 int kRow = toIndex_[iSet];
794 if (kRow < 0) {
795 int iColumn = keyVariable_[iSet];
796 if (iColumn < maximumGubColumns_) {
797 // key is not treated as basic
798 double b = 0.0;
799 // key is structural - where is slack
800 iStatus = getStatus(iSet);
801 assert (iStatus != ClpSimplex::basic);
802 if (iStatus == ClpSimplex::atLowerBound)
803 b = lowerSet_[iSet];
804 else
805 b = upperSet_[iSet];
806 if (b) {
807 objectiveOffset += b * cost_[iColumn];
808 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
809 int iRow = row_[j];
810 rhsOffset_[iRow] -= element_[j] * b;
811 }
812 }
813 }
814 }
815 }
816 }
817 model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
818#ifdef CLP_DEBUG
819 if (saveE) {
820 int iRow;
821 for (iRow = 0; iRow < numberRows; iRow++) {
822 if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3)
823 printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]);
824 }
825 delete [] saveE;
826 }
827#endif
828 lastRefresh_ = model->numberIterations();
829 }
830 }
831 return rhsOffset_;
832}
833/*
834 update information for a pivot (and effective rhs)
835*/
836int
837ClpDynamicMatrix::updatePivot(ClpSimplex * model, double oldInValue, double oldOutValue)
838{
839
840 // now update working model
841 int sequenceIn = model->sequenceIn();
842 int sequenceOut = model->sequenceOut();
843 int numberColumns = model->numberColumns();
844 if (sequenceIn != sequenceOut && sequenceIn < numberColumns)
845 backToPivotRow_[sequenceIn] = model->pivotRow();
846 if (sequenceIn >= firstDynamic_ && sequenceIn < numberColumns) {
847 int bigSequence = id_[sequenceIn-firstDynamic_];
848 if (getDynamicStatus(bigSequence) != inSmall) {
849 firstAvailable_++;
850 setDynamicStatus(bigSequence, inSmall);
851 }
852 }
853 // make sure slack is synchronized
854 if (sequenceIn >= numberColumns + numberStaticRows_) {
855 int iDynamic = sequenceIn - numberColumns - numberStaticRows_;
856 int iSet = fromIndex_[iDynamic];
857 setStatus(iSet, model->getStatus(sequenceIn));
858 }
859 if (sequenceOut >= numberColumns + numberStaticRows_) {
860 int iDynamic = sequenceOut - numberColumns - numberStaticRows_;
861 int iSet = fromIndex_[iDynamic];
862 // out may have gone through barrier - so check
863 double valueOut = model->lowerRegion()[sequenceOut];
864 if (fabs(valueOut - static_cast<double> (lowerSet_[iSet])) <
865 fabs(valueOut - static_cast<double> (upperSet_[iSet])))
866 setStatus(iSet, ClpSimplex::atLowerBound);
867 else
868 setStatus(iSet, ClpSimplex::atUpperBound);
869 if (lowerSet_[iSet] == upperSet_[iSet])
870 setStatus(iSet, ClpSimplex::isFixed);
871#if 0
872 if (getStatus(iSet) != model->getStatus(sequenceOut))
873 printf("** set %d status %d, var status %d\n", iSet,
874 getStatus(iSet), model->getStatus(sequenceOut));
875#endif
876 }
877 ClpMatrixBase::updatePivot(model, oldInValue, oldOutValue);
878#ifdef CLP_DEBUG
879 char * inSmall = new char [numberGubColumns_];
880 memset(inSmall, 0, numberGubColumns_);
881 const double * solution = model->solutionRegion();
882 for (int i = 0; i < numberGubColumns_; i++)
883 if (getDynamicStatus(i) == ClpDynamicMatrix::inSmall)
884 inSmall[i] = 1;
885 for (int i = firstDynamic_; i < firstAvailable_; i++) {
886 int k = id_[i-firstDynamic_];
887 inSmall[k] = 0;
888 //if (k>=23289&&k<23357&&solution[i])
889 //printf("var %d (in small %d) has value %g\n",k,i,solution[i]);
890 }
891 for (int i = 0; i < numberGubColumns_; i++)
892 assert (!inSmall[i]);
893 delete [] inSmall;
894#ifndef NDEBUG
895 for (int i = 0; i < numberActiveSets_; i++) {
896 int iSet = fromIndex_[i];
897 assert (toIndex_[iSet] == i);
898 //if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_))
899 //printf("*** set %d status %d, var status %d\n",iSet,
900 // getStatus(iSet),model->getRowStatus(i+numberStaticRows_));
901 //assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet));
902 //if (iSet==1035) {
903 //printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_],
904 // model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]);
905 //}
906 }
907#endif
908#endif
909 if (numberStaticRows_+numberActiveSets_<model->numberRows())
910 return 0;
911 else
912 return 1;
913}
914/*
915 utility dual function for dealing with dynamic constraints
916 mode=n see ClpGubMatrix.hpp for definition
917 Remember to update here when settled down
918*/
919void
920ClpDynamicMatrix::dualExpanded(ClpSimplex * model,
921 CoinIndexedVector * /*array*/,
922 double * /*other*/, int mode)
923{
924 switch (mode) {
925 // modify costs before transposeUpdate
926 case 0:
927 break;
928 // create duals for key variables (without check on dual infeasible)
929 case 1:
930 break;
931 // as 1 but check slacks and compute djs
932 case 2: {
933 // do pivots
934 int * pivotVariable = model->pivotVariable();
935 int numberRows = numberStaticRows_ + numberActiveSets_;
936 int numberColumns = model->numberColumns();
937 for (int iRow = 0; iRow < numberRows; iRow++) {
938 int iPivot = pivotVariable[iRow];
939 if (iPivot < numberColumns)
940 backToPivotRow_[iPivot] = iRow;
941 }
942 if (noCheck_ >= 0) {
943 if (infeasibilityWeight_ != model_->infeasibilityCost()) {
944 // don't bother checking
945 sumDualInfeasibilities_ = 100.0;
946 numberDualInfeasibilities_ = 1;
947 sumOfRelaxedDualInfeasibilities_ = 100.0;
948 return;
949 }
950 }
951 // In theory we should subtract out ones we have done but ....
952 // If key slack then dual 0.0
953 // If not then slack could be dual infeasible
954 // dj for key is zero so that defines dual on set
955 int i;
956 double * dual = model->dualRowSolution();
957 double dualTolerance = model->dualTolerance();
958 double relaxedTolerance = dualTolerance;
959 // we can't really trust infeasibilities if there is dual error
960 double error = CoinMin(1.0e-2, model->largestDualError());
961 // allow tolerance at least slightly bigger than standard
962 relaxedTolerance = relaxedTolerance + error;
963 // but we will be using difference
964 relaxedTolerance -= dualTolerance;
965 sumDualInfeasibilities_ = 0.0;
966 numberDualInfeasibilities_ = 0;
967 sumOfRelaxedDualInfeasibilities_ = 0.0;
968 for (i = 0; i < numberSets_; i++) {
969 double value = 0.0;
970 int gubRow = toIndex_[i];
971 if (gubRow < 0) {
972 int kColumn = keyVariable_[i];
973 if (kColumn < maximumGubColumns_) {
974 // dj without set
975 value = cost_[kColumn];
976 for (CoinBigIndex j = startColumn_[kColumn];
977 j < startColumn_[kColumn+1]; j++) {
978 int iRow = row_[j];
979 value -= dual[iRow] * element_[j];
980 }
981 double infeasibility = 0.0;
982 if (getStatus(i) == ClpSimplex::atLowerBound) {
983 if (-value > dualTolerance)
984 infeasibility = -value - dualTolerance;
985 } else if (getStatus(i) == ClpSimplex::atUpperBound) {
986 if (value > dualTolerance)
987 infeasibility = value - dualTolerance;
988 }
989 if (infeasibility > 0.0) {
990 sumDualInfeasibilities_ += infeasibility;
991 if (infeasibility > relaxedTolerance)
992 sumOfRelaxedDualInfeasibilities_ += infeasibility;
993 numberDualInfeasibilities_ ++;
994 }
995 }
996 } else {
997 value = dual[gubRow+numberStaticRows_];
998 }
999 // Now subtract out from all
1000 int k = startSet_[i];
1001 while (k >= 0) {
1002 if (getDynamicStatus(k) != inSmall) {
1003 double djValue = cost_[k] - value;
1004 for (CoinBigIndex j = startColumn_[k];
1005 j < startColumn_[k+1]; j++) {
1006 int iRow = row_[j];
1007 djValue -= dual[iRow] * element_[j];
1008 }
1009 double infeasibility = 0.0;
1010 if (getDynamicStatus(k) == atLowerBound) {
1011 if (djValue < -dualTolerance)
1012 infeasibility = -djValue - dualTolerance;
1013 } else if (getDynamicStatus(k) == atUpperBound) {
1014 // at upper bound
1015 if (djValue > dualTolerance)
1016 infeasibility = djValue - dualTolerance;
1017 }
1018 if (infeasibility > 0.0) {
1019 sumDualInfeasibilities_ += infeasibility;
1020 if (infeasibility > relaxedTolerance)
1021 sumOfRelaxedDualInfeasibilities_ += infeasibility;
1022 numberDualInfeasibilities_ ++;
1023 }
1024 }
1025 k = next_[k]; //onto next in set
1026 }
1027 }
1028 }
1029 infeasibilityWeight_ = -1.0;
1030 break;
1031 // Report on infeasibilities of key variables
1032 case 3: {
1033 model->setSumDualInfeasibilities(model->sumDualInfeasibilities() +
1034 sumDualInfeasibilities_);
1035 model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() +
1036 numberDualInfeasibilities_);
1037 model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() +
1038 sumOfRelaxedDualInfeasibilities_);
1039 }
1040 break;
1041 // modify costs before transposeUpdate for partial pricing
1042 case 4:
1043 break;
1044 }
1045}
1046/*
1047 general utility function for dealing with dynamic constraints
1048 mode=n see ClpGubMatrix.hpp for definition
1049 Remember to update here when settled down
1050*/
1051int
1052ClpDynamicMatrix::generalExpanded(ClpSimplex * model, int mode, int &number)
1053{
1054 int returnCode = 0;
1055#if 0 //ndef NDEBUG
1056 {
1057 int numberColumns = model->numberColumns();
1058 int numberRows = model->numberRows();
1059 int * pivotVariable = model->pivotVariable();
1060 if (pivotVariable&&model->numberIterations()) {
1061 for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
1062 assert (pivotVariable[i]==i+numberColumns);
1063 }
1064 }
1065 }
1066#endif
1067 switch (mode) {
1068 // Fill in pivotVariable
1069 case 0: {
1070 // If no effective rhs - form it
1071 if (!rhsOffset_) {
1072 rhsOffset_ = new double[model->numberRows()];
1073 rhsOffset(model, true);
1074 }
1075 int i;
1076 int numberBasic = number;
1077 int numberColumns = model->numberColumns();
1078 // Use different array so can build from true pivotVariable_
1079 //int * pivotVariable = model->pivotVariable();
1080 int * pivotVariable = model->rowArray(0)->getIndices();
1081 for (i = 0; i < numberColumns; i++) {
1082 if (model->getColumnStatus(i) == ClpSimplex::basic)
1083 pivotVariable[numberBasic++] = i;
1084 }
1085 number = numberBasic;
1086 }
1087 break;
1088 // Do initial extra rows + maximum basic
1089 case 2: {
1090 number = model->numberRows();
1091 }
1092 break;
1093 // Before normal replaceColumn
1094 case 3: {
1095 if (numberActiveSets_ + numberStaticRows_ == model_->numberRows()) {
1096 // no space - re-factorize
1097 returnCode = 4;
1098 number = -1; // say no need for normal replaceColumn
1099 }
1100 }
1101 break;
1102 // To see if can dual or primal
1103 case 4: {
1104 returnCode = 1;
1105 }
1106 break;
1107 // save status
1108 case 5: {
1109 memcpy(status_+numberSets_,status_,numberSets_);
1110 memcpy(status_+2*numberSets_,&numberActiveSets_,sizeof(int));
1111 memcpy(dynamicStatus_+maximumGubColumns_,
1112 dynamicStatus_,maximumGubColumns_);
1113 }
1114 break;
1115 // restore status
1116 case 6: {
1117 memcpy(status_,status_+numberSets_,numberSets_);
1118 memcpy(&numberActiveSets_,status_+2*numberSets_,sizeof(int));
1119 memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_,
1120 maximumGubColumns_);
1121 initialProblem();
1122 }
1123 break;
1124 // unflag all variables
1125 case 8: {
1126 for (int i = 0; i < numberGubColumns_; i++) {
1127 if (flagged(i)) {
1128 unsetFlagged(i);
1129 returnCode++;
1130 }
1131 }
1132 }
1133 break;
1134 // redo costs in primal
1135 case 9: {
1136 double * cost = model->costRegion();
1137 double * solution = model->solutionRegion();
1138 double * columnLower = model->lowerRegion();
1139 double * columnUpper = model->upperRegion();
1140 int i;
1141 bool doCosts = (number & 4) != 0;
1142 bool doBounds = (number & 1) != 0;
1143 for ( i = firstDynamic_; i < firstAvailable_; i++) {
1144 int jColumn = id_[i-firstDynamic_];
1145 if (doBounds) {
1146 if (!columnLower_ && !columnUpper_) {
1147 columnLower[i] = 0.0;
1148 columnUpper[i] = COIN_DBL_MAX;
1149 } else {
1150 if (columnLower_)
1151 columnLower[i] = columnLower_[jColumn];
1152 else
1153 columnLower[i] = 0.0;
1154 if (columnUpper_)
1155 columnUpper[i] = columnUpper_[jColumn];
1156 else
1157 columnUpper[i] = COIN_DBL_MAX;
1158 }
1159 }
1160 if (doCosts) {
1161 cost[i] = cost_[jColumn];
1162 // Original bounds
1163 if (model->nonLinearCost())
1164 model->nonLinearCost()->setOne(i, solution[i],
1165 this->columnLower(jColumn),
1166 this->columnUpper(jColumn), cost_[jColumn]);
1167 }
1168 }
1169 // and active sets
1170 for (i = 0; i < numberActiveSets_; i++ ) {
1171 int iSet = fromIndex_[i];
1172 int iSequence = lastDynamic_ + numberStaticRows_ + i;
1173 if (doBounds) {
1174 if (lowerSet_[iSet] > -1.0e20)
1175 columnLower[iSequence] = lowerSet_[iSet];
1176 else
1177 columnLower[iSequence] = -COIN_DBL_MAX;
1178 if (upperSet_[iSet] < 1.0e20)
1179 columnUpper[iSequence] = upperSet_[iSet];
1180 else
1181 columnUpper[iSequence] = COIN_DBL_MAX;
1182 }
1183 if (doCosts) {
1184 if (model->nonLinearCost()) {
1185 double trueLower;
1186 if (lowerSet_[iSet] > -1.0e20)
1187 trueLower = lowerSet_[iSet];
1188 else
1189 trueLower = -COIN_DBL_MAX;
1190 double trueUpper;
1191 if (upperSet_[iSet] < 1.0e20)
1192 trueUpper = upperSet_[iSet];
1193 else
1194 trueUpper = COIN_DBL_MAX;
1195 model->nonLinearCost()->setOne(iSequence, solution[iSequence],
1196 trueLower, trueUpper, 0.0);
1197 }
1198 }
1199 }
1200 }
1201 break;
1202 // return 1 if there may be changing bounds on variable (column generation)
1203 case 10: {
1204 // return 1 as bounds on rhs will change
1205 returnCode = 1;
1206 }
1207 break;
1208 // make sure set is clean
1209 case 7: {
1210 // first flag
1211 if (number >= firstDynamic_ && number < lastDynamic_) {
1212 int sequence = id_[number-firstDynamic_];
1213 setFlagged(sequence);
1214 //model->clearFlagged(number);
1215 } else if (number>=model_->numberColumns()+numberStaticRows_) {
1216 // slack
1217 int iSet = fromIndex_[number-model_->numberColumns()-
1218 numberStaticRows_];
1219 setFlaggedSlack(iSet);
1220 //model->clearFlagged(number);
1221 }
1222 }
1223 case 11: {
1224 //int sequenceIn = model->sequenceIn();
1225 if (number >= firstDynamic_ && number < lastDynamic_) {
1226 //assert (number == model->sequenceIn());
1227 // take out variable (but leave key)
1228 double * cost = model->costRegion();
1229 double * columnLower = model->lowerRegion();
1230 double * columnUpper = model->upperRegion();
1231 double * solution = model->solutionRegion();
1232 int * length = matrix_->getMutableVectorLengths();
1233#ifndef NDEBUG
1234 if (length[number]) {
1235 int * row = matrix_->getMutableIndices();
1236 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1237 int iRow = row[startColumn[number] + length[number] - 1];
1238 int which = iRow - numberStaticRows_;
1239 assert (which >= 0);
1240 int iSet = fromIndex_[which];
1241 assert (toIndex_[iSet] == which);
1242 }
1243#endif
1244 // no need firstAvailable_--;
1245 solution[firstAvailable_] = 0.0;
1246 cost[firstAvailable_] = 0.0;
1247 length[firstAvailable_] = 0;
1248 model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
1249 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1250 columnLower[firstAvailable_] = 0.0;
1251 columnUpper[firstAvailable_] = COIN_DBL_MAX;
1252
1253 // not really in small problem
1254 int iBig = id_[number-firstDynamic_];
1255 if (model->getStatus(number) == ClpSimplex::atLowerBound) {
1256 setDynamicStatus(iBig, atLowerBound);
1257 if (columnLower_)
1258 modifyOffset(number, columnLower_[iBig]);
1259 } else {
1260 setDynamicStatus(iBig, atUpperBound);
1261 modifyOffset(number, columnUpper_[iBig]);
1262 }
1263 } else if (number>=model_->numberColumns()+numberStaticRows_) {
1264 // slack
1265 int iSet = fromIndex_[number-model_->numberColumns()-
1266 numberStaticRows_];
1267 printf("what now - set %d\n",iSet);
1268 }
1269 }
1270 break;
1271 default:
1272 break;
1273 }
1274 return returnCode;
1275}
1276/* Purely for column generation and similar ideas. Allows
1277 matrix and any bounds or costs to be updated (sensibly).
1278 Returns non-zero if any changes.
1279*/
1280int
1281ClpDynamicMatrix::refresh(ClpSimplex * model)
1282{
1283 // If at beginning then create initial problem
1284 if (firstDynamic_ == firstAvailable_) {
1285 initialProblem();
1286 return 1;
1287 } else if (!model->nonLinearCost()) {
1288 // will be same as last time
1289 return 1;
1290 }
1291#ifndef NDEBUG
1292 {
1293 int numberColumns = model->numberColumns();
1294 int numberRows = model->numberRows();
1295 int * pivotVariable = model->pivotVariable();
1296 for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
1297 assert (pivotVariable[i]==i+numberColumns);
1298 }
1299 }
1300#endif
1301 // lookup array
1302 int * active = new int [numberActiveSets_];
1303 CoinZeroN(active, numberActiveSets_);
1304 int iColumn;
1305 int numberColumns = model->numberColumns();
1306 double * element = matrix_->getMutableElements();
1307 int * row = matrix_->getMutableIndices();
1308 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1309 int * length = matrix_->getMutableVectorLengths();
1310 double * cost = model->costRegion();
1311 double * columnLower = model->lowerRegion();
1312 double * columnUpper = model->upperRegion();
1313 CoinBigIndex numberElements = startColumn[firstDynamic_];
1314 // first just do lookup and basic stuff
1315 int currentNumber = firstAvailable_;
1316 firstAvailable_ = firstDynamic_;
1317 double objectiveChange = 0.0;
1318 double * solution = model->solutionRegion();
1319 int currentNumberActiveSets = numberActiveSets_;
1320 for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
1321 int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
1322 int which = iRow - numberStaticRows_;
1323 assert (which >= 0);
1324 int iSet = fromIndex_[which];
1325 assert (toIndex_[iSet] == which);
1326 if (model->getStatus(iColumn) == ClpSimplex::basic) {
1327 active[which]++;
1328 // may as well make key
1329 keyVariable_[iSet] = id_[iColumn-firstDynamic_];
1330 }
1331 }
1332 int i;
1333 numberActiveSets_ = 0;
1334 int numberDeleted = 0;
1335 for (i = 0; i < currentNumberActiveSets; i++) {
1336 int iRow = i + numberStaticRows_;
1337 int numberActive = active[i];
1338 int iSet = fromIndex_[i];
1339 if (model->getRowStatus(iRow) == ClpSimplex::basic) {
1340 numberActive++;
1341 // may as well make key
1342 keyVariable_[iSet] = maximumGubColumns_ + iSet;
1343 }
1344 if (numberActive > 1) {
1345 // keep
1346 active[i] = numberActiveSets_;
1347 numberActiveSets_++;
1348 } else {
1349 active[i] = -1;
1350 }
1351 }
1352
1353 for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
1354 int iRow = row[startColumn[iColumn] + length[iColumn] - 1];
1355 int which = iRow - numberStaticRows_;
1356 int jColumn = id_[iColumn-firstDynamic_];
1357 if (model->getStatus(iColumn) == ClpSimplex::atLowerBound ||
1358 model->getStatus(iColumn) == ClpSimplex::atUpperBound) {
1359 double value = solution[iColumn];
1360 bool toLowerBound = true;
1361 assert (jColumn>=0);assert (iColumn>=0);
1362 if (columnUpper_) {
1363 if (!columnLower_) {
1364 if (fabs(value - columnUpper_[jColumn]) < fabs(value))
1365 toLowerBound = false;
1366 } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) {
1367 toLowerBound = false;
1368 }
1369 }
1370 if (toLowerBound) {
1371 // throw out to lower bound
1372 if (columnLower_) {
1373 setDynamicStatus(jColumn, atLowerBound);
1374 // treat solution as if exactly at a bound
1375 double value = columnLower[iColumn];
1376 objectiveChange += cost[iColumn] * value;
1377 } else {
1378 setDynamicStatus(jColumn, atLowerBound);
1379 }
1380 } else {
1381 // throw out to upper bound
1382 setDynamicStatus(jColumn, atUpperBound);
1383 double value = columnUpper[iColumn];
1384 objectiveChange += cost[iColumn] * value;
1385 }
1386 numberDeleted++;
1387 } else {
1388 assert(model->getStatus(iColumn) != ClpSimplex::superBasic); // deal with later
1389 int iPut = active[which];
1390 if (iPut >= 0) {
1391 // move
1392 id_[firstAvailable_-firstDynamic_] = jColumn;
1393 int numberThis = startColumn_[jColumn+1] - startColumn_[jColumn] + 1;
1394 length[firstAvailable_] = numberThis;
1395 cost[firstAvailable_] = cost[iColumn];
1396 columnLower[firstAvailable_] = columnLower[iColumn];
1397 columnUpper[firstAvailable_] = columnUpper[iColumn];
1398 model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn], 0.0, COIN_DBL_MAX,
1399 cost_[jColumn]);
1400 CoinBigIndex base = startColumn_[jColumn];
1401 for (int j = 0; j < numberThis - 1; j++) {
1402 row[numberElements] = row_[base+j];
1403 element[numberElements++] = element_[base+j];
1404 }
1405 row[numberElements] = iPut + numberStaticRows_;
1406 element[numberElements++] = 1.0;
1407 model->setStatus(firstAvailable_, model->getStatus(iColumn));
1408 solution[firstAvailable_] = solution[iColumn];
1409 firstAvailable_++;
1410 startColumn[firstAvailable_] = numberElements;
1411 }
1412 }
1413 }
1414 // zero solution
1415 CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
1416 // zero cost
1417 CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
1418 // zero lengths
1419 CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
1420 for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
1421 model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
1422 model->setStatus(iColumn, ClpSimplex::atLowerBound);
1423 columnLower[iColumn] = 0.0;
1424 columnUpper[iColumn] = COIN_DBL_MAX;
1425 }
1426 // move rhs and set rest to infinite
1427 numberActiveSets_ = 0;
1428 for (i = 0; i < currentNumberActiveSets; i++) {
1429 int iSet = fromIndex_[i];
1430 assert (toIndex_[iSet] == i);
1431 int iRow = i + numberStaticRows_;
1432 int iPut = active[i];
1433 if (iPut >= 0) {
1434 int in = iRow + numberColumns;
1435 int out = iPut + numberColumns + numberStaticRows_;
1436 solution[out] = solution[in];
1437 columnLower[out] = columnLower[in];
1438 columnUpper[out] = columnUpper[in];
1439 cost[out] = cost[in];
1440 double trueLower;
1441 if (lowerSet_[iSet] > -1.0e20)
1442 trueLower = lowerSet_[iSet];
1443 else
1444 trueLower = -COIN_DBL_MAX;
1445 double trueUpper;
1446 if (upperSet_[iSet] < 1.0e20)
1447 trueUpper = upperSet_[iSet];
1448 else
1449 trueUpper = COIN_DBL_MAX;
1450 model->nonLinearCost()->setOne(out, solution[out],
1451 trueLower, trueUpper, 0.0);
1452 model->setStatus(out, model->getStatus(in));
1453 toIndex_[iSet] = numberActiveSets_;
1454 rhsOffset_[numberActiveSets_+numberStaticRows_] = rhsOffset_[i+numberStaticRows_];
1455 fromIndex_[numberActiveSets_++] = fromIndex_[i];
1456 } else {
1457 // adjust offset
1458 // put out as key
1459 int jColumn = keyVariable_[iSet];
1460 toIndex_[iSet] = -1;
1461 if (jColumn < maximumGubColumns_) {
1462 setDynamicStatus(jColumn, soloKey);
1463 double value = keyValue(iSet);
1464 objectiveChange += cost_[jColumn] * value;
1465 modifyOffset(jColumn, -value);
1466 }
1467 }
1468 }
1469 model->setObjectiveOffset(model->objectiveOffset() - objectiveChange);
1470 delete [] active;
1471 for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
1472 int iSequence = i + numberStaticRows_ + numberColumns;
1473 solution[iSequence] = 0.0;
1474 columnLower[iSequence] = -COIN_DBL_MAX;
1475 columnUpper[iSequence] = COIN_DBL_MAX;
1476 cost[iSequence] = 0.0;
1477 model->nonLinearCost()->setOne(iSequence, solution[iSequence],
1478 columnLower[iSequence],
1479 columnUpper[iSequence], 0.0);
1480 model->setStatus(iSequence, ClpSimplex::basic);
1481 rhsOffset_[i+numberStaticRows_] = 0.0;
1482 }
1483 if (currentNumberActiveSets != numberActiveSets_ || numberDeleted) {
1484 // deal with pivotVariable
1485 int * pivotVariable = model->pivotVariable();
1486 int sequence = firstDynamic_;
1487 int iRow = 0;
1488 int base1 = firstDynamic_;
1489 int base2 = lastDynamic_;
1490 int base3 = numberColumns + numberStaticRows_;
1491 int numberRows = numberStaticRows_ + currentNumberActiveSets;
1492 while (sequence < firstAvailable_) {
1493 int iPivot = pivotVariable[iRow];
1494 while (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
1495 iPivot = pivotVariable[++iRow];
1496 }
1497 pivotVariable[iRow++] = sequence;
1498 sequence++;
1499 }
1500 // move normal basic ones down
1501 int iPut = iRow;
1502 for (; iRow < numberRows; iRow++) {
1503 int iPivot = pivotVariable[iRow];
1504 if (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) {
1505 pivotVariable[iPut++] = iPivot;
1506 }
1507 }
1508 // and basic others
1509 for (i = 0; i < numberActiveSets_; i++) {
1510 if (model->getRowStatus(i + numberStaticRows_) == ClpSimplex::basic) {
1511 pivotVariable[iPut++] = i + base3;
1512 }
1513 }
1514 if (iPut<numberStaticRows_+numberActiveSets_) {
1515 printf("lost %d sets\n",
1516 numberStaticRows_+numberActiveSets_-iPut);
1517 iPut = numberStaticRows_+numberActiveSets_;
1518 }
1519 for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
1520 pivotVariable[iPut++] = i + base3;
1521 }
1522 //assert (iPut == numberRows);
1523 }
1524#ifdef CLP_DEBUG
1525#if 0
1526 printf("row for set 244 is %d, row status %d value %g ", toIndex_[244], status_[244],
1527 keyValue(244));
1528 int jj = startSet_[244];
1529 while (jj >= 0) {
1530 printf("var %d status %d ", jj, dynamicStatus_[jj]);
1531 jj = next_[jj];
1532 }
1533 printf("\n");
1534#endif
1535#ifdef CLP_DEBUG
1536 {
1537 // debug coding to analyze set
1538 int i;
1539 int kSet = -6;
1540 if (kSet >= 0) {
1541 int * back = new int [numberGubColumns_];
1542 for (i = 0; i < numberGubColumns_; i++)
1543 back[i] = -1;
1544 for (i = firstDynamic_; i < firstAvailable_; i++)
1545 back[id_[i-firstDynamic_]] = i;
1546 int sequence = startSet_[kSet];
1547 if (toIndex_[kSet] < 0) {
1548 printf("Not in - Set %d - slack status %d, key %d\n", kSet, status_[kSet], keyVariable_[kSet]);
1549 while (sequence >= 0) {
1550 printf("( %d status %d ) ", sequence, getDynamicStatus(sequence));
1551 sequence = next_[sequence];
1552 }
1553 } else {
1554 int iRow = numberStaticRows_ + toIndex_[kSet];
1555 printf("In - Set %d - slack status %d, key %d offset %g slack %g\n", kSet, status_[kSet], keyVariable_[kSet],
1556 rhsOffset_[iRow], model->solutionRegion(0)[iRow]);
1557 while (sequence >= 0) {
1558 int iBack = back[sequence];
1559 printf("( %d status %d value %g) ", sequence, getDynamicStatus(sequence), model->solutionRegion()[iBack]);
1560 sequence = next_[sequence];
1561 }
1562 }
1563 printf("\n");
1564 delete [] back;
1565 }
1566 }
1567#endif
1568 int n = numberActiveSets_;
1569 for (i = 0; i < numberSets_; i++) {
1570 if (toIndex_[i] < 0) {
1571 //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]);
1572 n++;
1573 }
1574 int k=0;
1575 for (int j=startSet_[i];j<startSet_[i+1];j++) {
1576 if (getDynamicStatus(j)==soloKey)
1577 k++;
1578 }
1579 assert (k<2);
1580 }
1581 assert (n == numberSets_);
1582#endif
1583 return 1;
1584}
1585void
1586ClpDynamicMatrix::times(double scalar,
1587 const double * x, double * y) const
1588{
1589 if (model_->specialOptions() != 16) {
1590 ClpPackedMatrix::times(scalar, x, y);
1591 } else {
1592 int iRow;
1593 const double * element = matrix_->getElements();
1594 const int * row = matrix_->getIndices();
1595 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1596 const int * length = matrix_->getVectorLengths();
1597 int * pivotVariable = model_->pivotVariable();
1598 for (iRow = 0; iRow < numberStaticRows_ + numberActiveSets_; iRow++) {
1599 y[iRow] -= scalar * rhsOffset_[iRow];
1600 int iColumn = pivotVariable[iRow];
1601 if (iColumn < lastDynamic_) {
1602 CoinBigIndex j;
1603 double value = scalar * x[iColumn];
1604 if (value) {
1605 for (j = startColumn[iColumn];
1606 j < startColumn[iColumn] + length[iColumn]; j++) {
1607 int jRow = row[j];
1608 y[jRow] += value * element[j];
1609 }
1610 }
1611 }
1612 }
1613 }
1614}
1615// Modifies rhs offset
1616void
1617ClpDynamicMatrix::modifyOffset(int sequence, double amount)
1618{
1619 if (amount) {
1620 assert (rhsOffset_);
1621 CoinBigIndex j;
1622 for (j = startColumn_[sequence]; j < startColumn_[sequence+1]; j++) {
1623 int iRow = row_[j];
1624 rhsOffset_[iRow] += amount * element_[j];
1625 }
1626 }
1627}
1628// Gets key value when none in small
1629double
1630ClpDynamicMatrix::keyValue(int iSet) const
1631{
1632 double value = 0.0;
1633 if (toIndex_[iSet] < 0) {
1634 int key = keyVariable_[iSet];
1635 if (key < maximumGubColumns_) {
1636 if (getStatus(iSet) == ClpSimplex::atLowerBound)
1637 value = lowerSet_[iSet];
1638 else
1639 value = upperSet_[iSet];
1640 int numberKey = 0;
1641 int j = startSet_[iSet];
1642 while (j >= 0) {
1643 DynamicStatus status = getDynamicStatus(j);
1644 assert (status != inSmall);
1645 if (status == soloKey) {
1646 numberKey++;
1647 } else if (status == atUpperBound) {
1648 value -= columnUpper_[j];
1649 } else if (columnLower_) {
1650 value -= columnLower_[j];
1651 }
1652 j = next_[j]; //onto next in set
1653 }
1654 assert (numberKey == 1);
1655 } else {
1656 int j = startSet_[iSet];
1657 while (j >= 0) {
1658 DynamicStatus status = getDynamicStatus(j);
1659 assert (status != inSmall);
1660 assert (status != soloKey);
1661 if (status == atUpperBound) {
1662 value += columnUpper_[j];
1663 } else if (columnLower_) {
1664 value += columnLower_[j];
1665 }
1666 j = next_[j]; //onto next in set
1667 }
1668#if 0
1669 // slack must be feasible
1670 double oldValue=value;
1671 value = CoinMax(value,lowerSet_[iSet]);
1672 value = CoinMin(value,upperSet_[iSet]);
1673 if (value!=oldValue)
1674 printf("using %g (not %g) for slack on set %d (%g,%g)\n",
1675 value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]);
1676#endif
1677 }
1678 }
1679 return value;
1680}
1681// Switches off dj checking each factorization (for BIG models)
1682void
1683ClpDynamicMatrix::switchOffCheck()
1684{
1685 noCheck_ = 0;
1686 infeasibilityWeight_ = 0.0;
1687}
1688/* Creates a variable. This is called after partial pricing and may modify matrix.
1689 May update bestSequence.
1690*/
1691void
1692ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence)
1693{
1694 int numberRows = model->numberRows();
1695 int slackOffset = lastDynamic_ + numberRows;
1696 int structuralOffset = slackOffset + numberSets_;
1697 int bestSequence2 = savedBestSequence_ - structuralOffset;
1698 if (bestSequence >= slackOffset) {
1699 double * columnLower = model->lowerRegion();
1700 double * columnUpper = model->upperRegion();
1701 double * solution = model->solutionRegion();
1702 double * reducedCost = model->djRegion();
1703 const double * duals = model->dualRowSolution();
1704 if (toIndex_[savedBestSet_] < 0) {
1705 // need to put key into basis
1706 int newRow = numberActiveSets_ + numberStaticRows_;
1707 model->dualRowSolution()[newRow] = savedBestGubDual_;
1708 double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set
1709 toIndex_[savedBestSet_] = numberActiveSets_;
1710 fromIndex_[numberActiveSets_++] = savedBestSet_;
1711 int iSequence = lastDynamic_ + newRow;
1712 // we need to get lower and upper correct
1713 double shift = 0.0;
1714 int j = startSet_[savedBestSet_];
1715 while (j >= 0) {
1716 if (getDynamicStatus(j) == atUpperBound)
1717 shift += columnUpper_[j];
1718 else if (getDynamicStatus(j) == atLowerBound && columnLower_)
1719 shift += columnLower_[j];
1720 j = next_[j]; //onto next in set
1721 }
1722 if (lowerSet_[savedBestSet_] > -1.0e20)
1723 columnLower[iSequence] = lowerSet_[savedBestSet_];
1724 else
1725 columnLower[iSequence] = -COIN_DBL_MAX;
1726 if (upperSet_[savedBestSet_] < 1.0e20)
1727 columnUpper[iSequence] = upperSet_[savedBestSet_];
1728 else
1729 columnUpper[iSequence] = COIN_DBL_MAX;
1730#ifdef CLP_DEBUG
1731 if (model_->logLevel() == 63) {
1732 printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n",
1733 bestSequence2, savedBestSet_, keyVariable_[savedBestSet_],
1734 columnLower[iSequence], columnUpper[iSequence], valueOfKey);
1735 int j = startSet_[savedBestSet_];
1736 while (j >= 0) {
1737 if (getDynamicStatus(j) == atUpperBound)
1738 printf("%d atup ", j);
1739 else if (getDynamicStatus(j) == atLowerBound)
1740 printf("%d atlo ", j);
1741 else if (getDynamicStatus(j) == soloKey)
1742 printf("%d solo ", j);
1743 else
1744 abort();
1745 j = next_[j]; //onto next in set
1746 }
1747 printf("\n");
1748 }
1749#endif
1750 if (keyVariable_[savedBestSet_] < maximumGubColumns_) {
1751 // slack not key
1752 model_->pivotVariable()[newRow] = firstAvailable_;
1753 backToPivotRow_[firstAvailable_] = newRow;
1754 model->setStatus(iSequence, getStatus(savedBestSet_));
1755 model->djRegion()[iSequence] = savedBestGubDual_;
1756 solution[iSequence] = valueOfKey;
1757 // create variable and pivot in
1758 int key = keyVariable_[savedBestSet_];
1759 setDynamicStatus(key, inSmall);
1760 double * element = matrix_->getMutableElements();
1761 int * row = matrix_->getMutableIndices();
1762 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1763 int * length = matrix_->getMutableVectorLengths();
1764 CoinBigIndex numberElements = startColumn[firstAvailable_];
1765 int numberThis = startColumn_[key+1] - startColumn_[key] + 1;
1766 if (numberElements + numberThis > numberElements_) {
1767 // need to redo
1768 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1769 matrix_->reserve(lastDynamic_, numberElements_);
1770 element = matrix_->getMutableElements();
1771 row = matrix_->getMutableIndices();
1772 // these probably okay but be safe
1773 startColumn = matrix_->getMutableVectorStarts();
1774 length = matrix_->getMutableVectorLengths();
1775 }
1776 // already set startColumn[firstAvailable_]=numberElements;
1777 length[firstAvailable_] = numberThis;
1778 model->costRegion()[firstAvailable_] = cost_[key];
1779 CoinBigIndex base = startColumn_[key];
1780 for (int j = 0; j < numberThis - 1; j++) {
1781 row[numberElements] = row_[base+j];
1782 element[numberElements++] = element_[base+j];
1783 }
1784 row[numberElements] = newRow;
1785 element[numberElements++] = 1.0;
1786 id_[firstAvailable_-firstDynamic_] = key;
1787 model->setObjectiveOffset(model->objectiveOffset() + cost_[key]*
1788 valueOfKey);
1789 model->solutionRegion()[firstAvailable_] = valueOfKey;
1790 model->setStatus(firstAvailable_, ClpSimplex::basic);
1791 // ***** need to adjust effective rhs
1792 if (!columnLower_ && !columnUpper_) {
1793 columnLower[firstAvailable_] = 0.0;
1794 columnUpper[firstAvailable_] = COIN_DBL_MAX;
1795 } else {
1796 if (columnLower_)
1797 columnLower[firstAvailable_] = columnLower_[key];
1798 else
1799 columnLower[firstAvailable_] = 0.0;
1800 if (columnUpper_)
1801 columnUpper[firstAvailable_] = columnUpper_[key];
1802 else
1803 columnUpper[firstAvailable_] = COIN_DBL_MAX;
1804 }
1805 model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
1806 columnLower[firstAvailable_],
1807 columnUpper[firstAvailable_], cost_[key]);
1808 startColumn[firstAvailable_+1] = numberElements;
1809 reducedCost[firstAvailable_] = 0.0;
1810 modifyOffset(key, valueOfKey);
1811 rhsOffset_[newRow] = -shift; // sign?
1812#ifdef CLP_DEBUG
1813 model->rowArray(1)->checkClear();
1814#endif
1815 // now pivot in
1816 unpack(model, model->rowArray(1), firstAvailable_);
1817 model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(1));
1818 double alpha = model->rowArray(1)->denseVector()[newRow];
1819 int updateStatus =
1820 model->factorization()->replaceColumn(model,
1821 model->rowArray(2),
1822 model->rowArray(1),
1823 newRow, alpha);
1824 model->rowArray(1)->clear();
1825 if (updateStatus) {
1826 if (updateStatus == 3) {
1827 // out of memory
1828 // increase space if not many iterations
1829 if (model->factorization()->pivots() <
1830 0.5 * model->factorization()->maximumPivots() &&
1831 model->factorization()->pivots() < 400)
1832 model->factorization()->areaFactor(
1833 model->factorization()->areaFactor() * 1.1);
1834 } else {
1835 printf("Bad returncode %d from replaceColumn\n", updateStatus);
1836 }
1837 bestSequence = -1;
1838 return;
1839 }
1840 // firstAvailable_ only finally updated if good pivot (in updatePivot)
1841 // otherwise it reverts to firstAvailableBefore_
1842 firstAvailable_++;
1843 } else {
1844 // slack key
1845 model->setStatus(iSequence, ClpSimplex::basic);
1846 model->djRegion()[iSequence] = 0.0;
1847 solution[iSequence] = valueOfKey+shift;
1848 rhsOffset_[newRow] = -shift; // sign?
1849 }
1850 // correct slack
1851 model->costRegion()[iSequence] = 0.0;
1852 model->nonLinearCost()->setOne(iSequence, solution[iSequence], columnLower[iSequence],
1853 columnUpper[iSequence], 0.0);
1854 }
1855 if (savedBestSequence_ >= structuralOffset) {
1856 // recompute dj and create
1857 double value = cost_[bestSequence2] - savedBestGubDual_;
1858 for (CoinBigIndex jBigIndex = startColumn_[bestSequence2];
1859 jBigIndex < startColumn_[bestSequence2+1]; jBigIndex++) {
1860 int jRow = row_[jBigIndex];
1861 value -= duals[jRow] * element_[jBigIndex];
1862 }
1863 int gubRow = toIndex_[savedBestSet_] + numberStaticRows_;
1864 double * element = matrix_->getMutableElements();
1865 int * row = matrix_->getMutableIndices();
1866 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1867 int * length = matrix_->getMutableVectorLengths();
1868 CoinBigIndex numberElements = startColumn[firstAvailable_];
1869 int numberThis = startColumn_[bestSequence2+1] - startColumn_[bestSequence2] + 1;
1870 if (numberElements + numberThis > numberElements_) {
1871 // need to redo
1872 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1873 matrix_->reserve(lastDynamic_, numberElements_);
1874 element = matrix_->getMutableElements();
1875 row = matrix_->getMutableIndices();
1876 // these probably okay but be safe
1877 startColumn = matrix_->getMutableVectorStarts();
1878 length = matrix_->getMutableVectorLengths();
1879 }
1880 // already set startColumn[firstAvailable_]=numberElements;
1881 length[firstAvailable_] = numberThis;
1882 model->costRegion()[firstAvailable_] = cost_[bestSequence2];
1883 CoinBigIndex base = startColumn_[bestSequence2];
1884 for (int j = 0; j < numberThis - 1; j++) {
1885 row[numberElements] = row_[base+j];
1886 element[numberElements++] = element_[base+j];
1887 }
1888 row[numberElements] = gubRow;
1889 element[numberElements++] = 1.0;
1890 id_[firstAvailable_-firstDynamic_] = bestSequence2;
1891 //printf("best %d\n",bestSequence2);
1892 model->solutionRegion()[firstAvailable_] = 0.0;
1893 model->clearFlagged(firstAvailable_);
1894 if (!columnLower_ && !columnUpper_) {
1895 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1896 columnLower[firstAvailable_] = 0.0;
1897 columnUpper[firstAvailable_] = COIN_DBL_MAX;
1898 } else {
1899 DynamicStatus status = getDynamicStatus(bestSequence2);
1900 if (columnLower_)
1901 columnLower[firstAvailable_] = columnLower_[bestSequence2];
1902 else
1903 columnLower[firstAvailable_] = 0.0;
1904 if (columnUpper_)
1905 columnUpper[firstAvailable_] = columnUpper_[bestSequence2];
1906 else
1907 columnUpper[firstAvailable_] = COIN_DBL_MAX;
1908 if (status == atLowerBound) {
1909 solution[firstAvailable_] = columnLower[firstAvailable_];
1910 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
1911 } else {
1912 solution[firstAvailable_] = columnUpper[firstAvailable_];
1913 model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
1914 }
1915 }
1916 model->setObjectiveOffset(model->objectiveOffset() + cost_[bestSequence2]*
1917 solution[firstAvailable_]);
1918 model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
1919 columnLower[firstAvailable_],
1920 columnUpper[firstAvailable_], cost_[bestSequence2]);
1921 bestSequence = firstAvailable_;
1922 // firstAvailable_ only updated if good pivot (in updatePivot)
1923 startColumn[firstAvailable_+1] = numberElements;
1924 //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_);
1925 reducedCost[bestSequence] = value;
1926 } else {
1927 // slack - row must just have been created
1928 assert (toIndex_[savedBestSet_] == numberActiveSets_ - 1);
1929 int newRow = numberStaticRows_ + numberActiveSets_ - 1;
1930 bestSequence = lastDynamic_ + newRow;
1931 reducedCost[bestSequence] = savedBestGubDual_;
1932 }
1933 }
1934 // clear for next iteration
1935 savedBestSequence_ = -1;
1936}
1937// Returns reduced cost of a variable
1938double
1939ClpDynamicMatrix::reducedCost(ClpSimplex * model, int sequence) const
1940{
1941 int numberRows = model->numberRows();
1942 int slackOffset = lastDynamic_ + numberRows;
1943 if (sequence < slackOffset)
1944 return model->djRegion()[sequence];
1945 else
1946 return savedBestDj_;
1947}
1948// Does gub crash
1949void
1950ClpDynamicMatrix::gubCrash()
1951{
1952 // Do basis - cheapest or slack if feasible
1953 int longestSet = 0;
1954 int iSet;
1955 for (iSet = 0; iSet < numberSets_; iSet++) {
1956 int n = 0;
1957 int j = startSet_[iSet];
1958 while (j >= 0) {
1959 n++;
1960 j = next_[j];
1961 }
1962 longestSet = CoinMax(longestSet, n);
1963 }
1964 double * upper = new double[longestSet+1];
1965 double * cost = new double[longestSet+1];
1966 double * lower = new double[longestSet+1];
1967 double * solution = new double[longestSet+1];
1968 int * back = new int[longestSet+1];
1969 double tolerance = model_->primalTolerance();
1970 double objectiveOffset = 0.0;
1971 for (iSet = 0; iSet < numberSets_; iSet++) {
1972 int iBasic = -1;
1973 double value = 0.0;
1974 // find cheapest
1975 int numberInSet = 0;
1976 int j = startSet_[iSet];
1977 while (j >= 0) {
1978 if (!columnLower_)
1979 lower[numberInSet] = 0.0;
1980 else
1981 lower[numberInSet] = columnLower_[j];
1982 if (!columnUpper_)
1983 upper[numberInSet] = COIN_DBL_MAX;
1984 else
1985 upper[numberInSet] = columnUpper_[j];
1986 back[numberInSet++] = j;
1987 j = next_[j];
1988 }
1989 CoinFillN(solution, numberInSet, 0.0);
1990 // and slack
1991 iBasic = numberInSet;
1992 solution[iBasic] = -value;
1993 lower[iBasic] = -upperSet_[iSet];
1994 upper[iBasic] = -lowerSet_[iSet];
1995 int kphase;
1996 if (value < lowerSet_[iSet] - tolerance || value > upperSet_[iSet] + tolerance) {
1997 // infeasible
1998 kphase = 0;
1999 // remember bounds are flipped so opposite to natural
2000 if (value < lowerSet_[iSet] - tolerance)
2001 cost[iBasic] = 1.0;
2002 else
2003 cost[iBasic] = -1.0;
2004 CoinZeroN(cost, numberInSet);
2005 double dualTolerance = model_->dualTolerance();
2006 for (int iphase = kphase; iphase < 2; iphase++) {
2007 if (iphase) {
2008 cost[numberInSet] = 0.0;
2009 for (int j = 0; j < numberInSet; j++)
2010 cost[j] = cost_[back[j]];
2011 }
2012 // now do one row lp
2013 bool improve = true;
2014 while (improve) {
2015 improve = false;
2016 double dual = cost[iBasic];
2017 int chosen = -1;
2018 double best = dualTolerance;
2019 int way = 0;
2020 for (int i = 0; i <= numberInSet; i++) {
2021 double dj = cost[i] - dual;
2022 double improvement = 0.0;
2023 if (iphase || i < numberInSet)
2024 assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
2025 if (dj > dualTolerance)
2026 improvement = dj * (solution[i] - lower[i]);
2027 else if (dj < -dualTolerance)
2028 improvement = dj * (solution[i] - upper[i]);
2029 if (improvement > best) {
2030 best = improvement;
2031 chosen = i;
2032 if (dj < 0.0) {
2033 way = 1;
2034 } else {
2035 way = -1;
2036 }
2037 }
2038 }
2039 if (chosen >= 0) {
2040 improve = true;
2041 // now see how far
2042 if (way > 0) {
2043 // incoming increasing so basic decreasing
2044 // if phase 0 then go to nearest bound
2045 double distance = upper[chosen] - solution[chosen];
2046 double basicDistance;
2047 if (!iphase) {
2048 assert (iBasic == numberInSet);
2049 assert (solution[iBasic] > upper[iBasic]);
2050 basicDistance = solution[iBasic] - upper[iBasic];
2051 } else {
2052 basicDistance = solution[iBasic] - lower[iBasic];
2053 }
2054 // need extra coding for unbounded
2055 assert (CoinMin(distance, basicDistance) < 1.0e20);
2056 if (distance > basicDistance) {
2057 // incoming becomes basic
2058 solution[chosen] += basicDistance;
2059 if (!iphase)
2060 solution[iBasic] = upper[iBasic];
2061 else
2062 solution[iBasic] = lower[iBasic];
2063 iBasic = chosen;
2064 } else {
2065 // flip
2066 solution[chosen] = upper[chosen];
2067 solution[iBasic] -= distance;
2068 }
2069 } else {
2070 // incoming decreasing so basic increasing
2071 // if phase 0 then go to nearest bound
2072 double distance = solution[chosen] - lower[chosen];
2073 double basicDistance;
2074 if (!iphase) {
2075 assert (iBasic == numberInSet);
2076 assert (solution[iBasic] < lower[iBasic]);
2077 basicDistance = lower[iBasic] - solution[iBasic];
2078 } else {
2079 basicDistance = upper[iBasic] - solution[iBasic];
2080 }
2081 // need extra coding for unbounded - for now just exit
2082 if (CoinMin(distance, basicDistance) > 1.0e20) {
2083 printf("unbounded on set %d\n", iSet);
2084 iphase = 1;
2085 iBasic = numberInSet;
2086 break;
2087 }
2088 if (distance > basicDistance) {
2089 // incoming becomes basic
2090 solution[chosen] -= basicDistance;
2091 if (!iphase)
2092 solution[iBasic] = lower[iBasic];
2093 else
2094 solution[iBasic] = upper[iBasic];
2095 iBasic = chosen;
2096 } else {
2097 // flip
2098 solution[chosen] = lower[chosen];
2099 solution[iBasic] += distance;
2100 }
2101 }
2102 if (!iphase) {
2103 if(iBasic < numberInSet)
2104 break; // feasible
2105 else if (solution[iBasic] >= lower[iBasic] &&
2106 solution[iBasic] <= upper[iBasic])
2107 break; // feasible (on flip)
2108 }
2109 }
2110 }
2111 }
2112 }
2113 // do solution i.e. bounds
2114 if (columnLower_ || columnUpper_) {
2115 for (int j = 0; j < numberInSet; j++) {
2116 if (j != iBasic) {
2117 objectiveOffset += solution[j] * cost[j];
2118 if (columnLower_ && columnUpper_) {
2119 if (fabs(solution[j] - columnLower_[back[j]]) >
2120 fabs(solution[j] - columnUpper_[back[j]]))
2121 setDynamicStatus(back[j], atUpperBound);
2122 } else if (columnUpper_ && solution[j] > 0.0) {
2123 setDynamicStatus(back[j], atUpperBound);
2124 } else {
2125 setDynamicStatus(back[j], atLowerBound);
2126 assert(!solution[j]);
2127 }
2128 }
2129 }
2130 }
2131 // convert iBasic back and do bounds
2132 if (iBasic == numberInSet) {
2133 // slack basic
2134 setStatus(iSet, ClpSimplex::basic);
2135 iBasic = iSet + maximumGubColumns_;
2136 } else {
2137 iBasic = back[iBasic];
2138 setDynamicStatus(iBasic, soloKey);
2139 // remember bounds flipped
2140 if (upper[numberInSet] == lower[numberInSet])
2141 setStatus(iSet, ClpSimplex::isFixed);
2142 else if (solution[numberInSet] == upper[numberInSet])
2143 setStatus(iSet, ClpSimplex::atLowerBound);
2144 else if (solution[numberInSet] == lower[numberInSet])
2145 setStatus(iSet, ClpSimplex::atUpperBound);
2146 else
2147 abort();
2148 }
2149 keyVariable_[iSet] = iBasic;
2150 }
2151 model_->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
2152 delete [] lower;
2153 delete [] solution;
2154 delete [] upper;
2155 delete [] cost;
2156 delete [] back;
2157 // make sure matrix is in good shape
2158 matrix_->orderMatrix();
2159}
2160// Populates initial matrix from dynamic status
2161void
2162ClpDynamicMatrix::initialProblem()
2163{
2164 int iSet;
2165 double * element = matrix_->getMutableElements();
2166 int * row = matrix_->getMutableIndices();
2167 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
2168 int * length = matrix_->getMutableVectorLengths();
2169 double * cost = model_->objective();
2170 double * solution = model_->primalColumnSolution();
2171 double * columnLower = model_->columnLower();
2172 double * columnUpper = model_->columnUpper();
2173 double * rowSolution = model_->primalRowSolution();
2174 double * rowLower = model_->rowLower();
2175 double * rowUpper = model_->rowUpper();
2176 CoinBigIndex numberElements = startColumn[firstDynamic_];
2177
2178 firstAvailable_ = firstDynamic_;
2179 numberActiveSets_ = 0;
2180 for (iSet = 0; iSet < numberSets_; iSet++) {
2181 toIndex_[iSet] = -1;
2182 int numberActive = 0;
2183 int whichKey = -1;
2184 if (getStatus(iSet) == ClpSimplex::basic) {
2185 whichKey = maximumGubColumns_ + iSet;
2186 numberActive = 1;
2187 } else {
2188 whichKey = -1;
2189 }
2190 int j = startSet_[iSet];
2191 while (j >= 0) {
2192 assert (getDynamicStatus(j) != soloKey || whichKey == -1);
2193 if (getDynamicStatus(j) == inSmall) {
2194 numberActive++;
2195 } else if (getDynamicStatus(j) == soloKey) {
2196 whichKey = j;
2197 numberActive++;
2198 }
2199 j = next_[j]; //onto next in set
2200 }
2201 if (numberActive > 1) {
2202 int iRow = numberActiveSets_ + numberStaticRows_;
2203 rowSolution[iRow] = 0.0;
2204 double lowerValue;
2205 if (lowerSet_[iSet] > -1.0e20)
2206 lowerValue = lowerSet_[iSet];
2207 else
2208 lowerValue = -COIN_DBL_MAX;
2209 double upperValue;
2210 if (upperSet_[iSet] < 1.0e20)
2211 upperValue = upperSet_[iSet];
2212 else
2213 upperValue = COIN_DBL_MAX;
2214 rowLower[iRow] = lowerValue;
2215 rowUpper[iRow] = upperValue;
2216 if (getStatus(iSet) == ClpSimplex::basic) {
2217 model_->setRowStatus(iRow, ClpSimplex::basic);
2218 rowSolution[iRow] = 0.0;
2219 } else if (getStatus(iSet) == ClpSimplex::atLowerBound) {
2220 model_->setRowStatus(iRow, ClpSimplex::atLowerBound);
2221 rowSolution[iRow] = lowerValue;
2222 } else {
2223 model_->setRowStatus(iRow, ClpSimplex::atUpperBound);
2224 rowSolution[iRow] = upperValue;
2225 }
2226 j = startSet_[iSet];
2227 while (j >= 0) {
2228 DynamicStatus status = getDynamicStatus(j);
2229 if (status == inSmall) {
2230 int numberThis = startColumn_[j+1] - startColumn_[j] + 1;
2231 if (numberElements + numberThis > numberElements_) {
2232 // need to redo
2233 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
2234 matrix_->reserve(lastDynamic_, numberElements_);
2235 element = matrix_->getMutableElements();
2236 row = matrix_->getMutableIndices();
2237 // these probably okay but be safe
2238 startColumn = matrix_->getMutableVectorStarts();
2239 length = matrix_->getMutableVectorLengths();
2240 }
2241 length[firstAvailable_] = numberThis;
2242 cost[firstAvailable_] = cost_[j];
2243 CoinBigIndex base = startColumn_[j];
2244 for (int k = 0; k < numberThis - 1; k++) {
2245 row[numberElements] = row_[base+k];
2246 element[numberElements++] = element_[base+k];
2247 }
2248 row[numberElements] = iRow;
2249 element[numberElements++] = 1.0;
2250 id_[firstAvailable_-firstDynamic_] = j;
2251 solution[firstAvailable_] = 0.0;
2252 model_->setStatus(firstAvailable_, ClpSimplex::basic);
2253 if (!columnLower_ && !columnUpper_) {
2254 columnLower[firstAvailable_] = 0.0;
2255 columnUpper[firstAvailable_] = COIN_DBL_MAX;
2256 } else {
2257 if (columnLower_)
2258 columnLower[firstAvailable_] = columnLower_[j];
2259 else
2260 columnLower[firstAvailable_] = 0.0;
2261 if (columnUpper_)
2262 columnUpper[firstAvailable_] = columnUpper_[j];
2263 else
2264 columnUpper[firstAvailable_] = COIN_DBL_MAX;
2265 if (status != atUpperBound) {
2266 solution[firstAvailable_] = columnLower[firstAvailable_];
2267 } else {
2268 solution[firstAvailable_] = columnUpper[firstAvailable_];
2269 }
2270 }
2271 firstAvailable_++;
2272 startColumn[firstAvailable_] = numberElements;
2273 }
2274 j = next_[j]; //onto next in set
2275 }
2276 model_->setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet));
2277 toIndex_[iSet] = numberActiveSets_;
2278 fromIndex_[numberActiveSets_++] = iSet;
2279 } else {
2280 // solo key
2281 bool needKey=false;
2282 if (numberActive) {
2283 if (whichKey<maximumGubColumns_) {
2284 // structural - assume ok
2285 needKey = false;
2286 } else {
2287 // slack
2288 keyVariable_[iSet] = maximumGubColumns_ + iSet;
2289 double value = keyValue(iSet);
2290 if (value<lowerSet_[iSet]-1.0e-8||
2291 value>upperSet_[iSet]+1.0e-8)
2292 needKey=true;
2293 }
2294 } else {
2295 needKey = true;
2296 }
2297 if (needKey) {
2298 // all to lb then up some (slack/null if possible)
2299 int length=99999999;
2300 int which=-1;
2301 double sum=0.0;
2302 for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) {
2303 setDynamicStatus(iColumn,atLowerBound);
2304 sum += columnLower_[iColumn];
2305 if (length>startColumn_[iColumn+1]-startColumn_[iColumn]) {
2306 which=iColumn;
2307 length=startColumn_[iColumn+1]-startColumn_[iColumn];
2308 }
2309 }
2310 if (sum>lowerSet_[iSet]-1.0e-8) {
2311 // slack can be basic
2312 setStatus(iSet,ClpSimplex::basic);
2313 keyVariable_[iSet] = maximumGubColumns_ + iSet;
2314 } else {
2315 // use shortest
2316 setDynamicStatus(which,soloKey);
2317 keyVariable_[iSet] = which;
2318 setStatus(iSet,ClpSimplex::atLowerBound);
2319 }
2320 }
2321 }
2322 assert (toIndex_[iSet] >= 0 || whichKey >= 0);
2323 keyVariable_[iSet] = whichKey;
2324 }
2325 // clean up pivotVariable
2326 int numberColumns = model_->numberColumns();
2327 int numberRows = model_->numberRows();
2328 int * pivotVariable = model_->pivotVariable();
2329 if (pivotVariable) {
2330 for (int i=0; i<numberStaticRows_+numberActiveSets_;i++) {
2331 if (model_->getRowStatus(i)!=ClpSimplex::basic)
2332 pivotVariable[i]=-1;
2333 else
2334 pivotVariable[i]=numberColumns+i;
2335 }
2336 for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
2337 pivotVariable[i]=i+numberColumns;
2338 }
2339 int put=-1;
2340 for (int i=0;i<numberColumns;i++) {
2341 if (model_->getColumnStatus(i)==ClpSimplex::basic) {
2342 while(put<numberRows) {
2343 put++;
2344 if (pivotVariable[put]==-1) {
2345 pivotVariable[put]=i;
2346 break;
2347 }
2348 }
2349 }
2350 }
2351 for (int i=CoinMax(put,0);i<numberRows;i++) {
2352 if (pivotVariable[i]==-1)
2353 pivotVariable[i]=i+numberColumns;
2354 }
2355 }
2356 if (rhsOffset_) {
2357 double * cost = model_->costRegion();
2358 double * columnLower = model_->lowerRegion();
2359 double * columnUpper = model_->upperRegion();
2360 double * solution = model_->solutionRegion();
2361 int numberRows = model_->numberRows();
2362 for (int i = numberActiveSets_; i < numberRows-numberStaticRows_; i++) {
2363 int iSequence = i + numberStaticRows_ + numberColumns;
2364 solution[iSequence] = 0.0;
2365 columnLower[iSequence] = -COIN_DBL_MAX;
2366 columnUpper[iSequence] = COIN_DBL_MAX;
2367 cost[iSequence] = 0.0;
2368 model_->nonLinearCost()->setOne(iSequence, solution[iSequence],
2369 columnLower[iSequence],
2370 columnUpper[iSequence], 0.0);
2371 model_->setStatus(iSequence, ClpSimplex::basic);
2372 rhsOffset_[i+numberStaticRows_] = 0.0;
2373 }
2374#if 0
2375 for (int i=0;i<numberStaticRows_;i++)
2376 printf("%d offset %g\n",
2377 i,rhsOffset_[i]);
2378#endif
2379 }
2380 numberActiveColumns_ = firstAvailable_;
2381#if 0
2382 for (iSet = 0; iSet < numberSets_; iSet++) {
2383 for (int j=startSet_[iSet];j<startSet_[iSet+1];j++) {
2384 if (getDynamicStatus(j)==soloKey)
2385 printf("%d in set %d is solo key\n",j,iSet);
2386 else if (getDynamicStatus(j)==inSmall)
2387 printf("%d in set %d is in small\n",j,iSet);
2388 }
2389 }
2390#endif
2391 return;
2392}
2393// Writes out model (without names)
2394void
2395ClpDynamicMatrix::writeMps(const char * name)
2396{
2397 int numberTotalRows = numberStaticRows_+numberSets_;
2398 int numberTotalColumns = firstDynamic_+numberGubColumns_;
2399 // over estimate
2400 int numberElements = getNumElements()+startColumn_[numberGubColumns_]
2401 + numberGubColumns_;
2402 double * columnLower = new double [numberTotalColumns];
2403 double * columnUpper = new double [numberTotalColumns];
2404 double * cost = new double [numberTotalColumns];
2405 double * rowLower = new double [numberTotalRows];
2406 double * rowUpper = new double [numberTotalRows];
2407 CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1];
2408 int * row = new int [numberElements];
2409 double * element = new double [numberElements];
2410 // Fill in
2411 const CoinBigIndex * startA = getVectorStarts();
2412 const int * lengthA = getVectorLengths();
2413 const int * rowA = getIndices();
2414 const double * elementA = getElements();
2415 const double * columnLowerA = model_->columnLower();
2416 const double * columnUpperA = model_->columnUpper();
2417 const double * costA = model_->objective();
2418 const double * rowLowerA = model_->rowLower();
2419 const double * rowUpperA = model_->rowUpper();
2420 start[0]=0;
2421 numberElements=0;
2422 for (int i=0;i<firstDynamic_;i++) {
2423 columnLower[i] = columnLowerA[i];
2424 columnUpper[i] = columnUpperA[i];
2425 cost[i] = costA[i];
2426 for (CoinBigIndex j = startA[i];j<startA[i]+lengthA[i];j++) {
2427 row[numberElements] = rowA[j];
2428 element[numberElements++]=elementA[j];
2429 }
2430 start[i+1]=numberElements;
2431 }
2432 for (int i=0;i<numberStaticRows_;i++) {
2433 rowLower[i] = rowLowerA[i];
2434 rowUpper[i] = rowUpperA[i];
2435 }
2436 int putC=firstDynamic_;
2437 int putR=numberStaticRows_;
2438 for (int i=0;i<numberSets_;i++) {
2439 rowLower[putR]=lowerSet_[i];
2440 rowUpper[putR]=upperSet_[i];
2441 for (CoinBigIndex k=startSet_[i];k<startSet_[i+1];k++) {
2442 columnLower[putC]=columnLower_[k];
2443 columnUpper[putC]=columnUpper_[k];
2444 cost[putC]=cost_[k];
2445 putC++;
2446 for (CoinBigIndex j = startColumn_[k];j<startColumn_[k+1];j++) {
2447 row[numberElements] = row_[j];
2448 element[numberElements++]=element_[j];
2449 }
2450 row[numberElements] = putR;
2451 element[numberElements++]=1.0;
2452 start[putC]=numberElements;
2453 }
2454 putR++;
2455 }
2456
2457 assert (putR==numberTotalRows);
2458 assert (putC==numberTotalColumns);
2459 ClpSimplex modelOut;
2460 modelOut.loadProblem(numberTotalColumns,numberTotalRows,
2461 start,row,element,
2462 columnLower,columnUpper,cost,
2463 rowLower,rowUpper);
2464 modelOut.writeMps(name);
2465 delete [] columnLower;
2466 delete [] columnUpper;
2467 delete [] cost;
2468 delete [] rowLower;
2469 delete [] rowUpper;
2470 delete [] start;
2471 delete [] row;
2472 delete [] element;
2473}
2474// Adds in a column to gub structure (called from descendant)
2475int
2476ClpDynamicMatrix::addColumn(int numberEntries, const int * row, const double * element,
2477 double cost, double lower, double upper, int iSet,
2478 DynamicStatus status)
2479{
2480 // check if already in
2481 int j = startSet_[iSet];
2482 while (j >= 0) {
2483 if (startColumn_[j+1] - startColumn_[j] == numberEntries) {
2484 const int * row2 = row_ + startColumn_[j];
2485 const double * element2 = element_ + startColumn_[j];
2486 bool same = true;
2487 for (int k = 0; k < numberEntries; k++) {
2488 if (row[k] != row2[k] || element[k] != element2[k]) {
2489 same = false;
2490 break;
2491 }
2492 }
2493 if (same) {
2494 bool odd = false;
2495 if (cost != cost_[j])
2496 odd = true;
2497 if (columnLower_ && lower != columnLower_[j])
2498 odd = true;
2499 if (columnUpper_ && upper != columnUpper_[j])
2500 odd = true;
2501 if (odd) {
2502 printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
2503 cost, lower, upper, cost_[j],
2504 columnLower_ ? columnLower_[j] : 0.0,
2505 columnUpper_ ? columnUpper_[j] : 1.0e100);
2506 } else {
2507 setDynamicStatus(j, status);
2508 return j;
2509 }
2510 }
2511 }
2512 j = next_[j];
2513 }
2514
2515 if (numberGubColumns_ == maximumGubColumns_ ||
2516 startColumn_[numberGubColumns_] + numberEntries > maximumElements_) {
2517 CoinBigIndex j;
2518 int i;
2519 int put = 0;
2520 int numberElements = 0;
2521 CoinBigIndex start = 0;
2522 // compress - leave ones at ub and basic
2523 int * which = new int [numberGubColumns_];
2524 for (i = 0; i < numberGubColumns_; i++) {
2525 CoinBigIndex end = startColumn_[i+1];
2526 // what about ubs if column generation?
2527 if (getDynamicStatus(i) != atLowerBound) {
2528 // keep in
2529 for (j = start; j < end; j++) {
2530 row_[numberElements] = row_[j];
2531 element_[numberElements++] = element_[j];
2532 }
2533 startColumn_[put+1] = numberElements;
2534 cost_[put] = cost_[i];
2535 if (columnLower_)
2536 columnLower_[put] = columnLower_[i];
2537 if (columnUpper_)
2538 columnUpper_[put] = columnUpper_[i];
2539 dynamicStatus_[put] = dynamicStatus_[i];
2540 id_[put] = id_[i];
2541 which[i] = put;
2542 put++;
2543 } else {
2544 // out
2545 which[i] = -1;
2546 }
2547 start = end;
2548 }
2549 // now redo startSet_ and next_
2550 int * newNext = new int [maximumGubColumns_];
2551 for (int jSet = 0; jSet < numberSets_; jSet++) {
2552 int sequence = startSet_[jSet];
2553 while (which[sequence] < 0) {
2554 // out
2555 assert (next_[sequence] >= 0);
2556 sequence = next_[sequence];
2557 }
2558 startSet_[jSet] = which[sequence];
2559 int last = which[sequence];
2560 while (next_[sequence] >= 0) {
2561 sequence = next_[sequence];
2562 if(which[sequence] >= 0) {
2563 // keep
2564 newNext[last] = which[sequence];
2565 last = which[sequence];
2566 }
2567 }
2568 newNext[last] = -jSet - 1;
2569 }
2570 delete [] next_;
2571 next_ = newNext;
2572 delete [] which;
2573 abort();
2574 }
2575 CoinBigIndex start = startColumn_[numberGubColumns_];
2576 CoinMemcpyN(row, numberEntries, row_ + start);
2577 CoinMemcpyN(element, numberEntries, element_ + start);
2578 startColumn_[numberGubColumns_+1] = start + numberEntries;
2579 cost_[numberGubColumns_] = cost;
2580 if (columnLower_)
2581 columnLower_[numberGubColumns_] = lower;
2582 else
2583 assert (!lower);
2584 if (columnUpper_)
2585 columnUpper_[numberGubColumns_] = upper;
2586 else
2587 assert (upper > 1.0e20);
2588 setDynamicStatus(numberGubColumns_, status);
2589 // Do next_
2590 j = startSet_[iSet];
2591 startSet_[iSet] = numberGubColumns_;
2592 next_[numberGubColumns_] = j;
2593 numberGubColumns_++;
2594 return numberGubColumns_ - 1;
2595}
2596// Returns which set a variable is in
2597int
2598ClpDynamicMatrix::whichSet (int sequence) const
2599{
2600 while (next_[sequence] >= 0)
2601 sequence = next_[sequence];
2602 int iSet = - next_[sequence] - 1;
2603 return iSet;
2604}
2605