1/* $Id: ClpGubMatrix.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 "ClpGubMatrix.hpp"
19//#include "ClpGubDynamicMatrix.hpp"
20#include "ClpMessage.hpp"
21//#define CLP_DEBUG
22//#define CLP_DEBUG_PRINT
23//#############################################################################
24// Constructors / Destructor / Assignment
25//#############################################################################
26
27//-------------------------------------------------------------------
28// Default Constructor
29//-------------------------------------------------------------------
30ClpGubMatrix::ClpGubMatrix ()
31 : ClpPackedMatrix(),
32 sumDualInfeasibilities_(0.0),
33 sumPrimalInfeasibilities_(0.0),
34 sumOfRelaxedDualInfeasibilities_(0.0),
35 sumOfRelaxedPrimalInfeasibilities_(0.0),
36 infeasibilityWeight_(0.0),
37 start_(NULL),
38 end_(NULL),
39 lower_(NULL),
40 upper_(NULL),
41 status_(NULL),
42 saveStatus_(NULL),
43 savedKeyVariable_(NULL),
44 backward_(NULL),
45 backToPivotRow_(NULL),
46 changeCost_(NULL),
47 keyVariable_(NULL),
48 next_(NULL),
49 toIndex_(NULL),
50 fromIndex_(NULL),
51 model_(NULL),
52 numberDualInfeasibilities_(0),
53 numberPrimalInfeasibilities_(0),
54 noCheck_(-1),
55 numberSets_(0),
56 saveNumber_(0),
57 possiblePivotKey_(0),
58 gubSlackIn_(-1),
59 firstGub_(0),
60 lastGub_(0),
61 gubType_(0)
62{
63 setType(16);
64}
65
66//-------------------------------------------------------------------
67// Copy constructor
68//-------------------------------------------------------------------
69ClpGubMatrix::ClpGubMatrix (const ClpGubMatrix & rhs)
70 : ClpPackedMatrix(rhs)
71{
72 numberSets_ = rhs.numberSets_;
73 saveNumber_ = rhs.saveNumber_;
74 possiblePivotKey_ = rhs.possiblePivotKey_;
75 gubSlackIn_ = rhs.gubSlackIn_;
76 start_ = ClpCopyOfArray(rhs.start_, numberSets_);
77 end_ = ClpCopyOfArray(rhs.end_, numberSets_);
78 lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
79 upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
80 status_ = ClpCopyOfArray(rhs.status_, numberSets_);
81 saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
82 savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
83 int numberColumns = getNumCols();
84 backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
85 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
86 changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
87 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
88 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
89 // find longest set
90 int * longest = new int[numberSets_];
91 CoinZeroN(longest, numberSets_);
92 int j;
93 for (j = 0; j < numberColumns; j++) {
94 int iSet = backward_[j];
95 if (iSet >= 0)
96 longest[iSet]++;
97 }
98 int length = 0;
99 for (j = 0; j < numberSets_; j++)
100 length = CoinMax(length, longest[j]);
101 next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
102 toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
103 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
104 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
105 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
106 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
107 infeasibilityWeight_ = rhs.infeasibilityWeight_;
108 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
109 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
110 noCheck_ = rhs.noCheck_;
111 firstGub_ = rhs.firstGub_;
112 lastGub_ = rhs.lastGub_;
113 gubType_ = rhs.gubType_;
114 model_ = rhs.model_;
115}
116
117//-------------------------------------------------------------------
118// assign matrix (for space reasons)
119//-------------------------------------------------------------------
120ClpGubMatrix::ClpGubMatrix (CoinPackedMatrix * rhs)
121 : ClpPackedMatrix(rhs),
122 sumDualInfeasibilities_(0.0),
123 sumPrimalInfeasibilities_(0.0),
124 sumOfRelaxedDualInfeasibilities_(0.0),
125 sumOfRelaxedPrimalInfeasibilities_(0.0),
126 infeasibilityWeight_(0.0),
127 start_(NULL),
128 end_(NULL),
129 lower_(NULL),
130 upper_(NULL),
131 status_(NULL),
132 saveStatus_(NULL),
133 savedKeyVariable_(NULL),
134 backward_(NULL),
135 backToPivotRow_(NULL),
136 changeCost_(NULL),
137 keyVariable_(NULL),
138 next_(NULL),
139 toIndex_(NULL),
140 fromIndex_(NULL),
141 model_(NULL),
142 numberDualInfeasibilities_(0),
143 numberPrimalInfeasibilities_(0),
144 noCheck_(-1),
145 numberSets_(0),
146 saveNumber_(0),
147 possiblePivotKey_(0),
148 gubSlackIn_(-1),
149 firstGub_(0),
150 lastGub_(0),
151 gubType_(0)
152{
153 setType(16);
154}
155
156/* This takes over ownership (for space reasons) and is the
157 real constructor*/
158ClpGubMatrix::ClpGubMatrix(ClpPackedMatrix * matrix, int numberSets,
159 const int * start, const int * end,
160 const double * lower, const double * upper,
161 const unsigned char * status)
162 : ClpPackedMatrix(matrix->matrix()),
163 sumDualInfeasibilities_(0.0),
164 sumPrimalInfeasibilities_(0.0),
165 sumOfRelaxedDualInfeasibilities_(0.0),
166 sumOfRelaxedPrimalInfeasibilities_(0.0),
167 numberDualInfeasibilities_(0),
168 numberPrimalInfeasibilities_(0),
169 saveNumber_(0),
170 possiblePivotKey_(0),
171 gubSlackIn_(-1)
172{
173 model_ = NULL;
174 numberSets_ = numberSets;
175 start_ = ClpCopyOfArray(start, numberSets_);
176 end_ = ClpCopyOfArray(end, numberSets_);
177 lower_ = ClpCopyOfArray(lower, numberSets_);
178 upper_ = ClpCopyOfArray(upper, numberSets_);
179 // Check valid and ordered
180 int last = -1;
181 int numberColumns = matrix_->getNumCols();
182 int numberRows = matrix_->getNumRows();
183 backward_ = new int[numberColumns];
184 backToPivotRow_ = new int[numberColumns];
185 changeCost_ = new double [numberRows+numberSets_];
186 keyVariable_ = new int[numberSets_];
187 // signal to need new ordering
188 next_ = NULL;
189 for (int iColumn = 0; iColumn < numberColumns; iColumn++)
190 backward_[iColumn] = -1;
191
192 int iSet;
193 for (iSet = 0; iSet < numberSets_; iSet++) {
194 // set key variable as slack
195 keyVariable_[iSet] = iSet + numberColumns;
196 if (start_[iSet] < 0 || start_[iSet] >= numberColumns)
197 throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
198 if (end_[iSet] < 0 || end_[iSet] > numberColumns)
199 throw CoinError("Index out of range", "constructor", "ClpGubMatrix");
200 if (end_[iSet] <= start_[iSet])
201 throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix");
202 if (start_[iSet] < last)
203 throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix");
204 last = end_[iSet];
205 int j;
206 for (j = start_[iSet]; j < end_[iSet]; j++)
207 backward_[j] = iSet;
208 }
209 // Find type of gub
210 firstGub_ = numberColumns + 1;
211 lastGub_ = -1;
212 int i;
213 for (i = 0; i < numberColumns; i++) {
214 if (backward_[i] >= 0) {
215 firstGub_ = CoinMin(firstGub_, i);
216 lastGub_ = CoinMax(lastGub_, i);
217 }
218 }
219 gubType_ = 0;
220 // adjust lastGub_
221 if (lastGub_ > 0)
222 lastGub_++;
223 for (i = firstGub_; i < lastGub_; i++) {
224 if (backward_[i] < 0) {
225 gubType_ = 1;
226 printf("interior non gub %d\n", i);
227 break;
228 }
229 }
230 if (status) {
231 status_ = ClpCopyOfArray(status, numberSets_);
232 } else {
233 status_ = new unsigned char [numberSets_];
234 memset(status_, 0, numberSets_);
235 int i;
236 for (i = 0; i < numberSets_; i++) {
237 // make slack key
238 setStatus(i, ClpSimplex::basic);
239 }
240 }
241 saveStatus_ = new unsigned char [numberSets_];
242 memset(saveStatus_, 0, numberSets_);
243 savedKeyVariable_ = new int [numberSets_];
244 memset(savedKeyVariable_, 0, numberSets_ * sizeof(int));
245 noCheck_ = -1;
246 infeasibilityWeight_ = 0.0;
247}
248
249ClpGubMatrix::ClpGubMatrix (const CoinPackedMatrix & rhs)
250 : ClpPackedMatrix(rhs),
251 sumDualInfeasibilities_(0.0),
252 sumPrimalInfeasibilities_(0.0),
253 sumOfRelaxedDualInfeasibilities_(0.0),
254 sumOfRelaxedPrimalInfeasibilities_(0.0),
255 infeasibilityWeight_(0.0),
256 start_(NULL),
257 end_(NULL),
258 lower_(NULL),
259 upper_(NULL),
260 status_(NULL),
261 saveStatus_(NULL),
262 savedKeyVariable_(NULL),
263 backward_(NULL),
264 backToPivotRow_(NULL),
265 changeCost_(NULL),
266 keyVariable_(NULL),
267 next_(NULL),
268 toIndex_(NULL),
269 fromIndex_(NULL),
270 model_(NULL),
271 numberDualInfeasibilities_(0),
272 numberPrimalInfeasibilities_(0),
273 noCheck_(-1),
274 numberSets_(0),
275 saveNumber_(0),
276 possiblePivotKey_(0),
277 gubSlackIn_(-1),
278 firstGub_(0),
279 lastGub_(0),
280 gubType_(0)
281{
282 setType(16);
283
284}
285
286//-------------------------------------------------------------------
287// Destructor
288//-------------------------------------------------------------------
289ClpGubMatrix::~ClpGubMatrix ()
290{
291 delete [] start_;
292 delete [] end_;
293 delete [] lower_;
294 delete [] upper_;
295 delete [] status_;
296 delete [] saveStatus_;
297 delete [] savedKeyVariable_;
298 delete [] backward_;
299 delete [] backToPivotRow_;
300 delete [] changeCost_;
301 delete [] keyVariable_;
302 delete [] next_;
303 delete [] toIndex_;
304 delete [] fromIndex_;
305}
306
307//----------------------------------------------------------------
308// Assignment operator
309//-------------------------------------------------------------------
310ClpGubMatrix &
311ClpGubMatrix::operator=(const ClpGubMatrix& rhs)
312{
313 if (this != &rhs) {
314 ClpPackedMatrix::operator=(rhs);
315 delete [] start_;
316 delete [] end_;
317 delete [] lower_;
318 delete [] upper_;
319 delete [] status_;
320 delete [] saveStatus_;
321 delete [] savedKeyVariable_;
322 delete [] backward_;
323 delete [] backToPivotRow_;
324 delete [] changeCost_;
325 delete [] keyVariable_;
326 delete [] next_;
327 delete [] toIndex_;
328 delete [] fromIndex_;
329 numberSets_ = rhs.numberSets_;
330 saveNumber_ = rhs.saveNumber_;
331 possiblePivotKey_ = rhs.possiblePivotKey_;
332 gubSlackIn_ = rhs.gubSlackIn_;
333 start_ = ClpCopyOfArray(rhs.start_, numberSets_);
334 end_ = ClpCopyOfArray(rhs.end_, numberSets_);
335 lower_ = ClpCopyOfArray(rhs.lower_, numberSets_);
336 upper_ = ClpCopyOfArray(rhs.upper_, numberSets_);
337 status_ = ClpCopyOfArray(rhs.status_, numberSets_);
338 saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_);
339 savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_);
340 int numberColumns = getNumCols();
341 backward_ = ClpCopyOfArray(rhs.backward_, numberColumns);
342 backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns);
343 changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_);
344 fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1);
345 keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_);
346 // find longest set
347 int * longest = new int[numberSets_];
348 CoinZeroN(longest, numberSets_);
349 int j;
350 for (j = 0; j < numberColumns; j++) {
351 int iSet = backward_[j];
352 if (iSet >= 0)
353 longest[iSet]++;
354 }
355 int length = 0;
356 for (j = 0; j < numberSets_; j++)
357 length = CoinMax(length, longest[j]);
358 next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length);
359 toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_);
360 sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
361 sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
362 sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
363 sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
364 infeasibilityWeight_ = rhs.infeasibilityWeight_;
365 numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
366 numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
367 noCheck_ = rhs.noCheck_;
368 firstGub_ = rhs.firstGub_;
369 lastGub_ = rhs.lastGub_;
370 gubType_ = rhs.gubType_;
371 model_ = rhs.model_;
372 }
373 return *this;
374}
375//-------------------------------------------------------------------
376// Clone
377//-------------------------------------------------------------------
378ClpMatrixBase * ClpGubMatrix::clone() const
379{
380 return new ClpGubMatrix(*this);
381}
382/* Subset clone (without gaps). Duplicates are allowed
383 and order is as given */
384ClpMatrixBase *
385ClpGubMatrix::subsetClone (int numberRows, const int * whichRows,
386 int numberColumns,
387 const int * whichColumns) const
388{
389 return new ClpGubMatrix(*this, numberRows, whichRows,
390 numberColumns, whichColumns);
391}
392/* Returns a new matrix in reverse order without gaps
393 Is allowed to return NULL if doesn't want to have row copy */
394ClpMatrixBase *
395ClpGubMatrix::reverseOrderedCopy() const
396{
397 return NULL;
398}
399int
400ClpGubMatrix::hiddenRows() const
401{
402 return numberSets_;
403}
404/* Subset constructor (without gaps). Duplicates are allowed
405 and order is as given */
406ClpGubMatrix::ClpGubMatrix (
407 const ClpGubMatrix & rhs,
408 int numberRows, const int * whichRows,
409 int numberColumns, const int * whichColumns)
410 : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
411{
412 // Assuming no gub rows deleted
413 // We also assume all sets in same order
414 // Get array with backward pointers
415 int numberColumnsOld = rhs.matrix_->getNumCols();
416 int * array = new int [ numberColumnsOld];
417 int i;
418 for (i = 0; i < numberColumnsOld; i++)
419 array[i] = -1;
420 for (int iSet = 0; iSet < numberSets_; iSet++) {
421 for (int j = start_[iSet]; j < end_[iSet]; j++)
422 array[j] = iSet;
423 }
424 numberSets_ = -1;
425 int lastSet = -1;
426 bool inSet = false;
427 for (i = 0; i < numberColumns; i++) {
428 int iColumn = whichColumns[i];
429 int iSet = array[iColumn];
430 if (iSet < 0) {
431 inSet = false;
432 } else {
433 if (!inSet) {
434 // start of new set but check okay
435 if (iSet <= lastSet)
436 throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
437 lastSet = iSet;
438 numberSets_++;
439 start_[numberSets_] = i;
440 end_[numberSets_] = i + 1;
441 lower_[numberSets_] = lower_[iSet];
442 upper_[numberSets_] = upper_[iSet];
443 inSet = true;
444 } else {
445 if (iSet < lastSet) {
446 throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix");
447 } else if (iSet == lastSet) {
448 end_[numberSets_] = i + 1;
449 } else {
450 // new set
451 lastSet = iSet;
452 numberSets_++;
453 start_[numberSets_] = i;
454 end_[numberSets_] = i + 1;
455 lower_[numberSets_] = lower_[iSet];
456 upper_[numberSets_] = upper_[iSet];
457 }
458 }
459 }
460 }
461 delete [] array;
462 numberSets_++; // adjust
463 // Find type of gub
464 firstGub_ = numberColumns + 1;
465 lastGub_ = -1;
466 for (i = 0; i < numberColumns; i++) {
467 if (backward_[i] >= 0) {
468 firstGub_ = CoinMin(firstGub_, i);
469 lastGub_ = CoinMax(lastGub_, i);
470 }
471 }
472 if (lastGub_ > 0)
473 lastGub_++;
474 gubType_ = 0;
475 for (i = firstGub_; i < lastGub_; i++) {
476 if (backward_[i] < 0) {
477 gubType_ = 1;
478 break;
479 }
480 }
481
482 // Make sure key is feasible if only key in set
483}
484ClpGubMatrix::ClpGubMatrix (
485 const CoinPackedMatrix & rhs,
486 int numberRows, const int * whichRows,
487 int numberColumns, const int * whichColumns)
488 : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns),
489 sumDualInfeasibilities_(0.0),
490 sumPrimalInfeasibilities_(0.0),
491 sumOfRelaxedDualInfeasibilities_(0.0),
492 sumOfRelaxedPrimalInfeasibilities_(0.0),
493 start_(NULL),
494 end_(NULL),
495 lower_(NULL),
496 upper_(NULL),
497 backward_(NULL),
498 backToPivotRow_(NULL),
499 changeCost_(NULL),
500 keyVariable_(NULL),
501 next_(NULL),
502 toIndex_(NULL),
503 fromIndex_(NULL),
504 numberDualInfeasibilities_(0),
505 numberPrimalInfeasibilities_(0),
506 numberSets_(0),
507 saveNumber_(0),
508 possiblePivotKey_(0),
509 gubSlackIn_(-1),
510 firstGub_(0),
511 lastGub_(0),
512 gubType_(0)
513{
514 setType(16);
515}
516/* Return <code>x * A + y</code> in <code>z</code>.
517 Squashes small elements and knows about ClpSimplex */
518void
519ClpGubMatrix::transposeTimes(const ClpSimplex * model, double scalar,
520 const CoinIndexedVector * rowArray,
521 CoinIndexedVector * y,
522 CoinIndexedVector * columnArray) const
523{
524 columnArray->clear();
525 double * pi = rowArray->denseVector();
526 int numberNonZero = 0;
527 int * index = columnArray->getIndices();
528 double * array = columnArray->denseVector();
529 int numberInRowArray = rowArray->getNumElements();
530 // maybe I need one in OsiSimplex
531 double zeroTolerance = model->zeroTolerance();
532 int numberRows = model->numberRows();
533 ClpPackedMatrix* rowCopy =
534 dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
535 bool packed = rowArray->packedMode();
536 double factor = 0.3;
537 // We may not want to do by row if there may be cache problems
538 int numberColumns = model->numberColumns();
539 // It would be nice to find L2 cache size - for moment 512K
540 // Be slightly optimistic
541 if (numberColumns * sizeof(double) > 1000000) {
542 if (numberRows * 10 < numberColumns)
543 factor = 0.1;
544 else if (numberRows * 4 < numberColumns)
545 factor = 0.15;
546 else if (numberRows * 2 < numberColumns)
547 factor = 0.2;
548 //if (model->numberIterations()%50==0)
549 //printf("%d nonzero\n",numberInRowArray);
550 }
551 // reduce for gub
552 factor *= 0.5;
553 assert (!y->getNumElements());
554 if (numberInRowArray > factor * numberRows || !rowCopy) {
555 // do by column
556 int iColumn;
557 // get matrix data pointers
558 const int * row = matrix_->getIndices();
559 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
560 const int * columnLength = matrix_->getVectorLengths();
561 const double * elementByColumn = matrix_->getElements();
562 const double * rowScale = model->rowScale();
563 int numberColumns = model->numberColumns();
564 int iSet = -1;
565 double djMod = 0.0;
566 if (packed) {
567 // need to expand pi into y
568 assert(y->capacity() >= numberRows);
569 double * piOld = pi;
570 pi = y->denseVector();
571 const int * whichRow = rowArray->getIndices();
572 int i;
573 if (!rowScale) {
574 // modify pi so can collapse to one loop
575 for (i = 0; i < numberInRowArray; i++) {
576 int iRow = whichRow[i];
577 pi[iRow] = scalar * piOld[i];
578 }
579 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
580 if (backward_[iColumn] != iSet) {
581 // get pi on gub row
582 iSet = backward_[iColumn];
583 if (iSet >= 0) {
584 int iBasic = keyVariable_[iSet];
585 if (iBasic < numberColumns) {
586 // get dj without
587 assert (model->getStatus(iBasic) == ClpSimplex::basic);
588 djMod = 0.0;
589 for (CoinBigIndex j = columnStart[iBasic];
590 j < columnStart[iBasic] + columnLength[iBasic]; j++) {
591 int jRow = row[j];
592 djMod -= pi[jRow] * elementByColumn[j];
593 }
594 } else {
595 djMod = 0.0;
596 }
597 } else {
598 djMod = 0.0;
599 }
600 }
601 double value = -djMod;
602 CoinBigIndex j;
603 for (j = columnStart[iColumn];
604 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
605 int iRow = row[j];
606 value += pi[iRow] * elementByColumn[j];
607 }
608 if (fabs(value) > zeroTolerance) {
609 array[numberNonZero] = value;
610 index[numberNonZero++] = iColumn;
611 }
612 }
613 } else {
614 // scaled
615 // modify pi so can collapse to one loop
616 for (i = 0; i < numberInRowArray; i++) {
617 int iRow = whichRow[i];
618 pi[iRow] = scalar * piOld[i] * rowScale[iRow];
619 }
620 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
621 if (backward_[iColumn] != iSet) {
622 // get pi on gub row
623 iSet = backward_[iColumn];
624 if (iSet >= 0) {
625 int iBasic = keyVariable_[iSet];
626 if (iBasic < numberColumns) {
627 // get dj without
628 assert (model->getStatus(iBasic) == ClpSimplex::basic);
629 djMod = 0.0;
630 // scaled
631 for (CoinBigIndex j = columnStart[iBasic];
632 j < columnStart[iBasic] + columnLength[iBasic]; j++) {
633 int jRow = row[j];
634 djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow];
635 }
636 } else {
637 djMod = 0.0;
638 }
639 } else {
640 djMod = 0.0;
641 }
642 }
643 double value = -djMod;
644 CoinBigIndex j;
645 const double * columnScale = model->columnScale();
646 for (j = columnStart[iColumn];
647 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
648 int iRow = row[j];
649 value += pi[iRow] * elementByColumn[j];
650 }
651 value *= columnScale[iColumn];
652 if (fabs(value) > zeroTolerance) {
653 array[numberNonZero] = value;
654 index[numberNonZero++] = iColumn;
655 }
656 }
657 }
658 // zero out
659 for (i = 0; i < numberInRowArray; i++) {
660 int iRow = whichRow[i];
661 pi[iRow] = 0.0;
662 }
663 } else {
664 // code later
665 assert (packed);
666 if (!rowScale) {
667 if (scalar == -1.0) {
668 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
669 double value = 0.0;
670 CoinBigIndex j;
671 for (j = columnStart[iColumn];
672 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
673 int iRow = row[j];
674 value += pi[iRow] * elementByColumn[j];
675 }
676 if (fabs(value) > zeroTolerance) {
677 index[numberNonZero++] = iColumn;
678 array[iColumn] = -value;
679 }
680 }
681 } else if (scalar == 1.0) {
682 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
683 double value = 0.0;
684 CoinBigIndex j;
685 for (j = columnStart[iColumn];
686 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
687 int iRow = row[j];
688 value += pi[iRow] * elementByColumn[j];
689 }
690 if (fabs(value) > zeroTolerance) {
691 index[numberNonZero++] = iColumn;
692 array[iColumn] = value;
693 }
694 }
695 } else {
696 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
697 double value = 0.0;
698 CoinBigIndex j;
699 for (j = columnStart[iColumn];
700 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
701 int iRow = row[j];
702 value += pi[iRow] * elementByColumn[j];
703 }
704 value *= scalar;
705 if (fabs(value) > zeroTolerance) {
706 index[numberNonZero++] = iColumn;
707 array[iColumn] = value;
708 }
709 }
710 }
711 } else {
712 // scaled
713 if (scalar == -1.0) {
714 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
715 double value = 0.0;
716 CoinBigIndex j;
717 const double * columnScale = model->columnScale();
718 for (j = columnStart[iColumn];
719 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
720 int iRow = row[j];
721 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
722 }
723 value *= columnScale[iColumn];
724 if (fabs(value) > zeroTolerance) {
725 index[numberNonZero++] = iColumn;
726 array[iColumn] = -value;
727 }
728 }
729 } else if (scalar == 1.0) {
730 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
731 double value = 0.0;
732 CoinBigIndex j;
733 const double * columnScale = model->columnScale();
734 for (j = columnStart[iColumn];
735 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
736 int iRow = row[j];
737 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
738 }
739 value *= columnScale[iColumn];
740 if (fabs(value) > zeroTolerance) {
741 index[numberNonZero++] = iColumn;
742 array[iColumn] = value;
743 }
744 }
745 } else {
746 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
747 double value = 0.0;
748 CoinBigIndex j;
749 const double * columnScale = model->columnScale();
750 for (j = columnStart[iColumn];
751 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
752 int iRow = row[j];
753 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
754 }
755 value *= scalar * columnScale[iColumn];
756 if (fabs(value) > zeroTolerance) {
757 index[numberNonZero++] = iColumn;
758 array[iColumn] = value;
759 }
760 }
761 }
762 }
763 }
764 columnArray->setNumElements(numberNonZero);
765 y->setNumElements(0);
766 } else {
767 // do by row
768 transposeTimesByRow(model, scalar, rowArray, y, columnArray);
769 }
770 if (packed)
771 columnArray->setPackedMode(true);
772 if (0) {
773 columnArray->checkClean();
774 int numberNonZero = columnArray->getNumElements();
775 int * index = columnArray->getIndices();
776 double * array = columnArray->denseVector();
777 int i;
778 for (i = 0; i < numberNonZero; i++) {
779 int j = index[i];
780 double value;
781 if (packed)
782 value = array[i];
783 else
784 value = array[j];
785 printf("Ti %d %d %g\n", i, j, value);
786 }
787 }
788}
789/* Return <code>x * A + y</code> in <code>z</code>.
790 Squashes small elements and knows about ClpSimplex */
791void
792ClpGubMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
793 const CoinIndexedVector * rowArray,
794 CoinIndexedVector * y,
795 CoinIndexedVector * columnArray) const
796{
797 // Do packed part
798 ClpPackedMatrix::transposeTimesByRow(model, scalar, rowArray, y, columnArray);
799 if (numberSets_) {
800 /* what we need to do is do by row as normal but get list of sets touched
801 and then update those ones */
802 abort();
803 }
804}
805/* Return <code>x *A in <code>z</code> but
806 just for indices in y. */
807void
808ClpGubMatrix::subsetTransposeTimes(const ClpSimplex * model,
809 const CoinIndexedVector * rowArray,
810 const CoinIndexedVector * y,
811 CoinIndexedVector * columnArray) const
812{
813 columnArray->clear();
814 double * pi = rowArray->denseVector();
815 double * array = columnArray->denseVector();
816 int jColumn;
817 // get matrix data pointers
818 const int * row = matrix_->getIndices();
819 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
820 const int * columnLength = matrix_->getVectorLengths();
821 const double * elementByColumn = matrix_->getElements();
822 const double * rowScale = model->rowScale();
823 int numberToDo = y->getNumElements();
824 const int * which = y->getIndices();
825 assert (!rowArray->packedMode());
826 columnArray->setPacked();
827 int numberTouched = 0;
828 if (!rowScale) {
829 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
830 int iColumn = which[jColumn];
831 double value = 0.0;
832 CoinBigIndex j;
833 for (j = columnStart[iColumn];
834 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
835 int iRow = row[j];
836 value += pi[iRow] * elementByColumn[j];
837 }
838 array[jColumn] = value;
839 if (value) {
840 int iSet = backward_[iColumn];
841 if (iSet >= 0) {
842 int iBasic = keyVariable_[iSet];
843 if (iBasic == iColumn) {
844 toIndex_[iSet] = jColumn;
845 fromIndex_[numberTouched++] = iSet;
846 }
847 }
848 }
849 }
850 } else {
851 // scaled
852 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
853 int iColumn = which[jColumn];
854 double value = 0.0;
855 CoinBigIndex j;
856 const double * columnScale = model->columnScale();
857 for (j = columnStart[iColumn];
858 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
859 int iRow = row[j];
860 value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
861 }
862 value *= columnScale[iColumn];
863 array[jColumn] = value;
864 if (value) {
865 int iSet = backward_[iColumn];
866 if (iSet >= 0) {
867 int iBasic = keyVariable_[iSet];
868 if (iBasic == iColumn) {
869 toIndex_[iSet] = jColumn;
870 fromIndex_[numberTouched++] = iSet;
871 }
872 }
873 }
874 }
875 }
876 // adjust djs
877 for (jColumn = 0; jColumn < numberToDo; jColumn++) {
878 int iColumn = which[jColumn];
879 int iSet = backward_[iColumn];
880 if (iSet >= 0) {
881 int kColumn = toIndex_[iSet];
882 if (kColumn >= 0)
883 array[jColumn] -= array[kColumn];
884 }
885 }
886 // and clear basic
887 for (int j = 0; j < numberTouched; j++) {
888 int iSet = fromIndex_[j];
889 int kColumn = toIndex_[iSet];
890 toIndex_[iSet] = -1;
891 array[kColumn] = 0.0;
892 }
893}
894/// returns number of elements in column part of basis,
895CoinBigIndex
896ClpGubMatrix::countBasis(const int * whichColumn,
897 int & numberColumnBasic)
898{
899 int i;
900 int numberColumns = getNumCols();
901 const int * columnLength = matrix_->getVectorLengths();
902 int numberRows = getNumRows();
903 int numberBasic = 0;
904 CoinBigIndex numberElements = 0;
905 int lastSet = -1;
906 int key = -1;
907 int keyLength = -1;
908 double * work = new double[numberRows];
909 CoinZeroN(work, numberRows);
910 char * mark = new char[numberRows];
911 CoinZeroN(mark, numberRows);
912 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
913 const int * row = matrix_->getIndices();
914 const double * elementByColumn = matrix_->getElements();
915 //ClpGubDynamicMatrix* gubx =
916 //dynamic_cast< ClpGubDynamicMatrix*>(this);
917 //int * id = gubx->id();
918 // just count
919 for (i = 0; i < numberColumnBasic; i++) {
920 int iColumn = whichColumn[i];
921 int iSet = backward_[iColumn];
922 int length = columnLength[iColumn];
923 if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
924 numberElements += length;
925 numberBasic++;
926 //printf("non gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,length);
927 } else {
928 // in gub set
929 if (iColumn != keyVariable_[iSet]) {
930 numberBasic++;
931 CoinBigIndex j;
932 // not key
933 if (lastSet < iSet) {
934 // erase work
935 if (key >= 0) {
936 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
937 work[row[j]] = 0.0;
938 }
939 key = keyVariable_[iSet];
940 lastSet = iSet;
941 keyLength = columnLength[key];
942 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++)
943 work[row[j]] = elementByColumn[j];
944 }
945 int extra = keyLength;
946 for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
947 int iRow = row[j];
948 double keyValue = work[iRow];
949 double value = elementByColumn[j];
950 if (!keyValue) {
951 if (fabs(value) > 1.0e-20)
952 extra++;
953 } else {
954 value -= keyValue;
955 if (fabs(value) <= 1.0e-20)
956 extra--;
957 }
958 }
959 numberElements += extra;
960 //printf("gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,extra);
961 }
962 }
963 }
964 delete [] work;
965 delete [] mark;
966 // update number of column basic
967 numberColumnBasic = numberBasic;
968 return numberElements;
969}
970void
971ClpGubMatrix::fillBasis(ClpSimplex * model,
972 const int * whichColumn,
973 int & numberColumnBasic,
974 int * indexRowU, int * start,
975 int * rowCount, int * columnCount,
976 CoinFactorizationDouble * elementU)
977{
978 int i;
979 int numberColumns = getNumCols();
980 const int * columnLength = matrix_->getVectorLengths();
981 int numberRows = getNumRows();
982 assert (next_ || !elementU) ;
983 CoinBigIndex numberElements = start[0];
984 int lastSet = -1;
985 int key = -1;
986 int keyLength = -1;
987 double * work = new double[numberRows];
988 CoinZeroN(work, numberRows);
989 char * mark = new char[numberRows];
990 CoinZeroN(mark, numberRows);
991 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
992 const int * row = matrix_->getIndices();
993 const double * elementByColumn = matrix_->getElements();
994 const double * rowScale = model->rowScale();
995 int numberBasic = 0;
996 if (0) {
997 printf("%d basiccolumns\n", numberColumnBasic);
998 int i;
999 for (i = 0; i < numberSets_; i++) {
1000 int k = keyVariable_[i];
1001 if (k < numberColumns) {
1002 printf("key %d on set %d, %d elements\n", k, i, columnStart[k+1] - columnStart[k]);
1003 for (int j = columnStart[k]; j < columnStart[k+1]; j++)
1004 printf("row %d el %g\n", row[j], elementByColumn[j]);
1005 } else {
1006 printf("slack key on set %d\n", i);
1007 }
1008 }
1009 }
1010 // fill
1011 if (!rowScale) {
1012 // no scaling
1013 for (i = 0; i < numberColumnBasic; i++) {
1014 int iColumn = whichColumn[i];
1015 int iSet = backward_[iColumn];
1016 int length = columnLength[iColumn];
1017 if (0) {
1018 int k = iColumn;
1019 printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k+1] - columnStart[k]);
1020 for (int j = columnStart[k]; j < columnStart[k+1]; j++)
1021 printf("row %d el %g\n", row[j], elementByColumn[j]);
1022 }
1023 CoinBigIndex j;
1024 if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
1025 for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1026 double value = elementByColumn[j];
1027 if (fabs(value) > 1.0e-20) {
1028 int iRow = row[j];
1029 indexRowU[numberElements] = iRow;
1030 rowCount[iRow]++;
1031 elementU[numberElements++] = value;
1032 }
1033 }
1034 // end of column
1035 columnCount[numberBasic] = numberElements - start[numberBasic];
1036 numberBasic++;
1037 start[numberBasic] = numberElements;
1038 } else {
1039 // in gub set
1040 if (iColumn != keyVariable_[iSet]) {
1041 // not key
1042 if (lastSet != iSet) {
1043 // erase work
1044 if (key >= 0) {
1045 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1046 int iRow = row[j];
1047 work[iRow] = 0.0;
1048 mark[iRow] = 0;
1049 }
1050 }
1051 key = keyVariable_[iSet];
1052 lastSet = iSet;
1053 keyLength = columnLength[key];
1054 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1055 int iRow = row[j];
1056 work[iRow] = elementByColumn[j];
1057 mark[iRow] = 1;
1058 }
1059 }
1060 for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
1061 int iRow = row[j];
1062 double value = elementByColumn[j];
1063 if (mark[iRow]) {
1064 mark[iRow] = 0;
1065 double keyValue = work[iRow];
1066 value -= keyValue;
1067 }
1068 if (fabs(value) > 1.0e-20) {
1069 indexRowU[numberElements] = iRow;
1070 rowCount[iRow]++;
1071 elementU[numberElements++] = value;
1072 }
1073 }
1074 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1075 int iRow = row[j];
1076 if (mark[iRow]) {
1077 double value = -work[iRow];
1078 if (fabs(value) > 1.0e-20) {
1079 indexRowU[numberElements] = iRow;
1080 rowCount[iRow]++;
1081 elementU[numberElements++] = value;
1082 }
1083 } else {
1084 // just put back mark
1085 mark[iRow] = 1;
1086 }
1087 }
1088 // end of column
1089 columnCount[numberBasic] = numberElements - start[numberBasic];
1090 numberBasic++;
1091 start[numberBasic] = numberElements;
1092 }
1093 }
1094 }
1095 } else {
1096 // scaling
1097 const double * columnScale = model->columnScale();
1098 for (i = 0; i < numberColumnBasic; i++) {
1099 int iColumn = whichColumn[i];
1100 int iSet = backward_[iColumn];
1101 int length = columnLength[iColumn];
1102 CoinBigIndex j;
1103 if (iSet < 0 || keyVariable_[iSet] >= numberColumns) {
1104 double scale = columnScale[iColumn];
1105 for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1106 int iRow = row[j];
1107 double value = elementByColumn[j] * scale * rowScale[iRow];
1108 if (fabs(value) > 1.0e-20) {
1109 indexRowU[numberElements] = iRow;
1110 rowCount[iRow]++;
1111 elementU[numberElements++] = value;
1112 }
1113 }
1114 // end of column
1115 columnCount[numberBasic] = numberElements - start[numberBasic];
1116 numberBasic++;
1117 start[numberBasic] = numberElements;
1118 } else {
1119 // in gub set
1120 if (iColumn != keyVariable_[iSet]) {
1121 double scale = columnScale[iColumn];
1122 // not key
1123 if (lastSet < iSet) {
1124 // erase work
1125 if (key >= 0) {
1126 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1127 int iRow = row[j];
1128 work[iRow] = 0.0;
1129 mark[iRow] = 0;
1130 }
1131 }
1132 key = keyVariable_[iSet];
1133 lastSet = iSet;
1134 keyLength = columnLength[key];
1135 double scale = columnScale[key];
1136 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1137 int iRow = row[j];
1138 work[iRow] = elementByColumn[j] * scale * rowScale[iRow];
1139 mark[iRow] = 1;
1140 }
1141 }
1142 for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) {
1143 int iRow = row[j];
1144 double value = elementByColumn[j] * scale * rowScale[iRow];
1145 if (mark[iRow]) {
1146 mark[iRow] = 0;
1147 double keyValue = work[iRow];
1148 value -= keyValue;
1149 }
1150 if (fabs(value) > 1.0e-20) {
1151 indexRowU[numberElements] = iRow;
1152 rowCount[iRow]++;
1153 elementU[numberElements++] = value;
1154 }
1155 }
1156 for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) {
1157 int iRow = row[j];
1158 if (mark[iRow]) {
1159 double value = -work[iRow];
1160 if (fabs(value) > 1.0e-20) {
1161 indexRowU[numberElements] = iRow;
1162 rowCount[iRow]++;
1163 elementU[numberElements++] = value;
1164 }
1165 } else {
1166 // just put back mark
1167 mark[iRow] = 1;
1168 }
1169 }
1170 // end of column
1171 columnCount[numberBasic] = numberElements - start[numberBasic];
1172 numberBasic++;
1173 start[numberBasic] = numberElements;
1174 }
1175 }
1176 }
1177 }
1178 delete [] work;
1179 delete [] mark;
1180 // update number of column basic
1181 numberColumnBasic = numberBasic;
1182}
1183/* Unpacks a column into an CoinIndexedvector
1184 */
1185void
1186ClpGubMatrix::unpack(const ClpSimplex * model, CoinIndexedVector * rowArray,
1187 int iColumn) const
1188{
1189 assert (iColumn < model->numberColumns());
1190 // Do packed part
1191 ClpPackedMatrix::unpack(model, rowArray, iColumn);
1192 int iSet = backward_[iColumn];
1193 if (iSet >= 0) {
1194 int iBasic = keyVariable_[iSet];
1195 if (iBasic < model->numberColumns()) {
1196 add(model, rowArray, iBasic, -1.0);
1197 }
1198 }
1199}
1200/* Unpacks a column into a CoinIndexedVector
1201** in packed format
1202Note that model is NOT const. Bounds and objective could
1203be modified if doing column generation (just for this variable) */
1204void
1205ClpGubMatrix::unpackPacked(ClpSimplex * model,
1206 CoinIndexedVector * rowArray,
1207 int iColumn) const
1208{
1209 int numberColumns = model->numberColumns();
1210 if (iColumn < numberColumns) {
1211 // Do packed part
1212 ClpPackedMatrix::unpackPacked(model, rowArray, iColumn);
1213 int iSet = backward_[iColumn];
1214 if (iSet >= 0) {
1215 // columns are in order
1216 int iBasic = keyVariable_[iSet];
1217 if (iBasic < numberColumns) {
1218 int number = rowArray->getNumElements();
1219 const double * rowScale = model->rowScale();
1220 const int * row = matrix_->getIndices();
1221 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1222 const int * columnLength = matrix_->getVectorLengths();
1223 const double * elementByColumn = matrix_->getElements();
1224 double * array = rowArray->denseVector();
1225 int * index = rowArray->getIndices();
1226 CoinBigIndex i;
1227 int numberOld = number;
1228 int lastIndex = 0;
1229 int next = index[lastIndex];
1230 if (!rowScale) {
1231 for (i = columnStart[iBasic];
1232 i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1233 int iRow = row[i];
1234 while (iRow > next) {
1235 lastIndex++;
1236 if (lastIndex == numberOld)
1237 next = matrix_->getNumRows();
1238 else
1239 next = index[lastIndex];
1240 }
1241 if (iRow < next) {
1242 array[number] = -elementByColumn[i];
1243 index[number++] = iRow;
1244 } else {
1245 assert (iRow == next);
1246 array[lastIndex] -= elementByColumn[i];
1247 if (!array[lastIndex])
1248 array[lastIndex] = 1.0e-100;
1249 }
1250 }
1251 } else {
1252 // apply scaling
1253 double scale = model->columnScale()[iBasic];
1254 for (i = columnStart[iBasic];
1255 i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1256 int iRow = row[i];
1257 while (iRow > next) {
1258 lastIndex++;
1259 if (lastIndex == numberOld)
1260 next = matrix_->getNumRows();
1261 else
1262 next = index[lastIndex];
1263 }
1264 if (iRow < next) {
1265 array[number] = -elementByColumn[i] * scale * rowScale[iRow];
1266 index[number++] = iRow;
1267 } else {
1268 assert (iRow == next);
1269 array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow];
1270 if (!array[lastIndex])
1271 array[lastIndex] = 1.0e-100;
1272 }
1273 }
1274 }
1275 rowArray->setNumElements(number);
1276 }
1277 }
1278 } else {
1279 // key slack entering
1280 int iBasic = keyVariable_[gubSlackIn_];
1281 assert (iBasic < numberColumns);
1282 int number = 0;
1283 const double * rowScale = model->rowScale();
1284 const int * row = matrix_->getIndices();
1285 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1286 const int * columnLength = matrix_->getVectorLengths();
1287 const double * elementByColumn = matrix_->getElements();
1288 double * array = rowArray->denseVector();
1289 int * index = rowArray->getIndices();
1290 CoinBigIndex i;
1291 if (!rowScale) {
1292 for (i = columnStart[iBasic];
1293 i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1294 int iRow = row[i];
1295 array[number] = elementByColumn[i];
1296 index[number++] = iRow;
1297 }
1298 } else {
1299 // apply scaling
1300 double scale = model->columnScale()[iBasic];
1301 for (i = columnStart[iBasic];
1302 i < columnStart[iBasic] + columnLength[iBasic]; i++) {
1303 int iRow = row[i];
1304 array[number] = elementByColumn[i] * scale * rowScale[iRow];
1305 index[number++] = iRow;
1306 }
1307 }
1308 rowArray->setNumElements(number);
1309 rowArray->setPacked();
1310 }
1311}
1312/* Adds multiple of a column into an CoinIndexedvector
1313 You can use quickAdd to add to vector */
1314void
1315ClpGubMatrix::add(const ClpSimplex * model, CoinIndexedVector * rowArray,
1316 int iColumn, double multiplier) const
1317{
1318 assert (iColumn < model->numberColumns());
1319 // Do packed part
1320 ClpPackedMatrix::add(model, rowArray, iColumn, multiplier);
1321 int iSet = backward_[iColumn];
1322 if (iSet >= 0 && iColumn != keyVariable_[iSet]) {
1323 ClpPackedMatrix::add(model, rowArray, keyVariable_[iSet], -multiplier);
1324 }
1325}
1326/* Adds multiple of a column into an array */
1327void
1328ClpGubMatrix::add(const ClpSimplex * model, double * array,
1329 int iColumn, double multiplier) const
1330{
1331 assert (iColumn < model->numberColumns());
1332 // Do packed part
1333 ClpPackedMatrix::add(model, array, iColumn, multiplier);
1334 if (iColumn < model->numberColumns()) {
1335 int iSet = backward_[iColumn];
1336 if (iSet >= 0 && iColumn != keyVariable_[iSet] && keyVariable_[iSet] < model->numberColumns()) {
1337 ClpPackedMatrix::add(model, array, keyVariable_[iSet], -multiplier);
1338 }
1339 }
1340}
1341// Partial pricing
1342void
1343ClpGubMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1344 int & bestSequence, int & numberWanted)
1345{
1346 numberWanted = currentWanted_;
1347 if (numberSets_) {
1348 // Do packed part before gub
1349 int numberColumns = matrix_->getNumCols();
1350 double ratio = static_cast<double> (firstGub_) /
1351 static_cast<double> (numberColumns);
1352 ClpPackedMatrix::partialPricing(model, startFraction * ratio,
1353 endFraction * ratio, bestSequence, numberWanted);
1354 if (numberWanted || minimumGoodReducedCosts_ < -1) {
1355 // do gub
1356 const double * element = matrix_->getElements();
1357 const int * row = matrix_->getIndices();
1358 const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1359 const int * length = matrix_->getVectorLengths();
1360 const double * rowScale = model->rowScale();
1361 const double * columnScale = model->columnScale();
1362 int iSequence;
1363 CoinBigIndex j;
1364 double tolerance = model->currentDualTolerance();
1365 double * reducedCost = model->djRegion();
1366 const double * duals = model->dualRowSolution();
1367 const double * cost = model->costRegion();
1368 double bestDj;
1369 int numberColumns = model->numberColumns();
1370 int numberRows = model->numberRows();
1371 if (bestSequence >= 0)
1372 bestDj = fabs(this->reducedCost(model, bestSequence));
1373 else
1374 bestDj = tolerance;
1375 int sequenceOut = model->sequenceOut();
1376 int saveSequence = bestSequence;
1377 int startG = firstGub_ + static_cast<int> (startFraction * (lastGub_ - firstGub_));
1378 int endG = firstGub_ + static_cast<int> (endFraction * (lastGub_ - firstGub_));
1379 endG = CoinMin(lastGub_, endG + 1);
1380 // If nothing found yet can go all the way to end
1381 int endAll = endG;
1382 if (bestSequence < 0 && !startG)
1383 endAll = lastGub_;
1384 int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_;
1385 int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_;
1386 int nSets = 0;
1387 int iSet = -1;
1388 double djMod = 0.0;
1389 double infeasibilityCost = model->infeasibilityCost();
1390 if (rowScale) {
1391 double bestDjMod = 0.0;
1392 // scaled
1393 for (iSequence = startG; iSequence < endAll; iSequence++) {
1394 if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
1395 // give up
1396 numberWanted = 0;
1397 break;
1398 } else if (iSequence == endG && bestSequence >= 0) {
1399 break;
1400 }
1401 if (backward_[iSequence] != iSet) {
1402 // get pi on gub row
1403 iSet = backward_[iSequence];
1404 if (iSet >= 0) {
1405 nSets++;
1406 int iBasic = keyVariable_[iSet];
1407 if (iBasic >= numberColumns) {
1408 djMod = - weight(iSet) * infeasibilityCost;
1409 } else {
1410 // get dj without
1411 assert (model->getStatus(iBasic) == ClpSimplex::basic);
1412 djMod = 0.0;
1413 // scaled
1414 for (j = startColumn[iBasic];
1415 j < startColumn[iBasic] + length[iBasic]; j++) {
1416 int jRow = row[j];
1417 djMod -= duals[jRow] * element[j] * rowScale[jRow];
1418 }
1419 // allow for scaling
1420 djMod += cost[iBasic] / columnScale[iBasic];
1421 // See if gub slack possible - dj is djMod
1422 if (getStatus(iSet) == ClpSimplex::atLowerBound) {
1423 double value = -djMod;
1424 if (value > tolerance) {
1425 numberWanted--;
1426 if (value > bestDj) {
1427 // check flagged variable and correct dj
1428 if (!flagged(iSet)) {
1429 bestDj = value;
1430 bestSequence = numberRows + numberColumns + iSet;
1431 bestDjMod = djMod;
1432 } else {
1433 // just to make sure we don't exit before got something
1434 numberWanted++;
1435 abort();
1436 }
1437 }
1438 }
1439 } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
1440 double value = djMod;
1441 if (value > tolerance) {
1442 numberWanted--;
1443 if (value > bestDj) {
1444 // check flagged variable and correct dj
1445 if (!flagged(iSet)) {
1446 bestDj = value;
1447 bestSequence = numberRows + numberColumns + iSet;
1448 bestDjMod = djMod;
1449 } else {
1450 // just to make sure we don't exit before got something
1451 numberWanted++;
1452 abort();
1453 }
1454 }
1455 }
1456 }
1457 }
1458 } else {
1459 // not in set
1460 djMod = 0.0;
1461 }
1462 }
1463 if (iSequence != sequenceOut) {
1464 double value;
1465 ClpSimplex::Status status = model->getStatus(iSequence);
1466
1467 switch(status) {
1468
1469 case ClpSimplex::basic:
1470 case ClpSimplex::isFixed:
1471 break;
1472 case ClpSimplex::isFree:
1473 case ClpSimplex::superBasic:
1474 value = -djMod;
1475 // scaled
1476 for (j = startColumn[iSequence];
1477 j < startColumn[iSequence] + length[iSequence]; j++) {
1478 int jRow = row[j];
1479 value -= duals[jRow] * element[j] * rowScale[jRow];
1480 }
1481 value = fabs(cost[iSequence] + value * columnScale[iSequence]);
1482 if (value > FREE_ACCEPT * tolerance) {
1483 numberWanted--;
1484 // we are going to bias towards free (but only if reasonable)
1485 value *= FREE_BIAS;
1486 if (value > bestDj) {
1487 // check flagged variable and correct dj
1488 if (!model->flagged(iSequence)) {
1489 bestDj = value;
1490 bestSequence = iSequence;
1491 bestDjMod = djMod;
1492 } else {
1493 // just to make sure we don't exit before got something
1494 numberWanted++;
1495 }
1496 }
1497 }
1498 break;
1499 case ClpSimplex::atUpperBound:
1500 value = -djMod;
1501 // scaled
1502 for (j = startColumn[iSequence];
1503 j < startColumn[iSequence] + length[iSequence]; j++) {
1504 int jRow = row[j];
1505 value -= duals[jRow] * element[j] * rowScale[jRow];
1506 }
1507 value = cost[iSequence] + value * columnScale[iSequence];
1508 if (value > tolerance) {
1509 numberWanted--;
1510 if (value > bestDj) {
1511 // check flagged variable and correct dj
1512 if (!model->flagged(iSequence)) {
1513 bestDj = value;
1514 bestSequence = iSequence;
1515 bestDjMod = djMod;
1516 } else {
1517 // just to make sure we don't exit before got something
1518 numberWanted++;
1519 }
1520 }
1521 }
1522 break;
1523 case ClpSimplex::atLowerBound:
1524 value = -djMod;
1525 // scaled
1526 for (j = startColumn[iSequence];
1527 j < startColumn[iSequence] + length[iSequence]; j++) {
1528 int jRow = row[j];
1529 value -= duals[jRow] * element[j] * rowScale[jRow];
1530 }
1531 value = -(cost[iSequence] + value * columnScale[iSequence]);
1532 if (value > tolerance) {
1533 numberWanted--;
1534 if (value > bestDj) {
1535 // check flagged variable and correct dj
1536 if (!model->flagged(iSequence)) {
1537 bestDj = value;
1538 bestSequence = iSequence;
1539 bestDjMod = djMod;
1540 } else {
1541 // just to make sure we don't exit before got something
1542 numberWanted++;
1543 }
1544 }
1545 }
1546 break;
1547 }
1548 }
1549 if (!numberWanted)
1550 break;
1551 }
1552 if (bestSequence != saveSequence) {
1553 if (bestSequence < numberRows + numberColumns) {
1554 // recompute dj
1555 double value = bestDjMod;
1556 // scaled
1557 for (j = startColumn[bestSequence];
1558 j < startColumn[bestSequence] + length[bestSequence]; j++) {
1559 int jRow = row[j];
1560 value -= duals[jRow] * element[j] * rowScale[jRow];
1561 }
1562 reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
1563 gubSlackIn_ = -1;
1564 } else {
1565 // slack - make last column
1566 gubSlackIn_ = bestSequence - numberRows - numberColumns;
1567 bestSequence = numberColumns + 2 * numberRows;
1568 reducedCost[bestSequence] = bestDjMod;
1569 model->setStatus(bestSequence, getStatus(gubSlackIn_));
1570 if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
1571 model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
1572 else
1573 model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
1574 model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
1575 model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
1576 model->costRegion()[bestSequence] = 0.0;
1577 }
1578 savedBestSequence_ = bestSequence;
1579 savedBestDj_ = reducedCost[savedBestSequence_];
1580 }
1581 } else {
1582 double bestDjMod = 0.0;
1583 //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
1584 // startG,endG,numberWanted);
1585 for (iSequence = startG; iSequence < endG; iSequence++) {
1586 if (numberWanted + minNeg < originalWanted_ && nSets > minSet) {
1587 // give up
1588 numberWanted = 0;
1589 break;
1590 } else if (iSequence == endG && bestSequence >= 0) {
1591 break;
1592 }
1593 if (backward_[iSequence] != iSet) {
1594 // get pi on gub row
1595 iSet = backward_[iSequence];
1596 if (iSet >= 0) {
1597 nSets++;
1598 int iBasic = keyVariable_[iSet];
1599 if (iBasic >= numberColumns) {
1600 djMod = - weight(iSet) * infeasibilityCost;
1601 } else {
1602 // get dj without
1603 assert (model->getStatus(iBasic) == ClpSimplex::basic);
1604 djMod = 0.0;
1605
1606 for (j = startColumn[iBasic];
1607 j < startColumn[iBasic] + length[iBasic]; j++) {
1608 int jRow = row[j];
1609 djMod -= duals[jRow] * element[j];
1610 }
1611 djMod += cost[iBasic];
1612 // See if gub slack possible - dj is djMod
1613 if (getStatus(iSet) == ClpSimplex::atLowerBound) {
1614 double value = -djMod;
1615 if (value > tolerance) {
1616 numberWanted--;
1617 if (value > bestDj) {
1618 // check flagged variable and correct dj
1619 if (!flagged(iSet)) {
1620 bestDj = value;
1621 bestSequence = numberRows + numberColumns + iSet;
1622 bestDjMod = djMod;
1623 } else {
1624 // just to make sure we don't exit before got something
1625 numberWanted++;
1626 abort();
1627 }
1628 }
1629 }
1630 } else if (getStatus(iSet) == ClpSimplex::atUpperBound) {
1631 double value = djMod;
1632 if (value > tolerance) {
1633 numberWanted--;
1634 if (value > bestDj) {
1635 // check flagged variable and correct dj
1636 if (!flagged(iSet)) {
1637 bestDj = value;
1638 bestSequence = numberRows + numberColumns + iSet;
1639 bestDjMod = djMod;
1640 } else {
1641 // just to make sure we don't exit before got something
1642 numberWanted++;
1643 abort();
1644 }
1645 }
1646 }
1647 }
1648 }
1649 } else {
1650 // not in set
1651 djMod = 0.0;
1652 }
1653 }
1654 if (iSequence != sequenceOut) {
1655 double value;
1656 ClpSimplex::Status status = model->getStatus(iSequence);
1657
1658 switch(status) {
1659
1660 case ClpSimplex::basic:
1661 case ClpSimplex::isFixed:
1662 break;
1663 case ClpSimplex::isFree:
1664 case ClpSimplex::superBasic:
1665 value = cost[iSequence] - djMod;
1666 for (j = startColumn[iSequence];
1667 j < startColumn[iSequence] + length[iSequence]; j++) {
1668 int jRow = row[j];
1669 value -= duals[jRow] * element[j];
1670 }
1671 value = fabs(value);
1672 if (value > FREE_ACCEPT * tolerance) {
1673 numberWanted--;
1674 // we are going to bias towards free (but only if reasonable)
1675 value *= FREE_BIAS;
1676 if (value > bestDj) {
1677 // check flagged variable and correct dj
1678 if (!model->flagged(iSequence)) {
1679 bestDj = value;
1680 bestSequence = iSequence;
1681 bestDjMod = djMod;
1682 } else {
1683 // just to make sure we don't exit before got something
1684 numberWanted++;
1685 }
1686 }
1687 }
1688 break;
1689 case ClpSimplex::atUpperBound:
1690 value = cost[iSequence] - djMod;
1691 for (j = startColumn[iSequence];
1692 j < startColumn[iSequence] + length[iSequence]; j++) {
1693 int jRow = row[j];
1694 value -= duals[jRow] * element[j];
1695 }
1696 if (value > tolerance) {
1697 numberWanted--;
1698 if (value > bestDj) {
1699 // check flagged variable and correct dj
1700 if (!model->flagged(iSequence)) {
1701 bestDj = value;
1702 bestSequence = iSequence;
1703 bestDjMod = djMod;
1704 } else {
1705 // just to make sure we don't exit before got something
1706 numberWanted++;
1707 }
1708 }
1709 }
1710 break;
1711 case ClpSimplex::atLowerBound:
1712 value = cost[iSequence] - djMod;
1713 for (j = startColumn[iSequence];
1714 j < startColumn[iSequence] + length[iSequence]; j++) {
1715 int jRow = row[j];
1716 value -= duals[jRow] * element[j];
1717 }
1718 value = -value;
1719 if (value > tolerance) {
1720 numberWanted--;
1721 if (value > bestDj) {
1722 // check flagged variable and correct dj
1723 if (!model->flagged(iSequence)) {
1724 bestDj = value;
1725 bestSequence = iSequence;
1726 bestDjMod = djMod;
1727 } else {
1728 // just to make sure we don't exit before got something
1729 numberWanted++;
1730 }
1731 }
1732 }
1733 break;
1734 }
1735 }
1736 if (!numberWanted)
1737 break;
1738 }
1739 if (bestSequence != saveSequence) {
1740 if (bestSequence < numberRows + numberColumns) {
1741 // recompute dj
1742 double value = cost[bestSequence] - bestDjMod;
1743 for (j = startColumn[bestSequence];
1744 j < startColumn[bestSequence] + length[bestSequence]; j++) {
1745 int jRow = row[j];
1746 value -= duals[jRow] * element[j];
1747 }
1748 //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
1749 reducedCost[bestSequence] = value;
1750 gubSlackIn_ = -1;
1751 } else {
1752 // slack - make last column
1753 gubSlackIn_ = bestSequence - numberRows - numberColumns;
1754 bestSequence = numberColumns + 2 * numberRows;
1755 reducedCost[bestSequence] = bestDjMod;
1756 //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
1757 model->setStatus(bestSequence, getStatus(gubSlackIn_));
1758 if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound)
1759 model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
1760 else
1761 model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
1762 model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
1763 model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
1764 model->costRegion()[bestSequence] = 0.0;
1765 }
1766 }
1767 }
1768 // See if may be finished
1769 if (startG == firstGub_ && bestSequence < 0)
1770 infeasibilityWeight_ = model_->infeasibilityCost();
1771 else if (bestSequence >= 0)
1772 infeasibilityWeight_ = -1.0;
1773 }
1774 if (numberWanted) {
1775 // Do packed part after gub
1776 double offset = static_cast<double> (lastGub_) /
1777 static_cast<double> (numberColumns);
1778 double ratio = static_cast<double> (numberColumns) /
1779 static_cast<double> (numberColumns) - offset;
1780 double start2 = offset + ratio * startFraction;
1781 double end2 = CoinMin(1.0, offset + ratio * endFraction + 1.0e-6);
1782 ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted);
1783 }
1784 } else {
1785 // no gub
1786 ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted);
1787 }
1788 if (bestSequence >= 0)
1789 infeasibilityWeight_ = -1.0; // not optimal
1790 currentWanted_ = numberWanted;
1791}
1792/* expands an updated column to allow for extra rows which the main
1793 solver does not know about and returns number added.
1794*/
1795int
1796ClpGubMatrix::extendUpdated(ClpSimplex * model, CoinIndexedVector * update, int mode)
1797{
1798 // I think we only need to bother about sets with two in basis or incoming set
1799 int number = update->getNumElements();
1800 double * array = update->denseVector();
1801 int * index = update->getIndices();
1802 int i;
1803 assert (!number || update->packedMode());
1804 int * pivotVariable = model->pivotVariable();
1805 int numberRows = model->numberRows();
1806 int numberColumns = model->numberColumns();
1807 int numberTotal = numberRows + numberColumns;
1808 int sequenceIn = model->sequenceIn();
1809 int returnCode = 0;
1810 int iSetIn;
1811 if (sequenceIn < numberColumns) {
1812 iSetIn = backward_[sequenceIn];
1813 gubSlackIn_ = -1; // in case set
1814 } else if (sequenceIn < numberRows + numberColumns) {
1815 iSetIn = -1;
1816 gubSlackIn_ = -1; // in case set
1817 } else {
1818 iSetIn = gubSlackIn_;
1819 }
1820 double * lower = model->lowerRegion();
1821 double * upper = model->upperRegion();
1822 double * cost = model->costRegion();
1823 double * solution = model->solutionRegion();
1824 int number2 = number;
1825 if (!mode) {
1826 double primalTolerance = model->primalTolerance();
1827 double infeasibilityCost = model->infeasibilityCost();
1828 // extend
1829 saveNumber_ = number;
1830 for (i = 0; i < number; i++) {
1831 int iRow = index[i];
1832 int iPivot = pivotVariable[iRow];
1833 if (iPivot < numberColumns) {
1834 int iSet = backward_[iPivot];
1835 if (iSet >= 0) {
1836 // two (or more) in set
1837 int iIndex = toIndex_[iSet];
1838 double otherValue = array[i];
1839 double value;
1840 if (iIndex < 0) {
1841 toIndex_[iSet] = number2;
1842 int iNew = number2 - number;
1843 fromIndex_[number2-number] = iSet;
1844 iIndex = number2;
1845 index[number2] = numberRows + iNew;
1846 // do key stuff
1847 int iKey = keyVariable_[iSet];
1848 if (iKey < numberColumns) {
1849 // Save current cost of key
1850 changeCost_[number2-number] = cost[iKey];
1851 if (iSet != iSetIn)
1852 value = 0.0;
1853 else if (iSetIn != gubSlackIn_)
1854 value = 1.0;
1855 else
1856 value = -1.0;
1857 pivotVariable[numberRows+iNew] = iKey;
1858 // Do I need to recompute?
1859 double sol;
1860 assert (getStatus(iSet) != ClpSimplex::basic);
1861 if (getStatus(iSet) == ClpSimplex::atLowerBound)
1862 sol = lower_[iSet];
1863 else
1864 sol = upper_[iSet];
1865 if ((gubType_ & 8) != 0) {
1866 int iColumn = next_[iKey];
1867 // sum all non-key variables
1868 while(iColumn >= 0) {
1869 sol -= solution[iColumn];
1870 iColumn = next_[iColumn];
1871 }
1872 } else {
1873 int stop = -(iKey + 1);
1874 int iColumn = next_[iKey];
1875 // sum all non-key variables
1876 while(iColumn != stop) {
1877 if (iColumn < 0)
1878 iColumn = -iColumn - 1;
1879 sol -= solution[iColumn];
1880 iColumn = next_[iColumn];
1881 }
1882 }
1883 solution[iKey] = sol;
1884 if (model->algorithm() > 0)
1885 model->nonLinearCost()->setOne(iKey, sol);
1886 //assert (fabs(sol-solution[iKey])<1.0e-3);
1887 } else {
1888 // gub slack is basic
1889 // Save current cost of key
1890 changeCost_[number2-number] = -weight(iSet) * infeasibilityCost;
1891 otherValue = - otherValue; //allow for - sign on slack
1892 if (iSet != iSetIn)
1893 value = 0.0;
1894 else
1895 value = -1.0;
1896 pivotVariable[numberRows+iNew] = iNew + numberTotal;
1897 model->djRegion()[iNew+numberTotal] = 0.0;
1898 double sol = 0.0;
1899 if ((gubType_ & 8) != 0) {
1900 int iColumn = next_[iKey];
1901 // sum all non-key variables
1902 while(iColumn >= 0) {
1903 sol += solution[iColumn];
1904 iColumn = next_[iColumn];
1905 }
1906 } else {
1907 int stop = -(iKey + 1);
1908 int iColumn = next_[iKey];
1909 // sum all non-key variables
1910 while(iColumn != stop) {
1911 if (iColumn < 0)
1912 iColumn = -iColumn - 1;
1913 sol += solution[iColumn];
1914 iColumn = next_[iColumn];
1915 }
1916 }
1917 solution[iNew+numberTotal] = sol;
1918 // and do cost in nonLinearCost
1919 if (model->algorithm() > 0)
1920 model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]);
1921 if (sol > upper_[iSet] + primalTolerance) {
1922 setAbove(iSet);
1923 lower[iNew+numberTotal] = upper_[iSet];
1924 upper[iNew+numberTotal] = COIN_DBL_MAX;
1925 } else if (sol < lower_[iSet] - primalTolerance) {
1926 setBelow(iSet);
1927 lower[iNew+numberTotal] = -COIN_DBL_MAX;
1928 upper[iNew+numberTotal] = lower_[iSet];
1929 } else {
1930 setFeasible(iSet);
1931 lower[iNew+numberTotal] = lower_[iSet];
1932 upper[iNew+numberTotal] = upper_[iSet];
1933 }
1934 cost[iNew+numberTotal] = weight(iSet) * infeasibilityCost;
1935 }
1936 number2++;
1937 } else {
1938 value = array[iIndex];
1939 int iKey = keyVariable_[iSet];
1940 if (iKey >= numberColumns)
1941 otherValue = - otherValue; //allow for - sign on slack
1942 }
1943 value -= otherValue;
1944 array[iIndex] = value;
1945 }
1946 }
1947 }
1948 if (iSetIn >= 0 && toIndex_[iSetIn] < 0) {
1949 // Do incoming
1950 update->setPacked(); // just in case no elements
1951 toIndex_[iSetIn] = number2;
1952 int iNew = number2 - number;
1953 fromIndex_[number2-number] = iSetIn;
1954 // Save current cost of key
1955 double currentCost;
1956 int key = keyVariable_[iSetIn];
1957 if (key < numberColumns)
1958 currentCost = cost[key];
1959 else
1960 currentCost = -weight(iSetIn) * infeasibilityCost;
1961 changeCost_[number2-number] = currentCost;
1962 index[number2] = numberRows + iNew;
1963 // do key stuff
1964 int iKey = keyVariable_[iSetIn];
1965 if (iKey < numberColumns) {
1966 if (gubSlackIn_ < 0)
1967 array[number2] = 1.0;
1968 else
1969 array[number2] = -1.0;
1970 pivotVariable[numberRows+iNew] = iKey;
1971 // Do I need to recompute?
1972 double sol;
1973 assert (getStatus(iSetIn) != ClpSimplex::basic);
1974 if (getStatus(iSetIn) == ClpSimplex::atLowerBound)
1975 sol = lower_[iSetIn];
1976 else
1977 sol = upper_[iSetIn];
1978 if ((gubType_ & 8) != 0) {
1979 int iColumn = next_[iKey];
1980 // sum all non-key variables
1981 while(iColumn >= 0) {
1982 sol -= solution[iColumn];
1983 iColumn = next_[iColumn];
1984 }
1985 } else {
1986 // bounds exist - sum over all except key
1987 int stop = -(iKey + 1);
1988 int iColumn = next_[iKey];
1989 // sum all non-key variables
1990 while(iColumn != stop) {
1991 if (iColumn < 0)
1992 iColumn = -iColumn - 1;
1993 sol -= solution[iColumn];
1994 iColumn = next_[iColumn];
1995 }
1996 }
1997 solution[iKey] = sol;
1998 if (model->algorithm() > 0)
1999 model->nonLinearCost()->setOne(iKey, sol);
2000 //assert (fabs(sol-solution[iKey])<1.0e-3);
2001 } else {
2002 // gub slack is basic
2003 array[number2] = -1.0;
2004 pivotVariable[numberRows+iNew] = iNew + numberTotal;
2005 model->djRegion()[iNew+numberTotal] = 0.0;
2006 double sol = 0.0;
2007 if ((gubType_ & 8) != 0) {
2008 int iColumn = next_[iKey];
2009 // sum all non-key variables
2010 while(iColumn >= 0) {
2011 sol += solution[iColumn];
2012 iColumn = next_[iColumn];
2013 }
2014 } else {
2015 // bounds exist - sum over all except key
2016 int stop = -(iKey + 1);
2017 int iColumn = next_[iKey];
2018 // sum all non-key variables
2019 while(iColumn != stop) {
2020 if (iColumn < 0)
2021 iColumn = -iColumn - 1;
2022 sol += solution[iColumn];
2023 iColumn = next_[iColumn];
2024 }
2025 }
2026 solution[iNew+numberTotal] = sol;
2027 // and do cost in nonLinearCost
2028 if (model->algorithm() > 0)
2029 model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]);
2030 if (sol > upper_[iSetIn] + primalTolerance) {
2031 setAbove(iSetIn);
2032 lower[iNew+numberTotal] = upper_[iSetIn];
2033 upper[iNew+numberTotal] = COIN_DBL_MAX;
2034 } else if (sol < lower_[iSetIn] - primalTolerance) {
2035 setBelow(iSetIn);
2036 lower[iNew+numberTotal] = -COIN_DBL_MAX;
2037 upper[iNew+numberTotal] = lower_[iSetIn];
2038 } else {
2039 setFeasible(iSetIn);
2040 lower[iNew+numberTotal] = lower_[iSetIn];
2041 upper[iNew+numberTotal] = upper_[iSetIn];
2042 }
2043 cost[iNew+numberTotal] = weight(iSetIn) * infeasibilityCost;
2044 }
2045 number2++;
2046 }
2047 // mark end
2048 fromIndex_[number2-number] = -1;
2049 returnCode = number2 - number;
2050 // make sure lower_ upper_ adjusted
2051 synchronize(model, 9);
2052 } else {
2053 // take off?
2054 if (number > saveNumber_) {
2055 // clear
2056 double theta = model->theta();
2057 double * solution = model->solutionRegion();
2058 for (i = saveNumber_; i < number; i++) {
2059 int iRow = index[i];
2060 int iColumn = pivotVariable[iRow];
2061#ifdef CLP_DEBUG_PRINT
2062 printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
2063 iColumn, fromIndex_[i-saveNumber_], lower[iColumn], upper[iColumn], array[i],
2064 solution[iColumn], solution[iColumn] - model->theta()*array[i], model->theta());
2065#endif
2066 double value = array[i];
2067 array[i] = 0.0;
2068 int iSet = fromIndex_[i-saveNumber_];
2069 toIndex_[iSet] = -1;
2070 if (iSet == iSetIn && iColumn < numberColumns) {
2071 // update as may need value
2072 solution[iColumn] -= theta * value;
2073 }
2074 }
2075 }
2076#ifdef CLP_DEBUG
2077 for (i = 0; i < numberSets_; i++)
2078 assert(toIndex_[i] == -1);
2079#endif
2080 number2 = saveNumber_;
2081 }
2082 update->setNumElements(number2);
2083 return returnCode;
2084}
2085/*
2086 utility primal function for dealing with dynamic constraints
2087 mode=n see ClpGubMatrix.hpp for definition
2088 Remember to update here when settled down
2089*/
2090void
2091ClpGubMatrix::primalExpanded(ClpSimplex * model, int mode)
2092{
2093 int numberColumns = model->numberColumns();
2094 switch (mode) {
2095 // If key variable then slot in gub rhs so will get correct contribution
2096 case 0: {
2097 int i;
2098 double * solution = model->solutionRegion();
2099 ClpSimplex::Status iStatus;
2100 for (i = 0; i < numberSets_; i++) {
2101 int iColumn = keyVariable_[i];
2102 if (iColumn < numberColumns) {
2103 // key is structural - where is slack
2104 iStatus = getStatus(i);
2105 assert (iStatus != ClpSimplex::basic);
2106 if (iStatus == ClpSimplex::atLowerBound)
2107 solution[iColumn] = lower_[i];
2108 else
2109 solution[iColumn] = upper_[i];
2110 }
2111 }
2112 }
2113 break;
2114 // Compute values of key variables
2115 case 1: {
2116 int i;
2117 double * solution = model->solutionRegion();
2118 ClpSimplex::Status iStatus;
2119 //const int * columnLength = matrix_->getVectorLengths();
2120 //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2121 //const int * row = matrix_->getIndices();
2122 //const double * elementByColumn = matrix_->getElements();
2123 //int * pivotVariable = model->pivotVariable();
2124 sumPrimalInfeasibilities_ = 0.0;
2125 numberPrimalInfeasibilities_ = 0;
2126 double primalTolerance = model->primalTolerance();
2127 double relaxedTolerance = primalTolerance;
2128 // we can't really trust infeasibilities if there is primal error
2129 double error = CoinMin(1.0e-2, model->largestPrimalError());
2130 // allow tolerance at least slightly bigger than standard
2131 relaxedTolerance = relaxedTolerance + error;
2132 // but we will be using difference
2133 relaxedTolerance -= primalTolerance;
2134 sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2135 for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds)
2136 int kColumn = keyVariable_[i];
2137 double value = 0.0;
2138 if ((gubType_ & 8) != 0) {
2139 int iColumn = next_[kColumn];
2140 // sum all non-key variables
2141 while(iColumn >= 0) {
2142 value += solution[iColumn];
2143 iColumn = next_[iColumn];
2144 }
2145 } else {
2146 // bounds exist - sum over all except key
2147 int stop = -(kColumn + 1);
2148 int iColumn = next_[kColumn];
2149 // sum all non-key variables
2150 while(iColumn != stop) {
2151 if (iColumn < 0)
2152 iColumn = -iColumn - 1;
2153 value += solution[iColumn];
2154 iColumn = next_[iColumn];
2155 }
2156 }
2157 if (kColumn < numberColumns) {
2158 // make sure key is basic - so will be skipped in values pass
2159 model->setStatus(kColumn, ClpSimplex::basic);
2160 // feasibility will be done later
2161 assert (getStatus(i) != ClpSimplex::basic);
2162 if (getStatus(i) == ClpSimplex::atUpperBound)
2163 solution[kColumn] = upper_[i] - value;
2164 else
2165 solution[kColumn] = lower_[i] - value;
2166 //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
2167 } else {
2168 // slack is key
2169 iStatus = getStatus(i);
2170 assert (iStatus == ClpSimplex::basic);
2171 double infeasibility = 0.0;
2172 if (value > upper_[i] + primalTolerance) {
2173 infeasibility = value - upper_[i] - primalTolerance;
2174 setAbove(i);
2175 } else if (value < lower_[i] - primalTolerance) {
2176 infeasibility = lower_[i] - value - primalTolerance ;
2177 setBelow(i);
2178 } else {
2179 setFeasible(i);
2180 }
2181 //printf("Value of key slack for set %d is %g\n",i,value);
2182 if (infeasibility > 0.0) {
2183 sumPrimalInfeasibilities_ += infeasibility;
2184 if (infeasibility > relaxedTolerance)
2185 sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
2186 numberPrimalInfeasibilities_ ++;
2187 }
2188 }
2189 }
2190 }
2191 break;
2192 // Report on infeasibilities of key variables
2193 case 2: {
2194 model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() +
2195 sumPrimalInfeasibilities_);
2196 model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() +
2197 numberPrimalInfeasibilities_);
2198 model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() +
2199 sumOfRelaxedPrimalInfeasibilities_);
2200 }
2201 break;
2202 }
2203}
2204/*
2205 utility dual function for dealing with dynamic constraints
2206 mode=n see ClpGubMatrix.hpp for definition
2207 Remember to update here when settled down
2208*/
2209void
2210ClpGubMatrix::dualExpanded(ClpSimplex * model,
2211 CoinIndexedVector * array,
2212 double * /*other*/, int mode)
2213{
2214 switch (mode) {
2215 // modify costs before transposeUpdate
2216 case 0: {
2217 int i;
2218 double * cost = model->costRegion();
2219 ClpSimplex::Status iStatus;
2220 // not dual values yet
2221 //assert (!other);
2222 //double * work = array->denseVector();
2223 double infeasibilityCost = model->infeasibilityCost();
2224 int * pivotVariable = model->pivotVariable();
2225 int numberRows = model->numberRows();
2226 int numberColumns = model->numberColumns();
2227 for (i = 0; i < numberRows; i++) {
2228 int iPivot = pivotVariable[i];
2229 if (iPivot < numberColumns) {
2230 int iSet = backward_[iPivot];
2231 if (iSet >= 0) {
2232 int kColumn = keyVariable_[iSet];
2233 double costValue;
2234 if (kColumn < numberColumns) {
2235 // structural has cost
2236 costValue = cost[kColumn];
2237 } else {
2238 // slack is key
2239 iStatus = getStatus(iSet);
2240 assert (iStatus == ClpSimplex::basic);
2241 // negative as -1.0 for slack
2242 costValue = -weight(iSet) * infeasibilityCost;
2243 }
2244 array->add(i, -costValue); // was work[i]-costValue
2245 }
2246 }
2247 }
2248 }
2249 break;
2250 // create duals for key variables (without check on dual infeasible)
2251 case 1: {
2252 // If key slack then dual 0.0 (if feasible)
2253 // dj for key is zero so that defines dual on set
2254 int i;
2255 double * dj = model->djRegion();
2256 int numberColumns = model->numberColumns();
2257 double infeasibilityCost = model->infeasibilityCost();
2258 for (i = 0; i < numberSets_; i++) {
2259 int kColumn = keyVariable_[i];
2260 if (kColumn < numberColumns) {
2261 // dj without set
2262 double value = dj[kColumn];
2263 // Now subtract out from all
2264 dj[kColumn] = 0.0;
2265 int iColumn = next_[kColumn];
2266 // modify all non-key variables
2267 while(iColumn >= 0) {
2268 dj[iColumn] -= value;
2269 iColumn = next_[iColumn];
2270 }
2271 } else {
2272 // slack key - may not be feasible
2273 assert (getStatus(i) == ClpSimplex::basic);
2274 // negative as -1.0 for slack
2275 double value = -weight(i) * infeasibilityCost;
2276 if (value) {
2277 int iColumn = next_[kColumn];
2278 // modify all non-key variables basic
2279 while(iColumn >= 0) {
2280 dj[iColumn] -= value;
2281 iColumn = next_[iColumn];
2282 }
2283 }
2284 }
2285 }
2286 }
2287 break;
2288 // as 1 but check slacks and compute djs
2289 case 2: {
2290 // If key slack then dual 0.0
2291 // If not then slack could be dual infeasible
2292 // dj for key is zero so that defines dual on set
2293 int i;
2294 // make sure fromIndex will not confuse pricing
2295 fromIndex_[0] = -1;
2296 possiblePivotKey_ = -1;
2297 // Create array
2298 int numberColumns = model->numberColumns();
2299 int * pivotVariable = model->pivotVariable();
2300 int numberRows = model->numberRows();
2301 for (i = 0; i < numberRows; i++) {
2302 int iPivot = pivotVariable[i];
2303 if (iPivot < numberColumns)
2304 backToPivotRow_[iPivot] = i;
2305 }
2306 if (noCheck_ >= 0) {
2307 if (infeasibilityWeight_ != model->infeasibilityCost()) {
2308 // don't bother checking
2309 sumDualInfeasibilities_ = 100.0;
2310 numberDualInfeasibilities_ = 1;
2311 sumOfRelaxedDualInfeasibilities_ = 100.0;
2312 return;
2313 }
2314 }
2315 double * dj = model->djRegion();
2316 double * dual = model->dualRowSolution();
2317 double * cost = model->costRegion();
2318 ClpSimplex::Status iStatus;
2319 const int * columnLength = matrix_->getVectorLengths();
2320 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2321 const int * row = matrix_->getIndices();
2322 const double * elementByColumn = matrix_->getElements();
2323 double infeasibilityCost = model->infeasibilityCost();
2324 sumDualInfeasibilities_ = 0.0;
2325 numberDualInfeasibilities_ = 0;
2326 double dualTolerance = model->dualTolerance();
2327 double relaxedTolerance = dualTolerance;
2328 // we can't really trust infeasibilities if there is dual error
2329 double error = CoinMin(1.0e-2, model->largestDualError());
2330 // allow tolerance at least slightly bigger than standard
2331 relaxedTolerance = relaxedTolerance + error;
2332 // but we will be using difference
2333 relaxedTolerance -= dualTolerance;
2334 sumOfRelaxedDualInfeasibilities_ = 0.0;
2335 for (i = 0; i < numberSets_; i++) {
2336 int kColumn = keyVariable_[i];
2337 if (kColumn < numberColumns) {
2338 // dj without set
2339 double value = cost[kColumn];
2340 for (CoinBigIndex j = columnStart[kColumn];
2341 j < columnStart[kColumn] + columnLength[kColumn]; j++) {
2342 int iRow = row[j];
2343 value -= dual[iRow] * elementByColumn[j];
2344 }
2345 // Now subtract out from all
2346 dj[kColumn] -= value;
2347 int stop = -(kColumn + 1);
2348 kColumn = next_[kColumn];
2349 while (kColumn != stop) {
2350 if (kColumn < 0)
2351 kColumn = -kColumn - 1;
2352 double djValue = dj[kColumn] - value;
2353 dj[kColumn] = djValue;
2354 double infeasibility = 0.0;
2355 iStatus = model->getStatus(kColumn);
2356 if (iStatus == ClpSimplex::atLowerBound) {
2357 if (djValue < -dualTolerance)
2358 infeasibility = -djValue - dualTolerance;
2359 } else if (iStatus == ClpSimplex::atUpperBound) {
2360 // at upper bound
2361 if (djValue > dualTolerance)
2362 infeasibility = djValue - dualTolerance;
2363 }
2364 if (infeasibility > 0.0) {
2365 sumDualInfeasibilities_ += infeasibility;
2366 if (infeasibility > relaxedTolerance)
2367 sumOfRelaxedDualInfeasibilities_ += infeasibility;
2368 numberDualInfeasibilities_ ++;
2369 }
2370 kColumn = next_[kColumn];
2371 }
2372 // check slack
2373 iStatus = getStatus(i);
2374 assert (iStatus != ClpSimplex::basic);
2375 double infeasibility = 0.0;
2376 // dj of slack is -(-1.0)value
2377 if (iStatus == ClpSimplex::atLowerBound) {
2378 if (value < -dualTolerance)
2379 infeasibility = -value - dualTolerance;
2380 } else if (iStatus == ClpSimplex::atUpperBound) {
2381 // at upper bound
2382 if (value > dualTolerance)
2383 infeasibility = value - dualTolerance;
2384 }
2385 if (infeasibility > 0.0) {
2386 sumDualInfeasibilities_ += infeasibility;
2387 if (infeasibility > relaxedTolerance)
2388 sumOfRelaxedDualInfeasibilities_ += infeasibility;
2389 numberDualInfeasibilities_ ++;
2390 }
2391 } else {
2392 // slack key - may not be feasible
2393 assert (getStatus(i) == ClpSimplex::basic);
2394 // negative as -1.0 for slack
2395 double value = -weight(i) * infeasibilityCost;
2396 if (value) {
2397 // Now subtract out from all
2398 int kColumn = i + numberColumns;
2399 int stop = -(kColumn + 1);
2400 kColumn = next_[kColumn];
2401 while (kColumn != stop) {
2402 if (kColumn < 0)
2403 kColumn = -kColumn - 1;
2404 double djValue = dj[kColumn] - value;
2405 dj[kColumn] = djValue;
2406 double infeasibility = 0.0;
2407 iStatus = model->getStatus(kColumn);
2408 if (iStatus == ClpSimplex::atLowerBound) {
2409 if (djValue < -dualTolerance)
2410 infeasibility = -djValue - dualTolerance;
2411 } else if (iStatus == ClpSimplex::atUpperBound) {
2412 // at upper bound
2413 if (djValue > dualTolerance)
2414 infeasibility = djValue - dualTolerance;
2415 }
2416 if (infeasibility > 0.0) {
2417 sumDualInfeasibilities_ += infeasibility;
2418 if (infeasibility > relaxedTolerance)
2419 sumOfRelaxedDualInfeasibilities_ += infeasibility;
2420 numberDualInfeasibilities_ ++;
2421 }
2422 kColumn = next_[kColumn];
2423 }
2424 }
2425 }
2426 }
2427 // and get statistics for column generation
2428 synchronize(model, 4);
2429 infeasibilityWeight_ = -1.0;
2430 }
2431 break;
2432 // Report on infeasibilities of key variables
2433 case 3: {
2434 model->setSumDualInfeasibilities(model->sumDualInfeasibilities() +
2435 sumDualInfeasibilities_);
2436 model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() +
2437 numberDualInfeasibilities_);
2438 model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() +
2439 sumOfRelaxedDualInfeasibilities_);
2440 }
2441 break;
2442 // modify costs before transposeUpdate for partial pricing
2443 case 4: {
2444 // First compute new costs etc for interesting gubs
2445 int iLook = 0;
2446 int iSet = fromIndex_[0];
2447 double primalTolerance = model->primalTolerance();
2448 const double * cost = model->costRegion();
2449 double * solution = model->solutionRegion();
2450 double infeasibilityCost = model->infeasibilityCost();
2451 int numberColumns = model->numberColumns();
2452 int numberChanged = 0;
2453 int * pivotVariable = model->pivotVariable();
2454 while (iSet >= 0) {
2455 int key = keyVariable_[iSet];
2456 double value = 0.0;
2457 // sum over all except key
2458 if ((gubType_ & 8) != 0) {
2459 int iColumn = next_[key];
2460 // sum all non-key variables
2461 while(iColumn >= 0) {
2462 value += solution[iColumn];
2463 iColumn = next_[iColumn];
2464 }
2465 } else {
2466 // bounds exist - sum over all except key
2467 int stop = -(key + 1);
2468 int iColumn = next_[key];
2469 // sum all non-key variables
2470 while(iColumn != stop) {
2471 if (iColumn < 0)
2472 iColumn = -iColumn - 1;
2473 value += solution[iColumn];
2474 iColumn = next_[iColumn];
2475 }
2476 }
2477 double costChange;
2478 double oldCost = changeCost_[iLook];
2479 if (key < numberColumns) {
2480 assert (getStatus(iSet) != ClpSimplex::basic);
2481 double sol;
2482 if (getStatus(iSet) == ClpSimplex::atUpperBound)
2483 sol = upper_[iSet] - value;
2484 else
2485 sol = lower_[iSet] - value;
2486 solution[key] = sol;
2487 // fix up cost
2488 model->nonLinearCost()->setOne(key, sol);
2489#ifdef CLP_DEBUG_PRINT
2490 printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol,
2491 cost[key], oldCost);
2492#endif
2493 costChange = cost[key] - oldCost;
2494 } else {
2495 // slack is key
2496 if (value > upper_[iSet] + primalTolerance) {
2497 setAbove(iSet);
2498 } else if (value < lower_[iSet] - primalTolerance) {
2499 setBelow(iSet);
2500 } else {
2501 setFeasible(iSet);
2502 }
2503 // negative as -1.0 for slack
2504 costChange = -weight(iSet) * infeasibilityCost - oldCost;
2505#ifdef CLP_DEBUG_PRINT
2506 printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value,
2507 weight(iSet)*infeasibilityCost, oldCost);
2508#endif
2509 }
2510 if (costChange) {
2511 fromIndex_[numberChanged] = iSet;
2512 toIndex_[iSet] = numberChanged;
2513 changeCost_[numberChanged++] = costChange;
2514 }
2515 iSet = fromIndex_[++iLook];
2516 }
2517 if (numberChanged || possiblePivotKey_ >= 0) {
2518 // first do those in list already
2519 int number = array->getNumElements();
2520 array->setPacked();
2521 int i;
2522 double * work = array->denseVector();
2523 int * which = array->getIndices();
2524 for (i = 0; i < number; i++) {
2525 int iRow = which[i];
2526 int iPivot = pivotVariable[iRow];
2527 if (iPivot < numberColumns) {
2528 int iSet = backward_[iPivot];
2529 if (iSet >= 0 && toIndex_[iSet] >= 0) {
2530 double newValue = work[i] + changeCost_[toIndex_[iSet]];
2531 if (!newValue)
2532 newValue = 1.0e-100;
2533 work[i] = newValue;
2534 // mark as done
2535 backward_[iPivot] = -1;
2536 }
2537 }
2538 if (possiblePivotKey_ == iRow) {
2539 double newValue = work[i] - model->dualIn();
2540 if (!newValue)
2541 newValue = 1.0e-100;
2542 work[i] = newValue;
2543 possiblePivotKey_ = -1;
2544 }
2545 }
2546 // now do rest and clean up
2547 for (i = 0; i < numberChanged; i++) {
2548 int iSet = fromIndex_[i];
2549 int key = keyVariable_[iSet];
2550 int iColumn = next_[key];
2551 double change = changeCost_[i];
2552 while (iColumn >= 0) {
2553 if (backward_[iColumn] >= 0) {
2554 int iRow = backToPivotRow_[iColumn];
2555 assert (iRow >= 0);
2556 work[number] = change;
2557 if (possiblePivotKey_ == iRow) {
2558 double newValue = work[number] - model->dualIn();
2559 if (!newValue)
2560 newValue = 1.0e-100;
2561 work[number] = newValue;
2562 possiblePivotKey_ = -1;
2563 }
2564 which[number++] = iRow;
2565 } else {
2566 // reset
2567 backward_[iColumn] = iSet;
2568 }
2569 iColumn = next_[iColumn];
2570 }
2571 toIndex_[iSet] = -1;
2572 }
2573 if (possiblePivotKey_ >= 0) {
2574 work[number] = -model->dualIn();
2575 which[number++] = possiblePivotKey_;
2576 possiblePivotKey_ = -1;
2577 }
2578 fromIndex_[0] = -1;
2579 array->setNumElements(number);
2580 }
2581 }
2582 break;
2583 }
2584}
2585// This is local to Gub to allow synchronization when status is good
2586int
2587ClpGubMatrix::synchronize(ClpSimplex *, int)
2588{
2589 return 0;
2590}
2591/*
2592 general utility function for dealing with dynamic constraints
2593 mode=n see ClpGubMatrix.hpp for definition
2594 Remember to update here when settled down
2595*/
2596int
2597ClpGubMatrix::generalExpanded(ClpSimplex * model, int mode, int &number)
2598{
2599 int returnCode = 0;
2600 int numberColumns = model->numberColumns();
2601 switch (mode) {
2602 // Fill in pivotVariable but not for key variables
2603 case 0: {
2604 if (!next_ ) {
2605 // do ordering
2606 assert (!rhsOffset_);
2607 // create and do gub crash
2608 useEffectiveRhs(model, false);
2609 }
2610 int i;
2611 int numberBasic = number;
2612 // Use different array so can build from true pivotVariable_
2613 //int * pivotVariable = model->pivotVariable();
2614 int * pivotVariable = model->rowArray(0)->getIndices();
2615 for (i = 0; i < numberColumns; i++) {
2616 if (model->getColumnStatus(i) == ClpSimplex::basic) {
2617 int iSet = backward_[i];
2618 if (iSet < 0 || i != keyVariable_[iSet])
2619 pivotVariable[numberBasic++] = i;
2620 }
2621 }
2622 number = numberBasic;
2623 if (model->numberIterations())
2624 assert (number == model->numberRows());
2625 }
2626 break;
2627 // Make all key variables basic
2628 case 1: {
2629 int i;
2630 for (i = 0; i < numberSets_; i++) {
2631 int iColumn = keyVariable_[i];
2632 if (iColumn < numberColumns)
2633 model->setColumnStatus(iColumn, ClpSimplex::basic);
2634 }
2635 }
2636 break;
2637 // Do initial extra rows + maximum basic
2638 case 2: {
2639 returnCode = getNumRows() + 1;
2640 number = model->numberRows() + numberSets_;
2641 }
2642 break;
2643 // Before normal replaceColumn
2644 case 3: {
2645 int sequenceIn = model->sequenceIn();
2646 int sequenceOut = model->sequenceOut();
2647 int numberColumns = model->numberColumns();
2648 int numberRows = model->numberRows();
2649 int pivotRow = model->pivotRow();
2650 if (gubSlackIn_ >= 0)
2651 assert (sequenceIn > numberRows + numberColumns);
2652 if (sequenceIn == sequenceOut)
2653 return -1;
2654 int iSetIn = -1;
2655 int iSetOut = -1;
2656 if (sequenceOut < numberColumns) {
2657 iSetOut = backward_[sequenceOut];
2658 } else if (sequenceOut >= numberRows + numberColumns) {
2659 assert (pivotRow >= numberRows);
2660 int iExtra = pivotRow - numberRows;
2661 assert (iExtra >= 0);
2662 if (iSetOut < 0)
2663 iSetOut = fromIndex_[iExtra];
2664 else
2665 assert(iSetOut == fromIndex_[iExtra]);
2666 }
2667 if (sequenceIn < numberColumns) {
2668 iSetIn = backward_[sequenceIn];
2669 } else if (gubSlackIn_ >= 0) {
2670 iSetIn = gubSlackIn_;
2671 }
2672 possiblePivotKey_ = -1;
2673 number = 0; // say do ordinary
2674 int * pivotVariable = model->pivotVariable();
2675 if (pivotRow >= numberRows) {
2676 int iExtra = pivotRow - numberRows;
2677 //const int * length = matrix_->getVectorLengths();
2678
2679 assert (sequenceOut >= numberRows + numberColumns ||
2680 sequenceOut == keyVariable_[iSetOut]);
2681 int incomingColumn = sequenceIn; // to be used in updates
2682 if (iSetIn != iSetOut) {
2683 // We need to find a possible pivot for incoming
2684 // look through rowArray_[1]
2685 int n = model->rowArray(1)->getNumElements();
2686 int * which = model->rowArray(1)->getIndices();
2687 double * array = model->rowArray(1)->denseVector();
2688 double bestAlpha = 1.0e-5;
2689 //int shortest=numberRows+1;
2690 for (int i = 0; i < n; i++) {
2691 int iRow = which[i];
2692 int iPivot = pivotVariable[iRow];
2693 if (iPivot < numberColumns && backward_[iPivot] == iSetOut) {
2694 if (fabs(array[i]) > fabs(bestAlpha)) {
2695 bestAlpha = array[i];
2696 possiblePivotKey_ = iRow;
2697 }
2698 }
2699 }
2700 assert (possiblePivotKey_ >= 0); // could set returnCode=4
2701 number = 1;
2702 if (sequenceIn >= numberRows + numberColumns) {
2703 number = 3;
2704 // need swap as gub slack in and must become key
2705 // is this best way
2706 int key = keyVariable_[iSetIn];
2707 assert (key < numberColumns);
2708 // check other basic
2709 int iColumn = next_[key];
2710 // set new key to be used by unpack
2711 keyVariable_[iSetIn] = iSetIn + numberColumns;
2712 // change cost in changeCost
2713 {
2714 int iLook = 0;
2715 int iSet = fromIndex_[0];
2716 while (iSet >= 0) {
2717 if (iSet == iSetIn) {
2718 changeCost_[iLook] = 0.0;
2719 break;
2720 }
2721 iSet = fromIndex_[++iLook];
2722 }
2723 }
2724 while (iColumn >= 0) {
2725 if (iColumn != sequenceOut) {
2726 // need partial ftran and skip accuracy check in replaceColumn
2727#ifdef CLP_DEBUG_PRINT
2728 printf("TTTTTry 5\n");
2729#endif
2730 int iRow = backToPivotRow_[iColumn];
2731 assert (iRow >= 0);
2732 unpack(model, model->rowArray(3), iColumn);
2733 model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2734 double alpha = model->rowArray(3)->denseVector()[iRow];
2735 //if (!alpha)
2736 //printf("zero alpha a\n");
2737 int updateStatus = model->factorization()->replaceColumn(model,
2738 model->rowArray(2),
2739 model->rowArray(3),
2740 iRow, alpha);
2741 returnCode = CoinMax(updateStatus, returnCode);
2742 model->rowArray(3)->clear();
2743 if (returnCode)
2744 break;
2745 }
2746 iColumn = next_[iColumn];
2747 }
2748 if (!returnCode) {
2749 // now factorization looks as if key is out
2750 // pivot back in
2751#ifdef CLP_DEBUG_PRINT
2752 printf("TTTTTry 6\n");
2753#endif
2754 unpack(model, model->rowArray(3), key);
2755 model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2756 pivotRow = possiblePivotKey_;
2757 double alpha = model->rowArray(3)->denseVector()[pivotRow];
2758 //if (!alpha)
2759 //printf("zero alpha b\n");
2760 int updateStatus = model->factorization()->replaceColumn(model,
2761 model->rowArray(2),
2762 model->rowArray(3),
2763 pivotRow, alpha);
2764 returnCode = CoinMax(updateStatus, returnCode);
2765 model->rowArray(3)->clear();
2766 }
2767 // restore key
2768 keyVariable_[iSetIn] = key;
2769 // now alternate column can replace key on out
2770 incomingColumn = pivotVariable[possiblePivotKey_];
2771 } else {
2772#ifdef CLP_DEBUG_PRINT
2773 printf("TTTTTTry 4 %d\n", possiblePivotKey_);
2774#endif
2775 int updateStatus = model->factorization()->replaceColumn(model,
2776 model->rowArray(2),
2777 model->rowArray(1),
2778 possiblePivotKey_,
2779 bestAlpha);
2780 returnCode = CoinMax(updateStatus, returnCode);
2781 incomingColumn = pivotVariable[possiblePivotKey_];
2782 }
2783
2784 //returnCode=4; // need swap
2785 } else {
2786 // key swap
2787 number = -1;
2788 }
2789 int key = keyVariable_[iSetOut];
2790 if (key < numberColumns)
2791 assert(key == sequenceOut);
2792 // check if any other basic
2793 int iColumn = next_[key];
2794 if (returnCode)
2795 iColumn = -1; // skip if error on previous
2796 // set new key to be used by unpack
2797 if (incomingColumn < numberColumns)
2798 keyVariable_[iSetOut] = incomingColumn;
2799 else
2800 keyVariable_[iSetOut] = iSetIn + numberColumns;
2801 double * cost = model->costRegion();
2802 if (possiblePivotKey_ < 0) {
2803 double dj = model->djRegion()[sequenceIn] - cost[sequenceIn];
2804 changeCost_[iExtra] = -dj;
2805#ifdef CLP_DEBUG_PRINT
2806 printf("modifying changeCost %d by %g - cost %g\n", iExtra, dj, cost[sequenceIn]);
2807#endif
2808 }
2809 while (iColumn >= 0) {
2810 if (iColumn != incomingColumn) {
2811 number = -2;
2812 // need partial ftran and skip accuracy check in replaceColumn
2813#ifdef CLP_DEBUG_PRINT
2814 printf("TTTTTTry 1\n");
2815#endif
2816 int iRow = backToPivotRow_[iColumn];
2817 assert (iRow >= 0 && iRow < numberRows);
2818 unpack(model, model->rowArray(3), iColumn);
2819 model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2820 double * array = model->rowArray(3)->denseVector();
2821 double alpha = array[iRow];
2822 //if (!alpha)
2823 //printf("zero alpha d\n");
2824 int updateStatus = model->factorization()->replaceColumn(model,
2825 model->rowArray(2),
2826 model->rowArray(3),
2827 iRow, alpha);
2828 returnCode = CoinMax(updateStatus, returnCode);
2829 model->rowArray(3)->clear();
2830 if (returnCode)
2831 break;
2832 }
2833 iColumn = next_[iColumn];
2834 }
2835 // restore key
2836 keyVariable_[iSetOut] = key;
2837 } else if (sequenceIn >= numberRows + numberColumns) {
2838 number = 2;
2839 //returnCode=4;
2840 // need swap as gub slack in and must become key
2841 // is this best way
2842 int key = keyVariable_[iSetIn];
2843 assert (key < numberColumns);
2844 // check other basic
2845 int iColumn = next_[key];
2846 // set new key to be used by unpack
2847 keyVariable_[iSetIn] = iSetIn + numberColumns;
2848 // change cost in changeCost
2849 {
2850 int iLook = 0;
2851 int iSet = fromIndex_[0];
2852 while (iSet >= 0) {
2853 if (iSet == iSetIn) {
2854 changeCost_[iLook] = 0.0;
2855 break;
2856 }
2857 iSet = fromIndex_[++iLook];
2858 }
2859 }
2860 while (iColumn >= 0) {
2861 if (iColumn != sequenceOut) {
2862 // need partial ftran and skip accuracy check in replaceColumn
2863#ifdef CLP_DEBUG_PRINT
2864 printf("TTTTTry 2\n");
2865#endif
2866 int iRow = backToPivotRow_[iColumn];
2867 assert (iRow >= 0);
2868 unpack(model, model->rowArray(3), iColumn);
2869 model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2870 double alpha = model->rowArray(3)->denseVector()[iRow];
2871 //if (!alpha)
2872 //printf("zero alpha e\n");
2873 int updateStatus = model->factorization()->replaceColumn(model,
2874 model->rowArray(2),
2875 model->rowArray(3),
2876 iRow, alpha);
2877 returnCode = CoinMax(updateStatus, returnCode);
2878 model->rowArray(3)->clear();
2879 if (returnCode)
2880 break;
2881 }
2882 iColumn = next_[iColumn];
2883 }
2884 if (!returnCode) {
2885 // now factorization looks as if key is out
2886 // pivot back in
2887#ifdef CLP_DEBUG_PRINT
2888 printf("TTTTTry 3\n");
2889#endif
2890 unpack(model, model->rowArray(3), key);
2891 model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3));
2892 double alpha = model->rowArray(3)->denseVector()[pivotRow];
2893 //if (!alpha)
2894 //printf("zero alpha f\n");
2895 int updateStatus = model->factorization()->replaceColumn(model,
2896 model->rowArray(2),
2897 model->rowArray(3),
2898 pivotRow, alpha);
2899 returnCode = CoinMax(updateStatus, returnCode);
2900 model->rowArray(3)->clear();
2901 }
2902 // restore key
2903 keyVariable_[iSetIn] = key;
2904 } else {
2905 // normal - but might as well do here
2906 returnCode = model->factorization()->replaceColumn(model,
2907 model->rowArray(2),
2908 model->rowArray(1),
2909 model->pivotRow(),
2910 model->alpha());
2911 }
2912 }
2913#ifdef CLP_DEBUG_PRINT
2914 printf("Update type after %d - status %d - pivot row %d\n",
2915 number, returnCode, model->pivotRow());
2916#endif
2917 // see if column generation says time to re-factorize
2918 returnCode = CoinMax(returnCode, synchronize(model, 5));
2919 number = -1; // say no need for normal replaceColumn
2920 break;
2921 // To see if can dual or primal
2922 case 4: {
2923 returnCode = 1;
2924 }
2925 break;
2926 // save status
2927 case 5: {
2928 synchronize(model, 0);
2929 CoinMemcpyN(status_, numberSets_, saveStatus_);
2930 CoinMemcpyN(keyVariable_, numberSets_, savedKeyVariable_);
2931 }
2932 break;
2933 // restore status
2934 case 6: {
2935 CoinMemcpyN(saveStatus_, numberSets_, status_);
2936 CoinMemcpyN(savedKeyVariable_, numberSets_, keyVariable_);
2937 // restore firstAvailable_
2938 synchronize(model, 7);
2939 // redo next_
2940 int i;
2941 int * last = new int[numberSets_];
2942 for (i = 0; i < numberSets_; i++) {
2943 int iKey = keyVariable_[i];
2944 assert(iKey >= numberColumns || backward_[iKey] == i);
2945 last[i] = iKey;
2946 // make sure basic
2947 //if (iKey<numberColumns)
2948 //model->setStatus(iKey,ClpSimplex::basic);
2949 }
2950 for (i = 0; i < numberColumns; i++) {
2951 int iSet = backward_[i];
2952 if (iSet >= 0) {
2953 next_[last[iSet]] = i;
2954 last[iSet] = i;
2955 }
2956 }
2957 for (i = 0; i < numberSets_; i++) {
2958 next_[last[i]] = -(keyVariable_[i] + 1);
2959 redoSet(model, keyVariable_[i], keyVariable_[i], i);
2960 }
2961 delete [] last;
2962 // redo pivotVariable
2963 int * pivotVariable = model->pivotVariable();
2964 int iRow;
2965 int numberBasic = 0;
2966 int numberRows = model->numberRows();
2967 for (iRow = 0; iRow < numberRows; iRow++) {
2968 if (model->getRowStatus(iRow) == ClpSimplex::basic) {
2969 numberBasic++;
2970 pivotVariable[iRow] = iRow + numberColumns;
2971 } else {
2972 pivotVariable[iRow] = -1;
2973 }
2974 }
2975 i = 0;
2976 int iColumn;
2977 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2978 if (model->getStatus(iColumn) == ClpSimplex::basic) {
2979 int iSet = backward_[iColumn];
2980 if (iSet < 0 || keyVariable_[iSet] != iColumn) {
2981 while (pivotVariable[i] >= 0) {
2982 i++;
2983 assert (i < numberRows);
2984 }
2985 pivotVariable[i] = iColumn;
2986 backToPivotRow_[iColumn] = i;
2987 numberBasic++;
2988 }
2989 }
2990 }
2991 assert (numberBasic == numberRows);
2992 rhsOffset(model, true);
2993 }
2994 break;
2995 // flag a variable
2996 case 7: {
2997 assert (number == model->sequenceIn());
2998 synchronize(model, 1);
2999 synchronize(model, 8);
3000 }
3001 break;
3002 // unflag all variables
3003 case 8: {
3004 returnCode = synchronize(model, 2);
3005 }
3006 break;
3007 // redo costs in primal
3008 case 9: {
3009 returnCode = synchronize(model, 3);
3010 }
3011 break;
3012 // return 1 if there may be changing bounds on variable (column generation)
3013 case 10: {
3014 returnCode = synchronize(model, 6);
3015 }
3016 break;
3017 // make sure set is clean
3018 case 11: {
3019 assert (number == model->sequenceIn());
3020 returnCode = synchronize(model, 8);
3021 }
3022 break;
3023 default:
3024 break;
3025 }
3026 return returnCode;
3027}
3028// Sets up an effective RHS and does gub crash if needed
3029void
3030ClpGubMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
3031{
3032 // Do basis - cheapest or slack if feasible (unless cheapest set)
3033 int longestSet = 0;
3034 int iSet;
3035 for (iSet = 0; iSet < numberSets_; iSet++)
3036 longestSet = CoinMax(longestSet, end_[iSet] - start_[iSet]);
3037
3038 double * upper = new double[longestSet+1];
3039 double * cost = new double[longestSet+1];
3040 double * lower = new double[longestSet+1];
3041 double * solution = new double[longestSet+1];
3042 assert (!next_);
3043 int numberColumns = getNumCols();
3044 const int * columnLength = matrix_->getVectorLengths();
3045 const double * columnLower = model->lowerRegion();
3046 const double * columnUpper = model->upperRegion();
3047 double * columnSolution = model->solutionRegion();
3048 const double * objective = model->costRegion();
3049 int numberRows = getNumRows();
3050 toIndex_ = new int[numberSets_];
3051 for (iSet = 0; iSet < numberSets_; iSet++)
3052 toIndex_[iSet] = -1;
3053 fromIndex_ = new int [getNumRows()+1];
3054 double tolerance = model->primalTolerance();
3055 bool noNormalBounds = true;
3056 gubType_ &= ~8;
3057 bool gotBasis = false;
3058 for (iSet = 0; iSet < numberSets_; iSet++) {
3059 if (keyVariable_[iSet] < numberColumns)
3060 gotBasis = true;
3061 CoinBigIndex j;
3062 CoinBigIndex iStart = start_[iSet];
3063 CoinBigIndex iEnd = end_[iSet];
3064 for (j = iStart; j < iEnd; j++) {
3065 if (columnLower[j] && columnLower[j] > -1.0e20)
3066 noNormalBounds = false;
3067 if (columnUpper[j] && columnUpper[j] < 1.0e20)
3068 noNormalBounds = false;
3069 }
3070 }
3071 if (noNormalBounds)
3072 gubType_ |= 8;
3073 if (!gotBasis) {
3074 for (iSet = 0; iSet < numberSets_; iSet++) {
3075 CoinBigIndex j;
3076 int numberBasic = 0;
3077 int iBasic = -1;
3078 CoinBigIndex iStart = start_[iSet];
3079 CoinBigIndex iEnd = end_[iSet];
3080 // find one with smallest length
3081 int smallest = numberRows + 1;
3082 double value = 0.0;
3083 for (j = iStart; j < iEnd; j++) {
3084 if (model->getStatus(j) == ClpSimplex::basic) {
3085 if (columnLength[j] < smallest) {
3086 smallest = columnLength[j];
3087 iBasic = j;
3088 }
3089 numberBasic++;
3090 }
3091 value += columnSolution[j];
3092 }
3093 bool done = false;
3094 if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) {
3095 if (getStatus(iSet) == ClpSimplex::basic)
3096 iBasic = iSet + numberColumns; // slack key - use
3097 done = true;
3098 } else if (numberBasic == 1) {
3099 // see if can be key
3100 double thisSolution = columnSolution[iBasic];
3101 if (thisSolution > columnUpper[iBasic]) {
3102 value -= thisSolution - columnUpper[iBasic];
3103 thisSolution = columnUpper[iBasic];
3104 columnSolution[iBasic] = thisSolution;
3105 }
3106 if (thisSolution < columnLower[iBasic]) {
3107 value -= thisSolution - columnLower[iBasic];
3108 thisSolution = columnLower[iBasic];
3109 columnSolution[iBasic] = thisSolution;
3110 }
3111 // try setting slack to a bound
3112 assert (upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20);
3113 double cost1 = COIN_DBL_MAX;
3114 int whichBound = -1;
3115 if (upper_[iSet] < 1.0e20) {
3116 // try slack at ub
3117 double newBasic = thisSolution + upper_[iSet] - value;
3118 if (newBasic >= columnLower[iBasic] - tolerance &&
3119 newBasic <= columnUpper[iBasic] + tolerance) {
3120 // can go
3121 whichBound = 1;
3122 cost1 = newBasic * objective[iBasic];
3123 // But if exact then may be good solution
3124 if (fabs(upper_[iSet] - value) < tolerance)
3125 cost1 = -COIN_DBL_MAX;
3126 }
3127 }
3128 if (lower_[iSet] > -1.0e20) {
3129 // try slack at lb
3130 double newBasic = thisSolution + lower_[iSet] - value;
3131 if (newBasic >= columnLower[iBasic] - tolerance &&
3132 newBasic <= columnUpper[iBasic] + tolerance) {
3133 // can go but is it cheaper
3134 double cost2 = newBasic * objective[iBasic];
3135 // But if exact then may be good solution
3136 if (fabs(lower_[iSet] - value) < tolerance)
3137 cost2 = -COIN_DBL_MAX;
3138 if (cost2 < cost1)
3139 whichBound = 0;
3140 }
3141 }
3142 if (whichBound != -1) {
3143 // key
3144 done = true;
3145 if (whichBound) {
3146 // slack to upper
3147 columnSolution[iBasic] = thisSolution + upper_[iSet] - value;
3148 setStatus(iSet, ClpSimplex::atUpperBound);
3149 } else {
3150 // slack to lower
3151 columnSolution[iBasic] = thisSolution + lower_[iSet] - value;
3152 setStatus(iSet, ClpSimplex::atLowerBound);
3153 }
3154 }
3155 }
3156 if (!done) {
3157 if (!cheapest) {
3158 // see if slack can be key
3159 if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
3160 done = true;
3161 setStatus(iSet, ClpSimplex::basic);
3162 iBasic = iSet + numberColumns;
3163 }
3164 }
3165 if (!done) {
3166 // set non basic if there was one
3167 if (iBasic >= 0)
3168 model->setStatus(iBasic, ClpSimplex::atLowerBound);
3169 // find cheapest
3170 int numberInSet = iEnd - iStart;
3171 CoinMemcpyN(columnLower + iStart, numberInSet, lower);
3172 CoinMemcpyN(columnUpper + iStart, numberInSet, upper);
3173 CoinMemcpyN(columnSolution + iStart, numberInSet, solution);
3174 // and slack
3175 iBasic = numberInSet;
3176 solution[iBasic] = -value;
3177 lower[iBasic] = -upper_[iSet];
3178 upper[iBasic] = -lower_[iSet];
3179 int kphase;
3180 if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) {
3181 // feasible
3182 kphase = 1;
3183 cost[iBasic] = 0.0;
3184 CoinMemcpyN(objective + iStart, numberInSet, cost);
3185 } else {
3186 // infeasible
3187 kphase = 0;
3188 // remember bounds are flipped so opposite to natural
3189 if (value < lower_[iSet] - tolerance)
3190 cost[iBasic] = 1.0;
3191 else
3192 cost[iBasic] = -1.0;
3193 CoinZeroN(cost, numberInSet);
3194 }
3195 double dualTolerance = model->dualTolerance();
3196 for (int iphase = kphase; iphase < 2; iphase++) {
3197 if (iphase) {
3198 cost[numberInSet] = 0.0;
3199 CoinMemcpyN(objective + iStart, numberInSet, cost);
3200 }
3201 // now do one row lp
3202 bool improve = true;
3203 while (improve) {
3204 improve = false;
3205 double dual = cost[iBasic];
3206 int chosen = -1;
3207 double best = dualTolerance;
3208 int way = 0;
3209 for (int i = 0; i <= numberInSet; i++) {
3210 double dj = cost[i] - dual;
3211 double improvement = 0.0;
3212 if (iphase || i < numberInSet)
3213 assert (solution[i] >= lower[i] && solution[i] <= upper[i]);
3214 if (dj > dualTolerance)
3215 improvement = dj * (solution[i] - lower[i]);
3216 else if (dj < -dualTolerance)
3217 improvement = dj * (solution[i] - upper[i]);
3218 if (improvement > best) {
3219 best = improvement;
3220 chosen = i;
3221 if (dj < 0.0) {
3222 way = 1;
3223 } else {
3224 way = -1;
3225 }
3226 }
3227 }
3228 if (chosen >= 0) {
3229 improve = true;
3230 // now see how far
3231 if (way > 0) {
3232 // incoming increasing so basic decreasing
3233 // if phase 0 then go to nearest bound
3234 double distance = upper[chosen] - solution[chosen];
3235 double basicDistance;
3236 if (!iphase) {
3237 assert (iBasic == numberInSet);
3238 assert (solution[iBasic] > upper[iBasic]);
3239 basicDistance = solution[iBasic] - upper[iBasic];
3240 } else {
3241 basicDistance = solution[iBasic] - lower[iBasic];
3242 }
3243 // need extra coding for unbounded
3244 assert (CoinMin(distance, basicDistance) < 1.0e20);
3245 if (distance > basicDistance) {
3246 // incoming becomes basic
3247 solution[chosen] += basicDistance;
3248 if (!iphase)
3249 solution[iBasic] = upper[iBasic];
3250 else
3251 solution[iBasic] = lower[iBasic];
3252 iBasic = chosen;
3253 } else {
3254 // flip
3255 solution[chosen] = upper[chosen];
3256 solution[iBasic] -= distance;
3257 }
3258 } else {
3259 // incoming decreasing so basic increasing
3260 // if phase 0 then go to nearest bound
3261 double distance = solution[chosen] - lower[chosen];
3262 double basicDistance;
3263 if (!iphase) {
3264 assert (iBasic == numberInSet);
3265 assert (solution[iBasic] < lower[iBasic]);
3266 basicDistance = lower[iBasic] - solution[iBasic];
3267 } else {
3268 basicDistance = upper[iBasic] - solution[iBasic];
3269 }
3270 // need extra coding for unbounded - for now just exit
3271 if (CoinMin(distance, basicDistance) > 1.0e20) {
3272 printf("unbounded on set %d\n", iSet);
3273 iphase = 1;
3274 iBasic = numberInSet;
3275 break;
3276 }
3277 if (distance > basicDistance) {
3278 // incoming becomes basic
3279 solution[chosen] -= basicDistance;
3280 if (!iphase)
3281 solution[iBasic] = lower[iBasic];
3282 else
3283 solution[iBasic] = upper[iBasic];
3284 iBasic = chosen;
3285 } else {
3286 // flip
3287 solution[chosen] = lower[chosen];
3288 solution[iBasic] += distance;
3289 }
3290 }
3291 if (!iphase) {
3292 if(iBasic < numberInSet)
3293 break; // feasible
3294 else if (solution[iBasic] >= lower[iBasic] &&
3295 solution[iBasic] <= upper[iBasic])
3296 break; // feasible (on flip)
3297 }
3298 }
3299 }
3300 }
3301 // convert iBasic back and do bounds
3302 if (iBasic == numberInSet) {
3303 // slack basic
3304 setStatus(iSet, ClpSimplex::basic);
3305 iBasic = iSet + numberColumns;
3306 } else {
3307 iBasic += start_[iSet];
3308 model->setStatus(iBasic, ClpSimplex::basic);
3309 // remember bounds flipped
3310 if (upper[numberInSet] == lower[numberInSet])
3311 setStatus(iSet, ClpSimplex::isFixed);
3312 else if (solution[numberInSet] == upper[numberInSet])
3313 setStatus(iSet, ClpSimplex::atLowerBound);
3314 else if (solution[numberInSet] == lower[numberInSet])
3315 setStatus(iSet, ClpSimplex::atUpperBound);
3316 else
3317 abort();
3318 }
3319 for (j = iStart; j < iEnd; j++) {
3320 if (model->getStatus(j) != ClpSimplex::basic) {
3321 int inSet = j - iStart;
3322 columnSolution[j] = solution[inSet];
3323 if (upper[inSet] == lower[inSet])
3324 model->setStatus(j, ClpSimplex::isFixed);
3325 else if (solution[inSet] == upper[inSet])
3326 model->setStatus(j, ClpSimplex::atUpperBound);
3327 else if (solution[inSet] == lower[inSet])
3328 model->setStatus(j, ClpSimplex::atLowerBound);
3329 }
3330 }
3331 }
3332 }
3333 keyVariable_[iSet] = iBasic;
3334 }
3335 }
3336 delete [] lower;
3337 delete [] solution;
3338 delete [] upper;
3339 delete [] cost;
3340 // make sure matrix is in good shape
3341 matrix_->orderMatrix();
3342 // create effective rhs
3343 delete [] rhsOffset_;
3344 rhsOffset_ = new double[numberRows];
3345 delete [] next_;
3346 next_ = new int[numberColumns+numberSets_+2*longestSet];
3347 char * mark = new char[numberColumns];
3348 memset(mark, 0, numberColumns);
3349 for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3350 next_[iColumn] = COIN_INT_MAX;
3351 int i;
3352 int * keys = new int[numberSets_];
3353 for (i = 0; i < numberSets_; i++)
3354 keys[i] = COIN_INT_MAX;
3355 // set up chains
3356 for (i = 0; i < numberColumns; i++) {
3357 if (model->getStatus(i) == ClpSimplex::basic)
3358 mark[i] = 1;
3359 int iSet = backward_[i];
3360 if (iSet >= 0) {
3361 int iNext = keys[iSet];
3362 next_[i] = iNext;
3363 keys[iSet] = i;
3364 }
3365 }
3366 for (i = 0; i < numberSets_; i++) {
3367 int j;
3368 if (getStatus(i) != ClpSimplex::basic) {
3369 // make sure fixed if it is
3370 if (upper_[i] == lower_[i])
3371 setStatus(i, ClpSimplex::isFixed);
3372 // slack not key - choose one with smallest length
3373 int smallest = numberRows + 1;
3374 int key = -1;
3375 j = keys[i];
3376 if (j != COIN_INT_MAX) {
3377 while (1) {
3378 if (mark[j] && columnLength[j] < smallest && !gotBasis) {
3379 key = j;
3380 smallest = columnLength[j];
3381 }
3382 if (next_[j] != COIN_INT_MAX) {
3383 j = next_[j];
3384 } else {
3385 // correct end
3386 next_[j] = -(keys[i] + 1);
3387 break;
3388 }
3389 }
3390 } else {
3391 next_[i+numberColumns] = -(numberColumns + i + 1);
3392 }
3393 if (gotBasis)
3394 key = keyVariable_[i];
3395 if (key >= 0) {
3396 keyVariable_[i] = key;
3397 } else {
3398 // nothing basic - make slack key
3399 //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
3400 // fudge to avoid const problem
3401 status_[i] = 1;
3402 }
3403 } else {
3404 // slack key
3405 keyVariable_[i] = numberColumns + i;
3406 int j;
3407 double sol = 0.0;
3408 j = keys[i];
3409 if (j != COIN_INT_MAX) {
3410 while (1) {
3411 sol += columnSolution[j];
3412 if (next_[j] != COIN_INT_MAX) {
3413 j = next_[j];
3414 } else {
3415 // correct end
3416 next_[j] = -(keys[i] + 1);
3417 break;
3418 }
3419 }
3420 } else {
3421 next_[i+numberColumns] = -(numberColumns + i + 1);
3422 }
3423 if (sol > upper_[i] + tolerance) {
3424 setAbove(i);
3425 } else if (sol < lower_[i] - tolerance) {
3426 setBelow(i);
3427 } else {
3428 setFeasible(i);
3429 }
3430 }
3431 // Create next_
3432 int key = keyVariable_[i];
3433 redoSet(model, key, keys[i], i);
3434 }
3435 delete [] keys;
3436 delete [] mark;
3437 rhsOffset(model, true);
3438}
3439// redoes next_ for a set.
3440void
3441ClpGubMatrix::redoSet(ClpSimplex * model, int newKey, int oldKey, int iSet)
3442{
3443 int numberColumns = model->numberColumns();
3444 int * save = next_ + numberColumns + numberSets_;
3445 int number = 0;
3446 int stop = -(oldKey + 1);
3447 int j = next_[oldKey];
3448 while (j != stop) {
3449 if (j < 0)
3450 j = -j - 1;
3451 if (j != newKey)
3452 save[number++] = j;
3453 j = next_[j];
3454 }
3455 // and add oldkey
3456 if (newKey != oldKey)
3457 save[number++] = oldKey;
3458 // now do basic
3459 int lastMarker = -(newKey + 1);
3460 keyVariable_[iSet] = newKey;
3461 next_[newKey] = lastMarker;
3462 int last = newKey;
3463 for ( j = 0; j < number; j++) {
3464 int iColumn = save[j];
3465 if (iColumn < numberColumns) {
3466 if (model->getStatus(iColumn) == ClpSimplex::basic) {
3467 next_[last] = iColumn;
3468 next_[iColumn] = lastMarker;
3469 last = iColumn;
3470 }
3471 }
3472 }
3473 // now add in non-basic
3474 for ( j = 0; j < number; j++) {
3475 int iColumn = save[j];
3476 if (iColumn < numberColumns) {
3477 if (model->getStatus(iColumn) != ClpSimplex::basic) {
3478 next_[last] = -(iColumn + 1);
3479 next_[iColumn] = lastMarker;
3480 last = iColumn;
3481 }
3482 }
3483 }
3484
3485}
3486/* Returns effective RHS if it is being used. This is used for long problems
3487 or big gub or anywhere where going through full columns is
3488 expensive. This may re-compute */
3489double *
3490ClpGubMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh, bool
3491#ifdef CLP_DEBUG
3492 check
3493#endif
3494 )
3495{
3496 //forceRefresh=true;
3497 if (rhsOffset_) {
3498#ifdef CLP_DEBUG
3499 if (check) {
3500 // no need - but check anyway
3501 // zero out basic
3502 int numberRows = model->numberRows();
3503 int numberColumns = model->numberColumns();
3504 double * solution = new double [numberColumns];
3505 double * rhs = new double[numberRows];
3506 CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
3507 CoinZeroN(rhs, numberRows);
3508 int iRow;
3509 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3510 if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
3511 solution[iColumn] = 0.0;
3512 }
3513 for (int iSet = 0; iSet < numberSets_; iSet++) {
3514 int iColumn = keyVariable_[iSet];
3515 if (iColumn < numberColumns)
3516 solution[iColumn] = 0.0;
3517 }
3518 times(-1.0, solution, rhs);
3519 delete [] solution;
3520 const double * columnSolution = model->solutionRegion();
3521 // and now subtract out non basic
3522 ClpSimplex::Status iStatus;
3523 for (int iSet = 0; iSet < numberSets_; iSet++) {
3524 int iColumn = keyVariable_[iSet];
3525 if (iColumn < numberColumns) {
3526 double b = 0.0;
3527 // key is structural - where is slack
3528 iStatus = getStatus(iSet);
3529 assert (iStatus != ClpSimplex::basic);
3530 if (iStatus == ClpSimplex::atLowerBound)
3531 b = lower_[iSet];
3532 else
3533 b = upper_[iSet];
3534 // subtract out others at bounds
3535 if ((gubType_ & 8) == 0) {
3536 int stop = -(iColumn + 1);
3537 int jColumn = next_[iColumn];
3538 // sum all non-basic variables - first skip basic
3539 while(jColumn >= 0)
3540 jColumn = next_[jColumn];
3541 while(jColumn != stop) {
3542 assert (jColumn < 0);
3543 jColumn = -jColumn - 1;
3544 b -= columnSolution[jColumn];
3545 jColumn = next_[jColumn];
3546 }
3547 }
3548 // subtract out
3549 ClpPackedMatrix::add(model, rhs, iColumn, -b);
3550 }
3551 }
3552 for (iRow = 0; iRow < numberRows; iRow++) {
3553 if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3)
3554 printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]);
3555 }
3556 delete [] rhs;
3557 }
3558#endif
3559 if (forceRefresh || (refreshFrequency_ && model->numberIterations() >=
3560 lastRefresh_ + refreshFrequency_)) {
3561 // zero out basic
3562 int numberRows = model->numberRows();
3563 int numberColumns = model->numberColumns();
3564 double * solution = new double [numberColumns];
3565 CoinMemcpyN(model->solutionRegion(), numberColumns, solution);
3566 CoinZeroN(rhsOffset_, numberRows);
3567 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
3568 if (model->getColumnStatus(iColumn) == ClpSimplex::basic)
3569 solution[iColumn] = 0.0;
3570 }
3571 int iSet;
3572 for ( iSet = 0; iSet < numberSets_; iSet++) {
3573 int iColumn = keyVariable_[iSet];
3574 if (iColumn < numberColumns)
3575 solution[iColumn] = 0.0;
3576 }
3577 times(-1.0, solution, rhsOffset_);
3578 delete [] solution;
3579 lastRefresh_ = model->numberIterations();
3580 const double * columnSolution = model->solutionRegion();
3581 // and now subtract out non basic
3582 ClpSimplex::Status iStatus;
3583 for ( iSet = 0; iSet < numberSets_; iSet++) {
3584 int iColumn = keyVariable_[iSet];
3585 if (iColumn < numberColumns) {
3586 double b = 0.0;
3587 // key is structural - where is slack
3588 iStatus = getStatus(iSet);
3589 assert (iStatus != ClpSimplex::basic);
3590 if (iStatus == ClpSimplex::atLowerBound)
3591 b = lower_[iSet];
3592 else
3593 b = upper_[iSet];
3594 // subtract out others at bounds
3595 if ((gubType_ & 8) == 0) {
3596 int stop = -(iColumn + 1);
3597 int jColumn = next_[iColumn];
3598 // sum all non-basic variables - first skip basic
3599 while(jColumn >= 0)
3600 jColumn = next_[jColumn];
3601 while(jColumn != stop) {
3602 assert (jColumn < 0);
3603 jColumn = -jColumn - 1;
3604 b -= columnSolution[jColumn];
3605 jColumn = next_[jColumn];
3606 }
3607 }
3608 // subtract out
3609 if (b)
3610 ClpPackedMatrix::add(model, rhsOffset_, iColumn, -b);
3611 }
3612 }
3613 }
3614 }
3615 return rhsOffset_;
3616}
3617/*
3618 update information for a pivot (and effective rhs)
3619*/
3620int
3621ClpGubMatrix::updatePivot(ClpSimplex * model, double oldInValue, double /*oldOutValue*/)
3622{
3623 int sequenceIn = model->sequenceIn();
3624 int sequenceOut = model->sequenceOut();
3625 double * solution = model->solutionRegion();
3626 int numberColumns = model->numberColumns();
3627 int numberRows = model->numberRows();
3628 int pivotRow = model->pivotRow();
3629 int iSetIn;
3630 // Correct sequence in
3631 trueSequenceIn_ = sequenceIn;
3632 if (sequenceIn < numberColumns) {
3633 iSetIn = backward_[sequenceIn];
3634 } else if (sequenceIn < numberColumns + numberRows) {
3635 iSetIn = -1;
3636 } else {
3637 iSetIn = gubSlackIn_;
3638 trueSequenceIn_ = numberColumns + numberRows + iSetIn;
3639 }
3640 int iSetOut = -1;
3641 trueSequenceOut_ = sequenceOut;
3642 if (sequenceOut < numberColumns) {
3643 iSetOut = backward_[sequenceOut];
3644 } else if (sequenceOut >= numberRows + numberColumns) {
3645 assert (pivotRow >= numberRows);
3646 int iExtra = pivotRow - numberRows;
3647 assert (iExtra >= 0);
3648 if (iSetOut < 0)
3649 iSetOut = fromIndex_[iExtra];
3650 else
3651 assert(iSetOut == fromIndex_[iExtra]);
3652 trueSequenceOut_ = numberColumns + numberRows + iSetOut;
3653 }
3654 if (rhsOffset_) {
3655 // update effective rhs
3656 if (sequenceIn == sequenceOut) {
3657 assert (sequenceIn < numberRows + numberColumns); // should be easy to deal with
3658 if (sequenceIn < numberColumns)
3659 add(model, rhsOffset_, sequenceIn, oldInValue - solution[sequenceIn]);
3660 } else {
3661 if (sequenceIn < numberColumns) {
3662 // we need to test if WILL be key
3663 ClpPackedMatrix::add(model, rhsOffset_, sequenceIn, oldInValue);
3664 if (iSetIn >= 0) {
3665 // old contribution to rhsOffset_
3666 int key = keyVariable_[iSetIn];
3667 if (key < numberColumns) {
3668 double oldB = 0.0;
3669 ClpSimplex::Status iStatus = getStatus(iSetIn);
3670 if (iStatus == ClpSimplex::atLowerBound)
3671 oldB = lower_[iSetIn];
3672 else
3673 oldB = upper_[iSetIn];
3674 // subtract out others at bounds
3675 if ((gubType_ & 8) == 0) {
3676 int stop = -(key + 1);
3677 int iColumn = next_[key];
3678 // skip basic
3679 while (iColumn >= 0)
3680 iColumn = next_[iColumn];
3681 // sum all non-key variables
3682 while(iColumn != stop) {
3683 assert (iColumn < 0);
3684 iColumn = -iColumn - 1;
3685 if (iColumn == sequenceIn)
3686 oldB -= oldInValue;
3687 else if ( iColumn != sequenceOut )
3688 oldB -= solution[iColumn];
3689 iColumn = next_[iColumn];
3690 }
3691 }
3692 if (oldB)
3693 ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3694 }
3695 }
3696 } else if (sequenceIn < numberRows + numberColumns) {
3697 //rhsOffset_[sequenceIn-numberColumns] -= oldInValue;
3698 } else {
3699#ifdef CLP_DEBUG_PRINT
3700 printf("** in is key slack %d\n", sequenceIn);
3701#endif
3702 // old contribution to rhsOffset_
3703 int key = keyVariable_[iSetIn];
3704 if (key < numberColumns) {
3705 double oldB = 0.0;
3706 ClpSimplex::Status iStatus = getStatus(iSetIn);
3707 if (iStatus == ClpSimplex::atLowerBound)
3708 oldB = lower_[iSetIn];
3709 else
3710 oldB = upper_[iSetIn];
3711 // subtract out others at bounds
3712 if ((gubType_ & 8) == 0) {
3713 int stop = -(key + 1);
3714 int iColumn = next_[key];
3715 // skip basic
3716 while (iColumn >= 0)
3717 iColumn = next_[iColumn];
3718 // sum all non-key variables
3719 while(iColumn != stop) {
3720 assert (iColumn < 0);
3721 iColumn = -iColumn - 1;
3722 if ( iColumn != sequenceOut )
3723 oldB -= solution[iColumn];
3724 iColumn = next_[iColumn];
3725 }
3726 }
3727 if (oldB)
3728 ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3729 }
3730 }
3731 if (sequenceOut < numberColumns) {
3732 ClpPackedMatrix::add(model, rhsOffset_, sequenceOut, -solution[sequenceOut]);
3733 if (iSetOut >= 0) {
3734 // old contribution to rhsOffset_
3735 int key = keyVariable_[iSetOut];
3736 if (key < numberColumns && iSetIn != iSetOut) {
3737 double oldB = 0.0;
3738 ClpSimplex::Status iStatus = getStatus(iSetOut);
3739 if (iStatus == ClpSimplex::atLowerBound)
3740 oldB = lower_[iSetOut];
3741 else
3742 oldB = upper_[iSetOut];
3743 // subtract out others at bounds
3744 if ((gubType_ & 8) == 0) {
3745 int stop = -(key + 1);
3746 int iColumn = next_[key];
3747 // skip basic
3748 while (iColumn >= 0)
3749 iColumn = next_[iColumn];
3750 // sum all non-key variables
3751 while(iColumn != stop) {
3752 assert (iColumn < 0);
3753 iColumn = -iColumn - 1;
3754 if (iColumn == sequenceIn)
3755 oldB -= oldInValue;
3756 else if ( iColumn != sequenceOut )
3757 oldB -= solution[iColumn];
3758 iColumn = next_[iColumn];
3759 }
3760 }
3761 if (oldB)
3762 ClpPackedMatrix::add(model, rhsOffset_, key, oldB);
3763 }
3764 }
3765 } else if (sequenceOut < numberRows + numberColumns) {
3766 //rhsOffset_[sequenceOut-numberColumns] -= -solution[sequenceOut];
3767 } else {
3768#ifdef CLP_DEBUG_PRINT
3769 printf("** out is key slack %d\n", sequenceOut);
3770#endif
3771 assert (pivotRow >= numberRows);
3772 }
3773 }
3774 }
3775 int * pivotVariable = model->pivotVariable();
3776 // may need to deal with key
3777 // Also need coding to mark/allow key slack entering
3778 if (pivotRow >= numberRows) {
3779 assert (sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]);
3780#ifdef CLP_DEBUG_PRINT
3781 if (sequenceIn >= numberRows + numberColumns)
3782 printf("key slack %d in, set out %d\n", gubSlackIn_, iSetOut);
3783 printf("** danger - key out for set %d in %d (set %d)\n", iSetOut, sequenceIn,
3784 iSetIn);
3785#endif
3786 // if slack out mark correctly
3787 if (sequenceOut >= numberRows + numberColumns) {
3788 double value = model->valueOut();
3789 if (value == upper_[iSetOut]) {
3790 setStatus(iSetOut, ClpSimplex::atUpperBound);
3791 } else if (value == lower_[iSetOut]) {
3792 setStatus(iSetOut, ClpSimplex::atLowerBound);
3793 } else {
3794 if (fabs(value - upper_[iSetOut]) <
3795 fabs(value - lower_[iSetOut])) {
3796 setStatus(iSetOut, ClpSimplex::atUpperBound);
3797 } else {
3798 setStatus(iSetOut, ClpSimplex::atLowerBound);
3799 }
3800 }
3801 if (upper_[iSetOut] == lower_[iSetOut])
3802 setStatus(iSetOut, ClpSimplex::isFixed);
3803 setFeasible(iSetOut);
3804 }
3805 if (iSetOut == iSetIn) {
3806 // key swap
3807 int key;
3808 if (sequenceIn >= numberRows + numberColumns) {
3809 key = numberColumns + iSetIn;
3810 setStatus(iSetIn, ClpSimplex::basic);
3811 } else {
3812 key = sequenceIn;
3813 }
3814 redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3815 } else {
3816 // key was chosen
3817 assert (possiblePivotKey_ >= 0 && possiblePivotKey_ < numberRows);
3818 int key = pivotVariable[possiblePivotKey_];
3819 // and set incoming here
3820 if (sequenceIn >= numberRows + numberColumns) {
3821 // slack in - so use old key
3822 sequenceIn = keyVariable_[iSetIn];
3823 model->setStatus(sequenceIn, ClpSimplex::basic);
3824 setStatus(iSetIn, ClpSimplex::basic);
3825 redoSet(model, iSetIn + numberColumns, keyVariable_[iSetIn], iSetIn);
3826 }
3827 //? do not do if iSetIn<0 ? as will be done later
3828 pivotVariable[possiblePivotKey_] = sequenceIn;
3829 if (sequenceIn < numberColumns)
3830 backToPivotRow_[sequenceIn] = possiblePivotKey_;
3831 redoSet(model, key, keyVariable_[iSetOut], iSetOut);
3832 }
3833 } else {
3834 if (sequenceOut < numberColumns) {
3835 if (iSetIn >= 0 && iSetOut == iSetIn) {
3836 // key not out - only problem is if slack in
3837 int key;
3838 if (sequenceIn >= numberRows + numberColumns) {
3839 key = numberColumns + iSetIn;
3840 setStatus(iSetIn, ClpSimplex::basic);
3841 assert (pivotRow < numberRows);
3842 // must swap with current key
3843 int key = keyVariable_[iSetIn];
3844 model->setStatus(key, ClpSimplex::basic);
3845 pivotVariable[pivotRow] = key;
3846 backToPivotRow_[key] = pivotRow;
3847 } else {
3848 key = keyVariable_[iSetIn];
3849 }
3850 redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3851 } else if (iSetOut >= 0) {
3852 // just redo set
3853 int key = keyVariable_[iSetOut];
3854 redoSet(model, key, keyVariable_[iSetOut], iSetOut);
3855 }
3856 }
3857 }
3858 if (iSetIn >= 0 && iSetIn != iSetOut) {
3859 int key = keyVariable_[iSetIn];
3860 if (sequenceIn == numberColumns + 2 * numberRows) {
3861 // key slack in
3862 assert (pivotRow < numberRows);
3863 // must swap with current key
3864 model->setStatus(key, ClpSimplex::basic);
3865 pivotVariable[pivotRow] = key;
3866 backToPivotRow_[key] = pivotRow;
3867 setStatus(iSetIn, ClpSimplex::basic);
3868 key = iSetIn + numberColumns;
3869 }
3870 // redo set to allow for new one
3871 redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3872 }
3873 // update pivot
3874 if (sequenceIn < numberColumns) {
3875 if (pivotRow < numberRows) {
3876 backToPivotRow_[sequenceIn] = pivotRow;
3877 } else {
3878 if (possiblePivotKey_ >= 0) {
3879 assert (possiblePivotKey_ < numberRows);
3880 backToPivotRow_[sequenceIn] = possiblePivotKey_;
3881 pivotVariable[possiblePivotKey_] = sequenceIn;
3882 }
3883 }
3884 } else if (sequenceIn >= numberRows + numberColumns) {
3885 // key in - something should have been done before
3886 int key = keyVariable_[iSetIn];
3887 assert (key == numberColumns + iSetIn);
3888 //pivotVariable[pivotRow]=key;
3889 //backToPivotRow_[key]=pivotRow;
3890 //model->setStatus(key,ClpSimplex::basic);
3891 //key=numberColumns+iSetIn;
3892 setStatus(iSetIn, ClpSimplex::basic);
3893 redoSet(model, key, keyVariable_[iSetIn], iSetIn);
3894 }
3895#ifdef CLP_DEBUG
3896 {
3897 char * xx = new char[numberColumns+numberRows];
3898 memset(xx, 0, numberRows + numberColumns);
3899 for (int i = 0; i < numberRows; i++) {
3900 int iPivot = pivotVariable[i];
3901 assert (iPivot < numberRows + numberColumns);
3902 assert (!xx[iPivot]);
3903 xx[iPivot] = 1;
3904 if (iPivot < numberColumns) {
3905 int iBack = backToPivotRow_[iPivot];
3906 assert (i == iBack);
3907 }
3908 }
3909 delete [] xx;
3910 }
3911#endif
3912 if (rhsOffset_) {
3913 // update effective rhs
3914 if (sequenceIn != sequenceOut) {
3915 if (sequenceIn < numberColumns) {
3916 if (iSetIn >= 0) {
3917 // new contribution to rhsOffset_
3918 int key = keyVariable_[iSetIn];
3919 if (key < numberColumns) {
3920 double newB = 0.0;
3921 ClpSimplex::Status iStatus = getStatus(iSetIn);
3922 if (iStatus == ClpSimplex::atLowerBound)
3923 newB = lower_[iSetIn];
3924 else
3925 newB = upper_[iSetIn];
3926 // subtract out others at bounds
3927 if ((gubType_ & 8) == 0) {
3928 int stop = -(key + 1);
3929 int iColumn = next_[key];
3930 // skip basic
3931 while (iColumn >= 0)
3932 iColumn = next_[iColumn];
3933 // sum all non-key variables
3934 while(iColumn != stop) {
3935 assert (iColumn < 0);
3936 iColumn = -iColumn - 1;
3937 newB -= solution[iColumn];
3938 iColumn = next_[iColumn];
3939 }
3940 }
3941 if (newB)
3942 ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
3943 }
3944 }
3945 }
3946 if (iSetOut >= 0) {
3947 // new contribution to rhsOffset_
3948 int key = keyVariable_[iSetOut];
3949 if (key < numberColumns && iSetIn != iSetOut) {
3950 double newB = 0.0;
3951 ClpSimplex::Status iStatus = getStatus(iSetOut);
3952 if (iStatus == ClpSimplex::atLowerBound)
3953 newB = lower_[iSetOut];
3954 else
3955 newB = upper_[iSetOut];
3956 // subtract out others at bounds
3957 if ((gubType_ & 8) == 0) {
3958 int stop = -(key + 1);
3959 int iColumn = next_[key];
3960 // skip basic
3961 while (iColumn >= 0)
3962 iColumn = next_[iColumn];
3963 // sum all non-key variables
3964 while(iColumn != stop) {
3965 assert (iColumn < 0);
3966 iColumn = -iColumn - 1;
3967 newB -= solution[iColumn];
3968 iColumn = next_[iColumn];
3969 }
3970 }
3971 if (newB)
3972 ClpPackedMatrix::add(model, rhsOffset_, key, -newB);
3973 }
3974 }
3975 }
3976 }
3977#ifdef CLP_DEBUG
3978 // debug
3979 {
3980 int i;
3981 char * xxxx = new char[numberColumns];
3982 memset(xxxx, 0, numberColumns);
3983 for (i = 0; i < numberRows; i++) {
3984 int iPivot = pivotVariable[i];
3985 assert (model->getStatus(iPivot) == ClpSimplex::basic);
3986 if (iPivot < numberColumns && backward_[iPivot] >= 0)
3987 xxxx[iPivot] = 1;
3988 }
3989 double primalTolerance = model->primalTolerance();
3990 for (i = 0; i < numberSets_; i++) {
3991 int key = keyVariable_[i];
3992 double value = 0.0;
3993 // sum over all except key
3994 int iColumn = next_[key];
3995 // sum all non-key variables
3996 int k = 0;
3997 int stop = -(key + 1);
3998 while (iColumn != stop) {
3999 if (iColumn < 0)
4000 iColumn = -iColumn - 1;
4001 value += solution[iColumn];
4002 k++;
4003 assert (k < 100);
4004 assert (backward_[iColumn] == i);
4005 iColumn = next_[iColumn];
4006 }
4007 iColumn = next_[key];
4008 if (key < numberColumns) {
4009 // feasibility will be done later
4010 assert (getStatus(i) != ClpSimplex::basic);
4011 double sol;
4012 if (getStatus(i) == ClpSimplex::atUpperBound)
4013 sol = upper_[i] - value;
4014 else
4015 sol = lower_[i] - value;
4016 //printf("xx Value of key structural %d for set %d is %g - cost %g\n",key,i,sol,
4017 // cost[key]);
4018 //if (fabs(sol-solution[key])>1.0e-3)
4019 //printf("** stored value was %g\n",solution[key]);
4020 } else {
4021 // slack is key
4022 double infeasibility = 0.0;
4023 if (value > upper_[i] + primalTolerance) {
4024 infeasibility = value - upper_[i] - primalTolerance;
4025 //setAbove(i);
4026 } else if (value < lower_[i] - primalTolerance) {
4027 infeasibility = lower_[i] - value - primalTolerance ;
4028 //setBelow(i);
4029 } else {
4030 //setFeasible(i);
4031 }
4032 //printf("xx Value of key slack for set %d is %g\n",i,value);
4033 }
4034 while (iColumn >= 0) {
4035 assert (xxxx[iColumn]);
4036 xxxx[iColumn] = 0;
4037 iColumn = next_[iColumn];
4038 }
4039 }
4040 for (i = 0; i < numberColumns; i++) {
4041 if (i < numberColumns && backward_[i] >= 0) {
4042 assert (!xxxx[i] || i == keyVariable_[backward_[i]]);
4043 }
4044 }
4045 delete [] xxxx;
4046 }
4047#endif
4048 return 0;
4049}
4050// Switches off dj checking each factorization (for BIG models)
4051void
4052ClpGubMatrix::switchOffCheck()
4053{
4054 noCheck_ = 0;
4055 infeasibilityWeight_ = 0.0;
4056}
4057// Correct sequence in and out to give true value
4058void
4059ClpGubMatrix::correctSequence(const ClpSimplex * /*model*/, int & sequenceIn, int & sequenceOut)
4060{
4061 if (sequenceIn != -999) {
4062 sequenceIn = trueSequenceIn_;
4063 sequenceOut = trueSequenceOut_;
4064 }
4065}
4066