1/* $Id: ClpGubDynamicMatrix.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2002, 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 "ClpGubDynamicMatrix.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//-------------------------------------------------------------------
29ClpGubDynamicMatrix::ClpGubDynamicMatrix ()
30 : ClpGubMatrix(),
31 objectiveOffset_(0.0),
32 startColumn_(NULL),
33 row_(NULL),
34 element_(NULL),
35 cost_(NULL),
36 fullStart_(NULL),
37 id_(NULL),
38 dynamicStatus_(NULL),
39 lowerColumn_(NULL),
40 upperColumn_(NULL),
41 lowerSet_(NULL),
42 upperSet_(NULL),
43 numberGubColumns_(0),
44 firstAvailable_(0),
45 savedFirstAvailable_(0),
46 firstDynamic_(0),
47 lastDynamic_(0),
48 numberElements_(0)
49{
50 setType(13);
51}
52
53//-------------------------------------------------------------------
54// Copy constructor
55//-------------------------------------------------------------------
56ClpGubDynamicMatrix::ClpGubDynamicMatrix (const ClpGubDynamicMatrix & rhs)
57 : ClpGubMatrix(rhs)
58{
59 objectiveOffset_ = rhs.objectiveOffset_;
60 numberGubColumns_ = rhs.numberGubColumns_;
61 firstAvailable_ = rhs.firstAvailable_;
62 savedFirstAvailable_ = rhs.savedFirstAvailable_;
63 firstDynamic_ = rhs.firstDynamic_;
64 lastDynamic_ = rhs.lastDynamic_;
65 numberElements_ = rhs.numberElements_;
66 startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
67 CoinBigIndex numberElements = startColumn_[numberGubColumns_];
68 row_ = ClpCopyOfArray(rhs.row_, numberElements);
69 element_ = ClpCopyOfArray(rhs.element_, numberElements);
70 cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
71 fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
72 id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
73 lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
74 upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
75 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
76 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
77 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
78}
79
80/* This is the real constructor*/
81ClpGubDynamicMatrix::ClpGubDynamicMatrix(ClpSimplex * model, int numberSets,
82 int numberGubColumns, const int * starts,
83 const double * lower, const double * upper,
84 const CoinBigIndex * startColumn, const int * row,
85 const double * element, const double * cost,
86 const double * lowerColumn, const double * upperColumn,
87 const unsigned char * status)
88 : ClpGubMatrix()
89{
90 objectiveOffset_ = model->objectiveOffset();
91 model_ = model;
92 numberSets_ = numberSets;
93 numberGubColumns_ = numberGubColumns;
94 fullStart_ = ClpCopyOfArray(starts, numberSets_ + 1);
95 lower_ = ClpCopyOfArray(lower, numberSets_);
96 upper_ = ClpCopyOfArray(upper, numberSets_);
97 int numberColumns = model->numberColumns();
98 int numberRows = model->numberRows();
99 // Number of columns needed
100 int numberGubInSmall = numberSets_ + numberRows + 2 * model->factorizationFrequency() + 2;
101 // for small problems this could be too big
102 //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
103 int numberNeeded = numberGubInSmall + numberColumns;
104 firstAvailable_ = numberColumns;
105 savedFirstAvailable_ = numberColumns;
106 firstDynamic_ = numberColumns;
107 lastDynamic_ = numberNeeded;
108 startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1);
109 CoinBigIndex numberElements = startColumn_[numberGubColumns_];
110 row_ = ClpCopyOfArray(row, numberElements);
111 element_ = new double[numberElements];
112 CoinBigIndex i;
113 for (i = 0; i < numberElements; i++)
114 element_[i] = element[i];
115 cost_ = new double[numberGubColumns_];
116 for (i = 0; i < numberGubColumns_; i++) {
117 cost_[i] = cost[i];
118 // need sorted
119 CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i+1], element_ + startColumn_[i]);
120 }
121 if (lowerColumn) {
122 lowerColumn_ = new double[numberGubColumns_];
123 for (i = 0; i < numberGubColumns_; i++)
124 lowerColumn_[i] = lowerColumn[i];
125 } else {
126 lowerColumn_ = NULL;
127 }
128 if (upperColumn) {
129 upperColumn_ = new double[numberGubColumns_];
130 for (i = 0; i < numberGubColumns_; i++)
131 upperColumn_[i] = upperColumn[i];
132 } else {
133 upperColumn_ = NULL;
134 }
135 if (upperColumn || lowerColumn) {
136 lowerSet_ = new double[numberSets_];
137 for (i = 0; i < numberSets_; i++) {
138 if (lower[i] > -1.0e20)
139 lowerSet_[i] = lower[i];
140 else
141 lowerSet_[i] = -1.0e30;
142 }
143 upperSet_ = new double[numberSets_];
144 for (i = 0; i < numberSets_; i++) {
145 if (upper[i] < 1.0e20)
146 upperSet_[i] = upper[i];
147 else
148 upperSet_[i] = 1.0e30;
149 }
150 } else {
151 lowerSet_ = NULL;
152 upperSet_ = NULL;
153 }
154 start_ = NULL;
155 end_ = NULL;
156 dynamicStatus_ = NULL;
157 id_ = new int[numberGubInSmall];
158 for (i = 0; i < numberGubInSmall; i++)
159 id_[i] = -1;
160 ClpPackedMatrix* originalMatrixA =
161 dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
162 assert (originalMatrixA);
163 CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
164 originalMatrixA->setMatrixNull(); // so can be deleted safely
165 // guess how much space needed
166 double guess = originalMatrix->getNumElements() + 10;
167 guess /= static_cast<double> (numberColumns);
168 guess *= 2 * numberGubColumns_;
169 numberElements_ = static_cast<int> (CoinMin(guess, 10000000.0));
170 numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements();
171 matrix_ = originalMatrix;
172 flags_ &= ~1;
173 // resize model (matrix stays same)
174 model->resize(numberRows, numberNeeded);
175 if (upperColumn_) {
176 // set all upper bounds so we have enough space
177 double * columnUpper = model->columnUpper();
178 for(i = firstDynamic_; i < lastDynamic_; i++)
179 columnUpper[i] = 1.0e10;
180 }
181 // resize matrix
182 // extra 1 is so can keep number of elements handy
183 originalMatrix->reserve(numberNeeded, numberElements_, true);
184 originalMatrix->reserve(numberNeeded + 1, numberElements_, false);
185 originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements();
186 // redo number of columns
187 numberColumns = matrix_->getNumCols();
188 backward_ = new int[numberNeeded];
189 backToPivotRow_ = new int[numberNeeded];
190 // We know a bit better
191 delete [] changeCost_;
192 changeCost_ = new double [numberRows+numberSets_];
193 keyVariable_ = new int[numberSets_];
194 // signal to need new ordering
195 next_ = NULL;
196 for (int iColumn = 0; iColumn < numberNeeded; iColumn++)
197 backward_[iColumn] = -1;
198
199 firstGub_ = firstDynamic_;
200 lastGub_ = lastDynamic_;
201 if (!lowerColumn_ && !upperColumn_)
202 gubType_ = 8;
203 if (status) {
204 status_ = ClpCopyOfArray(status, numberSets_);
205 } else {
206 status_ = new unsigned char [numberSets_];
207 memset(status_, 0, numberSets_);
208 int i;
209 for (i = 0; i < numberSets_; i++) {
210 // make slack key
211 setStatus(i, ClpSimplex::basic);
212 }
213 }
214 saveStatus_ = new unsigned char [numberSets_];
215 memset(saveStatus_, 0, numberSets_);
216 savedKeyVariable_ = new int [numberSets_];
217 memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
218}
219
220//-------------------------------------------------------------------
221// Destructor
222//-------------------------------------------------------------------
223ClpGubDynamicMatrix::~ClpGubDynamicMatrix ()
224{
225 delete [] startColumn_;
226 delete [] row_;
227 delete [] element_;
228 delete [] cost_;
229 delete [] fullStart_;
230 delete [] id_;
231 delete [] dynamicStatus_;
232 delete [] lowerColumn_;
233 delete [] upperColumn_;
234 delete [] lowerSet_;
235 delete [] upperSet_;
236}
237
238//----------------------------------------------------------------
239// Assignment operator
240//-------------------------------------------------------------------
241ClpGubDynamicMatrix &
242ClpGubDynamicMatrix::operator=(const ClpGubDynamicMatrix& rhs)
243{
244 if (this != &rhs) {
245 ClpGubMatrix::operator=(rhs);
246 delete [] startColumn_;
247 delete [] row_;
248 delete [] element_;
249 delete [] cost_;
250 delete [] fullStart_;
251 delete [] id_;
252 delete [] dynamicStatus_;
253 delete [] lowerColumn_;
254 delete [] upperColumn_;
255 delete [] lowerSet_;
256 delete [] upperSet_;
257 objectiveOffset_ = rhs.objectiveOffset_;
258 numberGubColumns_ = rhs.numberGubColumns_;
259 firstAvailable_ = rhs.firstAvailable_;
260 savedFirstAvailable_ = rhs.savedFirstAvailable_;
261 firstDynamic_ = rhs.firstDynamic_;
262 lastDynamic_ = rhs.lastDynamic_;
263 numberElements_ = rhs.numberElements_;
264 startColumn_ = ClpCopyOfArray(rhs.startColumn_, numberGubColumns_ + 1);
265 int numberElements = startColumn_[numberGubColumns_];
266 row_ = ClpCopyOfArray(rhs.row_, numberElements);
267 element_ = ClpCopyOfArray(rhs.element_, numberElements);
268 cost_ = ClpCopyOfArray(rhs.cost_, numberGubColumns_);
269 fullStart_ = ClpCopyOfArray(rhs.fullStart_, numberSets_ + 1);
270 id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_);
271 lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_, numberGubColumns_);
272 upperColumn_ = ClpCopyOfArray(rhs.upperColumn_, numberGubColumns_);
273 dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, numberGubColumns_);
274 lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
275 upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
276 }
277 return *this;
278}
279//-------------------------------------------------------------------
280// Clone
281//-------------------------------------------------------------------
282ClpMatrixBase * ClpGubDynamicMatrix::clone() const
283{
284 return new ClpGubDynamicMatrix(*this);
285}
286// Partial pricing
287void
288ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
289 int & bestSequence, int & numberWanted)
290{
291 assert(!model->rowScale());
292 numberWanted = currentWanted_;
293 if (!numberSets_) {
294 // no gub
295 ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
296 return;
297 } else {
298 // and do some proportion of full set
299 int startG2 = static_cast<int> (startFraction * numberSets_);
300 int endG2 = static_cast<int> (endFraction * numberSets_ + 0.1);
301 endG2 = CoinMin(endG2, numberSets_);
302 //printf("gub price - set start %d end %d\n",
303 // startG2,endG2);
304 double tolerance = model->currentDualTolerance();
305 double * reducedCost = model->djRegion();
306 const double * duals = model->dualRowSolution();
307 double * cost = model->costRegion();
308 double bestDj;
309 int numberRows = model->numberRows();
310 int numberColumns = lastDynamic_;
311 // If nothing found yet can go all the way to end
312 int endAll = endG2;
313 if (bestSequence < 0 && !startG2)
314 endAll = numberSets_;
315 if (bestSequence >= 0)
316 bestDj = fabs(reducedCost[bestSequence]);
317 else
318 bestDj = tolerance;
319 int saveSequence = bestSequence;
320 double djMod = 0.0;
321 double infeasibilityCost = model->infeasibilityCost();
322 double bestDjMod = 0.0;
323 //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
324 // startG2,endG2,numberWanted);
325 int bestType = -1;
326 int bestSet = -1;
327 const double * element = matrix_->getElements();
328 const int * row = matrix_->getIndices();
329 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
330 int * length = matrix_->getMutableVectorLengths();
331#if 0
332 // make sure first available is clean (in case last iteration rejected)
333 cost[firstAvailable_] = 0.0;
334 length[firstAvailable_] = 0;
335 model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0);
336 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
337 {
338 for (int i = firstAvailable_; i < lastDynamic_; i++)
339 assert(!cost[i]);
340 }
341#endif
342#ifdef CLP_DEBUG
343 {
344 for (int i = firstDynamic_; i < firstAvailable_; i++) {
345 assert (getDynamicStatus(id_[i-firstDynamic_]) == inSmall);
346 }
347 }
348#endif
349 int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
350 int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_;
351 for (int iSet = startG2; iSet < endAll; iSet++) {
352 if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) {
353 // give up
354 numberWanted = 0;
355 break;
356 } else if (iSet == endG2 && bestSequence >= 0) {
357 break;
358 }
359 CoinBigIndex j;
360 int iBasic = keyVariable_[iSet];
361 if (iBasic >= numberColumns) {
362 djMod = - weight(iSet) * infeasibilityCost;
363 } else {
364 // get dj without
365 assert (model->getStatus(iBasic) == ClpSimplex::basic);
366 djMod = 0.0;
367
368 for (j = startColumn[iBasic];
369 j < startColumn[iBasic] + length[iBasic]; j++) {
370 int jRow = row[j];
371 djMod -= duals[jRow] * element[j];
372 }
373 djMod += cost[iBasic];
374 // See if gub slack possible - dj is djMod
375 if (getStatus(iSet) == ClpSimplex::atLowerBound) {
376 double value = -djMod;
377 if (value > tolerance) {
378 numberWanted--;
379 if (value > bestDj) {
380 // check flagged variable and correct dj
381 if (!flagged(iSet)) {
382 bestDj = value;
383 bestSequence = numberRows + numberColumns + iSet;
384 bestDjMod = djMod;
385 bestType = 0;
386 bestSet = iSet;
387 } else {
388 // just to make sure we don't exit before got something
389 numberWanted++;
390 abort();
391 }
392 }
393 }
394 } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
395 double value = djMod;
396 if (value > tolerance) {
397 numberWanted--;
398 if (value > bestDj) {
399 // check flagged variable and correct dj
400 if (!flagged(iSet)) {
401 bestDj = value;
402 bestSequence = numberRows + numberColumns + iSet;
403 bestDjMod = djMod;
404 bestType = 0;
405 bestSet = iSet;
406 } else {
407 // just to make sure we don't exit before got something
408 numberWanted++;
409 abort();
410 }
411 }
412 }
413 }
414 }
415 for (int iSequence = fullStart_[iSet]; iSequence < fullStart_[iSet+1]; iSequence++) {
416 DynamicStatus status = getDynamicStatus(iSequence);
417 if (status != inSmall) {
418 double value = cost_[iSequence] - djMod;
419 for (j = startColumn_[iSequence];
420 j < startColumn_[iSequence+1]; j++) {
421 int jRow = row_[j];
422 value -= duals[jRow] * element_[j];
423 }
424 // change sign if at lower bound
425 if (status == atLowerBound)
426 value = -value;
427 if (value > tolerance) {
428 numberWanted--;
429 if (value > bestDj) {
430 // check flagged variable and correct dj
431 if (!flagged(iSequence)) {
432 bestDj = value;
433 bestSequence = iSequence;
434 bestDjMod = djMod;
435 bestType = 1;
436 bestSet = iSet;
437 } else {
438 // just to make sure we don't exit before got something
439 numberWanted++;
440 }
441 }
442 }
443 }
444 }
445 if (numberWanted <= 0) {
446 numberWanted = 0;
447 break;
448 }
449 }
450 // Do packed part before gub and small gub - but lightly
451 int saveMinNeg = minimumGoodReducedCosts_;
452 int saveSequence2 = bestSequence;
453 if (bestSequence >= 0)
454 minimumGoodReducedCosts_ = -2;
455 int saveLast = lastGub_;
456 lastGub_ = firstAvailable_;
457 currentWanted_ = numberWanted;
458 ClpGubMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
459 minimumGoodReducedCosts_ = saveMinNeg;
460 lastGub_ = saveLast;
461 if (bestSequence != saveSequence2) {
462 bestType = -1; // in normal or small gub part
463 saveSequence = bestSequence;
464 }
465 if (bestSequence != saveSequence || bestType >= 0) {
466 double * lowerColumn = model->lowerRegion();
467 double * upperColumn = model->upperRegion();
468 double * solution = model->solutionRegion();
469 if (bestType > 0) {
470 // recompute dj and create
471 double value = cost_[bestSequence] - bestDjMod;
472 for (CoinBigIndex jBigIndex = startColumn_[bestSequence];
473 jBigIndex < startColumn_[bestSequence+1]; jBigIndex++) {
474 int jRow = row_[jBigIndex];
475 value -= duals[jRow] * element_[jBigIndex];
476 }
477 double * element = matrix_->getMutableElements();
478 int * row = matrix_->getMutableIndices();
479 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
480 int * length = matrix_->getMutableVectorLengths();
481 CoinBigIndex numberElements = startColumn[firstAvailable_];
482 int numberThis = startColumn_[bestSequence+1] - startColumn_[bestSequence];
483 if (numberElements + numberThis > numberElements_) {
484 // need to redo
485 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
486 matrix_->reserve(numberColumns, numberElements_);
487 element = matrix_->getMutableElements();
488 row = matrix_->getMutableIndices();
489 // these probably okay but be safe
490 startColumn = matrix_->getMutableVectorStarts();
491 length = matrix_->getMutableVectorLengths();
492 }
493 // already set startColumn[firstAvailable_]=numberElements;
494 length[firstAvailable_] = numberThis;
495 model->costRegion()[firstAvailable_] = cost_[bestSequence];
496 CoinBigIndex base = startColumn_[bestSequence];
497 for (int j = 0; j < numberThis; j++) {
498 row[numberElements] = row_[base+j];
499 element[numberElements++] = element_[base+j];
500 }
501 id_[firstAvailable_-firstDynamic_] = bestSequence;
502 //printf("best %d\n",bestSequence);
503 backward_[firstAvailable_] = bestSet;
504 model->solutionRegion()[firstAvailable_] = 0.0;
505 if (!lowerColumn_ && !upperColumn_) {
506 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
507 lowerColumn[firstAvailable_] = 0.0;
508 upperColumn[firstAvailable_] = COIN_DBL_MAX;
509 } else {
510 DynamicStatus status = getDynamicStatus(bestSequence);
511 if (lowerColumn_)
512 lowerColumn[firstAvailable_] = lowerColumn_[bestSequence];
513 else
514 lowerColumn[firstAvailable_] = 0.0;
515 if (upperColumn_)
516 upperColumn[firstAvailable_] = upperColumn_[bestSequence];
517 else
518 upperColumn[firstAvailable_] = COIN_DBL_MAX;
519 if (status == atLowerBound) {
520 solution[firstAvailable_] = lowerColumn[firstAvailable_];
521 model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
522 } else {
523 solution[firstAvailable_] = upperColumn[firstAvailable_];
524 model->setStatus(firstAvailable_, ClpSimplex::atUpperBound);
525 }
526 }
527 model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_],
528 lowerColumn[firstAvailable_],
529 upperColumn[firstAvailable_], cost_[bestSequence]);
530 bestSequence = firstAvailable_;
531 // firstAvailable_ only updated if good pivot (in updatePivot)
532 startColumn[firstAvailable_+1] = numberElements;
533 //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
534 reducedCost[bestSequence] = value;
535 gubSlackIn_ = -1;
536 } else {
537 // slack - make last column
538 gubSlackIn_ = bestSequence - numberRows - numberColumns;
539 bestSequence = numberColumns + 2 * numberRows;
540 reducedCost[bestSequence] = bestDjMod;
541 //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
542 model->setStatus(bestSequence, getStatus(gubSlackIn_));
543 if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
544 solution[bestSequence] = upper_[gubSlackIn_];
545 else
546 solution[bestSequence] = lower_[gubSlackIn_];
547 lowerColumn[bestSequence] = lower_[gubSlackIn_];
548 upperColumn[bestSequence] = upper_[gubSlackIn_];
549 model->costRegion()[bestSequence] = 0.0;
550 model->nonLinearCost()->setOne(bestSequence, solution[bestSequence], lowerColumn[bestSequence],
551 upperColumn[bestSequence], 0.0);
552 }
553 savedBestSequence_ = bestSequence;
554 savedBestDj_ = reducedCost[savedBestSequence_];
555 }
556 // See if may be finished
557 if (!startG2 && bestSequence < 0)
558 infeasibilityWeight_ = model_->infeasibilityCost();
559 else if (bestSequence >= 0)
560 infeasibilityWeight_ = -1.0;
561 }
562 currentWanted_ = numberWanted;
563}
564// This is local to Gub to allow synchronization when status is good
565int
566ClpGubDynamicMatrix::synchronize(ClpSimplex * model, int mode)
567{
568 int returnNumber = 0;
569 switch (mode) {
570 case 0: {
571#ifdef CLP_DEBUG
572 {
573 for (int i = 0; i < numberSets_; i++)
574 assert(toIndex_[i] == -1);
575 }
576#endif
577 // lookup array
578 int * lookup = new int[lastDynamic_];
579 int iColumn;
580 int numberColumns = model->numberColumns();
581 double * element = matrix_->getMutableElements();
582 int * row = matrix_->getMutableIndices();
583 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
584 int * length = matrix_->getMutableVectorLengths();
585 double * cost = model->costRegion();
586 double * lowerColumn = model->lowerRegion();
587 double * upperColumn = model->upperRegion();
588 int * pivotVariable = model->pivotVariable();
589 CoinBigIndex numberElements = startColumn[firstDynamic_];
590 // first just do lookup and basic stuff
591 int currentNumber = firstAvailable_;
592 firstAvailable_ = firstDynamic_;
593 int numberToDo = 0;
594 double objectiveChange = 0.0;
595 double * solution = model->solutionRegion();
596 for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
597 int iSet = backward_[iColumn];
598 if (toIndex_[iSet] < 0) {
599 toIndex_[iSet] = 0;
600 fromIndex_[numberToDo++] = iSet;
601 }
602 if (model->getStatus(iColumn) == ClpSimplex::basic || iColumn == keyVariable_[iSet]) {
603 lookup[iColumn] = firstAvailable_;
604 if (iColumn != keyVariable_[iSet]) {
605 int iPivot = backToPivotRow_[iColumn];
606 backToPivotRow_[firstAvailable_] = iPivot;
607 pivotVariable[iPivot] = firstAvailable_;
608 }
609 firstAvailable_++;
610 } else {
611 int jColumn = id_[iColumn-firstDynamic_];
612 setDynamicStatus(jColumn, atLowerBound);
613 if (lowerColumn_ || upperColumn_) {
614 if (model->getStatus(iColumn) == ClpSimplex::atUpperBound)
615 setDynamicStatus(jColumn, atUpperBound);
616 // treat solution as if exactly at a bound
617 double value = solution[iColumn];
618 if (fabs(value - lowerColumn[iColumn]) < fabs(value - upperColumn[iColumn]))
619 value = lowerColumn[iColumn];
620 else
621 value = upperColumn[iColumn];
622 objectiveChange += cost[iColumn] * value;
623 // redo lower and upper on sets
624 double shift = value;
625 if (lowerSet_[iSet] > -1.0e20)
626 lower_[iSet] = lowerSet_[iSet] - shift;
627 if (upperSet_[iSet] < 1.0e20)
628 upper_[iSet] = upperSet_[iSet] - shift;
629 }
630 lookup[iColumn] = -1;
631 }
632 }
633 model->setObjectiveOffset(model->objectiveOffset() + objectiveChange);
634 firstAvailable_ = firstDynamic_;
635 for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) {
636 if (lookup[iColumn] >= 0) {
637 // move
638 int jColumn = id_[iColumn-firstDynamic_];
639 id_[firstAvailable_-firstDynamic_] = jColumn;
640 int numberThis = startColumn_[jColumn+1] - startColumn_[jColumn];
641 length[firstAvailable_] = numberThis;
642 cost[firstAvailable_] = cost[iColumn];
643 lowerColumn[firstAvailable_] = lowerColumn[iColumn];
644 upperColumn[firstAvailable_] = upperColumn[iColumn];
645 double originalLower = lowerColumn_ ? lowerColumn_[jColumn] : 0.0;
646 double originalUpper = upperColumn_ ? upperColumn_[jColumn] : COIN_DBL_MAX;
647 if (originalUpper > 1.0e30)
648 originalUpper = COIN_DBL_MAX;
649 model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn],
650 originalLower, originalUpper,
651 cost_[jColumn]);
652 CoinBigIndex base = startColumn_[jColumn];
653 for (int j = 0; j < numberThis; j++) {
654 row[numberElements] = row_[base+j];
655 element[numberElements++] = element_[base+j];
656 }
657 model->setStatus(firstAvailable_, model->getStatus(iColumn));
658 backward_[firstAvailable_] = backward_[iColumn];
659 solution[firstAvailable_] = solution[iColumn];
660 firstAvailable_++;
661 startColumn[firstAvailable_] = numberElements;
662 }
663 }
664 // clean up next_
665 int * temp = new int [firstAvailable_];
666 for (int jSet = 0; jSet < numberToDo; jSet++) {
667 int iSet = fromIndex_[jSet];
668 toIndex_[iSet] = -1;
669 int last = keyVariable_[iSet];
670 int j = next_[last];
671 bool setTemp = true;
672 if (last < lastDynamic_) {
673 last = lookup[last];
674 assert (last >= 0);
675 keyVariable_[iSet] = last;
676 } else if (j >= 0) {
677 int newJ = lookup[j];
678 assert (newJ >= 0);
679 j = next_[j];
680 next_[last] = newJ;
681 last = newJ;
682 } else {
683 next_[last] = -(iSet + numberColumns + 1);
684 setTemp = false;
685 }
686 while (j >= 0) {
687 int newJ = lookup[j];
688 assert (newJ >= 0);
689 temp[last] = newJ;
690 last = newJ;
691 j = next_[j];
692 }
693 if (setTemp)
694 temp[last] = -(keyVariable_[iSet] + 1);
695 if (lowerSet_) {
696 // we only need to get lower_ and upper_ correct
697 double shift = 0.0;
698 for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
699 if (getDynamicStatus(j) == atUpperBound)
700 shift += upperColumn_[j];
701 else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
702 shift += lowerColumn_[j];
703 if (lowerSet_[iSet] > -1.0e20)
704 lower_[iSet] = lowerSet_[iSet] - shift;
705 if (upperSet_[iSet] < 1.0e20)
706 upper_[iSet] = upperSet_[iSet] - shift;
707 }
708 }
709 // move to next_
710 CoinMemcpyN(temp + firstDynamic_, (firstAvailable_ - firstDynamic_), next_ + firstDynamic_);
711 // if odd iterations may be one out so adjust currentNumber
712 currentNumber = CoinMin(currentNumber + 1, lastDynamic_);
713 // zero solution
714 CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
715 // zero cost
716 CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
717 // zero lengths
718 CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
719 for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
720 model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
721 model->setStatus(iColumn, ClpSimplex::atLowerBound);
722 backward_[iColumn] = -1;
723 }
724 delete [] lookup;
725 delete [] temp;
726 // make sure fromIndex clean
727 fromIndex_[0] = -1;
728 //#define CLP_DEBUG
729#ifdef CLP_DEBUG
730 // debug
731 {
732 int i;
733 int numberRows = model->numberRows();
734 char * xxxx = new char[numberColumns];
735 memset(xxxx, 0, numberColumns);
736 for (i = 0; i < numberRows; i++) {
737 int iPivot = pivotVariable[i];
738 assert (model->getStatus(iPivot) == ClpSimplex::basic);
739 if (iPivot < numberColumns && backward_[iPivot] >= 0)
740 xxxx[iPivot] = 1;
741 }
742 for (i = 0; i < numberSets_; i++) {
743 int key = keyVariable_[i];
744 int iColumn = next_[key];
745 int k = 0;
746 while(iColumn >= 0) {
747 k++;
748 assert (k < 100);
749 assert (backward_[iColumn] == i);
750 iColumn = next_[iColumn];
751 }
752 int stop = -(key + 1);
753 while (iColumn != stop) {
754 assert (iColumn < 0);
755 iColumn = -iColumn - 1;
756 k++;
757 assert (k < 100);
758 assert (backward_[iColumn] == i);
759 iColumn = next_[iColumn];
760 }
761 iColumn = next_[key];
762 while (iColumn >= 0) {
763 assert (xxxx[iColumn]);
764 xxxx[iColumn] = 0;
765 iColumn = next_[iColumn];
766 }
767 }
768 for (i = 0; i < numberColumns; i++) {
769 if (i < numberColumns && backward_[i] >= 0) {
770 assert (!xxxx[i] || i == keyVariable_[backward_[i]]);
771 }
772 }
773 delete [] xxxx;
774 }
775 {
776 for (int i = 0; i < numberSets_; i++)
777 assert(toIndex_[i] == -1);
778 }
779#endif
780 savedFirstAvailable_ = firstAvailable_;
781 }
782 break;
783 // flag a variable
784 case 1: {
785 // id will be sitting at firstAvailable
786 int sequence = id_[firstAvailable_-firstDynamic_];
787 assert (!flagged(sequence));
788 setFlagged(sequence);
789 model->clearFlagged(firstAvailable_);
790 }
791 break;
792 // unflag all variables
793 case 2: {
794 for (int i = 0; i < numberGubColumns_; i++) {
795 if (flagged(i)) {
796 unsetFlagged(i);
797 returnNumber++;
798 }
799 }
800 }
801 break;
802 // just reset costs and bounds (primal)
803 case 3: {
804 double * cost = model->costRegion();
805 double * solution = model->solutionRegion();
806 double * lowerColumn = model->columnLower();
807 double * upperColumn = model->columnUpper();
808 for (int i = firstDynamic_; i < firstAvailable_; i++) {
809 int jColumn = id_[i-firstDynamic_];
810 cost[i] = cost_[jColumn];
811 if (!lowerColumn_ && !upperColumn_) {
812 lowerColumn[i] = 0.0;
813 upperColumn[i] = COIN_DBL_MAX;
814 } else {
815 if (lowerColumn_)
816 lowerColumn[i] = lowerColumn_[jColumn];
817 else
818 lowerColumn[i] = 0.0;
819 if (upperColumn_)
820 upperColumn[i] = upperColumn_[jColumn];
821 else
822 upperColumn[i] = COIN_DBL_MAX;
823 }
824 if (model->nonLinearCost())
825 model->nonLinearCost()->setOne(i, solution[i],
826 lowerColumn[i],
827 upperColumn[i], cost_[jColumn]);
828 }
829 if (!model->numberIterations() && rhsOffset_) {
830 lastRefresh_ = - refreshFrequency_; // force refresh
831 }
832 }
833 break;
834 // and get statistics for column generation
835 case 4: {
836 // In theory we should subtract out ones we have done but ....
837 // If key slack then dual 0.0
838 // If not then slack could be dual infeasible
839 // dj for key is zero so that defines dual on set
840 int i;
841 int numberColumns = model->numberColumns();
842 double * dual = model->dualRowSolution();
843 double infeasibilityCost = model->infeasibilityCost();
844 double dualTolerance = model->dualTolerance();
845 double relaxedTolerance = dualTolerance;
846 // we can't really trust infeasibilities if there is dual error
847 double error = CoinMin(1.0e-2, model->largestDualError());
848 // allow tolerance at least slightly bigger than standard
849 relaxedTolerance = relaxedTolerance + error;
850 // but we will be using difference
851 relaxedTolerance -= dualTolerance;
852 double objectiveOffset = 0.0;
853 for (i = 0; i < numberSets_; i++) {
854 int kColumn = keyVariable_[i];
855 double value = 0.0;
856 if (kColumn < numberColumns) {
857 kColumn = id_[kColumn-firstDynamic_];
858 // dj without set
859 value = cost_[kColumn];
860 for (CoinBigIndex j = startColumn_[kColumn];
861 j < startColumn_[kColumn+1]; j++) {
862 int iRow = row_[j];
863 value -= dual[iRow] * element_[j];
864 }
865 double infeasibility = 0.0;
866 if (getStatus(i) == ClpSimplex::atLowerBound) {
867 if (-value > dualTolerance)
868 infeasibility = -value - dualTolerance;
869 } else if (getStatus(i) == ClpSimplex::atUpperBound) {
870 if (value > dualTolerance)
871 infeasibility = value - dualTolerance;
872 }
873 if (infeasibility > 0.0) {
874 sumDualInfeasibilities_ += infeasibility;
875 if (infeasibility > relaxedTolerance)
876 sumOfRelaxedDualInfeasibilities_ += infeasibility;
877 numberDualInfeasibilities_ ++;
878 }
879 } else {
880 // slack key - may not be feasible
881 assert (getStatus(i) == ClpSimplex::basic);
882 // negative as -1.0 for slack
883 value = -weight(i) * infeasibilityCost;
884 }
885 // Now subtract out from all
886 for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i+1]; k++) {
887 if (getDynamicStatus(k) != inSmall) {
888 double djValue = cost_[k] - value;
889 for (CoinBigIndex j = startColumn_[k];
890 j < startColumn_[k+1]; j++) {
891 int iRow = row_[j];
892 djValue -= dual[iRow] * element_[j];
893 }
894 double infeasibility = 0.0;
895 double shift = 0.0;
896 if (getDynamicStatus(k) == atLowerBound) {
897 if (lowerColumn_)
898 shift = lowerColumn_[k];
899 if (djValue < -dualTolerance)
900 infeasibility = -djValue - dualTolerance;
901 } else {
902 // at upper bound
903 shift = upperColumn_[k];
904 if (djValue > dualTolerance)
905 infeasibility = djValue - dualTolerance;
906 }
907 objectiveOffset += shift * cost_[k];
908 if (infeasibility > 0.0) {
909 sumDualInfeasibilities_ += infeasibility;
910 if (infeasibility > relaxedTolerance)
911 sumOfRelaxedDualInfeasibilities_ += infeasibility;
912 numberDualInfeasibilities_ ++;
913 }
914 }
915 }
916 }
917 model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
918 }
919 break;
920 // see if time to re-factorize
921 case 5: {
922 if (firstAvailable_ > numberSets_ + model->numberRows() + model->factorizationFrequency())
923 returnNumber = 4;
924 }
925 break;
926 // return 1 if there may be changing bounds on variable (column generation)
927 case 6: {
928 returnNumber = (lowerColumn_ != NULL || upperColumn_ != NULL) ? 1 : 0;
929#if 0
930 if (!returnNumber) {
931 // may be gub slacks
932 for (int i = 0; i < numberSets_; i++) {
933 if (upper_[i] > lower_[i]) {
934 returnNumber = 1;
935 break;
936 }
937 }
938 }
939#endif
940 }
941 break;
942 // restore firstAvailable_
943 case 7: {
944 int iColumn;
945 int * length = matrix_->getMutableVectorLengths();
946 double * cost = model->costRegion();
947 double * solution = model->solutionRegion();
948 int currentNumber = firstAvailable_;
949 firstAvailable_ = savedFirstAvailable_;
950 // zero solution
951 CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_);
952 // zero cost
953 CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_);
954 // zero lengths
955 CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_);
956 for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) {
957 model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0);
958 model->setStatus(iColumn, ClpSimplex::atLowerBound);
959 backward_[iColumn] = -1;
960 }
961 }
962 break;
963 // make sure set is clean
964 case 8: {
965 int sequenceIn = model->sequenceIn();
966 if (sequenceIn < model->numberColumns()) {
967 int iSet = backward_[sequenceIn];
968 if (iSet >= 0 && lowerSet_) {
969 // we only need to get lower_ and upper_ correct
970 double shift = 0.0;
971 for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
972 if (getDynamicStatus(j) == atUpperBound)
973 shift += upperColumn_[j];
974 else if (getDynamicStatus(j) == atLowerBound && lowerColumn_)
975 shift += lowerColumn_[j];
976 if (lowerSet_[iSet] > -1.0e20)
977 lower_[iSet] = lowerSet_[iSet] - shift;
978 if (upperSet_[iSet] < 1.0e20)
979 upper_[iSet] = upperSet_[iSet] - shift;
980 }
981 if (sequenceIn == firstAvailable_) {
982 // not really in small problem
983 int iBig = id_[sequenceIn-firstDynamic_];
984 if (model->getStatus(sequenceIn) == ClpSimplex::atLowerBound)
985 setDynamicStatus(iBig, atLowerBound);
986 else
987 setDynamicStatus(iBig, atUpperBound);
988 }
989 }
990 }
991 break;
992 // adjust lower,upper
993 case 9: {
994 int sequenceIn = model->sequenceIn();
995 if (sequenceIn >= firstDynamic_ && sequenceIn < lastDynamic_ && lowerSet_) {
996 int iSet = backward_[sequenceIn];
997 assert (iSet >= 0);
998 int inBig = id_[sequenceIn-firstDynamic_];
999 const double * solution = model->solutionRegion();
1000 setDynamicStatus(inBig, inSmall);
1001 if (lowerSet_[iSet] > -1.0e20)
1002 lower_[iSet] += solution[sequenceIn];
1003 if (upperSet_[iSet] < 1.0e20)
1004 upper_[iSet] += solution[sequenceIn];
1005 model->setObjectiveOffset(model->objectiveOffset() -
1006 solution[sequenceIn]*cost_[inBig]);
1007 }
1008 }
1009 }
1010 return returnNumber;
1011}
1012// Add a new variable to a set
1013void
1014ClpGubDynamicMatrix::insertNonBasic(int sequence, int iSet)
1015{
1016 int last = keyVariable_[iSet];
1017 int j = next_[last];
1018 while (j >= 0) {
1019 last = j;
1020 j = next_[j];
1021 }
1022 next_[last] = -(sequence + 1);
1023 next_[sequence] = j;
1024}
1025// Sets up an effective RHS and does gub crash if needed
1026void
1027ClpGubDynamicMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
1028{
1029 // Do basis - cheapest or slack if feasible (unless cheapest set)
1030 int longestSet = 0;
1031 int iSet;
1032 for (iSet = 0; iSet < numberSets_; iSet++)
1033 longestSet = CoinMax(longestSet, fullStart_[iSet+1] - fullStart_[iSet]);
1034
1035 double * upper = new double[longestSet+1];
1036 double * cost = new double[longestSet+1];
1037 double * lower = new double[longestSet+1];
1038 double * solution = new double[longestSet+1];
1039 assert (!next_);
1040 delete [] next_;
1041 int numberColumns = model->numberColumns();
1042 next_ = new int[numberColumns+numberSets_+CoinMax(2*longestSet, lastDynamic_-firstDynamic_)];
1043 char * mark = new char[numberColumns];
1044 memset(mark, 0, numberColumns);
1045 for (int iColumn = 0; iColumn < numberColumns; iColumn++)
1046 next_[iColumn] = COIN_INT_MAX;
1047 int i;
1048 int * keys = new int[numberSets_];
1049 int * back = new int[numberGubColumns_];
1050 CoinFillN(back, numberGubColumns_, -1);
1051 for (i = 0; i < numberSets_; i++)
1052 keys[i] = COIN_INT_MAX;
1053 delete [] dynamicStatus_;
1054 dynamicStatus_ = new unsigned char [numberGubColumns_];
1055 memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
1056 for (i = 0; i < numberGubColumns_; i++)
1057 setDynamicStatus(i, atLowerBound);
1058 // set up chains
1059 for (i = firstDynamic_; i < lastDynamic_; i++) {
1060 if (id_[i-firstDynamic_] >= 0) {
1061 if (model->getStatus(i) == ClpSimplex::basic)
1062 mark[i] = 1;
1063 int iSet = backward_[i];
1064 assert (iSet >= 0);
1065 int iNext = keys[iSet];
1066 next_[i] = iNext;
1067 keys[iSet] = i;
1068 back[id_[i-firstDynamic_]] = i;
1069 } else {
1070 model->setStatus(i, ClpSimplex::atLowerBound);
1071 backward_[i] = -1;
1072 }
1073 }
1074 double * columnSolution = model->solutionRegion();
1075 int numberRows = getNumRows();
1076 toIndex_ = new int[numberSets_];
1077 for (iSet = 0; iSet < numberSets_; iSet++)
1078 toIndex_[iSet] = -1;
1079 fromIndex_ = new int [numberRows+numberSets_];
1080 double tolerance = model->primalTolerance();
1081 double * element = matrix_->getMutableElements();
1082 int * row = matrix_->getMutableIndices();
1083 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1084 int * length = matrix_->getMutableVectorLengths();
1085 double objectiveOffset = 0.0;
1086 for (iSet = 0; iSet < numberSets_; iSet++) {
1087 int j;
1088 int numberBasic = 0;
1089 int iBasic = -1;
1090 int iStart = fullStart_[iSet];
1091 int iEnd = fullStart_[iSet+1];
1092 // find one with smallest length
1093 int smallest = numberRows + 1;
1094 double value = 0.0;
1095 j = keys[iSet];
1096 while (j != COIN_INT_MAX) {
1097 if (model->getStatus(j) == ClpSimplex::basic) {
1098 if (length[j] < smallest) {
1099 smallest = length[j];
1100 iBasic = j;
1101 }
1102 numberBasic++;
1103 }
1104 value += columnSolution[j];
1105 j = next_[j];
1106 }
1107 bool done = false;
1108 if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
1109 if (getStatus(iSet) == ClpSimplex::basic)
1110 iBasic = iSet + numberColumns; // slack key - use
1111 done = true;
1112 } else if (numberBasic == 1) {
1113 // see if can be key
1114 double thisSolution = columnSolution[iBasic];
1115 if (thisSolution < 0.0) {
1116 value -= thisSolution;
1117 thisSolution = 0.0;
1118 columnSolution[iBasic] = thisSolution;
1119 }
1120 // try setting slack to a bound
1121 assert (upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
1122 double cost1 = COIN_DBL_MAX;
1123 int whichBound = -1;
1124 if (upper_[iSet] < 1.0e20) {
1125 // try slack at ub
1126 double newBasic = thisSolution + upper_[iSet] - value;
1127 if (newBasic >= -tolerance) {
1128 // can go
1129 whichBound = 1;
1130 cost1 = newBasic * cost_[iBasic];
1131 // But if exact then may be good solution
1132 if (fabs(upper_[iSet] - value) < tolerance)
1133 cost1 = -COIN_DBL_MAX;
1134 }
1135 }
1136 if (lower_[iSet] > -1.0e20) {
1137 // try slack at lb
1138 double newBasic = thisSolution + lower_[iSet] - value;
1139 if (newBasic >= -tolerance) {
1140 // can go but is it cheaper
1141 double cost2 = newBasic * cost_[iBasic];
1142 // But if exact then may be good solution
1143 if (fabs(lower_[iSet] - value) < tolerance)
1144 cost2 = -COIN_DBL_MAX;
1145 if (cost2 < cost1)
1146 whichBound = 0;
1147 }
1148 }
1149 if (whichBound != -1) {
1150 // key
1151 done = true;
1152 if (whichBound) {
1153 // slack to upper
1154 columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
1155 setStatus(iSet, ClpSimplex::atUpperBound);
1156 } else {
1157 // slack to lower
1158 columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
1159 setStatus(iSet, ClpSimplex::atLowerBound);
1160 }
1161 }
1162 }
1163 if (!done) {
1164 if (!cheapest) {
1165 // see if slack can be key
1166 if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
1167 done = true;
1168 setStatus(iSet, ClpSimplex::basic);
1169 iBasic = iSet + numberColumns;
1170 }
1171 }
1172 if (!done) {
1173 // set non basic if there was one
1174 if (iBasic >= 0)
1175 model->setStatus(iBasic, ClpSimplex::atLowerBound);
1176 // find cheapest
1177 int numberInSet = iEnd - iStart;
1178 if (!lowerColumn_) {
1179 CoinZeroN(lower, numberInSet);
1180 } else {
1181 for (int j = 0; j < numberInSet; j++)
1182 lower[j] = lowerColumn_[j+iStart];
1183 }
1184 if (!upperColumn_) {
1185 CoinFillN(upper, numberInSet, COIN_DBL_MAX);
1186 } else {
1187 for (int j = 0; j < numberInSet; j++)
1188 upper[j] = upperColumn_[j+iStart];
1189 }
1190 CoinFillN(solution, numberInSet, 0.0);
1191 // and slack
1192 iBasic = numberInSet;
1193 solution[iBasic] = -value;
1194 lower[iBasic] = -upper_[iSet];
1195 upper[iBasic] = -lower_[iSet];
1196 int kphase;
1197 if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
1198 // feasible
1199 kphase = 1;
1200 cost[iBasic] = 0.0;
1201 for (int j = 0; j < numberInSet; j++)
1202 cost[j] = cost_[j+iStart];
1203 } else {
1204 // infeasible
1205 kphase = 0;
1206 // remember bounds are flipped so opposite to natural
1207 if (value < lower_[iSet] - tolerance)
1208 cost[iBasic] = 1.0;
1209 else
1210 cost[iBasic] = -1.0;
1211 CoinZeroN(cost, numberInSet);
1212 }
1213 double dualTolerance = model->dualTolerance();
1214 for (int iphase = kphase; iphase < 2; iphase++) {
1215 if (iphase) {
1216 cost[numberInSet] = 0.0;
1217 for (int j = 0; j < numberInSet; j++)
1218 cost[j] = cost_[j+iStart];
1219 }
1220 // now do one row lp
1221 bool improve = true;
1222 while (improve) {
1223 improve = false;
1224 double dual = cost[iBasic];
1225 int chosen = -1;
1226 double best = dualTolerance;
1227 int way = 0;
1228 for (int i = 0; i <= numberInSet; i++) {
1229 double dj = cost[i] - dual;
1230 double improvement = 0.0;
1231 if (iphase || i < numberInSet)
1232 assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
1233 if (dj > dualTolerance)
1234 improvement = dj * (solution[i] - lower[i]);
1235 else if (dj < -dualTolerance)
1236 improvement = dj * (solution[i] - upper[i]);
1237 if (improvement > best) {
1238 best = improvement;
1239 chosen = i;
1240 if (dj < 0.0) {
1241 way = 1;
1242 } else {
1243 way = -1;
1244 }
1245 }
1246 }
1247 if (chosen >= 0) {
1248 improve = true;
1249 // now see how far
1250 if (way > 0) {
1251 // incoming increasing so basic decreasing
1252 // if phase 0 then go to nearest bound
1253 double distance = upper[chosen] - solution[chosen];
1254 double basicDistance;
1255 if (!iphase) {
1256 assert (iBasic == numberInSet);
1257 assert (solution[iBasic] > upper[iBasic]);
1258 basicDistance = solution[iBasic] - upper[iBasic];
1259 } else {
1260 basicDistance = solution[iBasic] - lower[iBasic];
1261 }
1262 // need extra coding for unbounded
1263 assert (CoinMin(distance, basicDistance) < 1.0e20);
1264 if (distance > basicDistance) {
1265 // incoming becomes basic
1266 solution[chosen] += basicDistance;
1267 if (!iphase)
1268 solution[iBasic] = upper[iBasic];
1269 else
1270 solution[iBasic] = lower[iBasic];
1271 iBasic = chosen;
1272 } else {
1273 // flip
1274 solution[chosen] = upper[chosen];
1275 solution[iBasic] -= distance;
1276 }
1277 } else {
1278 // incoming decreasing so basic increasing
1279 // if phase 0 then go to nearest bound
1280 double distance = solution[chosen] - lower[chosen];
1281 double basicDistance;
1282 if (!iphase) {
1283 assert (iBasic == numberInSet);
1284 assert (solution[iBasic] < lower[iBasic]);
1285 basicDistance = lower[iBasic] - solution[iBasic];
1286 } else {
1287 basicDistance = upper[iBasic] - solution[iBasic];
1288 }
1289 // need extra coding for unbounded - for now just exit
1290 if (CoinMin(distance, basicDistance) > 1.0e20) {
1291 printf("unbounded on set %d\n", iSet);
1292 iphase = 1;
1293 iBasic = numberInSet;
1294 break;
1295 }
1296 if (distance > basicDistance) {
1297 // incoming becomes basic
1298 solution[chosen] -= basicDistance;
1299 if (!iphase)
1300 solution[iBasic] = lower[iBasic];
1301 else
1302 solution[iBasic] = upper[iBasic];
1303 iBasic = chosen;
1304 } else {
1305 // flip
1306 solution[chosen] = lower[chosen];
1307 solution[iBasic] += distance;
1308 }
1309 }
1310 if (!iphase) {
1311 if(iBasic < numberInSet)
1312 break; // feasible
1313 else if (solution[iBasic] >= lower[iBasic] &&
1314 solution[iBasic] <= upper[iBasic])
1315 break; // feasible (on flip)
1316 }
1317 }
1318 }
1319 }
1320 // do solution i.e. bounds
1321 if (lowerColumn_ || upperColumn_) {
1322 for (int j = 0; j < numberInSet; j++) {
1323 if (j != iBasic) {
1324 objectiveOffset += solution[j] * cost[j];
1325 if (lowerColumn_ && upperColumn_) {
1326 if (fabs(solution[j] - lowerColumn_[j+iStart]) >
1327 fabs(solution[j] - upperColumn_[j+iStart]))
1328 setDynamicStatus(j + iStart, atUpperBound);
1329 } else if (upperColumn_ && solution[j] > 0.0) {
1330 setDynamicStatus(j + iStart, atUpperBound);
1331 } else {
1332 setDynamicStatus(j + iStart, atLowerBound);
1333 }
1334 }
1335 }
1336 }
1337 // convert iBasic back and do bounds
1338 if (iBasic == numberInSet) {
1339 // slack basic
1340 setStatus(iSet, ClpSimplex::basic);
1341 iBasic = iSet + numberColumns;
1342 } else {
1343 iBasic += fullStart_[iSet];
1344 if (back[iBasic] >= 0) {
1345 // exists
1346 iBasic = back[iBasic];
1347 } else {
1348 // create
1349 CoinBigIndex numberElements = startColumn[firstAvailable_];
1350 int numberThis = startColumn_[iBasic+1] - startColumn_[iBasic];
1351 if (numberElements + numberThis > numberElements_) {
1352 // need to redo
1353 numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis);
1354 matrix_->reserve(numberColumns, numberElements_);
1355 element = matrix_->getMutableElements();
1356 row = matrix_->getMutableIndices();
1357 // these probably okay but be safe
1358 startColumn = matrix_->getMutableVectorStarts();
1359 length = matrix_->getMutableVectorLengths();
1360 }
1361 length[firstAvailable_] = numberThis;
1362 model->costRegion()[firstAvailable_] = cost_[iBasic];
1363 if (lowerColumn_)
1364 model->lowerRegion()[firstAvailable_] = lowerColumn_[iBasic];
1365 else
1366 model->lowerRegion()[firstAvailable_] = 0.0;
1367 if (upperColumn_)
1368 model->upperRegion()[firstAvailable_] = upperColumn_[iBasic];
1369 else
1370 model->upperRegion()[firstAvailable_] = COIN_DBL_MAX;
1371 columnSolution[firstAvailable_] = solution[iBasic-fullStart_[iSet]];
1372 CoinBigIndex base = startColumn_[iBasic];
1373 for (int j = 0; j < numberThis; j++) {
1374 row[numberElements] = row_[base+j];
1375 element[numberElements++] = element_[base+j];
1376 }
1377 // already set startColumn[firstAvailable_]=numberElements;
1378 id_[firstAvailable_-firstDynamic_] = iBasic;
1379 setDynamicStatus(iBasic, inSmall);
1380 backward_[firstAvailable_] = iSet;
1381 iBasic = firstAvailable_;
1382 firstAvailable_++;
1383 startColumn[firstAvailable_] = numberElements;
1384 }
1385 model->setStatus(iBasic, ClpSimplex::basic);
1386 // remember bounds flipped
1387 if (upper[numberInSet] == lower[numberInSet])
1388 setStatus(iSet, ClpSimplex::isFixed);
1389 else if (solution[numberInSet] == upper[numberInSet])
1390 setStatus(iSet, ClpSimplex::atLowerBound);
1391 else if (solution[numberInSet] == lower[numberInSet])
1392 setStatus(iSet, ClpSimplex::atUpperBound);
1393 else
1394 abort();
1395 }
1396 for (j = iStart; j < iEnd; j++) {
1397 int iBack = back[j];
1398 if (iBack >= 0) {
1399 if (model->getStatus(iBack) != ClpSimplex::basic) {
1400 int inSet = j - iStart;
1401 columnSolution[iBack] = solution[inSet];
1402 if (upper[inSet] == lower[inSet])
1403 model->setStatus(iBack, ClpSimplex::isFixed);
1404 else if (solution[inSet] == upper[inSet])
1405 model->setStatus(iBack, ClpSimplex::atUpperBound);
1406 else if (solution[inSet] == lower[inSet])
1407 model->setStatus(iBack, ClpSimplex::atLowerBound);
1408 }
1409 }
1410 }
1411 }
1412 }
1413 keyVariable_[iSet] = iBasic;
1414 }
1415 model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
1416 delete [] lower;
1417 delete [] solution;
1418 delete [] upper;
1419 delete [] cost;
1420 // make sure matrix is in good shape
1421 matrix_->orderMatrix();
1422 // create effective rhs
1423 delete [] rhsOffset_;
1424 rhsOffset_ = new double[numberRows];
1425 // and redo chains
1426 memset(mark, 0, numberColumns);
1427 for (int iColumnX = 0; iColumnX < firstAvailable_; iColumnX++)
1428 next_[iColumnX] = COIN_INT_MAX;
1429 for (i = 0; i < numberSets_; i++) {
1430 keys[i] = COIN_INT_MAX;
1431 int iKey = keyVariable_[i];
1432 if (iKey < numberColumns)
1433 model->setStatus(iKey, ClpSimplex::basic);
1434 }
1435 // set up chains
1436 for (i = 0; i < firstAvailable_; i++) {
1437 if (model->getStatus(i) == ClpSimplex::basic)
1438 mark[i] = 1;
1439 int iSet = backward_[i];
1440 if (iSet >= 0) {
1441 int iNext = keys[iSet];
1442 next_[i] = iNext;
1443 keys[iSet] = i;
1444 }
1445 }
1446 for (i = 0; i < numberSets_; i++) {
1447 if (keys[i] != COIN_INT_MAX) {
1448 // something in set
1449 int j;
1450 if (getStatus(i) != ClpSimplex::basic) {
1451 // make sure fixed if it is
1452 if (upper_[i] == lower_[i])
1453 setStatus(i, ClpSimplex::isFixed);
1454 // slack not key - choose one with smallest length
1455 int smallest = numberRows + 1;
1456 int key = -1;
1457 j = keys[i];
1458 while (1) {
1459 if (mark[j] && length[j] < smallest) {
1460 key = j;
1461 smallest = length[j];
1462 }
1463 if (next_[j] != COIN_INT_MAX) {
1464 j = next_[j];
1465 } else {
1466 // correct end
1467 next_[j] = -(keys[i] + 1);
1468 break;
1469 }
1470 }
1471 if (key >= 0) {
1472 keyVariable_[i] = key;
1473 } else {
1474 // nothing basic - make slack key
1475 //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
1476 // fudge to avoid const problem
1477 status_[i] = 1;
1478 }
1479 } else {
1480 // slack key
1481 keyVariable_[i] = numberColumns + i;
1482 int j;
1483 double sol = 0.0;
1484 j = keys[i];
1485 while (1) {
1486 sol += columnSolution[j];
1487 if (next_[j] != COIN_INT_MAX) {
1488 j = next_[j];
1489 } else {
1490 // correct end
1491 next_[j] = -(keys[i] + 1);
1492 break;
1493 }
1494 }
1495 if (sol > upper_[i] + tolerance) {
1496 setAbove(i);
1497 } else if (sol < lower_[i] - tolerance) {
1498 setBelow(i);
1499 } else {
1500 setFeasible(i);
1501 }
1502 }
1503 // Create next_
1504 int key = keyVariable_[i];
1505 redoSet(model, key, keys[i], i);
1506 } else {
1507 // nothing in set!
1508 next_[i+numberColumns] = -(i + numberColumns + 1);
1509 keyVariable_[i] = numberColumns + i;
1510 double sol = 0.0;
1511 if (sol > upper_[i] + tolerance) {
1512 setAbove(i);
1513 } else if (sol < lower_[i] - tolerance) {
1514 setBelow(i);
1515 } else {
1516 setFeasible(i);
1517 }
1518 }
1519 }
1520 delete [] keys;
1521 delete [] mark;
1522 delete [] back;
1523 rhsOffset(model, true);
1524}
1525/* Returns effective RHS if it is being used. This is used for long problems
1526 or big gub or anywhere where going through full columns is
1527 expensive. This may re-compute */
1528double *
1529ClpGubDynamicMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh,
1530 bool
1531#ifdef CLP_DEBUG
1532 check
1533#endif
1534 )
1535{
1536 //forceRefresh=true;
1537 //check=false;
1538#ifdef CLP_DEBUG
1539 double * saveE = NULL;
1540 if (rhsOffset_ && check) {
1541 int numberRows = model->numberRows();
1542 saveE = new double[numberRows];
1543 }
1544#endif
1545 if (rhsOffset_) {
1546#ifdef CLP_DEBUG
1547 if (check) {
1548 // no need - but check anyway
1549 int numberRows = model->numberRows();
1550 double * rhs = new double[numberRows];
1551 int numberColumns = model->numberColumns();
1552 int iRow;
1553 CoinZeroN(rhs, numberRows);
1554 // do ones at bounds before gub
1555 const double * smallSolution = model->solutionRegion();
1556 const double * element = matrix_->getElements();
1557 const int * row = matrix_->getIndices();
1558 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1559 const int * length = matrix_->getVectorLengths();
1560 int iColumn;
1561 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
1562 if (model->getStatus(iColumn) != ClpSimplex::basic) {
1563 double value = smallSolution[iColumn];
1564 for (CoinBigIndex j = startColumn[iColumn];
1565 j < startColumn[iColumn] + length[iColumn]; j++) {
1566 int jRow = row[j];
1567 rhs[jRow] -= value * element[j];
1568 }
1569 }
1570 }
1571 if (lowerColumn_ || upperColumn_) {
1572 double * solution = new double [numberGubColumns_];
1573 for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1574 double value = 0.0;
1575 if(getDynamicStatus(iColumn) == atUpperBound)
1576 value = upperColumn_[iColumn];
1577 else if (lowerColumn_)
1578 value = lowerColumn_[iColumn];
1579 solution[iColumn] = value;
1580 }
1581 // ones at bounds in small and gub
1582 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
1583 int jFull = id_[iColumn-firstDynamic_];
1584 solution[jFull] = smallSolution[iColumn];
1585 }
1586 // zero all basic in small model
1587 int * pivotVariable = model->pivotVariable();
1588 for (iRow = 0; iRow < numberRows; iRow++) {
1589 int iColumn = pivotVariable[iRow];
1590 if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
1591 int iSequence = id_[iColumn-firstDynamic_];
1592 solution[iSequence] = 0.0;
1593 }
1594 }
1595 // and now compute value to use for key
1596 ClpSimplex::Status iStatus;
1597 for (int iSet = 0; iSet < numberSets_; iSet++) {
1598 iColumn = keyVariable_[iSet];
1599 if (iColumn < numberColumns) {
1600 int iSequence = id_[iColumn-firstDynamic_];
1601 solution[iSequence] = 0.0;
1602 double b = 0.0;
1603 // key is structural - where is slack
1604 iStatus = getStatus(iSet);
1605 assert (iStatus != ClpSimplex::basic);
1606 if (iStatus == ClpSimplex::atLowerBound)
1607 b = lowerSet_[iSet];
1608 else
1609 b = upperSet_[iSet];
1610 // subtract out others at bounds
1611 for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
1612 b -= solution[j];
1613 solution[iSequence] = b;
1614 }
1615 }
1616 for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1617 double value = solution[iColumn];
1618 if (value) {
1619 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
1620 int iRow = row_[j];
1621 rhs[iRow] -= element_[j] * value;
1622 }
1623 }
1624 }
1625 // now do lower and upper bounds on sets
1626 for (int iSet = 0; iSet < numberSets_; iSet++) {
1627 iColumn = keyVariable_[iSet];
1628 double shift = 0.0;
1629 for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++) {
1630 if (getDynamicStatus(j) != inSmall && j != iColumn) {
1631 if (getDynamicStatus(j) == atLowerBound) {
1632 if (lowerColumn_)
1633 shift += lowerColumn_[j];
1634 } else {
1635 shift += upperColumn_[j];
1636 }
1637 }
1638 }
1639 if (lowerSet_[iSet] > -1.0e20)
1640 assert(fabs(lower_[iSet] - (lowerSet_[iSet] - shift)) < 1.0e-3);
1641 if (upperSet_[iSet] < 1.0e20)
1642 assert(fabs(upper_[iSet] - ( upperSet_[iSet] - shift)) < 1.0e-3);
1643 }
1644 delete [] solution;
1645 } else {
1646 // no bounds
1647 ClpSimplex::Status iStatus;
1648 for (int iSet = 0; iSet < numberSets_; iSet++) {
1649 int iColumn = keyVariable_[iSet];
1650 if (iColumn < numberColumns) {
1651 int iSequence = id_[iColumn-firstDynamic_];
1652 double b = 0.0;
1653 // key is structural - where is slack
1654 iStatus = getStatus(iSet);
1655 assert (iStatus != ClpSimplex::basic);
1656 if (iStatus == ClpSimplex::atLowerBound)
1657 b = lower_[iSet];
1658 else
1659 b = upper_[iSet];
1660 if (b) {
1661 for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence+1]; j++) {
1662 int iRow = row_[j];
1663 rhs[iRow] -= element_[j] * b;
1664 }
1665 }
1666 }
1667 }
1668 }
1669 for (iRow = 0; iRow < numberRows; iRow++) {
1670 if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
1671 printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
1672 }
1673 CoinMemcpyN(rhs, numberRows, saveE);
1674 delete [] rhs;
1675 }
1676#endif
1677 if (forceRefresh || (refreshFrequency_ && model->numberIterations() >=
1678 lastRefresh_ + refreshFrequency_)) {
1679 int numberRows = model->numberRows();
1680 int numberColumns = model->numberColumns();
1681 int iRow;
1682 CoinZeroN(rhsOffset_, numberRows);
1683 // do ones at bounds before gub
1684 const double * smallSolution = model->solutionRegion();
1685 const double * element = matrix_->getElements();
1686 const int * row = matrix_->getIndices();
1687 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1688 const int * length = matrix_->getVectorLengths();
1689 int iColumn;
1690 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
1691 if (model->getStatus(iColumn) != ClpSimplex::basic) {
1692 double value = smallSolution[iColumn];
1693 for (CoinBigIndex j = startColumn[iColumn];
1694 j < startColumn[iColumn] + length[iColumn]; j++) {
1695 int jRow = row[j];
1696 rhsOffset_[jRow] -= value * element[j];
1697 }
1698 }
1699 }
1700 if (lowerColumn_ || upperColumn_) {
1701 double * solution = new double [numberGubColumns_];
1702 for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1703 double value = 0.0;
1704 if(getDynamicStatus(iColumn) == atUpperBound)
1705 value = upperColumn_[iColumn];
1706 else if (lowerColumn_)
1707 value = lowerColumn_[iColumn];
1708 solution[iColumn] = value;
1709 }
1710 // ones in gub and in small problem
1711 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
1712 int jFull = id_[iColumn-firstDynamic_];
1713 solution[jFull] = smallSolution[iColumn];
1714 }
1715 // zero all basic in small model
1716 int * pivotVariable = model->pivotVariable();
1717 for (iRow = 0; iRow < numberRows; iRow++) {
1718 int iColumn = pivotVariable[iRow];
1719 if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
1720 int iSequence = id_[iColumn-firstDynamic_];
1721 solution[iSequence] = 0.0;
1722 }
1723 }
1724 // and now compute value to use for key
1725 ClpSimplex::Status iStatus;
1726 int iSet;
1727 for ( iSet = 0; iSet < numberSets_; iSet++) {
1728 iColumn = keyVariable_[iSet];
1729 if (iColumn < numberColumns) {
1730 int iSequence = id_[iColumn-firstDynamic_];
1731 solution[iSequence] = 0.0;
1732 double b = 0.0;
1733 // key is structural - where is slack
1734 iStatus = getStatus(iSet);
1735 assert (iStatus != ClpSimplex::basic);
1736 if (iStatus == ClpSimplex::atLowerBound)
1737 b = lowerSet_[iSet];
1738 else
1739 b = upperSet_[iSet];
1740 // subtract out others at bounds
1741 for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
1742 b -= solution[j];
1743 solution[iSequence] = b;
1744 }
1745 }
1746 for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
1747 double value = solution[iColumn];
1748 if (value) {
1749 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
1750 int iRow = row_[j];
1751 rhsOffset_[iRow] -= element_[j] * value;
1752 }
1753 }
1754 }
1755 // now do lower and upper bounds on sets
1756 // and offset
1757 double objectiveOffset = 0.0;
1758 for ( iSet = 0; iSet < numberSets_; iSet++) {
1759 iColumn = keyVariable_[iSet];
1760 double shift = 0.0;
1761 for (CoinBigIndex j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++) {
1762 if (getDynamicStatus(j) != inSmall) {
1763 double value = 0.0;
1764 if (getDynamicStatus(j) == atLowerBound) {
1765 if (lowerColumn_)
1766 value = lowerColumn_[j];
1767 } else {
1768 value = upperColumn_[j];
1769 }
1770 if (j != iColumn)
1771 shift += value;
1772 objectiveOffset += value * cost_[j];
1773 }
1774 }
1775 if (lowerSet_[iSet] > -1.0e20)
1776 lower_[iSet] = lowerSet_[iSet] - shift;
1777 if (upperSet_[iSet] < 1.0e20)
1778 upper_[iSet] = upperSet_[iSet] - shift;
1779 }
1780 delete [] solution;
1781 model->setObjectiveOffset(objectiveOffset_ - objectiveOffset);
1782 } else {
1783 // no bounds
1784 ClpSimplex::Status iStatus;
1785 for (int iSet = 0; iSet < numberSets_; iSet++) {
1786 int iColumn = keyVariable_[iSet];
1787 if (iColumn < numberColumns) {
1788 int iSequence = id_[iColumn-firstDynamic_];
1789 double b = 0.0;
1790 // key is structural - where is slack
1791 iStatus = getStatus(iSet);
1792 assert (iStatus != ClpSimplex::basic);
1793 if (iStatus == ClpSimplex::atLowerBound)
1794 b = lower_[iSet];
1795 else
1796 b = upper_[iSet];
1797 if (b) {
1798 for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence+1]; j++) {
1799 int iRow = row_[j];
1800 rhsOffset_[iRow] -= element_[j] * b;
1801 }
1802 }
1803 }
1804 }
1805 }
1806#ifdef CLP_DEBUG
1807 if (saveE) {
1808 for (iRow = 0; iRow < numberRows; iRow++) {
1809 if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3)
1810 printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]);
1811 }
1812 delete [] saveE;
1813 }
1814#endif
1815 lastRefresh_ = model->numberIterations();
1816 }
1817 }
1818 return rhsOffset_;
1819}
1820/*
1821 update information for a pivot (and effective rhs)
1822*/
1823int
1824ClpGubDynamicMatrix::updatePivot(ClpSimplex * model, double oldInValue, double oldOutValue)
1825{
1826
1827 // now update working model
1828 int sequenceIn = model->sequenceIn();
1829 int sequenceOut = model->sequenceOut();
1830 bool doPrinting = (model->messageHandler()->logLevel() == 63);
1831 bool print = false;
1832 int iSet;
1833 int trueIn = -1;
1834 int trueOut = -1;
1835 int numberRows = model->numberRows();
1836 int numberColumns = model->numberColumns();
1837 if (sequenceIn == firstAvailable_) {
1838 if (doPrinting)
1839 printf("New variable ");
1840 if (sequenceIn != sequenceOut) {
1841 insertNonBasic(firstAvailable_, backward_[firstAvailable_]);
1842 setDynamicStatus(id_[sequenceIn-firstDynamic_], inSmall);
1843 firstAvailable_++;
1844 } else {
1845 int bigSequence = id_[sequenceIn-firstDynamic_];
1846 if (model->getStatus(sequenceIn) == ClpSimplex::atUpperBound)
1847 setDynamicStatus(bigSequence, atUpperBound);
1848 else
1849 setDynamicStatus(bigSequence, atLowerBound);
1850 }
1851 synchronize(model, 8);
1852 }
1853 if (sequenceIn < lastDynamic_) {
1854 iSet = backward_[sequenceIn];
1855 if (iSet >= 0) {
1856 int bigSequence = id_[sequenceIn-firstDynamic_];
1857 trueIn = bigSequence + numberRows + numberColumns + numberSets_;
1858 if (doPrinting)
1859 printf(" incoming set %d big seq %d", iSet, bigSequence);
1860 print = true;
1861 }
1862 } else if (sequenceIn >= numberRows + numberColumns) {
1863 trueIn = numberRows + numberColumns + gubSlackIn_;
1864 }
1865 if (sequenceOut < lastDynamic_) {
1866 iSet = backward_[sequenceOut];
1867 if (iSet >= 0) {
1868 int bigSequence = id_[sequenceOut-firstDynamic_];
1869 trueOut = bigSequence + firstDynamic_;
1870 if (getDynamicStatus(bigSequence) != inSmall) {
1871 if (model->getStatus(sequenceOut) == ClpSimplex::atUpperBound)
1872 setDynamicStatus(bigSequence, atUpperBound);
1873 else
1874 setDynamicStatus(bigSequence, atLowerBound);
1875 }
1876 if (doPrinting)
1877 printf(" ,outgoing set %d big seq %d,", iSet, bigSequence);
1878 print = true;
1879 model->setSequenceIn(sequenceOut);
1880 synchronize(model, 8);
1881 model->setSequenceIn(sequenceIn);
1882 }
1883 }
1884 if (print && doPrinting)
1885 printf("\n");
1886 ClpGubMatrix::updatePivot(model, oldInValue, oldOutValue);
1887 // Redo true in and out
1888 if (trueIn >= 0)
1889 trueSequenceIn_ = trueIn;
1890 if (trueOut >= 0)
1891 trueSequenceOut_ = trueOut;
1892 if (doPrinting && 0) {
1893 for (int i = 0; i < numberSets_; i++) {
1894 printf("set %d key %d lower %g upper %g\n", i, keyVariable_[i], lower_[i], upper_[i]);
1895 for (int j = fullStart_[i]; j < fullStart_[i+1]; j++)
1896 if (getDynamicStatus(j) == atUpperBound) {
1897 bool print = true;
1898 for (int k = firstDynamic_; k < firstAvailable_; k++) {
1899 if (id_[k-firstDynamic_] == j)
1900 print = false;
1901 if (id_[k-firstDynamic_] == j)
1902 assert(getDynamicStatus(j) == inSmall);
1903 }
1904 if (print)
1905 printf("variable %d at ub\n", j);
1906 }
1907 }
1908 }
1909#ifdef CLP_DEBUG
1910 char * inSmall = new char [numberGubColumns_];
1911 memset(inSmall, 0, numberGubColumns_);
1912 for (int i = 0; i < numberGubColumns_; i++)
1913 if (getDynamicStatus(i) == ClpGubDynamicMatrix::inSmall)
1914 inSmall[i] = 1;
1915 for (int i = firstDynamic_; i < firstAvailable_; i++) {
1916 int k = id_[i-firstDynamic_];
1917 inSmall[k] = 0;
1918 }
1919 for (int i = 0; i < numberGubColumns_; i++)
1920 assert (!inSmall[i]);
1921 delete [] inSmall;
1922#endif
1923 return 0;
1924}
1925void
1926ClpGubDynamicMatrix::times(double scalar,
1927 const double * x, double * y) const
1928{
1929 if (model_->specialOptions() != 16) {
1930 ClpPackedMatrix::times(scalar, x, y);
1931 } else {
1932 int iRow;
1933 int numberColumns = model_->numberColumns();
1934 int numberRows = model_->numberRows();
1935 const double * element = matrix_->getElements();
1936 const int * row = matrix_->getIndices();
1937 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1938 const int * length = matrix_->getVectorLengths();
1939 int * pivotVariable = model_->pivotVariable();
1940 int numberToDo = 0;
1941 for (iRow = 0; iRow < numberRows; iRow++) {
1942 y[iRow] -= scalar * rhsOffset_[iRow];
1943 int iColumn = pivotVariable[iRow];
1944 if (iColumn < numberColumns) {
1945 int iSet = backward_[iColumn];
1946 if (iSet >= 0 && toIndex_[iSet] < 0) {
1947 toIndex_[iSet] = 0;
1948 fromIndex_[numberToDo++] = iSet;
1949 }
1950 CoinBigIndex j;
1951 double value = scalar * x[iColumn];
1952 if (value) {
1953 for (j = startColumn[iColumn];
1954 j < startColumn[iColumn] + length[iColumn]; j++) {
1955 int jRow = row[j];
1956 y[jRow] += value * element[j];
1957 }
1958 }
1959 }
1960 }
1961 // and gubs which are interacting
1962 for (int jSet = 0; jSet < numberToDo; jSet++) {
1963 int iSet = fromIndex_[jSet];
1964 toIndex_[iSet] = -1;
1965 int iKey = keyVariable_[iSet];
1966 if (iKey < numberColumns) {
1967 double valueKey;
1968 if (getStatus(iSet) == ClpSimplex::atLowerBound)
1969 valueKey = lower_[iSet];
1970 else
1971 valueKey = upper_[iSet];
1972 double value = scalar * (x[iKey] - valueKey);
1973 if (value) {
1974 for (CoinBigIndex j = startColumn[iKey];
1975 j < startColumn[iKey] + length[iKey]; j++) {
1976 int jRow = row[j];
1977 y[jRow] += value * element[j];
1978 }
1979 }
1980 }
1981 }
1982 }
1983}
1984/* Just for debug - may be extended to other matrix types later.
1985 Returns number and sum of primal infeasibilities.
1986*/
1987int
1988ClpGubDynamicMatrix::checkFeasible(ClpSimplex * /*model*/, double & sum) const
1989{
1990 int numberRows = model_->numberRows();
1991 double * rhs = new double[numberRows];
1992 int numberColumns = model_->numberColumns();
1993 int iRow;
1994 CoinZeroN(rhs, numberRows);
1995 // do ones at bounds before gub
1996 const double * smallSolution = model_->solutionRegion();
1997 const double * element = matrix_->getElements();
1998 const int * row = matrix_->getIndices();
1999 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
2000 const int * length = matrix_->getVectorLengths();
2001 int iColumn;
2002 int numberInfeasible = 0;
2003 const double * rowLower = model_->rowLower();
2004 const double * rowUpper = model_->rowUpper();
2005 sum = 0.0;
2006 for (iRow = 0; iRow < numberRows; iRow++) {
2007 double value = smallSolution[numberColumns+iRow];
2008 if (value < rowLower[iRow] - 1.0e-5 ||
2009 value > rowUpper[iRow] + 1.0e-5) {
2010 //printf("row %d %g %g %g\n",
2011 // iRow,rowLower[iRow],value,rowUpper[iRow]);
2012 numberInfeasible++;
2013 sum += CoinMax(rowLower[iRow] - value, value - rowUpper[iRow]);
2014 }
2015 rhs[iRow] = value;
2016 }
2017 const double * columnLower = model_->columnLower();
2018 const double * columnUpper = model_->columnUpper();
2019 for (iColumn = 0; iColumn < firstDynamic_; iColumn++) {
2020 double value = smallSolution[iColumn];
2021 if (value < columnLower[iColumn] - 1.0e-5 ||
2022 value > columnUpper[iColumn] + 1.0e-5) {
2023 //printf("column %d %g %g %g\n",
2024 // iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
2025 numberInfeasible++;
2026 sum += CoinMax(columnLower[iColumn] - value, value - columnUpper[iColumn]);
2027 }
2028 for (CoinBigIndex j = startColumn[iColumn];
2029 j < startColumn[iColumn] + length[iColumn]; j++) {
2030 int jRow = row[j];
2031 rhs[jRow] -= value * element[j];
2032 }
2033 }
2034 double * solution = new double [numberGubColumns_];
2035 for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
2036 double value = 0.0;
2037 if(getDynamicStatus(iColumn) == atUpperBound)
2038 value = upperColumn_[iColumn];
2039 else if (lowerColumn_)
2040 value = lowerColumn_[iColumn];
2041 solution[iColumn] = value;
2042 }
2043 // ones in small and gub
2044 for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) {
2045 int jFull = id_[iColumn-firstDynamic_];
2046 solution[jFull] = smallSolution[iColumn];
2047 }
2048 // fill in all basic in small model
2049 int * pivotVariable = model_->pivotVariable();
2050 for (iRow = 0; iRow < numberRows; iRow++) {
2051 int iColumn = pivotVariable[iRow];
2052 if (iColumn >= firstDynamic_ && iColumn < lastDynamic_) {
2053 int iSequence = id_[iColumn-firstDynamic_];
2054 solution[iSequence] = smallSolution[iColumn];
2055 }
2056 }
2057 // and now compute value to use for key
2058 ClpSimplex::Status iStatus;
2059 for (int iSet = 0; iSet < numberSets_; iSet++) {
2060 iColumn = keyVariable_[iSet];
2061 if (iColumn < numberColumns) {
2062 int iSequence = id_[iColumn-firstDynamic_];
2063 solution[iSequence] = 0.0;
2064 double b = 0.0;
2065 // key is structural - where is slack
2066 iStatus = getStatus(iSet);
2067 assert (iStatus != ClpSimplex::basic);
2068 if (iStatus == ClpSimplex::atLowerBound)
2069 b = lower_[iSet];
2070 else
2071 b = upper_[iSet];
2072 // subtract out others at bounds
2073 for (int j = fullStart_[iSet]; j < fullStart_[iSet+1]; j++)
2074 b -= solution[j];
2075 solution[iSequence] = b;
2076 }
2077 }
2078 for (iColumn = 0; iColumn < numberGubColumns_; iColumn++) {
2079 double value = solution[iColumn];
2080 if ((lowerColumn_ && value < lowerColumn_[iColumn] - 1.0e-5) ||
2081 (!lowerColumn_ && value < -1.0e-5) ||
2082 (upperColumn_ && value > upperColumn_[iColumn] + 1.0e-5)) {
2083 //printf("column %d %g %g %g\n",
2084 // iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
2085 numberInfeasible++;
2086 }
2087 if (value) {
2088 for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) {
2089 int iRow = row_[j];
2090 rhs[iRow] -= element_[j] * value;
2091 }
2092 }
2093 }
2094 for (iRow = 0; iRow < numberRows; iRow++) {
2095 if (fabs(rhs[iRow]) > 1.0e-5)
2096 printf("rhs mismatch %d %g\n", iRow, rhs[iRow]);
2097 }
2098 delete [] solution;
2099 delete [] rhs;
2100 return numberInfeasible;
2101}
2102// Cleans data after setWarmStart
2103void
2104ClpGubDynamicMatrix::cleanData(ClpSimplex * model)
2105{
2106 // and redo chains
2107 int numberColumns = model->numberColumns();
2108 int iColumn;
2109 // do backward
2110 int * mark = new int [numberGubColumns_];
2111 for (iColumn = 0; iColumn < numberGubColumns_; iColumn++)
2112 mark[iColumn] = -1;
2113 int i;
2114 for (i = 0; i < firstDynamic_; i++) {
2115 assert (backward_[i] == -1);
2116 next_[i] = -1;
2117 }
2118 for (i = firstDynamic_; i < firstAvailable_; i++) {
2119 iColumn = id_[i-firstDynamic_];
2120 mark[iColumn] = i;
2121 }
2122 for (i = 0; i < numberSets_; i++) {
2123 int iKey = keyVariable_[i];
2124 int lastNext = -1;
2125 int firstNext = -1;
2126 for (CoinBigIndex k = fullStart_[i]; k < fullStart_[i+1]; k++) {
2127 iColumn = mark[k];
2128 if (iColumn >= 0) {
2129 if (iColumn != iKey) {
2130 if (lastNext >= 0)
2131 next_[lastNext] = iColumn;
2132 else
2133 firstNext = iColumn;
2134 lastNext = iColumn;
2135 }
2136 backward_[iColumn] = i;
2137 }
2138 }
2139 setFeasible(i);
2140 if (firstNext >= 0) {
2141 // others
2142 next_[iKey] = firstNext;
2143 next_[lastNext] = -(iKey + 1);
2144 } else if (iKey < numberColumns) {
2145 next_[iKey] = -(iKey + 1);
2146 }
2147 }
2148 delete [] mark;
2149 // fill matrix
2150 double * element = matrix_->getMutableElements();
2151 int * row = matrix_->getMutableIndices();
2152 CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
2153 int * length = matrix_->getMutableVectorLengths();
2154 CoinBigIndex numberElements = startColumn[firstDynamic_];
2155 for (i = firstDynamic_; i < firstAvailable_; i++) {
2156 int iColumn = id_[i-firstDynamic_];
2157 int numberThis = startColumn_[iColumn+1] - startColumn_[iColumn];
2158 length[i] = numberThis;
2159 for (CoinBigIndex jBigIndex = startColumn_[iColumn];
2160 jBigIndex < startColumn_[iColumn+1]; jBigIndex++) {
2161 row[numberElements] = row_[jBigIndex];
2162 element[numberElements++] = element_[jBigIndex];
2163 }
2164 startColumn[i+1] = numberElements;
2165 }
2166}
2167