1/* $Id: ClpSimplexOther.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7
8#include <math.h>
9
10#include "CoinHelperFunctions.hpp"
11#include "ClpSimplexOther.hpp"
12#include "ClpSimplexDual.hpp"
13#include "ClpSimplexPrimal.hpp"
14#include "ClpEventHandler.hpp"
15#include "ClpHelperFunctions.hpp"
16#include "ClpFactorization.hpp"
17#include "ClpDualRowDantzig.hpp"
18#include "ClpDynamicMatrix.hpp"
19#include "CoinPackedMatrix.hpp"
20#include "CoinIndexedVector.hpp"
21#include "CoinBuild.hpp"
22#include "CoinMpsIO.hpp"
23#include "CoinFloatEqual.hpp"
24#include "ClpMessage.hpp"
25#include <cfloat>
26#include <cassert>
27#include <string>
28#include <stdio.h>
29#include <iostream>
30/* Dual ranging.
31 This computes increase/decrease in cost for each given variable and corresponding
32 sequence numbers which would change basis. Sequence numbers are 0..numberColumns
33 and numberColumns.. for artificials/slacks.
34 For non-basic variables the sequence number will be that of the non-basic variables.
35
36 Up to user to provide correct length arrays.
37
38*/
39void ClpSimplexOther::dualRanging(int numberCheck, const int * which,
40 double * costIncreased, int * sequenceIncreased,
41 double * costDecreased, int * sequenceDecreased,
42 double * valueIncrease, double * valueDecrease)
43{
44 rowArray_[1]->clear();
45 columnArray_[1]->clear();
46 // long enough for rows+columns
47 assert(rowArray_[3]->capacity() >= numberRows_ + numberColumns_);
48 rowArray_[3]->clear();
49 int * backPivot = rowArray_[3]->getIndices();
50 int i;
51 for ( i = 0; i < numberRows_ + numberColumns_; i++) {
52 backPivot[i] = -1;
53 }
54 for (i = 0; i < numberRows_; i++) {
55 int iSequence = pivotVariable_[i];
56 backPivot[iSequence] = i;
57 }
58 // dualTolerance may be zero if from CBC. In fact use that fact
59 bool inCBC = !dualTolerance_;
60 if (inCBC)
61 assert (integerType_);
62 dualTolerance_ = dblParam_[ClpDualTolerance];
63 double * arrayX = rowArray_[0]->denseVector();
64 for ( i = 0; i < numberCheck; i++) {
65 rowArray_[0]->clear();
66 //rowArray_[0]->checkClear();
67 //rowArray_[1]->checkClear();
68 //columnArray_[1]->checkClear();
69 columnArray_[0]->clear();
70 //columnArray_[0]->checkClear();
71 int iSequence = which[i];
72 if (iSequence < 0) {
73 costIncreased[i] = 0.0;
74 sequenceIncreased[i] = -1;
75 costDecreased[i] = 0.0;
76 sequenceDecreased[i] = -1;
77 continue;
78 }
79 double costIncrease = COIN_DBL_MAX;
80 double costDecrease = COIN_DBL_MAX;
81 int sequenceIncrease = -1;
82 int sequenceDecrease = -1;
83 if (valueIncrease) {
84 assert (valueDecrease);
85 valueIncrease[i] = iSequence < numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_];
86 valueDecrease[i] = valueIncrease[i];
87 }
88
89 switch(getStatus(iSequence)) {
90
91 case basic: {
92 // non-trvial
93 // Get pivot row
94 int iRow = backPivot[iSequence];
95 assert (iRow >= 0);
96 double plusOne = 1.0;
97 rowArray_[0]->createPacked(1, &iRow, &plusOne);
98 factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
99 // put row of tableau in rowArray[0] and columnArray[0]
100 matrix_->transposeTimes(this, -1.0,
101 rowArray_[0], columnArray_[1], columnArray_[0]);
102 double alphaIncrease;
103 double alphaDecrease;
104 // do ratio test up and down
105 checkDualRatios(rowArray_[0], columnArray_[0], costIncrease, sequenceIncrease, alphaIncrease,
106 costDecrease, sequenceDecrease, alphaDecrease);
107 if (!inCBC) {
108 if (valueIncrease) {
109 if (sequenceIncrease >= 0)
110 valueIncrease[i] = primalRanging1(sequenceIncrease, iSequence);
111 if (sequenceDecrease >= 0)
112 valueDecrease[i] = primalRanging1(sequenceDecrease, iSequence);
113 }
114 } else {
115 int number = rowArray_[0]->getNumElements();
116 double scale2 = 0.0;
117 int j;
118 for (j = 0; j < number; j++) {
119 scale2 += arrayX[j] * arrayX[j];
120 }
121 scale2 = 1.0 / sqrt(scale2);
122 //valueIncrease[i] = scale2;
123 if (sequenceIncrease >= 0) {
124 double djValue = dj_[sequenceIncrease];
125 if (fabs(djValue) > 10.0 * dualTolerance_) {
126 // we are going to use for cutoff so be exact
127 costIncrease = fabs(djValue / alphaIncrease);
128 /* Not sure this is good idea as I don't think correct e.g.
129 suppose a continuous variable has dj slightly greater. */
130 if(false && sequenceIncrease < numberColumns_ && integerType_[sequenceIncrease]) {
131 // can improve
132 double movement = (columnScale_ == NULL) ? 1.0 :
133 rhsScale_ * inverseColumnScale_[sequenceIncrease];
134 costIncrease = CoinMax(fabs(djValue * movement), costIncrease);
135 }
136 } else {
137 costIncrease = 0.0;
138 }
139 }
140 if (sequenceDecrease >= 0) {
141 double djValue = dj_[sequenceDecrease];
142 if (fabs(djValue) > 10.0 * dualTolerance_) {
143 // we are going to use for cutoff so be exact
144 costDecrease = fabs(djValue / alphaDecrease);
145 if(sequenceDecrease < numberColumns_ && integerType_[sequenceDecrease]) {
146 // can improve
147 double movement = (columnScale_ == NULL) ? 1.0 :
148 rhsScale_ * inverseColumnScale_[sequenceDecrease];
149 costDecrease = CoinMax(fabs(djValue * movement), costDecrease);
150 }
151 } else {
152 costDecrease = 0.0;
153 }
154 }
155 costIncrease *= scale2;
156 costDecrease *= scale2;
157 }
158 }
159 break;
160 case isFixed:
161 break;
162 case isFree:
163 case superBasic:
164 costIncrease = 0.0;
165 costDecrease = 0.0;
166 sequenceIncrease = iSequence;
167 sequenceDecrease = iSequence;
168 break;
169 case atUpperBound:
170 costIncrease = CoinMax(0.0, -dj_[iSequence]);
171 sequenceIncrease = iSequence;
172 if (valueIncrease)
173 valueIncrease[i] = primalRanging1(iSequence, iSequence);
174 break;
175 case atLowerBound:
176 costDecrease = CoinMax(0.0, dj_[iSequence]);
177 sequenceDecrease = iSequence;
178 if (valueIncrease)
179 valueDecrease[i] = primalRanging1(iSequence, iSequence);
180 break;
181 }
182 double scaleFactor;
183 if (rowScale_) {
184 if (iSequence < numberColumns_)
185 scaleFactor = 1.0 / (objectiveScale_ * columnScale_[iSequence]);
186 else
187 scaleFactor = rowScale_[iSequence-numberColumns_] / objectiveScale_;
188 } else {
189 scaleFactor = 1.0 / objectiveScale_;
190 }
191 if (costIncrease < 1.0e30)
192 costIncrease *= scaleFactor;
193 if (costDecrease < 1.0e30)
194 costDecrease *= scaleFactor;
195 if (optimizationDirection_ == 1.0) {
196 costIncreased[i] = costIncrease;
197 sequenceIncreased[i] = sequenceIncrease;
198 costDecreased[i] = costDecrease;
199 sequenceDecreased[i] = sequenceDecrease;
200 } else if (optimizationDirection_ == -1.0) {
201 costIncreased[i] = costDecrease;
202 sequenceIncreased[i] = sequenceDecrease;
203 costDecreased[i] = costIncrease;
204 sequenceDecreased[i] = sequenceIncrease;
205 if (valueIncrease) {
206 double temp = valueIncrease[i];
207 valueIncrease[i] = valueDecrease[i];
208 valueDecrease[i] = temp;
209 }
210 } else if (optimizationDirection_ == 0.0) {
211 // !!!!!! ???
212 costIncreased[i] = COIN_DBL_MAX;
213 sequenceIncreased[i] = -1;
214 costDecreased[i] = COIN_DBL_MAX;
215 sequenceDecreased[i] = -1;
216 } else {
217 abort();
218 }
219 }
220 rowArray_[0]->clear();
221 //rowArray_[1]->clear();
222 //columnArray_[1]->clear();
223 columnArray_[0]->clear();
224 //rowArray_[3]->clear();
225 if (!optimizationDirection_)
226 printf("*** ????? Ranging with zero optimization costs\n");
227}
228/*
229 Row array has row part of pivot row
230 Column array has column part.
231 This is used in dual ranging
232*/
233void
234ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
235 CoinIndexedVector * columnArray,
236 double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
237 double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
238{
239 double acceptablePivot = 1.0e-9;
240 double * work;
241 int number;
242 int * which;
243 int iSection;
244
245 double thetaDown = 1.0e31;
246 double thetaUp = 1.0e31;
247 int sequenceDown = -1;
248 int sequenceUp = -1;
249 double alphaDown = 0.0;
250 double alphaUp = 0.0;
251
252 int addSequence;
253
254 for (iSection = 0; iSection < 2; iSection++) {
255
256 int i;
257 if (!iSection) {
258 work = rowArray->denseVector();
259 number = rowArray->getNumElements();
260 which = rowArray->getIndices();
261 addSequence = numberColumns_;
262 } else {
263 work = columnArray->denseVector();
264 number = columnArray->getNumElements();
265 which = columnArray->getIndices();
266 addSequence = 0;
267 }
268
269 for (i = 0; i < number; i++) {
270 int iSequence = which[i];
271 int iSequence2 = iSequence + addSequence;
272 double alpha = work[i];
273 if (fabs(alpha) < acceptablePivot)
274 continue;
275 double oldValue = dj_[iSequence2];
276
277 switch(getStatus(iSequence2)) {
278
279 case basic:
280 break;
281 case ClpSimplex::isFixed:
282 break;
283 case isFree:
284 case superBasic:
285 // treat dj as if zero
286 thetaDown = 0.0;
287 thetaUp = 0.0;
288 sequenceDown = iSequence2;
289 sequenceUp = iSequence2;
290 break;
291 case atUpperBound:
292 if (alpha > 0.0) {
293 // test up
294 if (oldValue + thetaUp * alpha > dualTolerance_) {
295 thetaUp = (dualTolerance_ - oldValue) / alpha;
296 sequenceUp = iSequence2;
297 alphaUp = alpha;
298 }
299 } else {
300 // test down
301 if (oldValue - thetaDown * alpha > dualTolerance_) {
302 thetaDown = -(dualTolerance_ - oldValue) / alpha;
303 sequenceDown = iSequence2;
304 alphaDown = alpha;
305 }
306 }
307 break;
308 case atLowerBound:
309 if (alpha < 0.0) {
310 // test up
311 if (oldValue + thetaUp * alpha < - dualTolerance_) {
312 thetaUp = -(dualTolerance_ + oldValue) / alpha;
313 sequenceUp = iSequence2;
314 alphaUp = alpha;
315 }
316 } else {
317 // test down
318 if (oldValue - thetaDown * alpha < -dualTolerance_) {
319 thetaDown = (dualTolerance_ + oldValue) / alpha;
320 sequenceDown = iSequence2;
321 alphaDown = alpha;
322 }
323 }
324 break;
325 }
326 }
327 }
328 if (sequenceUp >= 0) {
329 costIncrease = thetaUp;
330 sequenceIncrease = sequenceUp;
331 alphaIncrease = alphaUp;
332 }
333 if (sequenceDown >= 0) {
334 costDecrease = thetaDown;
335 sequenceDecrease = sequenceDown;
336 alphaDecrease = alphaDown;
337 }
338}
339/** Primal ranging.
340 This computes increase/decrease in value for each given variable and corresponding
341 sequence numbers which would change basis. Sequence numbers are 0..numberColumns
342 and numberColumns.. for artificials/slacks.
343 For basic variables the sequence number will be that of the basic variables.
344
345 Up to user to provide correct length arrays.
346
347 When here - guaranteed optimal
348*/
349void
350ClpSimplexOther::primalRanging(int numberCheck, const int * which,
351 double * valueIncreased, int * sequenceIncreased,
352 double * valueDecreased, int * sequenceDecreased)
353{
354 rowArray_[0]->clear();
355 rowArray_[1]->clear();
356 lowerIn_ = -COIN_DBL_MAX;
357 upperIn_ = COIN_DBL_MAX;
358 valueIn_ = 0.0;
359 for ( int i = 0; i < numberCheck; i++) {
360 int iSequence = which[i];
361 double valueIncrease = COIN_DBL_MAX;
362 double valueDecrease = COIN_DBL_MAX;
363 int sequenceIncrease = -1;
364 int sequenceDecrease = -1;
365
366 switch(getStatus(iSequence)) {
367
368 case basic:
369 case isFree:
370 case superBasic:
371 // Easy
372 valueDecrease = CoinMax(0.0, upper_[iSequence] - solution_[iSequence]);
373 valueIncrease = CoinMax(0.0, solution_[iSequence] - lower_[iSequence]);
374 sequenceDecrease = iSequence;
375 sequenceIncrease = iSequence;
376 break;
377 case isFixed:
378 case atUpperBound:
379 case atLowerBound: {
380 // Non trivial
381 // Other bound is ignored
382 unpackPacked(rowArray_[1], iSequence);
383 factorization_->updateColumn(rowArray_[2], rowArray_[1]);
384 // Get extra rows
385 matrix_->extendUpdated(this, rowArray_[1], 0);
386 // do ratio test
387 checkPrimalRatios(rowArray_[1], 1);
388 if (pivotRow_ >= 0) {
389 valueIncrease = theta_;
390 sequenceIncrease = pivotVariable_[pivotRow_];
391 }
392 checkPrimalRatios(rowArray_[1], -1);
393 if (pivotRow_ >= 0) {
394 valueDecrease = theta_;
395 sequenceDecrease = pivotVariable_[pivotRow_];
396 }
397 rowArray_[1]->clear();
398 }
399 break;
400 }
401 double scaleFactor;
402 if (rowScale_) {
403 if (iSequence < numberColumns_)
404 scaleFactor = columnScale_[iSequence] / rhsScale_;
405 else
406 scaleFactor = 1.0 / (rowScale_[iSequence-numberColumns_] * rhsScale_);
407 } else {
408 scaleFactor = 1.0 / rhsScale_;
409 }
410 if (valueIncrease < 1.0e30)
411 valueIncrease *= scaleFactor;
412 else
413 valueIncrease = COIN_DBL_MAX;
414 if (valueDecrease < 1.0e30)
415 valueDecrease *= scaleFactor;
416 else
417 valueDecrease = COIN_DBL_MAX;
418 valueIncreased[i] = valueIncrease;
419 sequenceIncreased[i] = sequenceIncrease;
420 valueDecreased[i] = valueDecrease;
421 sequenceDecreased[i] = sequenceDecrease;
422 }
423}
424// Returns new value of whichOther when whichIn enters basis
425double
426ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
427{
428 rowArray_[0]->clear();
429 rowArray_[1]->clear();
430 int iSequence = whichIn;
431 double newValue = solution_[whichOther];
432 double alphaOther = 0.0;
433 Status status = getStatus(iSequence);
434 assert (status == atLowerBound || status == atUpperBound);
435 int wayIn = (status == atLowerBound) ? 1 : -1;
436
437 switch(getStatus(iSequence)) {
438
439 case basic:
440 case isFree:
441 case superBasic:
442 assert (whichIn == whichOther);
443 // Easy
444 newValue = wayIn > 0 ? upper_[iSequence] : lower_[iSequence];
445 break;
446 case isFixed:
447 case atUpperBound:
448 case atLowerBound:
449 // Non trivial
450 {
451 // Other bound is ignored
452 unpackPacked(rowArray_[1], iSequence);
453 factorization_->updateColumn(rowArray_[2], rowArray_[1]);
454 // Get extra rows
455 matrix_->extendUpdated(this, rowArray_[1], 0);
456 // do ratio test
457 double acceptablePivot = 1.0e-7;
458 double * work = rowArray_[1]->denseVector();
459 int number = rowArray_[1]->getNumElements();
460 int * which = rowArray_[1]->getIndices();
461
462 // we may need to swap sign
463 double way = wayIn;
464 double theta = 1.0e30;
465 for (int iIndex = 0; iIndex < number; iIndex++) {
466
467 int iRow = which[iIndex];
468 double alpha = work[iIndex] * way;
469 int iPivot = pivotVariable_[iRow];
470 if (iPivot == whichOther) {
471 alphaOther = alpha;
472 continue;
473 }
474 double oldValue = solution_[iPivot];
475 if (fabs(alpha) > acceptablePivot) {
476 if (alpha > 0.0) {
477 // basic variable going towards lower bound
478 double bound = lower_[iPivot];
479 oldValue -= bound;
480 if (oldValue - theta * alpha < 0.0) {
481 theta = CoinMax(0.0, oldValue / alpha);
482 }
483 } else {
484 // basic variable going towards upper bound
485 double bound = upper_[iPivot];
486 oldValue = oldValue - bound;
487 if (oldValue - theta * alpha > 0.0) {
488 theta = CoinMax(0.0, oldValue / alpha);
489 }
490 }
491 }
492 }
493 if (whichIn != whichOther) {
494 if (theta < 1.0e30)
495 newValue -= theta * alphaOther;
496 else
497 newValue = alphaOther > 0.0 ? -1.0e30 : 1.0e30;
498 } else {
499 newValue += theta * wayIn;
500 }
501 }
502 rowArray_[1]->clear();
503 break;
504 }
505 double scaleFactor;
506 if (rowScale_) {
507 if (whichOther < numberColumns_)
508 scaleFactor = columnScale_[whichOther] / rhsScale_;
509 else
510 scaleFactor = 1.0 / (rowScale_[whichOther-numberColumns_] * rhsScale_);
511 } else {
512 scaleFactor = 1.0 / rhsScale_;
513 }
514 if (newValue < 1.0e29)
515 if (newValue > -1.0e29)
516 newValue *= scaleFactor;
517 else
518 newValue = -COIN_DBL_MAX;
519 else
520 newValue = COIN_DBL_MAX;
521 return newValue;
522}
523/*
524 Row array has pivot column
525 This is used in primal ranging
526*/
527void
528ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
529 int direction)
530{
531 // sequence stays as row number until end
532 pivotRow_ = -1;
533 double acceptablePivot = 1.0e-7;
534 double * work = rowArray->denseVector();
535 int number = rowArray->getNumElements();
536 int * which = rowArray->getIndices();
537
538 // we need to swap sign if going down
539 double way = direction;
540 theta_ = 1.0e30;
541 for (int iIndex = 0; iIndex < number; iIndex++) {
542
543 int iRow = which[iIndex];
544 double alpha = work[iIndex] * way;
545 int iPivot = pivotVariable_[iRow];
546 double oldValue = solution_[iPivot];
547 if (fabs(alpha) > acceptablePivot) {
548 if (alpha > 0.0) {
549 // basic variable going towards lower bound
550 double bound = lower_[iPivot];
551 oldValue -= bound;
552 if (oldValue - theta_ * alpha < 0.0) {
553 pivotRow_ = iRow;
554 theta_ = CoinMax(0.0, oldValue / alpha);
555 }
556 } else {
557 // basic variable going towards upper bound
558 double bound = upper_[iPivot];
559 oldValue = oldValue - bound;
560 if (oldValue - theta_ * alpha > 0.0) {
561 pivotRow_ = iRow;
562 theta_ = CoinMax(0.0, oldValue / alpha);
563 }
564 }
565 }
566 }
567}
568/* Write the basis in MPS format to the specified file.
569 If writeValues true writes values of structurals
570 (and adds VALUES to end of NAME card)
571
572 Row and column names may be null.
573 formatType is
574 <ul>
575 <li> 0 - normal
576 <li> 1 - extra accuracy
577 <li> 2 - IEEE hex (later)
578 </ul>
579
580 Returns non-zero on I/O error
581
582 This is based on code contributed by Thorsten Koch
583*/
584int
585ClpSimplexOther::writeBasis(const char *filename,
586 bool writeValues,
587 int formatType) const
588{
589 formatType = CoinMax(0, formatType);
590 formatType = CoinMin(2, formatType);
591 if (!writeValues)
592 formatType = 0;
593 // See if INTEL if IEEE
594 if (formatType == 2) {
595 // test intel here and add 1 if not intel
596 double value = 1.0;
597 char x[8];
598 memcpy(x, &value, 8);
599 if (x[0] == 63) {
600 formatType ++; // not intel
601 } else {
602 assert (x[0] == 0);
603 }
604 }
605
606 char number[20];
607 FILE * fp = fopen(filename, "w");
608 if (!fp)
609 return -1;
610
611 // NAME card
612
613 if (strcmp(strParam_[ClpProbName].c_str(), "") == 0) {
614 fprintf(fp, "NAME BLANK ");
615 } else {
616 fprintf(fp, "NAME %s ", strParam_[ClpProbName].c_str());
617 }
618 if (formatType >= 2)
619 fprintf(fp, "FREEIEEE");
620 else if (writeValues)
621 fprintf(fp, "VALUES");
622 // finish off name
623 fprintf(fp, "\n");
624 int iRow = 0;
625 for(int iColumn = 0; iColumn < numberColumns_; iColumn++) {
626 bool printit = false;
627 if( getColumnStatus(iColumn) == ClpSimplex::basic) {
628 printit = true;
629 // Find non basic row
630 for(; iRow < numberRows_; iRow++) {
631 if (getRowStatus(iRow) != ClpSimplex::basic)
632 break;
633 }
634 if (lengthNames_) {
635 if (iRow != numberRows_) {
636 fprintf(fp, " %s %-8s %s",
637 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
638 columnNames_[iColumn].c_str(),
639 rowNames_[iRow].c_str());
640 iRow++;
641 } else {
642 // Allow for too many basics!
643 fprintf(fp, " BS %-8s ",
644 columnNames_[iColumn].c_str());
645 // Dummy row name if values
646 if (writeValues)
647 fprintf(fp, " _dummy_");
648 }
649 } else {
650 // no names
651 if (iRow != numberRows_) {
652 fprintf(fp, " %s C%7.7d R%7.7d",
653 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
654 iColumn, iRow);
655 iRow++;
656 } else {
657 // Allow for too many basics!
658 fprintf(fp, " BS C%7.7d", iColumn);
659 // Dummy row name if values
660 if (writeValues)
661 fprintf(fp, " _dummy_");
662 }
663 }
664 } else {
665 if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
666 printit = true;
667 if (lengthNames_)
668 fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
669 else
670 fprintf(fp, " UL C%7.7d", iColumn);
671 // Dummy row name if values
672 if (writeValues)
673 fprintf(fp, " _dummy_");
674 }
675 }
676 if (printit && writeValues) {
677 // add value
678 CoinConvertDouble(0, formatType, columnActivity_[iColumn], number);
679 fprintf(fp, " %s", number);
680 }
681 if (printit)
682 fprintf(fp, "\n");
683 }
684 fprintf(fp, "ENDATA\n");
685 fclose(fp);
686 return 0;
687}
688// Read a basis from the given filename
689int
690ClpSimplexOther::readBasis(const char *fileName)
691{
692 int status = 0;
693 bool canOpen = false;
694 if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
695 // stdin
696 canOpen = true;
697 } else {
698 FILE *fp = fopen(fileName, "r");
699 if (fp) {
700 // can open - lets go for it
701 fclose(fp);
702 canOpen = true;
703 } else {
704 handler_->message(CLP_UNABLE_OPEN, messages_)
705 << fileName << CoinMessageEol;
706 return -1;
707 }
708 }
709 CoinMpsIO m;
710 m.passInMessageHandler(handler_);
711 *m.messagesPointer() = coinMessages();
712 bool savePrefix = m.messageHandler()->prefix();
713 m.messageHandler()->setPrefix(handler_->prefix());
714 status = m.readBasis(fileName, "", columnActivity_, status_ + numberColumns_,
715 status_,
716 columnNames_, numberColumns_,
717 rowNames_, numberRows_);
718 m.messageHandler()->setPrefix(savePrefix);
719 if (status >= 0) {
720 if (!status) {
721 // set values
722 int iColumn, iRow;
723 for (iRow = 0; iRow < numberRows_; iRow++) {
724 if (getRowStatus(iRow) == atLowerBound)
725 rowActivity_[iRow] = rowLower_[iRow];
726 else if (getRowStatus(iRow) == atUpperBound)
727 rowActivity_[iRow] = rowUpper_[iRow];
728 }
729 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
730 if (getColumnStatus(iColumn) == atLowerBound)
731 columnActivity_[iColumn] = columnLower_[iColumn];
732 else if (getColumnStatus(iColumn) == atUpperBound)
733 columnActivity_[iColumn] = columnUpper_[iColumn];
734 }
735 } else {
736 memset(rowActivity_, 0, numberRows_ * sizeof(double));
737 matrix_->times(-1.0, columnActivity_, rowActivity_);
738 }
739 } else {
740 // errors
741 handler_->message(CLP_IMPORT_ERRORS, messages_)
742 << status << fileName << CoinMessageEol;
743 }
744 return status;
745}
746/* Creates dual of a problem if looks plausible
747 (defaults will always create model)
748 fractionRowRanges is fraction of rows allowed to have ranges
749 fractionColumnRanges is fraction of columns allowed to have ranges
750*/
751ClpSimplex *
752ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const
753{
754 const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
755 bool changed = false;
756 int numberChanged = 0;
757 int iColumn;
758 // check if we need to change bounds to rows
759 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
760 if (columnUpper_[iColumn] < 1.0e20 &&
761 columnLower_[iColumn] > -1.0e20) {
762 changed = true;
763 numberChanged++;
764 }
765 }
766 int iRow;
767 int numberExtraRows = 0;
768 if (numberChanged <= fractionColumnRanges * numberColumns_) {
769 for (iRow = 0; iRow < numberRows_; iRow++) {
770 if (rowLower_[iRow] > -1.0e20 &&
771 rowUpper_[iRow] < 1.0e20) {
772 if (rowUpper_[iRow] != rowLower_[iRow])
773 numberExtraRows++;
774 }
775 }
776 if (numberExtraRows > fractionRowRanges * numberRows_)
777 return NULL;
778 } else {
779 return NULL;
780 }
781 if (changed) {
782 ClpSimplex * model3 = new ClpSimplex(*model2);
783 CoinBuild build;
784 double one = 1.0;
785 int numberColumns = model3->numberColumns();
786 const double * columnLower = model3->columnLower();
787 const double * columnUpper = model3->columnUpper();
788 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
789 if (columnUpper[iColumn] < 1.0e20 &&
790 columnLower[iColumn] > -1.0e20) {
791 if (fabs(columnLower[iColumn]) < fabs(columnUpper[iColumn])) {
792 double value = columnUpper[iColumn];
793 model3->setColumnUpper(iColumn, COIN_DBL_MAX);
794 build.addRow(1, &iColumn, &one, -COIN_DBL_MAX, value);
795 } else {
796 double value = columnLower[iColumn];
797 model3->setColumnLower(iColumn, -COIN_DBL_MAX);
798 build.addRow(1, &iColumn, &one, value, COIN_DBL_MAX);
799 }
800 }
801 }
802 model3->addRows(build);
803 model2 = model3;
804 }
805 int numberColumns = model2->numberColumns();
806 const double * columnLower = model2->columnLower();
807 const double * columnUpper = model2->columnUpper();
808 int numberRows = model2->numberRows();
809 double * rowLower = CoinCopyOfArray(model2->rowLower(), numberRows);
810 double * rowUpper = CoinCopyOfArray(model2->rowUpper(), numberRows);
811
812 const double * objective = model2->objective();
813 CoinPackedMatrix * matrix = model2->matrix();
814 // get transpose
815 CoinPackedMatrix rowCopy = *matrix;
816 const int * row = matrix->getIndices();
817 const int * columnLength = matrix->getVectorLengths();
818 const CoinBigIndex * columnStart = matrix->getVectorStarts();
819 const double * elementByColumn = matrix->getElements();
820 double objOffset = 0.0;
821 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
822 double offset = 0.0;
823 double objValue = optimizationDirection_ * objective[iColumn];
824 if (columnUpper[iColumn] > 1.0e20) {
825 if (columnLower[iColumn] > -1.0e20)
826 offset = columnLower[iColumn];
827 } else if (columnLower[iColumn] < -1.0e20) {
828 offset = columnUpper[iColumn];
829 } else {
830 // taken care of before
831 abort();
832 }
833 if (offset) {
834 objOffset += offset * objValue;
835 for (CoinBigIndex j = columnStart[iColumn];
836 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
837 int iRow = row[j];
838 if (rowLower[iRow] > -1.0e20)
839 rowLower[iRow] -= offset * elementByColumn[j];
840 if (rowUpper[iRow] < 1.0e20)
841 rowUpper[iRow] -= offset * elementByColumn[j];
842 }
843 }
844 }
845 int * which = new int[numberRows+numberExtraRows];
846 rowCopy.reverseOrdering();
847 rowCopy.transpose();
848 double * fromRowsLower = new double[numberRows+numberExtraRows];
849 double * fromRowsUpper = new double[numberRows+numberExtraRows];
850 double * newObjective = new double[numberRows+numberExtraRows];
851 double * fromColumnsLower = new double[numberColumns];
852 double * fromColumnsUpper = new double[numberColumns];
853 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
854 double objValue = optimizationDirection_ * objective[iColumn];
855 // Offset is already in
856 if (columnUpper[iColumn] > 1.0e20) {
857 if (columnLower[iColumn] > -1.0e20) {
858 fromColumnsLower[iColumn] = -COIN_DBL_MAX;
859 fromColumnsUpper[iColumn] = objValue;
860 } else {
861 // free
862 fromColumnsLower[iColumn] = objValue;
863 fromColumnsUpper[iColumn] = objValue;
864 }
865 } else if (columnLower[iColumn] < -1.0e20) {
866 fromColumnsLower[iColumn] = objValue;
867 fromColumnsUpper[iColumn] = COIN_DBL_MAX;
868 } else {
869 abort();
870 }
871 }
872 int kRow = 0;
873 int kExtraRow = numberRows;
874 for (iRow = 0; iRow < numberRows; iRow++) {
875 if (rowLower[iRow] < -1.0e20) {
876 assert (rowUpper[iRow] < 1.0e20);
877 newObjective[kRow] = -rowUpper[iRow];
878 fromRowsLower[kRow] = -COIN_DBL_MAX;
879 fromRowsUpper[kRow] = 0.0;
880 which[kRow] = iRow;
881 kRow++;
882 } else if (rowUpper[iRow] > 1.0e20) {
883 newObjective[kRow] = -rowLower[iRow];
884 fromRowsLower[kRow] = 0.0;
885 fromRowsUpper[kRow] = COIN_DBL_MAX;
886 which[kRow] = iRow;
887 kRow++;
888 } else {
889 if (rowUpper[iRow] == rowLower[iRow]) {
890 newObjective[kRow] = -rowLower[iRow];
891 fromRowsLower[kRow] = -COIN_DBL_MAX;
892 fromRowsUpper[kRow] = COIN_DBL_MAX;
893 which[kRow] = iRow;
894 kRow++;
895 } else {
896 // range
897 newObjective[kRow] = -rowUpper[iRow];
898 fromRowsLower[kRow] = -COIN_DBL_MAX;
899 fromRowsUpper[kRow] = 0.0;
900 which[kRow] = iRow;
901 kRow++;
902 newObjective[kExtraRow] = -rowLower[iRow];
903 fromRowsLower[kExtraRow] = 0.0;
904 fromRowsUpper[kExtraRow] = COIN_DBL_MAX;
905 which[kExtraRow] = iRow;
906 kExtraRow++;
907 }
908 }
909 }
910 if (numberExtraRows) {
911 CoinPackedMatrix newCopy;
912 newCopy.setExtraGap(0.0);
913 newCopy.setExtraMajor(0.0);
914 newCopy.submatrixOfWithDuplicates(rowCopy, kExtraRow, which);
915 rowCopy = newCopy;
916 }
917 ClpSimplex * modelDual = new ClpSimplex();
918 modelDual->loadProblem(rowCopy, fromRowsLower, fromRowsUpper, newObjective,
919 fromColumnsLower, fromColumnsUpper);
920 modelDual->setObjectiveOffset(objOffset);
921 modelDual->setDualBound(model2->dualBound());
922 modelDual->setInfeasibilityCost(model2->infeasibilityCost());
923 modelDual->setDualTolerance(model2->dualTolerance());
924 modelDual->setPrimalTolerance(model2->primalTolerance());
925 modelDual->setPerturbation(model2->perturbation());
926 modelDual->setSpecialOptions(model2->specialOptions());
927 modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
928 delete [] fromRowsLower;
929 delete [] fromRowsUpper;
930 delete [] fromColumnsLower;
931 delete [] fromColumnsUpper;
932 delete [] newObjective;
933 delete [] which;
934 delete [] rowLower;
935 delete [] rowUpper;
936 if (changed)
937 delete model2;
938 modelDual->createStatus();
939 return modelDual;
940}
941// Restores solution from dualized problem
942int
943ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem)
944{
945 int returnCode = 0;
946 createStatus();
947 // Number of rows in dual problem was original number of columns
948 assert (numberColumns_ == dualProblem->numberRows());
949 // If slack on d-row basic then column at bound otherwise column basic
950 // If d-column basic then rhs tight
951 int numberBasic = 0;
952 int iRow, iColumn = 0;
953 // Get number of extra rows from ranges
954 int numberExtraRows = 0;
955 for (iRow = 0; iRow < numberRows_; iRow++) {
956 if (rowLower_[iRow] > -1.0e20 &&
957 rowUpper_[iRow] < 1.0e20) {
958 if (rowUpper_[iRow] != rowLower_[iRow])
959 numberExtraRows++;
960 }
961 }
962 const double * objective = this->objective();
963 const double * dualDual = dualProblem->dualRowSolution();
964 const double * dualDj = dualProblem->dualColumnSolution();
965 const double * dualSol = dualProblem->primalColumnSolution();
966 const double * dualActs = dualProblem->primalRowSolution();
967#if 0
968 ClpSimplex thisCopy = *this;
969 thisCopy.dual(); // for testing
970 const double * primalDual = thisCopy.dualRowSolution();
971 const double * primalDj = thisCopy.dualColumnSolution();
972 const double * primalSol = thisCopy.primalColumnSolution();
973 const double * primalActs = thisCopy.primalRowSolution();
974 char ss[] = {'F', 'B', 'U', 'L', 'S', 'F'};
975 printf ("Dual problem row info %d rows\n", dualProblem->numberRows());
976 for (iRow = 0; iRow < dualProblem->numberRows(); iRow++)
977 printf("%d at %c primal %g dual %g\n",
978 iRow, ss[dualProblem->getRowStatus(iRow)],
979 dualActs[iRow], dualDual[iRow]);
980 printf ("Dual problem column info %d columns\n", dualProblem->numberColumns());
981 for (iColumn = 0; iColumn < dualProblem->numberColumns(); iColumn++)
982 printf("%d at %c primal %g dual %g\n",
983 iColumn, ss[dualProblem->getColumnStatus(iColumn)],
984 dualSol[iColumn], dualDj[iColumn]);
985 printf ("Primal problem row info %d rows\n", thisCopy.numberRows());
986 for (iRow = 0; iRow < thisCopy.numberRows(); iRow++)
987 printf("%d at %c primal %g dual %g\n",
988 iRow, ss[thisCopy.getRowStatus(iRow)],
989 primalActs[iRow], primalDual[iRow]);
990 printf ("Primal problem column info %d columns\n", thisCopy.numberColumns());
991 for (iColumn = 0; iColumn < thisCopy.numberColumns(); iColumn++)
992 printf("%d at %c primal %g dual %g\n",
993 iColumn, ss[thisCopy.getColumnStatus(iColumn)],
994 primalSol[iColumn], primalDj[iColumn]);
995#endif
996 // position at bound information
997 int jColumn = numberRows_;
998 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
999 double objValue = optimizationDirection_ * objective[iColumn];
1000 Status status = dualProblem->getRowStatus(iColumn);
1001 double otherValue = COIN_DBL_MAX;
1002 if (columnUpper_[iColumn] < 1.0e20 &&
1003 columnLower_[iColumn] > -1.0e20) {
1004 if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) {
1005 otherValue = columnUpper_[iColumn] + dualDj[jColumn];
1006 } else {
1007 otherValue = columnLower_[iColumn] + dualDj[jColumn];
1008 }
1009 jColumn++;
1010 }
1011 if (status == basic) {
1012 // column is at bound
1013 if (otherValue == COIN_DBL_MAX) {
1014 reducedCost_[iColumn] = objValue - dualActs[iColumn];
1015 if (columnUpper_[iColumn] > 1.0e20) {
1016 if (columnLower_[iColumn] > -1.0e20) {
1017 if (columnUpper_[iColumn] > columnLower_[iColumn])
1018 setColumnStatus(iColumn, atLowerBound);
1019 else
1020 setColumnStatus(iColumn, isFixed);
1021 columnActivity_[iColumn] = columnLower_[iColumn];
1022 } else {
1023 // free
1024 setColumnStatus(iColumn, isFree);
1025 columnActivity_[iColumn] = 0.0;
1026 }
1027 } else {
1028 setColumnStatus(iColumn, atUpperBound);
1029 columnActivity_[iColumn] = columnUpper_[iColumn];
1030 }
1031 } else {
1032 reducedCost_[iColumn] = objValue - dualActs[iColumn];
1033 //printf("other dual sol %g\n",otherValue);
1034 if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1035 if (columnUpper_[iColumn] > columnLower_[iColumn])
1036 setColumnStatus(iColumn, atLowerBound);
1037 else
1038 setColumnStatus(iColumn, isFixed);
1039 columnActivity_[iColumn] = columnLower_[iColumn];
1040 } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1041 if (columnUpper_[iColumn] > columnLower_[iColumn])
1042 setColumnStatus(iColumn, atUpperBound);
1043 else
1044 setColumnStatus(iColumn, isFixed);
1045 columnActivity_[iColumn] = columnUpper_[iColumn];
1046 } else {
1047 abort();
1048 }
1049 }
1050 } else {
1051 if (otherValue == COIN_DBL_MAX) {
1052 // column basic
1053 setColumnStatus(iColumn, basic);
1054 numberBasic++;
1055 if (columnLower_[iColumn] > -1.0e20) {
1056 columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn];
1057 } else if (columnUpper_[iColumn] < 1.0e20) {
1058 columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn];
1059 } else {
1060 columnActivity_[iColumn] = -dualDual[iColumn];
1061 }
1062 reducedCost_[iColumn] = 0.0;
1063 } else {
1064 // may be at other bound
1065 //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1066 if (dualProblem->getColumnStatus(jColumn - 1) != basic) {
1067 // column basic
1068 setColumnStatus(iColumn, basic);
1069 numberBasic++;
1070 //printf("Col %d otherV %g dualDual %g\n",iColumn,
1071 // otherValue,dualDual[iColumn]);
1072 columnActivity_[iColumn] = -dualDual[iColumn];
1073 columnActivity_[iColumn] = otherValue;
1074 reducedCost_[iColumn] = 0.0;
1075 } else {
1076 reducedCost_[iColumn] = objValue - dualActs[iColumn];
1077 if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1078 if (columnUpper_[iColumn] > columnLower_[iColumn])
1079 setColumnStatus(iColumn, atLowerBound);
1080 else
1081 setColumnStatus(iColumn, isFixed);
1082 columnActivity_[iColumn] = columnLower_[iColumn];
1083 } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1084 if (columnUpper_[iColumn] > columnLower_[iColumn])
1085 setColumnStatus(iColumn, atUpperBound);
1086 else
1087 setColumnStatus(iColumn, isFixed);
1088 columnActivity_[iColumn] = columnUpper_[iColumn];
1089 } else {
1090 abort();
1091 }
1092 }
1093 }
1094 }
1095 }
1096 // now rows
1097 int kExtraRow = jColumn;
1098 int numberRanges = 0;
1099 for (iRow = 0; iRow < numberRows_; iRow++) {
1100 Status status = dualProblem->getColumnStatus(iRow);
1101 if (status == basic) {
1102 // row is at bound
1103 dual_[iRow] = dualSol[iRow];
1104 } else {
1105 // row basic
1106 setRowStatus(iRow, basic);
1107 numberBasic++;
1108 dual_[iRow] = 0.0;
1109 }
1110 if (rowLower_[iRow] < -1.0e20) {
1111 if (status == basic) {
1112 rowActivity_[iRow] = rowUpper_[iRow];
1113 setRowStatus(iRow, atUpperBound);
1114 } else {
1115 assert (dualDj[iRow] < 1.0e-5);
1116 rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
1117 }
1118 } else if (rowUpper_[iRow] > 1.0e20) {
1119 if (status == basic) {
1120 rowActivity_[iRow] = rowLower_[iRow];
1121 setRowStatus(iRow, atLowerBound);
1122 } else {
1123 rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
1124 assert (dualDj[iRow] > -1.0e-5);
1125 }
1126 } else {
1127 if (rowUpper_[iRow] == rowLower_[iRow]) {
1128 rowActivity_[iRow] = rowLower_[iRow];
1129 if (status == basic) {
1130 setRowStatus(iRow, isFixed);
1131 }
1132 } else {
1133 // range
1134 numberRanges++;
1135 Status statusL = dualProblem->getColumnStatus(kExtraRow);
1136 //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
1137 // iRow,status,kExtraRow,statusL, dualSol[iRow],
1138 // dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
1139 if (status == basic) {
1140 assert (statusL != basic);
1141 rowActivity_[iRow] = rowUpper_[iRow];
1142 setRowStatus(iRow, atUpperBound);
1143 } else if (statusL == basic) {
1144 numberBasic--; // already counted
1145 rowActivity_[iRow] = rowLower_[iRow];
1146 setRowStatus(iRow, atLowerBound);
1147 dual_[iRow] = dualSol[kExtraRow];
1148 } else {
1149 rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
1150 assert (dualDj[iRow] < 1.0e-5);
1151 // row basic
1152 //setRowStatus(iRow,basic);
1153 //numberBasic++;
1154 dual_[iRow] = 0.0;
1155 }
1156 kExtraRow++;
1157 }
1158 }
1159 }
1160 if (numberBasic != numberRows_) {
1161 printf("Bad basis - ranges - coding needed\n");
1162 assert (numberRanges);
1163 abort();
1164 }
1165 if (optimizationDirection_ < 0.0) {
1166 for (iRow = 0; iRow < numberRows_; iRow++) {
1167 dual_[iRow] = -dual_[iRow];
1168 }
1169 }
1170 // redo row activities
1171 memset(rowActivity_, 0, numberRows_ * sizeof(double));
1172 matrix_->times(1.0, columnActivity_, rowActivity_);
1173 // redo reduced costs
1174 memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
1175 matrix_->transposeTimes(-1.0, dual_, reducedCost_);
1176 checkSolutionInternal();
1177 if (sumDualInfeasibilities_ > 1.0e-5 || sumPrimalInfeasibilities_ > 1.0e-5) {
1178 returnCode = 1;
1179#ifdef CLP_INVESTIGATE
1180 printf("There are %d dual infeasibilities summing to %g ",
1181 numberDualInfeasibilities_, sumDualInfeasibilities_);
1182 printf("and %d primal infeasibilities summing to %g\n",
1183 numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
1184#endif
1185 }
1186 // Below will go to ..DEBUG later
1187#if 1 //ndef NDEBUG
1188 // Check if correct
1189 double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
1190 double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
1191 double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
1192 double * dual = CoinCopyOfArray(dual_, numberRows_);
1193 this->dual(); //primal();
1194 CoinRelFltEq eq(1.0e-5);
1195 for (iRow = 0; iRow < numberRows_; iRow++) {
1196 assert(eq(dual[iRow], dual_[iRow]));
1197 }
1198 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1199 assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
1200 }
1201 for (iRow = 0; iRow < numberRows_; iRow++) {
1202 assert(eq(rowActivity[iRow], rowActivity_[iRow]));
1203 }
1204 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1205 assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
1206 }
1207 delete [] columnActivity;
1208 delete [] rowActivity;
1209 delete [] reducedCost;
1210 delete [] dual;
1211#endif
1212 return returnCode;
1213}
1214/* Does very cursory presolve.
1215 rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
1216*/
1217ClpSimplex *
1218ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
1219 int & nBound, bool moreBounds, bool tightenBounds)
1220{
1221 //#define CHECK_STATUS
1222#ifdef CHECK_STATUS
1223 {
1224 int n = 0;
1225 int i;
1226 for (i = 0; i < numberColumns_; i++)
1227 if (getColumnStatus(i) == ClpSimplex::basic)
1228 n++;
1229 for (i = 0; i < numberRows_; i++)
1230 if (getRowStatus(i) == ClpSimplex::basic)
1231 n++;
1232 assert (n == numberRows_);
1233 }
1234#endif
1235
1236 const double * element = matrix_->getElements();
1237 const int * row = matrix_->getIndices();
1238 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1239 const int * columnLength = matrix_->getVectorLengths();
1240
1241 CoinZeroN(rhs, numberRows_);
1242 int iColumn;
1243 int iRow;
1244 CoinZeroN(whichRow, numberRows_);
1245 int * backColumn = whichColumn + numberColumns_;
1246 int numberRows2 = 0;
1247 int numberColumns2 = 0;
1248 double offset = 0.0;
1249 const double * objective = this->objective();
1250 double * solution = columnActivity_;
1251 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1252 double lower = columnLower_[iColumn];
1253 double upper = columnUpper_[iColumn];
1254 if (upper > lower || getColumnStatus(iColumn) == ClpSimplex::basic) {
1255 backColumn[iColumn] = numberColumns2;
1256 whichColumn[numberColumns2++] = iColumn;
1257 for (CoinBigIndex j = columnStart[iColumn];
1258 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1259 int iRow = row[j];
1260 int n = whichRow[iRow];
1261 if (n == 0 && element[j])
1262 whichRow[iRow] = -iColumn - 1;
1263 else if (n < 0)
1264 whichRow[iRow] = 2;
1265 }
1266 } else {
1267 // fixed
1268 backColumn[iColumn] = -1;
1269 solution[iColumn] = upper;
1270 if (upper) {
1271 offset += objective[iColumn] * upper;
1272 for (CoinBigIndex j = columnStart[iColumn];
1273 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1274 int iRow = row[j];
1275 double value = element[j];
1276 rhs[iRow] += upper * value;
1277 }
1278 }
1279 }
1280 }
1281 int returnCode = 0;
1282 double tolerance = primalTolerance();
1283 nBound = 2 * numberRows_;
1284 for (iRow = 0; iRow < numberRows_; iRow++) {
1285 int n = whichRow[iRow];
1286 if (n > 0) {
1287 whichRow[numberRows2++] = iRow;
1288 } else if (n < 0) {
1289 //whichRow[numberRows2++]=iRow;
1290 //continue;
1291 // Can only do in certain circumstances as we don't know current value
1292 if (rowLower_[iRow] == rowUpper_[iRow] || getRowStatus(iRow) == ClpSimplex::basic) {
1293 // save row and column for bound
1294 whichRow[--nBound] = iRow;
1295 whichRow[nBound+numberRows_] = -n - 1;
1296 } else if (moreBounds) {
1297 // save row and column for bound
1298 whichRow[--nBound] = iRow;
1299 whichRow[nBound+numberRows_] = -n - 1;
1300 } else {
1301 whichRow[numberRows2++] = iRow;
1302 }
1303 } else {
1304 // empty
1305 double rhsValue = rhs[iRow];
1306 if (rhsValue < rowLower_[iRow] - tolerance || rhsValue > rowUpper_[iRow] + tolerance) {
1307 returnCode = 1; // infeasible
1308 }
1309 }
1310 }
1311 ClpSimplex * small = NULL;
1312 if (!returnCode) {
1313 //printf("CRUNCH from (%d,%d) to (%d,%d)\n",
1314 // numberRows_,numberColumns_,numberRows2,numberColumns2);
1315 small = new ClpSimplex(this, numberRows2, whichRow,
1316 numberColumns2, whichColumn, true, false);
1317#if 0
1318 ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
1319 if (rowCopy) {
1320 assert(!small->rowCopy());
1321 small->setNewRowCopy(new ClpPackedMatrix(*rowCopy, numberRows2, whichRow,
1322 numberColumns2, whichColumn));
1323 }
1324#endif
1325 // Set some stuff
1326 small->setDualBound(dualBound_);
1327 small->setInfeasibilityCost(infeasibilityCost_);
1328 small->setSpecialOptions(specialOptions_);
1329 small->setPerturbation(perturbation_);
1330 small->defaultFactorizationFrequency();
1331 small->setAlphaAccuracy(alphaAccuracy_);
1332 // If no rows left then no tightening!
1333 if (!numberRows2 || !numberColumns2)
1334 tightenBounds = false;
1335
1336 int numberElements = getNumElements();
1337 int numberElements2 = small->getNumElements();
1338 small->setObjectiveOffset(objectiveOffset() - offset);
1339 handler_->message(CLP_CRUNCH_STATS, messages_)
1340 << numberRows2 << -(numberRows_ - numberRows2)
1341 << numberColumns2 << -(numberColumns_ - numberColumns2)
1342 << numberElements2 << -(numberElements - numberElements2)
1343 << CoinMessageEol;
1344 // And set objective value to match
1345 small->setObjectiveValue(this->objectiveValue());
1346 double * rowLower2 = small->rowLower();
1347 double * rowUpper2 = small->rowUpper();
1348 int jRow;
1349 for (jRow = 0; jRow < numberRows2; jRow++) {
1350 iRow = whichRow[jRow];
1351 if (rowLower2[jRow] > -1.0e20)
1352 rowLower2[jRow] -= rhs[iRow];
1353 if (rowUpper2[jRow] < 1.0e20)
1354 rowUpper2[jRow] -= rhs[iRow];
1355 }
1356 // and bounds
1357 double * columnLower2 = small->columnLower();
1358 double * columnUpper2 = small->columnUpper();
1359 const char * integerInformation = integerType_;
1360 for (jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1361 iRow = whichRow[jRow];
1362 iColumn = whichRow[jRow+numberRows_];
1363 double lowerRow = rowLower_[iRow];
1364 if (lowerRow > -1.0e20)
1365 lowerRow -= rhs[iRow];
1366 double upperRow = rowUpper_[iRow];
1367 if (upperRow < 1.0e20)
1368 upperRow -= rhs[iRow];
1369 int jColumn = backColumn[iColumn];
1370 double lower = columnLower2[jColumn];
1371 double upper = columnUpper2[jColumn];
1372 double value = 0.0;
1373 for (CoinBigIndex j = columnStart[iColumn];
1374 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1375 if (iRow == row[j]) {
1376 value = element[j];
1377 break;
1378 }
1379 }
1380 assert (value);
1381 // convert rowLower and Upper to implied bounds on column
1382 double newLower = -COIN_DBL_MAX;
1383 double newUpper = COIN_DBL_MAX;
1384 if (value > 0.0) {
1385 if (lowerRow > -1.0e20)
1386 newLower = lowerRow / value;
1387 if (upperRow < 1.0e20)
1388 newUpper = upperRow / value;
1389 } else {
1390 if (upperRow < 1.0e20)
1391 newLower = upperRow / value;
1392 if (lowerRow > -1.0e20)
1393 newUpper = lowerRow / value;
1394 }
1395 if (integerInformation && integerInformation[iColumn]) {
1396 if (newLower - floor(newLower) < 10.0 * tolerance)
1397 newLower = floor(newLower);
1398 else
1399 newLower = ceil(newLower);
1400 if (ceil(newUpper) - newUpper < 10.0 * tolerance)
1401 newUpper = ceil(newUpper);
1402 else
1403 newUpper = floor(newUpper);
1404 }
1405 newLower = CoinMax(lower, newLower);
1406 newUpper = CoinMin(upper, newUpper);
1407 if (newLower > newUpper + tolerance) {
1408 //printf("XXYY inf on bound\n");
1409 returnCode = 1;
1410 }
1411 columnLower2[jColumn] = newLower;
1412 columnUpper2[jColumn] = CoinMax(newLower, newUpper);
1413 if (getRowStatus(iRow) != ClpSimplex::basic) {
1414 if (getColumnStatus(iColumn) == ClpSimplex::basic) {
1415 if (columnLower2[jColumn] == columnUpper2[jColumn]) {
1416 // can only get here if will be fixed
1417 small->setColumnStatus(jColumn, ClpSimplex::isFixed);
1418 } else {
1419 // solution is valid
1420 if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) <
1421 fabs(columnActivity_[iColumn] - columnUpper2[jColumn]))
1422 small->setColumnStatus(jColumn, ClpSimplex::atLowerBound);
1423 else
1424 small->setColumnStatus(jColumn, ClpSimplex::atUpperBound);
1425 }
1426 } else {
1427 //printf("what now neither basic\n");
1428 }
1429 }
1430 }
1431 if (returnCode) {
1432 delete small;
1433 small = NULL;
1434 } else if (tightenBounds && integerInformation) {
1435 // See if we can tighten any bounds
1436 // use rhs for upper and small duals for lower
1437 double * up = rhs;
1438 double * lo = small->dualRowSolution();
1439 const double * element = small->clpMatrix()->getElements();
1440 const int * row = small->clpMatrix()->getIndices();
1441 const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1442 //const int * columnLength = small->clpMatrix()->getVectorLengths();
1443 CoinZeroN(lo, numberRows2);
1444 CoinZeroN(up, numberRows2);
1445 for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1446 double upper = columnUpper2[iColumn];
1447 double lower = columnLower2[iColumn];
1448 //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1449 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1450 int iRow = row[j];
1451 double value = element[j];
1452 if (value > 0.0) {
1453 if (upper < 1.0e20)
1454 up[iRow] += upper * value;
1455 else
1456 up[iRow] = COIN_DBL_MAX;
1457 if (lower > -1.0e20)
1458 lo[iRow] += lower * value;
1459 else
1460 lo[iRow] = -COIN_DBL_MAX;
1461 } else {
1462 if (upper < 1.0e20)
1463 lo[iRow] += upper * value;
1464 else
1465 lo[iRow] = -COIN_DBL_MAX;
1466 if (lower > -1.0e20)
1467 up[iRow] += lower * value;
1468 else
1469 up[iRow] = COIN_DBL_MAX;
1470 }
1471 }
1472 }
1473 double * rowLower2 = small->rowLower();
1474 double * rowUpper2 = small->rowUpper();
1475 bool feasible = true;
1476 // make safer
1477 for (int iRow = 0; iRow < numberRows2; iRow++) {
1478 double lower = lo[iRow];
1479 if (lower > rowUpper2[iRow] + tolerance) {
1480 feasible = false;
1481 break;
1482 } else {
1483 lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance;
1484 }
1485 double upper = up[iRow];
1486 if (upper < rowLower2[iRow] - tolerance) {
1487 feasible = false;
1488 break;
1489 } else {
1490 up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
1491 }
1492 }
1493 if (!feasible) {
1494 delete small;
1495 small = NULL;
1496 } else {
1497 // and tighten
1498 for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1499 if (integerInformation[whichColumn[iColumn]]) {
1500 double upper = columnUpper2[iColumn];
1501 double lower = columnLower2[iColumn];
1502 double newUpper = upper;
1503 double newLower = lower;
1504 double difference = upper - lower;
1505 if (lower > -1000.0 && upper < 1000.0) {
1506 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1507 int iRow = row[j];
1508 double value = element[j];
1509 if (value > 0.0) {
1510 double upWithOut = up[iRow] - value * difference;
1511 if (upWithOut < 0.0) {
1512 newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1513 }
1514 double lowWithOut = lo[iRow] + value * difference;
1515 if (lowWithOut > 0.0) {
1516 newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1517 }
1518 } else {
1519 double upWithOut = up[iRow] + value * difference;
1520 if (upWithOut < 0.0) {
1521 newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1522 }
1523 double lowWithOut = lo[iRow] - value * difference;
1524 if (lowWithOut > 0.0) {
1525 newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1526 }
1527 }
1528 }
1529 if (newLower > lower || newUpper < upper) {
1530 if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1531 newUpper = floor(newUpper);
1532 else
1533 newUpper = floor(newUpper + 0.5);
1534 if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1535 newLower = ceil(newLower);
1536 else
1537 newLower = ceil(newLower - 0.5);
1538 // change may be too small - check
1539 if (newLower > lower || newUpper < upper) {
1540 if (newUpper >= newLower) {
1541 // Could also tighten in this
1542 //printf("%d bounds %g %g tightened to %g %g\n",
1543 // iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1544 // newLower,newUpper);
1545#if 1
1546 columnUpper2[iColumn] = newUpper;
1547 columnLower2[iColumn] = newLower;
1548 columnUpper_[whichColumn[iColumn]] = newUpper;
1549 columnLower_[whichColumn[iColumn]] = newLower;
1550#endif
1551 // and adjust bounds on rows
1552 newUpper -= upper;
1553 newLower -= lower;
1554 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1555 int iRow = row[j];
1556 double value = element[j];
1557 if (value > 0.0) {
1558 up[iRow] += newUpper * value;
1559 lo[iRow] += newLower * value;
1560 } else {
1561 lo[iRow] += newUpper * value;
1562 up[iRow] += newLower * value;
1563 }
1564 }
1565 } else {
1566 // infeasible
1567 //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1568 // iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1569 // newLower,newUpper);
1570#if 1
1571 delete small;
1572 small = NULL;
1573 break;
1574#endif
1575 }
1576 }
1577 }
1578 }
1579 }
1580 }
1581 }
1582 }
1583 }
1584#if 0
1585 if (small) {
1586 static int which = 0;
1587 which++;
1588 char xxxx[20];
1589 sprintf(xxxx, "bad%d.mps", which);
1590 small->writeMps(xxxx, 0, 1);
1591 sprintf(xxxx, "largebad%d.mps", which);
1592 writeMps(xxxx, 0, 1);
1593 printf("bad%d %x old size %d %d new %d %d\n", which, small,
1594 numberRows_, numberColumns_, small->numberRows(), small->numberColumns());
1595#if 0
1596 for (int i = 0; i < numberColumns_; i++)
1597 printf("Bound %d %g %g\n", i, columnLower_[i], columnUpper_[i]);
1598 for (int i = 0; i < numberRows_; i++)
1599 printf("Row bound %d %g %g\n", i, rowLower_[i], rowUpper_[i]);
1600#endif
1601 }
1602#endif
1603#ifdef CHECK_STATUS
1604 {
1605 int n = 0;
1606 int i;
1607 for (i = 0; i < small->numberColumns(); i++)
1608 if (small->getColumnStatus(i) == ClpSimplex::basic)
1609 n++;
1610 for (i = 0; i < small->numberRows(); i++)
1611 if (small->getRowStatus(i) == ClpSimplex::basic)
1612 n++;
1613 assert (n == small->numberRows());
1614 }
1615#endif
1616 return small;
1617}
1618/* After very cursory presolve.
1619 rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
1620*/
1621void
1622ClpSimplexOther::afterCrunch(const ClpSimplex & small,
1623 const int * whichRow,
1624 const int * whichColumn, int nBound)
1625{
1626#ifndef NDEBUG
1627 for (int i = 0; i < small.numberRows(); i++)
1628 assert (whichRow[i] >= 0 && whichRow[i] < numberRows_);
1629 for (int i = 0; i < small.numberColumns(); i++)
1630 assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_);
1631#endif
1632 getbackSolution(small, whichRow, whichColumn);
1633 // and deal with status for bounds
1634 const double * element = matrix_->getElements();
1635 const int * row = matrix_->getIndices();
1636 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1637 const int * columnLength = matrix_->getVectorLengths();
1638 double tolerance = primalTolerance();
1639 double djTolerance = dualTolerance();
1640 for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1641 int iRow = whichRow[jRow];
1642 int iColumn = whichRow[jRow+numberRows_];
1643 if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1644 double lower = columnLower_[iColumn];
1645 double upper = columnUpper_[iColumn];
1646 double value = columnActivity_[iColumn];
1647 double djValue = reducedCost_[iColumn];
1648 dual_[iRow] = 0.0;
1649 if (upper > lower) {
1650 if (value < lower + tolerance && djValue > -djTolerance) {
1651 setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1652 setRowStatus(iRow, ClpSimplex::basic);
1653 } else if (value > upper - tolerance && djValue < djTolerance) {
1654 setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1655 setRowStatus(iRow, ClpSimplex::basic);
1656 } else {
1657 // has to be basic
1658 setColumnStatus(iColumn, ClpSimplex::basic);
1659 reducedCost_[iColumn] = 0.0;
1660 double value = 0.0;
1661 for (CoinBigIndex j = columnStart[iColumn];
1662 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1663 if (iRow == row[j]) {
1664 value = element[j];
1665 break;
1666 }
1667 }
1668 dual_[iRow] = djValue / value;
1669 if (rowUpper_[iRow] > rowLower_[iRow]) {
1670 if (fabs(rowActivity_[iRow] - rowLower_[iRow]) <
1671 fabs(rowActivity_[iRow] - rowUpper_[iRow]))
1672 setRowStatus(iRow, ClpSimplex::atLowerBound);
1673 else
1674 setRowStatus(iRow, ClpSimplex::atUpperBound);
1675 } else {
1676 setRowStatus(iRow, ClpSimplex::isFixed);
1677 }
1678 }
1679 } else {
1680 // row can always be basic
1681 setRowStatus(iRow, ClpSimplex::basic);
1682 }
1683 } else {
1684 // row can always be basic
1685 setRowStatus(iRow, ClpSimplex::basic);
1686 }
1687 }
1688 //#ifndef NDEBUG
1689#if 0
1690 if (small.status() == 0) {
1691 int n = 0;
1692 int i;
1693 for (i = 0; i < numberColumns; i++)
1694 if (getColumnStatus(i) == ClpSimplex::basic)
1695 n++;
1696 for (i = 0; i < numberRows; i++)
1697 if (getRowStatus(i) == ClpSimplex::basic)
1698 n++;
1699 assert (n == numberRows);
1700 }
1701#endif
1702}
1703/* Tightens integer bounds - returns number tightened or -1 if infeasible
1704 */
1705int
1706ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1707{
1708 // See if we can tighten any bounds
1709 // use rhs for upper and small duals for lower
1710 double * up = rhsSpace;
1711 double * lo = dual_;
1712 const double * element = matrix_->getElements();
1713 const int * row = matrix_->getIndices();
1714 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1715 const int * columnLength = matrix_->getVectorLengths();
1716 CoinZeroN(lo, numberRows_);
1717 CoinZeroN(up, numberRows_);
1718 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1719 double upper = columnUpper_[iColumn];
1720 double lower = columnLower_[iColumn];
1721 //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1722 for (CoinBigIndex j = columnStart[iColumn];
1723 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1724 int iRow = row[j];
1725 double value = element[j];
1726 if (value > 0.0) {
1727 if (upper < 1.0e20)
1728 up[iRow] += upper * value;
1729 else
1730 up[iRow] = COIN_DBL_MAX;
1731 if (lower > -1.0e20)
1732 lo[iRow] += lower * value;
1733 else
1734 lo[iRow] = -COIN_DBL_MAX;
1735 } else {
1736 if (upper < 1.0e20)
1737 lo[iRow] += upper * value;
1738 else
1739 lo[iRow] = -COIN_DBL_MAX;
1740 if (lower > -1.0e20)
1741 up[iRow] += lower * value;
1742 else
1743 up[iRow] = COIN_DBL_MAX;
1744 }
1745 }
1746 }
1747 bool feasible = true;
1748 // make safer
1749 double tolerance = primalTolerance();
1750 for (int iRow = 0; iRow < numberRows_; iRow++) {
1751 double lower = lo[iRow];
1752 if (lower > rowUpper_[iRow] + tolerance) {
1753 feasible = false;
1754 break;
1755 } else {
1756 lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance;
1757 }
1758 double upper = up[iRow];
1759 if (upper < rowLower_[iRow] - tolerance) {
1760 feasible = false;
1761 break;
1762 } else {
1763 up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance;
1764 }
1765 }
1766 int numberTightened = 0;
1767 if (!feasible) {
1768 return -1;
1769 } else if (integerType_) {
1770 // and tighten
1771 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1772 if (integerType_[iColumn]) {
1773 double upper = columnUpper_[iColumn];
1774 double lower = columnLower_[iColumn];
1775 double newUpper = upper;
1776 double newLower = lower;
1777 double difference = upper - lower;
1778 if (lower > -1000.0 && upper < 1000.0) {
1779 for (CoinBigIndex j = columnStart[iColumn];
1780 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1781 int iRow = row[j];
1782 double value = element[j];
1783 if (value > 0.0) {
1784 double upWithOut = up[iRow] - value * difference;
1785 if (upWithOut < 0.0) {
1786 newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1787 }
1788 double lowWithOut = lo[iRow] + value * difference;
1789 if (lowWithOut > 0.0) {
1790 newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1791 }
1792 } else {
1793 double upWithOut = up[iRow] + value * difference;
1794 if (upWithOut < 0.0) {
1795 newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1796 }
1797 double lowWithOut = lo[iRow] - value * difference;
1798 if (lowWithOut > 0.0) {
1799 newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1800 }
1801 }
1802 }
1803 if (newLower > lower || newUpper < upper) {
1804 if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1805 newUpper = floor(newUpper);
1806 else
1807 newUpper = floor(newUpper + 0.5);
1808 if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1809 newLower = ceil(newLower);
1810 else
1811 newLower = ceil(newLower - 0.5);
1812 // change may be too small - check
1813 if (newLower > lower || newUpper < upper) {
1814 if (newUpper >= newLower) {
1815 numberTightened++;
1816 //printf("%d bounds %g %g tightened to %g %g\n",
1817 // iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1818 // newLower,newUpper);
1819 columnUpper_[iColumn] = newUpper;
1820 columnLower_[iColumn] = newLower;
1821 // and adjust bounds on rows
1822 newUpper -= upper;
1823 newLower -= lower;
1824 for (CoinBigIndex j = columnStart[iColumn];
1825 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1826 int iRow = row[j];
1827 double value = element[j];
1828 if (value > 0.0) {
1829 up[iRow] += newUpper * value;
1830 lo[iRow] += newLower * value;
1831 } else {
1832 lo[iRow] += newUpper * value;
1833 up[iRow] += newLower * value;
1834 }
1835 }
1836 } else {
1837 // infeasible
1838 //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1839 // iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1840 // newLower,newUpper);
1841 return -1;
1842 }
1843 }
1844 }
1845 }
1846 }
1847 }
1848 }
1849 return numberTightened;
1850}
1851/* Parametrics
1852 This is an initial slow version.
1853 The code uses current bounds + theta * change (if change array not NULL)
1854 and similarly for objective.
1855 It starts at startingTheta and returns ending theta in endingTheta.
1856 If reportIncrement 0.0 it will report on any movement
1857 If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1858 If it can not reach input endingTheta return code will be 1 for infeasible,
1859 2 for unbounded, if error on ranges -1, otherwise 0.
1860 Normal report is just theta and objective but
1861 if event handler exists it may do more
1862 On exit endingTheta is maximum reached (can be used for next startingTheta)
1863*/
1864int
1865ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement,
1866 const double * changeLowerBound, const double * changeUpperBound,
1867 const double * changeLowerRhs, const double * changeUpperRhs,
1868 const double * changeObjective)
1869{
1870 bool needToDoSomething = true;
1871 bool canTryQuick = (reportIncrement) ? true : false;
1872 // Save copy of model
1873 ClpSimplex copyModel = *this;
1874 int savePerturbation = perturbation_;
1875 perturbation_ = 102; // switch off
1876 while (needToDoSomething) {
1877 needToDoSomething = false;
1878 algorithm_ = -1;
1879
1880 // save data
1881 ClpDataSave data = saveData();
1882 int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
1883 int iRow, iColumn;
1884 double * chgUpper = NULL;
1885 double * chgLower = NULL;
1886 double * chgObjective = NULL;
1887
1888 // Dantzig (as will not be used) (out later)
1889 ClpDualRowPivot * savePivot = dualRowPivot_;
1890 dualRowPivot_ = new ClpDualRowDantzig();
1891
1892 if (!returnCode) {
1893 // Find theta when bounds will cross over and create arrays
1894 int numberTotal = numberRows_ + numberColumns_;
1895 chgLower = new double[numberTotal];
1896 memset(chgLower, 0, numberTotal * sizeof(double));
1897 chgUpper = new double[numberTotal];
1898 memset(chgUpper, 0, numberTotal * sizeof(double));
1899 chgObjective = new double[numberTotal];
1900 memset(chgObjective, 0, numberTotal * sizeof(double));
1901 assert (!rowScale_);
1902 double maxTheta = 1.0e50;
1903 if (changeLowerRhs || changeUpperRhs) {
1904 for (iRow = 0; iRow < numberRows_; iRow++) {
1905 double lower = rowLower_[iRow];
1906 double upper = rowUpper_[iRow];
1907 if (lower > upper) {
1908 maxTheta = -1.0;
1909 break;
1910 }
1911 double changeLower = (changeLowerRhs) ? changeLowerRhs[iRow] : 0.0;
1912 double changeUpper = (changeUpperRhs) ? changeUpperRhs[iRow] : 0.0;
1913 if (lower > -1.0e20 && upper < 1.0e20) {
1914 if (lower + maxTheta * changeLower > upper + maxTheta * changeUpper) {
1915 maxTheta = (upper - lower) / (changeLower - changeUpper);
1916 }
1917 }
1918 if (lower > -1.0e20) {
1919 lower_[numberColumns_+iRow] += startingTheta * changeLower;
1920 chgLower[numberColumns_+iRow] = changeLower;
1921 }
1922 if (upper < 1.0e20) {
1923 upper_[numberColumns_+iRow] += startingTheta * changeUpper;
1924 chgUpper[numberColumns_+iRow] = changeUpper;
1925 }
1926 }
1927 }
1928 if (maxTheta > 0.0) {
1929 if (changeLowerBound || changeUpperBound) {
1930 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1931 double lower = columnLower_[iColumn];
1932 double upper = columnUpper_[iColumn];
1933 if (lower > upper) {
1934 maxTheta = -1.0;
1935 break;
1936 }
1937 double changeLower = (changeLowerBound) ? changeLowerBound[iColumn] : 0.0;
1938 double changeUpper = (changeUpperBound) ? changeUpperBound[iColumn] : 0.0;
1939 if (lower > -1.0e20 && upper < 1.0e20) {
1940 if (lower + maxTheta * changeLower > upper + maxTheta * changeUpper) {
1941 maxTheta = (upper - lower) / (changeLower - changeUpper);
1942 }
1943 }
1944 if (lower > -1.0e20) {
1945 lower_[iColumn] += startingTheta * changeLower;
1946 chgLower[iColumn] = changeLower;
1947 }
1948 if (upper < 1.0e20) {
1949 upper_[iColumn] += startingTheta * changeUpper;
1950 chgUpper[iColumn] = changeUpper;
1951 }
1952 }
1953 }
1954 if (maxTheta == 1.0e50)
1955 maxTheta = COIN_DBL_MAX;
1956 }
1957 if (maxTheta < 0.0) {
1958 // bad ranges or initial
1959 returnCode = -1;
1960 }
1961 if (maxTheta < endingTheta) {
1962 char line[100];
1963 sprintf(line,"Crossover considerations reduce ending theta from %g to %g\n",
1964 endingTheta,maxTheta);
1965 handler_->message(CLP_GENERAL,messages_)
1966 << line << CoinMessageEol;
1967 endingTheta = maxTheta;
1968 }
1969 if (endingTheta < startingTheta) {
1970 // bad initial
1971 returnCode = -2;
1972 }
1973 }
1974 double saveEndingTheta = endingTheta;
1975 if (!returnCode) {
1976 if (changeObjective) {
1977 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1978 chgObjective[iColumn] = changeObjective[iColumn];
1979 cost_[iColumn] += startingTheta * changeObjective[iColumn];
1980 }
1981 }
1982 double * saveDuals = NULL;
1983 reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
1984 assert (!problemStatus_);
1985 // Now do parametrics
1986 handler_->message(CLP_PARAMETRICS_STATS, messages_)
1987 << startingTheta << objectiveValue() << CoinMessageEol;
1988 while (!returnCode) {
1989 //assert (reportIncrement);
1990 returnCode = parametricsLoop(startingTheta, endingTheta, reportIncrement,
1991 chgLower, chgUpper, chgObjective, data,
1992 canTryQuick);
1993 if (!returnCode) {
1994 //double change = endingTheta-startingTheta;
1995 startingTheta = endingTheta;
1996 endingTheta = saveEndingTheta;
1997 //for (int i=0;i<numberTotal;i++) {
1998 //lower_[i] += change*chgLower[i];
1999 //upper_[i] += change*chgUpper[i];
2000 //cost_[i] += change*chgObjective[i];
2001 //}
2002 handler_->message(CLP_PARAMETRICS_STATS, messages_)
2003 << startingTheta << objectiveValue() << CoinMessageEol;
2004 if (startingTheta >= endingTheta)
2005 break;
2006 } else if (returnCode == -1) {
2007 // trouble - do external solve
2008 needToDoSomething = true;
2009 } else if (problemStatus_==1) {
2010 // can't move any further
2011 if (!canTryQuick) {
2012 handler_->message(CLP_PARAMETRICS_STATS, messages_)
2013 << endingTheta << objectiveValue() << CoinMessageEol;
2014 problemStatus_=0;
2015 }
2016 } else {
2017 abort();
2018 }
2019 }
2020 }
2021 reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
2022
2023 delete dualRowPivot_;
2024 dualRowPivot_ = savePivot;
2025 // Restore any saved stuff
2026 restoreData(data);
2027 if (needToDoSomething) {
2028 double saveStartingTheta = startingTheta; // known to be feasible
2029 int cleanedUp = 1;
2030 while (cleanedUp) {
2031 // tweak
2032 if (cleanedUp == 1) {
2033 if (!reportIncrement)
2034 startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
2035 else
2036 startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
2037 } else {
2038 // restoring to go slowly
2039 startingTheta = saveStartingTheta;
2040 }
2041 // only works if not scaled
2042 int i;
2043 const double * obj1 = objective();
2044 double * obj2 = copyModel.objective();
2045 const double * lower1 = columnLower_;
2046 double * lower2 = copyModel.columnLower();
2047 const double * upper1 = columnUpper_;
2048 double * upper2 = copyModel.columnUpper();
2049 for (i = 0; i < numberColumns_; i++) {
2050 obj2[i] = obj1[i] + startingTheta * chgObjective[i];
2051 lower2[i] = lower1[i] + startingTheta * chgLower[i];
2052 upper2[i] = upper1[i] + startingTheta * chgUpper[i];
2053 }
2054 lower1 = rowLower_;
2055 lower2 = copyModel.rowLower();
2056 upper1 = rowUpper_;
2057 upper2 = copyModel.rowUpper();
2058 for (i = 0; i < numberRows_; i++) {
2059 lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
2060 upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
2061 }
2062 copyModel.dual();
2063 if (copyModel.problemStatus()) {
2064 char line[100];
2065 sprintf(line,"Can not get to theta of %g\n", startingTheta);
2066 handler_->message(CLP_GENERAL,messages_)
2067 << line << CoinMessageEol;
2068 canTryQuick = false; // do slowly to get exact amount
2069 // back to last known good
2070 if (cleanedUp == 1)
2071 cleanedUp = 2;
2072 else
2073 abort();
2074 } else {
2075 // and move stuff back
2076 int numberTotal = numberRows_ + numberColumns_;
2077 CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
2078 CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
2079 CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
2080 cleanedUp = 0;
2081 }
2082 }
2083 }
2084 delete [] chgLower;
2085 delete [] chgUpper;
2086 delete [] chgObjective;
2087 }
2088 perturbation_ = savePerturbation;
2089 char line[100];
2090 sprintf(line,"Ending theta %g\n", endingTheta);
2091 handler_->message(CLP_GENERAL,messages_)
2092 << line << CoinMessageEol;
2093 return problemStatus_;
2094}
2095/* Version of parametrics which reads from file
2096 See CbcClpParam.cpp for details of format
2097 Returns -2 if unable to open file */
2098int
2099ClpSimplexOther::parametrics(const char * dataFile)
2100{
2101 int returnCode=-2;
2102 FILE *fp = fopen(dataFile, "r");
2103 char line[200];
2104 if (!fp) {
2105 handler_->message(CLP_UNABLE_OPEN, messages_)
2106 << dataFile << CoinMessageEol;
2107 return -2;
2108 }
2109
2110 if (!fgets(line, 200, fp)) {
2111 sprintf(line,"Empty parametrics file %s?",dataFile);
2112 handler_->message(CLP_GENERAL,messages_)
2113 << line << CoinMessageEol;
2114 fclose(fp);
2115 return -2;
2116 }
2117 char * pos = line;
2118 char * put = line;
2119 while (*pos >= ' ' && *pos != '\n') {
2120 if (*pos != ' ' && *pos != '\t') {
2121 *put = static_cast<char>(tolower(*pos));
2122 put++;
2123 }
2124 pos++;
2125 }
2126 *put = '\0';
2127 pos = line;
2128 double startTheta=0.0;
2129 double endTheta=0.0;
2130 double intervalTheta=COIN_DBL_MAX;
2131 int detail=0;
2132 bool good = true;
2133 while (good) {
2134 good=false;
2135 // check ROWS
2136 char * comma = strchr(pos, ',');
2137 if (!comma)
2138 break;
2139 *comma = '\0';
2140 if (strcmp(pos,"rows"))
2141 break;
2142 *comma = ',';
2143 pos = comma+1;
2144 // check lower theta
2145 comma = strchr(pos, ',');
2146 if (!comma)
2147 break;
2148 *comma = '\0';
2149 startTheta = atof(pos);
2150 *comma = ',';
2151 pos = comma+1;
2152 // check upper theta
2153 comma = strchr(pos, ',');
2154 good=true;
2155 if (comma)
2156 *comma = '\0';
2157 endTheta = atof(pos);
2158 if (comma) {
2159 *comma = ',';
2160 pos = comma+1;
2161 comma = strchr(pos, ',');
2162 if (comma)
2163 *comma = '\0';
2164 intervalTheta = atof(pos);
2165 if (comma) {
2166 *comma = ',';
2167 pos = comma+1;
2168 comma = strchr(pos, ',');
2169 if (comma)
2170 *comma = '\0';
2171 detail = atoi(pos);
2172 if (comma)
2173 *comma = ',';
2174 }
2175 }
2176 break;
2177 }
2178 if (good) {
2179 if (startTheta<0.0||
2180 startTheta>endTheta||
2181 intervalTheta<0.0)
2182 good=false;
2183 if (detail<0||detail>1)
2184 good=false;
2185 }
2186 if (intervalTheta>=endTheta)
2187 intervalTheta=0.0;
2188 if (!good) {
2189 sprintf(line,"Odd first line %s on file %s?",line,dataFile);
2190 handler_->message(CLP_GENERAL,messages_)
2191 << line << CoinMessageEol;
2192 fclose(fp);
2193 return -2;
2194 }
2195 if (!fgets(line, 200, fp)) {
2196 sprintf(line,"Not enough records on parametrics file %s?",dataFile);
2197 handler_->message(CLP_GENERAL,messages_)
2198 << line << CoinMessageEol;
2199 fclose(fp);
2200 return -2;
2201 }
2202 double * lowerRowMove = NULL;
2203 double * upperRowMove = NULL;
2204 double * lowerColumnMove = NULL;
2205 double * upperColumnMove = NULL;
2206 double * objectiveMove = NULL;
2207 char saveLine[200];
2208 saveLine[0]='\0';
2209 std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
2210 int gotRow[] = { -1, -1, -1, -1, -1};
2211 int orderRow[5];
2212 assert(sizeof(gotRow) == sizeof(orderRow));
2213 int nAcross = 0;
2214 pos = line;
2215 put = line;
2216 while (*pos >= ' ' && *pos != '\n') {
2217 if (*pos != ' ' && *pos != '\t') {
2218 *put = static_cast<char>(tolower(*pos));
2219 put++;
2220 }
2221 pos++;
2222 }
2223 *put = '\0';
2224 pos = line;
2225 int i;
2226 good = true;
2227 if (strncmp(line,"column",6)) {
2228 while (pos) {
2229 char * comma = strchr(pos, ',');
2230 if (comma)
2231 *comma = '\0';
2232 for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
2233 if (headingsRow[i] == pos) {
2234 if (gotRow[i] < 0) {
2235 orderRow[nAcross] = i;
2236 gotRow[i] = nAcross++;
2237 } else {
2238 // duplicate
2239 good = false;
2240 }
2241 break;
2242 }
2243 }
2244 if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
2245 good = false;
2246 if (comma) {
2247 *comma = ',';
2248 pos = comma + 1;
2249 } else {
2250 break;
2251 }
2252 }
2253 if (gotRow[0] < 0 && gotRow[1] < 0)
2254 good = false;
2255 if (gotRow[0] >= 0 && gotRow[1] >= 0)
2256 good = false;
2257 if (gotRow[0] >= 0 && !lengthNames())
2258 good = false;
2259 if (gotRow[4]<0) {
2260 if (gotRow[2] < 0 && gotRow[3] >= 0)
2261 good = false;
2262 else if (gotRow[3] < 0 && gotRow[2] >= 0)
2263 good = false;
2264 } else if (gotRow[2]>=0||gotRow[3]>=0) {
2265 good = false;
2266 }
2267 if (good) {
2268 char ** rowNames = new char * [numberRows_];
2269 int iRow;
2270 for (iRow = 0; iRow < numberRows_; iRow++) {
2271 rowNames[iRow] =
2272 CoinStrdup(rowName(iRow).c_str());
2273 }
2274 lowerRowMove = new double [numberRows_];
2275 memset(lowerRowMove,0,numberRows_*sizeof(double));
2276 upperRowMove = new double [numberRows_];
2277 memset(upperRowMove,0,numberRows_*sizeof(double));
2278 int nLine = 0;
2279 int nBadLine = 0;
2280 int nBadName = 0;
2281 bool goodLine=false;
2282 while (fgets(line, 200, fp)) {
2283 goodLine=true;
2284 if (!strncmp(line, "ENDATA", 6)||
2285 !strncmp(line, "COLUMN",6))
2286 break;
2287 goodLine=false;
2288 nLine++;
2289 iRow = -1;
2290 double upper = 0.0;
2291 double lower = 0.0;
2292 char * pos = line;
2293 char * put = line;
2294 while (*pos >= ' ' && *pos != '\n') {
2295 if (*pos != ' ' && *pos != '\t') {
2296 *put = *pos;
2297 put++;
2298 }
2299 pos++;
2300 }
2301 *put = '\0';
2302 pos = line;
2303 for (int i = 0; i < nAcross; i++) {
2304 char * comma = strchr(pos, ',');
2305 if (comma) {
2306 *comma = '\0';
2307 } else if (i < nAcross - 1) {
2308 nBadLine++;
2309 break;
2310 }
2311 switch (orderRow[i]) {
2312 // name
2313 case 0:
2314 // For large problems this could be slow
2315 for (iRow = 0; iRow < numberRows_; iRow++) {
2316 if (!strcmp(rowNames[iRow], pos))
2317 break;
2318 }
2319 if (iRow == numberRows_)
2320 iRow = -1;
2321 break;
2322 // number
2323 case 1:
2324 iRow = atoi(pos);
2325 if (iRow < 0 || iRow >= numberRows_)
2326 iRow = -1;
2327 break;
2328 // lower
2329 case 2:
2330 upper = atof(pos);
2331 break;
2332 // upper
2333 case 3:
2334 lower = atof(pos);
2335 break;
2336 // rhs
2337 case 4:
2338 lower = atof(pos);
2339 upper = lower;
2340 break;
2341 }
2342 if (comma) {
2343 *comma = ',';
2344 pos = comma + 1;
2345 }
2346 }
2347 if (iRow >= 0) {
2348 if (rowLower_[iRow]>-1.0e20)
2349 lowerRowMove[iRow] = lower;
2350 else
2351 lowerRowMove[iRow]=0.0;
2352 if (rowUpper_[iRow]<1.0e20)
2353 upperRowMove[iRow] = upper;
2354 else
2355 upperRowMove[iRow] = lower;
2356 } else {
2357 nBadName++;
2358 if(saveLine[0]=='\0')
2359 strcpy(saveLine,line);
2360 }
2361 }
2362 sprintf(line,"%d Row fields and %d records", nAcross, nLine);
2363 handler_->message(CLP_GENERAL,messages_)
2364 << line << CoinMessageEol;
2365 if (nBadName) {
2366 sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2367 handler_->message(CLP_GENERAL,messages_)
2368 << line << CoinMessageEol;
2369 returnCode=-1;
2370 good=false;
2371 }
2372 for (iRow = 0; iRow < numberRows_; iRow++) {
2373 free(rowNames[iRow]);
2374 }
2375 delete [] rowNames;
2376 } else {
2377 sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2378 handler_->message(CLP_GENERAL,messages_)
2379 << line << CoinMessageEol;
2380 returnCode=-1;
2381 good=false;
2382 }
2383 }
2384 if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) {
2385 if (!fgets(line, 200, fp)) {
2386 sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
2387 handler_->message(CLP_GENERAL,messages_)
2388 << line << CoinMessageEol;
2389 fclose(fp);
2390 return -2;
2391 }
2392 std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
2393 saveLine[0]='\0';
2394 int gotColumn[] = { -1, -1, -1, -1, -1};
2395 int orderColumn[5];
2396 assert(sizeof(gotColumn) == sizeof(orderColumn));
2397 nAcross = 0;
2398 pos = line;
2399 put = line;
2400 while (*pos >= ' ' && *pos != '\n') {
2401 if (*pos != ' ' && *pos != '\t') {
2402 *put = static_cast<char>(tolower(*pos));
2403 put++;
2404 }
2405 pos++;
2406 }
2407 *put = '\0';
2408 pos = line;
2409 int i;
2410 if (strncmp(line,"endata",6)&&good) {
2411 while (pos) {
2412 char * comma = strchr(pos, ',');
2413 if (comma)
2414 *comma = '\0';
2415 for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
2416 if (headingsColumn[i] == pos) {
2417 if (gotColumn[i] < 0) {
2418 orderColumn[nAcross] = i;
2419 gotColumn[i] = nAcross++;
2420 } else {
2421 // duplicate
2422 good = false;
2423 }
2424 break;
2425 }
2426 }
2427 if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
2428 good = false;
2429 if (comma) {
2430 *comma = ',';
2431 pos = comma + 1;
2432 } else {
2433 break;
2434 }
2435 }
2436 if (gotColumn[0] < 0 && gotColumn[1] < 0)
2437 good = false;
2438 if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
2439 good = false;
2440 if (gotColumn[0] >= 0 && !lengthNames())
2441 good = false;
2442 if (good) {
2443 char ** columnNames = new char * [numberColumns_];
2444 int iColumn;
2445 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2446 columnNames[iColumn] =
2447 CoinStrdup(columnName(iColumn).c_str());
2448 }
2449 lowerColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2450 memset(lowerColumnMove,0,numberColumns_*sizeof(double));
2451 upperColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2452 memset(upperColumnMove,0,numberColumns_*sizeof(double));
2453 objectiveMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2454 memset(objectiveMove,0,numberColumns_*sizeof(double));
2455 int nLine = 0;
2456 int nBadLine = 0;
2457 int nBadName = 0;
2458 bool goodLine=false;
2459 while (fgets(line, 200, fp)) {
2460 goodLine=true;
2461 if (!strncmp(line, "ENDATA", 6))
2462 break;
2463 goodLine=false;
2464 nLine++;
2465 iColumn = -1;
2466 double upper = 0.0;
2467 double lower = 0.0;
2468 double obj =0.0;
2469 char * pos = line;
2470 char * put = line;
2471 while (*pos >= ' ' && *pos != '\n') {
2472 if (*pos != ' ' && *pos != '\t') {
2473 *put = *pos;
2474 put++;
2475 }
2476 pos++;
2477 }
2478 *put = '\0';
2479 pos = line;
2480 for (int i = 0; i < nAcross; i++) {
2481 char * comma = strchr(pos, ',');
2482 if (comma) {
2483 *comma = '\0';
2484 } else if (i < nAcross - 1) {
2485 nBadLine++;
2486 break;
2487 }
2488 switch (orderColumn[i]) {
2489 // name
2490 case 0:
2491 // For large problems this could be slow
2492 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2493 if (!strcmp(columnNames[iColumn], pos))
2494 break;
2495 }
2496 if (iColumn == numberColumns_)
2497 iColumn = -1;
2498 break;
2499 // number
2500 case 1:
2501 iColumn = atoi(pos);
2502 if (iColumn < 0 || iColumn >= numberColumns_)
2503 iColumn = -1;
2504 break;
2505 // lower
2506 case 2:
2507 upper = atof(pos);
2508 break;
2509 // upper
2510 case 3:
2511 lower = atof(pos);
2512 break;
2513 // objective
2514 case 4:
2515 obj = atof(pos);
2516 upper = lower;
2517 break;
2518 }
2519 if (comma) {
2520 *comma = ',';
2521 pos = comma + 1;
2522 }
2523 }
2524 if (iColumn >= 0) {
2525 if (columnLower_[iColumn]>-1.0e20)
2526 lowerColumnMove[iColumn] = lower;
2527 else
2528 lowerColumnMove[iColumn]=0.0;
2529 if (columnUpper_[iColumn]<1.0e20)
2530 upperColumnMove[iColumn] = upper;
2531 else
2532 upperColumnMove[iColumn] = lower;
2533 objectiveMove[iColumn] = obj;
2534 } else {
2535 nBadName++;
2536 if(saveLine[0]=='\0')
2537 strcpy(saveLine,line);
2538 }
2539 }
2540 sprintf(line,"%d Column fields and %d records", nAcross, nLine);
2541 handler_->message(CLP_GENERAL,messages_)
2542 << line << CoinMessageEol;
2543 if (nBadName) {
2544 sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2545 handler_->message(CLP_GENERAL,messages_)
2546 << line << CoinMessageEol;
2547 returnCode=-1;
2548 good=false;
2549 }
2550 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2551 free(columnNames[iColumn]);
2552 }
2553 delete [] columnNames;
2554 } else {
2555 sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2556 handler_->message(CLP_GENERAL,messages_)
2557 << line << CoinMessageEol;
2558 returnCode=-1;
2559 good=false;
2560 }
2561 }
2562 }
2563 returnCode=-1;
2564 if (good) {
2565 // clean arrays
2566 if (lowerRowMove) {
2567 bool empty=true;
2568 for (int i=0;i<numberRows_;i++) {
2569 if (lowerRowMove[i]) {
2570 empty=false;
2571 break;
2572 }
2573 }
2574 if (empty) {
2575 delete [] lowerRowMove;
2576 lowerRowMove=NULL;
2577 }
2578 }
2579 if (upperRowMove) {
2580 bool empty=true;
2581 for (int i=0;i<numberRows_;i++) {
2582 if (upperRowMove[i]) {
2583 empty=false;
2584 break;
2585 }
2586 }
2587 if (empty) {
2588 delete [] upperRowMove;
2589 upperRowMove=NULL;
2590 }
2591 }
2592 if (lowerColumnMove) {
2593 bool empty=true;
2594 for (int i=0;i<numberColumns_;i++) {
2595 if (lowerColumnMove[i]) {
2596 empty=false;
2597 break;
2598 }
2599 }
2600 if (empty) {
2601 delete [] lowerColumnMove;
2602 lowerColumnMove=NULL;
2603 }
2604 }
2605 if (upperColumnMove) {
2606 bool empty=true;
2607 for (int i=0;i<numberColumns_;i++) {
2608 if (upperColumnMove[i]) {
2609 empty=false;
2610 break;
2611 }
2612 }
2613 if (empty) {
2614 delete [] upperColumnMove;
2615 upperColumnMove=NULL;
2616 }
2617 }
2618 if (objectiveMove) {
2619 bool empty=true;
2620 for (int i=0;i<numberColumns_;i++) {
2621 if (objectiveMove[i]) {
2622 empty=false;
2623 break;
2624 }
2625 }
2626 if (empty) {
2627 delete [] objectiveMove;
2628 objectiveMove=NULL;
2629 }
2630 }
2631 int saveScaling = scalingFlag_;
2632 scalingFlag_ = 0;
2633 int saveLogLevel = handler_->logLevel();
2634 if (detail>0&&!intervalTheta)
2635 handler_->setLogLevel(3);
2636 else
2637 handler_->setLogLevel(1);
2638 returnCode = parametrics(startTheta,endTheta,intervalTheta,
2639 lowerColumnMove,upperColumnMove,
2640 lowerRowMove,upperRowMove,
2641 objectiveMove);
2642 scalingFlag_ = saveScaling;
2643 handler_->setLogLevel(saveLogLevel);
2644 }
2645 delete [] lowerRowMove;
2646 delete [] upperRowMove;
2647 delete [] lowerColumnMove;
2648 delete [] upperColumnMove;
2649 delete [] objectiveMove;
2650 fclose(fp);
2651 return returnCode;
2652}
2653int
2654ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta, double reportIncrement,
2655 const double * changeLower, const double * changeUpper,
2656 const double * changeObjective, ClpDataSave & data,
2657 bool canTryQuick)
2658{
2659 // stuff is already at starting
2660 // For this crude version just try and go to end
2661 double change = 0.0;
2662 if (reportIncrement && canTryQuick) {
2663 endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
2664 change = endingTheta - startingTheta;
2665 }
2666 int numberTotal = numberRows_ + numberColumns_;
2667 int i;
2668 for ( i = 0; i < numberTotal; i++) {
2669 lower_[i] += change * changeLower[i];
2670 upper_[i] += change * changeUpper[i];
2671 switch(getStatus(i)) {
2672
2673 case basic:
2674 case isFree:
2675 case superBasic:
2676 break;
2677 case isFixed:
2678 case atUpperBound:
2679 solution_[i] = upper_[i];
2680 break;
2681 case atLowerBound:
2682 solution_[i] = lower_[i];
2683 break;
2684 }
2685 cost_[i] += change * changeObjective[i];
2686 }
2687 problemStatus_ = -1;
2688
2689 // This says whether to restore things etc
2690 // startup will have factorized so can skip
2691 int factorType = 0;
2692 // Start check for cycles
2693 progress_.startCheck();
2694 // Say change made on first iteration
2695 changeMade_ = 1;
2696 /*
2697 Status of problem:
2698 0 - optimal
2699 1 - infeasible
2700 2 - unbounded
2701 -1 - iterating
2702 -2 - factorization wanted
2703 -3 - redo checking without factorization
2704 -4 - looks infeasible
2705 */
2706 while (problemStatus_ < 0) {
2707 int iRow, iColumn;
2708 // clear
2709 for (iRow = 0; iRow < 4; iRow++) {
2710 rowArray_[iRow]->clear();
2711 }
2712
2713 for (iColumn = 0; iColumn < 2; iColumn++) {
2714 columnArray_[iColumn]->clear();
2715 }
2716
2717 // give matrix (and model costs and bounds a chance to be
2718 // refreshed (normally null)
2719 matrix_->refresh(this);
2720 // may factorize, checks if problem finished
2721 statusOfProblemInParametrics(factorType, data);
2722 // Say good factorization
2723 factorType = 1;
2724 if (data.sparseThreshold_) {
2725 // use default at present
2726 factorization_->sparseThreshold(0);
2727 factorization_->goSparse();
2728 }
2729
2730 // exit if victory declared
2731 if (problemStatus_ >= 0 &&
2732 (canTryQuick || startingTheta>=endingTheta-1.0e-7) )
2733 break;
2734
2735 // test for maximum iterations
2736 if (hitMaximumIterations()) {
2737 problemStatus_ = 3;
2738 break;
2739 }
2740 // Check event
2741 {
2742 int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2743 if (status >= 0) {
2744 problemStatus_ = 5;
2745 secondaryStatus_ = ClpEventHandler::endOfFactorization;
2746 break;
2747 }
2748 }
2749 // Do iterations
2750 problemStatus_=-1;
2751 if (canTryQuick) {
2752 double * saveDuals = NULL;
2753 reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
2754 } else {
2755 whileIterating(startingTheta, endingTheta, reportIncrement,
2756 changeLower, changeUpper,
2757 changeObjective);
2758 startingTheta = endingTheta;
2759 }
2760 }
2761 if (!problemStatus_) {
2762 theta_ = change + startingTheta;
2763 eventHandler_->event(ClpEventHandler::theta);
2764 return 0;
2765 } else if (problemStatus_ == 10) {
2766 return -1;
2767 } else {
2768 return problemStatus_;
2769 }
2770}
2771/* Checks if finished. Updates status */
2772void
2773ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
2774{
2775 if (type == 2) {
2776 // trouble - go to recovery
2777 problemStatus_ = 10;
2778 return;
2779 }
2780 if (problemStatus_ > -3 || factorization_->pivots()) {
2781 // factorize
2782 // later on we will need to recover from singularities
2783 // also we could skip if first time
2784 if (type) {
2785 // is factorization okay?
2786 if (internalFactorize(1)) {
2787 // trouble - go to recovery
2788 problemStatus_ = 10;
2789 return;
2790 }
2791 }
2792 if (problemStatus_ != -4 || factorization_->pivots() > 10)
2793 problemStatus_ = -3;
2794 }
2795 // at this stage status is -3 or -4 if looks infeasible
2796 // get primal and dual solutions
2797 gutsOfSolution(NULL, NULL);
2798 double realDualInfeasibilities = sumDualInfeasibilities_;
2799 // If bad accuracy treat as singular
2800 if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
2801 // trouble - go to recovery
2802 problemStatus_ = 10;
2803 return;
2804 } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
2805 // Can reduce tolerance
2806 double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
2807 factorization_->pivotTolerance(newTolerance);
2808 }
2809 // Check if looping
2810 int loop;
2811 if (type != 2)
2812 loop = progress_.looping();
2813 else
2814 loop = -1;
2815 if (loop >= 0) {
2816 problemStatus_ = loop; //exit if in loop
2817 if (!problemStatus_) {
2818 // declaring victory
2819 numberPrimalInfeasibilities_ = 0;
2820 sumPrimalInfeasibilities_ = 0.0;
2821 } else {
2822 problemStatus_ = 10; // instead - try other algorithm
2823 }
2824 return;
2825 } else if (loop < -1) {
2826 // something may have changed
2827 gutsOfSolution(NULL, NULL);
2828 }
2829 progressFlag_ = 0; //reset progress flag
2830 if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
2831 handler_->message(CLP_SIMPLEX_STATUS, messages_)
2832 << numberIterations_ << objectiveValue();
2833 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
2834 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
2835 handler_->printing(sumDualInfeasibilities_ > 0.0)
2836 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
2837 handler_->printing(numberDualInfeasibilitiesWithoutFree_
2838 < numberDualInfeasibilities_)
2839 << numberDualInfeasibilitiesWithoutFree_;
2840 handler_->message() << CoinMessageEol;
2841 }
2842 /* If we are primal feasible and any dual infeasibilities are on
2843 free variables then it is better to go to primal */
2844 if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
2845 numberDualInfeasibilities_) {
2846 problemStatus_ = 10;
2847 return;
2848 }
2849
2850 // check optimal
2851 // give code benefit of doubt
2852 if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
2853 sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
2854 // say optimal (with these bounds etc)
2855 numberDualInfeasibilities_ = 0;
2856 sumDualInfeasibilities_ = 0.0;
2857 numberPrimalInfeasibilities_ = 0;
2858 sumPrimalInfeasibilities_ = 0.0;
2859 }
2860 if (dualFeasible() || problemStatus_ == -4) {
2861 progress_.modifyObjective(objectiveValue_
2862 - sumDualInfeasibilities_ * dualBound_);
2863 }
2864 if (numberPrimalInfeasibilities_) {
2865 if (problemStatus_ == -4 || problemStatus_ == -5) {
2866 problemStatus_ = 1; // infeasible
2867 }
2868 } else if (numberDualInfeasibilities_) {
2869 // clean up
2870 problemStatus_ = 10;
2871 } else {
2872 problemStatus_ = 0;
2873 }
2874 lastGoodIteration_ = numberIterations_;
2875 if (problemStatus_ < 0) {
2876 sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
2877 if (sumDualInfeasibilities_)
2878 numberDualInfeasibilities_ = 1;
2879 }
2880 // Allow matrices to be sorted etc
2881 int fake = -999; // signal sort
2882 matrix_->correctSequence(this, fake, fake);
2883}
2884/* This has the flow between re-factorizations
2885 Reasons to come out:
2886 -1 iterations etc
2887 -2 inaccuracy
2888 -3 slight inaccuracy (and done iterations)
2889 +0 looks optimal (might be unbounded - but we will investigate)
2890 +1 looks infeasible
2891 +3 max iterations
2892 +4 accuracy problems
2893*/
2894int
2895ClpSimplexOther::whileIterating(double startingTheta, double & endingTheta, double /*reportIncrement*/,
2896 const double * changeLower, const double * changeUpper,
2897 const double * changeObjective)
2898{
2899 {
2900 int i;
2901 for (i = 0; i < 4; i++) {
2902 rowArray_[i]->clear();
2903 }
2904 for (i = 0; i < 2; i++) {
2905 columnArray_[i]->clear();
2906 }
2907 }
2908 // if can't trust much and long way from optimal then relax
2909 if (largestPrimalError_ > 10.0)
2910 factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
2911 else
2912 factorization_->relaxAccuracyCheck(1.0);
2913 // status stays at -1 while iterating, >=0 finished, -2 to invert
2914 // status -3 to go to top without an invert
2915 int returnCode = -1;
2916 double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
2917 double lastTheta = startingTheta;
2918 double useTheta = startingTheta;
2919 int numberTotal = numberColumns_ + numberRows_;
2920 double * primalChange = new double[numberTotal];
2921 double * dualChange = new double[numberTotal];
2922 int iSequence;
2923 // See if bounds
2924 int type = 0;
2925 for (iSequence = 0; iSequence < numberTotal; iSequence++) {
2926 if (changeLower[iSequence] || changeUpper[iSequence]) {
2927 type = 1;
2928 break;
2929 }
2930 }
2931 // See if objective
2932 for (iSequence = 0; iSequence < numberTotal; iSequence++) {
2933 if (changeObjective[iSequence]) {
2934 type |= 2;
2935 break;
2936 }
2937 }
2938 assert (type);
2939 while (problemStatus_ == -1) {
2940 double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
2941
2942 // Get theta for bounds - we know can't crossover
2943 int pivotType = nextTheta(type, increaseTheta, primalChange, dualChange,
2944 changeLower, changeUpper, changeObjective);
2945 useTheta += theta_;
2946 double change = useTheta - lastTheta;
2947 for (int i = 0; i < numberTotal; i++) {
2948 lower_[i] += change * changeLower[i];
2949 upper_[i] += change * changeUpper[i];
2950 switch(getStatus(i)) {
2951
2952 case basic:
2953 case isFree:
2954 case superBasic:
2955 break;
2956 case isFixed:
2957 case atUpperBound:
2958 solution_[i] = upper_[i];
2959 break;
2960 case atLowerBound:
2961 solution_[i] = lower_[i];
2962 break;
2963 }
2964 cost_[i] += change * changeObjective[i];
2965 assert (solution_[i]>lower_[i]-1.0e-5&&
2966 solution_[i]<upper_[i]+1.0e-5);
2967 }
2968 sequenceIn_=-1;
2969 if (pivotType) {
2970 problemStatus_ = -2;
2971 endingTheta = useTheta;
2972 return 4;
2973 }
2974 // choose row to go out
2975 //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
2976 if (pivotRow_ >= 0) {
2977 // we found a pivot row
2978 if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
2979 handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
2980 << pivotRow_
2981 << CoinMessageEol;
2982 }
2983 // check accuracy of weights
2984 dualRowPivot_->checkAccuracy();
2985 // Get good size for pivot
2986 // Allow first few iterations to take tiny
2987 double acceptablePivot = 1.0e-9;
2988 if (numberIterations_ > 100)
2989 acceptablePivot = 1.0e-8;
2990 if (factorization_->pivots() > 10 ||
2991 (factorization_->pivots() && saveSumDual))
2992 acceptablePivot = 1.0e-5; // if we have iterated be more strict
2993 else if (factorization_->pivots() > 5)
2994 acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
2995 else if (factorization_->pivots())
2996 acceptablePivot = 1.0e-8; // relax
2997 double bestPossiblePivot = 1.0;
2998 // get sign for finding row of tableau
2999 // normal iteration
3000 // create as packed
3001 double direction = directionOut_;
3002 rowArray_[0]->createPacked(1, &pivotRow_, &direction);
3003 factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
3004 // put row of tableau in rowArray[0] and columnArray[0]
3005 matrix_->transposeTimes(this, -1.0,
3006 rowArray_[0], rowArray_[3], columnArray_[0]);
3007 // do ratio test for normal iteration
3008 bestPossiblePivot = reinterpret_cast<ClpSimplexDual *> ( this)->dualColumn(rowArray_[0],
3009 columnArray_[0], columnArray_[1],
3010 rowArray_[3], acceptablePivot, NULL);
3011 if (sequenceIn_ >= 0) {
3012 // normal iteration
3013 // update the incoming column
3014 double btranAlpha = -alpha_ * directionOut_; // for check
3015 unpackPacked(rowArray_[1]);
3016 // moved into updateWeights factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
3017 // and update dual weights (can do in parallel - with extra array)
3018 alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
3019 rowArray_[2],
3020 rowArray_[3],
3021 rowArray_[1]);
3022 // see if update stable
3023#ifdef CLP_DEBUG
3024 if ((handler_->logLevel() & 32))
3025 printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
3026#endif
3027 double checkValue = 1.0e-7;
3028 // if can't trust much and long way from optimal then relax
3029 if (largestPrimalError_ > 10.0)
3030 checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
3031 if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3032 fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
3033 handler_->message(CLP_DUAL_CHECK, messages_)
3034 << btranAlpha
3035 << alpha_
3036 << CoinMessageEol;
3037 if (factorization_->pivots()) {
3038 dualRowPivot_->unrollWeights();
3039 problemStatus_ = -2; // factorize now
3040 rowArray_[0]->clear();
3041 rowArray_[1]->clear();
3042 columnArray_[0]->clear();
3043 returnCode = -2;
3044 break;
3045 } else {
3046 // take on more relaxed criterion
3047 double test;
3048 if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
3049 test = 1.0e-1 * fabs(alpha_);
3050 else
3051 test = 1.0e-4 * (1.0 + fabs(alpha_));
3052 if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3053 fabs(btranAlpha - alpha_) > test) {
3054 dualRowPivot_->unrollWeights();
3055 // need to reject something
3056 char x = isColumn(sequenceOut_) ? 'C' : 'R';
3057 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3058 << x << sequenceWithin(sequenceOut_)
3059 << CoinMessageEol;
3060 setFlagged(sequenceOut_);
3061 progress_.clearBadTimes();
3062 lastBadIteration_ = numberIterations_; // say be more cautious
3063 rowArray_[0]->clear();
3064 rowArray_[1]->clear();
3065 columnArray_[0]->clear();
3066 if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
3067 //printf("I think should declare infeasible\n");
3068 problemStatus_ = 1;
3069 returnCode = 1;
3070 break;
3071 }
3072 continue;
3073 }
3074 }
3075 }
3076 // update duals BEFORE replaceColumn so can do updateColumn
3077 double objectiveChange = 0.0;
3078 // do duals first as variables may flip bounds
3079 // rowArray_[0] and columnArray_[0] may have flips
3080 // so use rowArray_[3] for work array from here on
3081 int nswapped = 0;
3082 //rowArray_[0]->cleanAndPackSafe(1.0e-60);
3083 //columnArray_[0]->cleanAndPackSafe(1.0e-60);
3084 nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
3085 rowArray_[2], theta_,
3086 objectiveChange, false);
3087
3088 // which will change basic solution
3089 if (nswapped) {
3090 factorization_->updateColumn(rowArray_[3], rowArray_[2]);
3091 dualRowPivot_->updatePrimalSolution(rowArray_[2],
3092 1.0, objectiveChange);
3093 // recompute dualOut_
3094 valueOut_ = solution_[sequenceOut_];
3095 if (directionOut_ < 0) {
3096 dualOut_ = valueOut_ - upperOut_;
3097 } else {
3098 dualOut_ = lowerOut_ - valueOut_;
3099 }
3100 }
3101 // amount primal will move
3102 double movement = -dualOut_ * directionOut_ / alpha_;
3103 // so objective should increase by fabs(dj)*movement
3104 // but we already have objective change - so check will be good
3105 if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
3106#ifdef CLP_DEBUG
3107 if (handler_->logLevel() & 32)
3108 printf("movement %g, swap change %g, rest %g * %g\n",
3109 objectiveChange + fabs(movement * dualIn_),
3110 objectiveChange, movement, dualIn_);
3111#endif
3112 if(factorization_->pivots()) {
3113 // going backwards - factorize
3114 dualRowPivot_->unrollWeights();
3115 problemStatus_ = -2; // factorize now
3116 returnCode = -2;
3117 break;
3118 }
3119 }
3120 CoinAssert(fabs(dualOut_) < 1.0e50);
3121 // if stable replace in basis
3122 int updateStatus = factorization_->replaceColumn(this,
3123 rowArray_[2],
3124 rowArray_[1],
3125 pivotRow_,
3126 alpha_);
3127 // if no pivots, bad update but reasonable alpha - take and invert
3128 if (updateStatus == 2 &&
3129 !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
3130 updateStatus = 4;
3131 if (updateStatus == 1 || updateStatus == 4) {
3132 // slight error
3133 if (factorization_->pivots() > 5 || updateStatus == 4) {
3134 problemStatus_ = -2; // factorize now
3135 returnCode = -3;
3136 }
3137 } else if (updateStatus == 2) {
3138 // major error
3139 dualRowPivot_->unrollWeights();
3140 // later we may need to unwind more e.g. fake bounds
3141 if (factorization_->pivots()) {
3142 problemStatus_ = -2; // factorize now
3143 returnCode = -2;
3144 break;
3145 } else {
3146 // need to reject something
3147 char x = isColumn(sequenceOut_) ? 'C' : 'R';
3148 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3149 << x << sequenceWithin(sequenceOut_)
3150 << CoinMessageEol;
3151 setFlagged(sequenceOut_);
3152 progress_.clearBadTimes();
3153 lastBadIteration_ = numberIterations_; // say be more cautious
3154 rowArray_[0]->clear();
3155 rowArray_[1]->clear();
3156 columnArray_[0]->clear();
3157 // make sure dual feasible
3158 // look at all rows and columns
3159 double objectiveChange = 0.0;
3160 reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
3161 0.0, objectiveChange, true);
3162 continue;
3163 }
3164 } else if (updateStatus == 3) {
3165 // out of memory
3166 // increase space if not many iterations
3167 if (factorization_->pivots() <
3168 0.5 * factorization_->maximumPivots() &&
3169 factorization_->pivots() < 200)
3170 factorization_->areaFactor(
3171 factorization_->areaFactor() * 1.1);
3172 problemStatus_ = -2; // factorize now
3173 } else if (updateStatus == 5) {
3174 problemStatus_ = -2; // factorize now
3175 }
3176 // update primal solution
3177 if (theta_ < 0.0) {
3178#ifdef CLP_DEBUG
3179 if (handler_->logLevel() & 32)
3180 printf("negative theta %g\n", theta_);
3181#endif
3182 theta_ = 0.0;
3183 }
3184 // do actual flips
3185 reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
3186 //rowArray_[1]->expand();
3187 dualRowPivot_->updatePrimalSolution(rowArray_[1],
3188 movement,
3189 objectiveChange);
3190 // modify dualout
3191 dualOut_ /= alpha_;
3192 dualOut_ *= -directionOut_;
3193 //setStatus(sequenceIn_,basic);
3194 dj_[sequenceIn_] = 0.0;
3195 //double oldValue = valueIn_;
3196 if (directionIn_ == -1) {
3197 // as if from upper bound
3198 valueIn_ = upperIn_ + dualOut_;
3199 } else {
3200 // as if from lower bound
3201 valueIn_ = lowerIn_ + dualOut_;
3202 }
3203 objectiveChange = 0.0;
3204 for (int i=0;i<numberTotal;i++)
3205 objectiveChange += solution_[i]*cost_[i];
3206 objectiveChange -= objectiveValue_;
3207 // outgoing
3208 // set dj to zero unless values pass
3209 if (directionOut_ > 0) {
3210 valueOut_ = lowerOut_;
3211 dj_[sequenceOut_] = theta_;
3212 } else {
3213 valueOut_ = upperOut_;
3214 dj_[sequenceOut_] = -theta_;
3215 }
3216 solution_[sequenceOut_] = valueOut_;
3217 int whatNext = housekeeping(objectiveChange);
3218 {
3219 char in[200],out[200];
3220 int iSequence=sequenceIn_;
3221 if (iSequence<numberColumns_) {
3222 if (lengthNames_)
3223 strcpy(in,columnNames_[iSequence].c_str());
3224 else
3225 sprintf(in,"C%7.7d",iSequence);
3226 } else {
3227 iSequence -= numberColumns_;
3228 if (lengthNames_)
3229 strcpy(in,rowNames_[iSequence].c_str());
3230 else
3231 sprintf(in,"R%7.7d",iSequence);
3232 }
3233 iSequence=sequenceOut_;
3234 if (iSequence<numberColumns_) {
3235 if (lengthNames_)
3236 strcpy(out,columnNames_[iSequence].c_str());
3237 else
3238 sprintf(out,"C%7.7d",iSequence);
3239 } else {
3240 iSequence -= numberColumns_;
3241 if (lengthNames_)
3242 strcpy(out,rowNames_[iSequence].c_str());
3243 else
3244 sprintf(out,"R%7.7d",iSequence);
3245 }
3246 handler_->message(CLP_PARAMETRICS_STATS2, messages_)
3247 << useTheta << objectiveValue()
3248 << in << out << CoinMessageEol;
3249 }
3250 if (useTheta>lastTheta+1.0e-9) {
3251 handler_->message(CLP_PARAMETRICS_STATS, messages_)
3252 << useTheta << objectiveValue() << CoinMessageEol;
3253 lastTheta = useTheta;
3254 }
3255 // and set bounds correctly
3256 reinterpret_cast<ClpSimplexDual *> ( this)->originalBound(sequenceIn_);
3257 reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
3258 if (whatNext == 1) {
3259 problemStatus_ = -2; // refactorize
3260 } else if (whatNext == 2) {
3261 // maximum iterations or equivalent
3262 problemStatus_ = 3;
3263 returnCode = 3;
3264 break;
3265 }
3266 // Check event
3267 {
3268 int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3269 if (status >= 0) {
3270 problemStatus_ = 5;
3271 secondaryStatus_ = ClpEventHandler::endOfIteration;
3272 returnCode = 4;
3273 break;
3274 }
3275 }
3276 } else {
3277 // no incoming column is valid
3278 pivotRow_ = -1;
3279#ifdef CLP_DEBUG
3280 if (handler_->logLevel() & 32)
3281 printf("** no column pivot\n");
3282#endif
3283 if (factorization_->pivots() < 5) {
3284 // If not in branch and bound etc save ray
3285 if ((specialOptions_&(1024 | 4096)) == 0) {
3286 // create ray anyway
3287 delete [] ray_;
3288 ray_ = new double [ numberRows_];
3289 rowArray_[0]->expand(); // in case packed
3290 ClpDisjointCopyN(rowArray_[0]->denseVector(), numberRows_, ray_);
3291 }
3292 // If we have just factorized and infeasibility reasonable say infeas
3293 if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
3294 if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3
3295 || (specialOptions_ & 64) == 0) {
3296 // say infeasible
3297 problemStatus_ = 1;
3298 // unless primal feasible!!!!
3299 //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
3300 // numberDualInfeasibilities_,sumDualInfeasibilities_);
3301 if (numberDualInfeasibilities_)
3302 problemStatus_ = 10;
3303 rowArray_[0]->clear();
3304 columnArray_[0]->clear();
3305 returnCode = 1;
3306 break;
3307 }
3308 }
3309 // If special option set - put off as long as possible
3310 if ((specialOptions_ & 64) == 0) {
3311 problemStatus_ = -4; //say looks infeasible
3312 } else {
3313 // flag
3314 char x = isColumn(sequenceOut_) ? 'C' : 'R';
3315 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3316 << x << sequenceWithin(sequenceOut_)
3317 << CoinMessageEol;
3318 setFlagged(sequenceOut_);
3319 if (!factorization_->pivots()) {
3320 rowArray_[0]->clear();
3321 columnArray_[0]->clear();
3322 continue;
3323 }
3324 }
3325 }
3326 rowArray_[0]->clear();
3327 columnArray_[0]->clear();
3328 returnCode = 1;
3329 break;
3330 }
3331 } else {
3332 // no pivot row
3333#ifdef CLP_DEBUG
3334 if (handler_->logLevel() & 32)
3335 printf("** no row pivot\n");
3336#endif
3337 int numberPivots = factorization_->pivots();
3338 bool specialCase;
3339 int useNumberFake;
3340 returnCode = 0;
3341 if (numberPivots < 20 &&
3342 (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
3343 && dualBound_ > 1.0e8) {
3344 specialCase = true;
3345 // as dual bound high - should be okay
3346 useNumberFake = 0;
3347 } else {
3348 specialCase = false;
3349 useNumberFake = numberFake_;
3350 }
3351 if (!numberPivots || specialCase) {
3352 // may have crept through - so may be optimal
3353 // check any flagged variables
3354 int iRow;
3355 for (iRow = 0; iRow < numberRows_; iRow++) {
3356 int iPivot = pivotVariable_[iRow];
3357 if (flagged(iPivot))
3358 break;
3359 }
3360 if (iRow < numberRows_ && numberPivots) {
3361 // try factorization
3362 returnCode = -2;
3363 }
3364
3365 if (useNumberFake || numberDualInfeasibilities_) {
3366 // may be dual infeasible
3367 problemStatus_ = -5;
3368 } else {
3369 if (iRow < numberRows_) {
3370 problemStatus_ = -5;
3371 } else {
3372 if (numberPivots) {
3373 // objective may be wrong
3374 objectiveValue_ = innerProduct(cost_,
3375 numberColumns_ + numberRows_,
3376 solution_);
3377 objectiveValue_ += objective_->nonlinearOffset();
3378 objectiveValue_ /= (objectiveScale_ * rhsScale_);
3379 if ((specialOptions_ & 16384) == 0) {
3380 // and dual_ may be wrong (i.e. for fixed or basic)
3381 CoinIndexedVector * arrayVector = rowArray_[1];
3382 arrayVector->clear();
3383 int iRow;
3384 double * array = arrayVector->denseVector();
3385 /* Use dual_ instead of array
3386 Even though dual_ is only numberRows_ long this is
3387 okay as gets permuted to longer rowArray_[2]
3388 */
3389 arrayVector->setDenseVector(dual_);
3390 int * index = arrayVector->getIndices();
3391 int number = 0;
3392 for (iRow = 0; iRow < numberRows_; iRow++) {
3393 int iPivot = pivotVariable_[iRow];
3394 double value = cost_[iPivot];
3395 dual_[iRow] = value;
3396 if (value) {
3397 index[number++] = iRow;
3398 }
3399 }
3400 arrayVector->setNumElements(number);
3401 // Extended duals before "updateTranspose"
3402 matrix_->dualExpanded(this, arrayVector, NULL, 0);
3403 // Btran basic costs
3404 rowArray_[2]->clear();
3405 factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
3406 // and return vector
3407 arrayVector->setDenseVector(array);
3408 }
3409 }
3410 problemStatus_ = 0;
3411 sumPrimalInfeasibilities_ = 0.0;
3412 if ((specialOptions_&(1024 + 16384)) != 0) {
3413 CoinIndexedVector * arrayVector = rowArray_[1];
3414 arrayVector->clear();
3415 double * rhs = arrayVector->denseVector();
3416 times(1.0, solution_, rhs);
3417 bool bad2 = false;
3418 int i;
3419 for ( i = 0; i < numberRows_; i++) {
3420 if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
3421 rhs[i] > rowUpperWork_[i] + primalTolerance_) {
3422 bad2 = true;
3423 } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
3424 }
3425 rhs[i] = 0.0;
3426 }
3427 for ( i = 0; i < numberColumns_; i++) {
3428 if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
3429 solution_[i] > columnUpperWork_[i] + primalTolerance_) {
3430 bad2 = true;
3431 }
3432 }
3433 if (bad2) {
3434 problemStatus_ = -3;
3435 returnCode = -2;
3436 // Force to re-factorize early next time
3437 int numberPivots = factorization_->pivots();
3438 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
3439 }
3440 }
3441 }
3442 }
3443 } else {
3444 problemStatus_ = -3;
3445 returnCode = -2;
3446 // Force to re-factorize early next time
3447 int numberPivots = factorization_->pivots();
3448 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
3449 }
3450 break;
3451 }
3452 }
3453 delete [] primalChange;
3454 delete [] dualChange;
3455 endingTheta = lastTheta;
3456 return returnCode;
3457}
3458// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
3459int
3460ClpSimplexOther::nextTheta(int type, double maxTheta, double * primalChange, double * /*dualChange*/,
3461 const double * changeLower, const double * changeUpper,
3462 const double * /*changeObjective*/)
3463{
3464 int numberTotal = numberColumns_ + numberRows_;
3465 int iSequence;
3466 int iRow;
3467 theta_ = maxTheta;
3468 bool toLower = false;
3469 if ((type & 1) != 0) {
3470 // get change
3471 for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3472 primalChange[iSequence] = 0.0;
3473 switch(getStatus(iSequence)) {
3474
3475 case basic:
3476 case isFree:
3477 case superBasic:
3478 break;
3479 case isFixed:
3480 case atUpperBound:
3481 primalChange[iSequence] = changeUpper[iSequence];
3482 break;
3483 case atLowerBound:
3484 primalChange[iSequence] = changeLower[iSequence];
3485 break;
3486 }
3487 }
3488 // use array
3489 double * array = rowArray_[1]->denseVector();
3490 // put slacks in
3491 for (int i=0;i<numberRows_;i++)
3492 array[i] = - primalChange[i+numberColumns_];
3493 times(1.0, primalChange, array);
3494 int * index = rowArray_[1]->getIndices();
3495 int number = 0;
3496 pivotRow_ = -1;
3497 for (iRow = 0; iRow < numberRows_; iRow++) {
3498 double value = array[iRow];
3499 if (value) {
3500 index[number++] = iRow;
3501 }
3502 }
3503 // ftran it
3504 rowArray_[1]->setNumElements(number);
3505 factorization_->updateColumn(rowArray_[0], rowArray_[1]);
3506 //number = rowArray_[1]->getNumElements();
3507 for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
3508 //int iPivot = index[iRow];
3509 iSequence = pivotVariable_[iPivot];
3510 // solution value will be sol - theta*alpha
3511 // bounds will be bounds + change *theta
3512 double currentSolution = solution_[iSequence];
3513 double currentLower = lower_[iSequence];
3514 double currentUpper = upper_[iSequence];
3515 double alpha = array[iPivot];
3516 assert (currentSolution >= currentLower - primalTolerance_);
3517 assert (currentSolution <= currentUpper + primalTolerance_);
3518 double thetaCoefficient;
3519 double hitsLower = COIN_DBL_MAX;
3520 thetaCoefficient = changeLower[iSequence] + alpha;
3521 if (thetaCoefficient > 1.0e-8)
3522 hitsLower = (currentSolution - currentLower) / thetaCoefficient;
3523 //if (hitsLower < 0.0) {
3524 // does not hit - but should we check further
3525 // hitsLower = COIN_DBL_MAX;
3526 //}
3527 double hitsUpper = COIN_DBL_MAX;
3528 thetaCoefficient = changeUpper[iSequence] + alpha;
3529 if (thetaCoefficient < -1.0e-8)
3530 hitsUpper = (currentSolution - currentUpper) / thetaCoefficient;
3531 //if (hitsUpper < 0.0) {
3532 // does not hit - but should we check further
3533 // hitsUpper = COIN_DBL_MAX;
3534 //}
3535 if (CoinMin(hitsLower, hitsUpper) < theta_) {
3536 theta_ = CoinMin(hitsLower, hitsUpper);
3537 toLower = hitsLower < hitsUpper;
3538 pivotRow_ = iPivot;
3539 }
3540 }
3541 }
3542 if ((type & 2) != 0) {
3543 abort();
3544 }
3545 theta_ = CoinMax(theta_,0.0);
3546 // update solution
3547 double * array = rowArray_[1]->denseVector();
3548 int * index = rowArray_[1]->getIndices();
3549 int number = rowArray_[1]->getNumElements();
3550 for (int iRow = 0; iRow < number; iRow++) {
3551 int iPivot = index[iRow];
3552 iSequence = pivotVariable_[iPivot];
3553 // solution value will be sol - theta*alpha
3554 double alpha = array[iPivot];
3555 solution_[iSequence] -= theta_ * alpha;
3556 }
3557 if (pivotRow_ >= 0) {
3558 sequenceOut_ = pivotVariable_[pivotRow_];
3559 valueOut_ = solution_[sequenceOut_];
3560 lowerOut_ = lower_[sequenceOut_]+theta_*changeLower[sequenceOut_];
3561 upperOut_ = upper_[sequenceOut_]+theta_*changeUpper[sequenceOut_];
3562 if (!toLower) {
3563 directionOut_ = -1;
3564 dualOut_ = valueOut_ - upperOut_;
3565 } else {
3566 directionOut_ = 1;
3567 dualOut_ = lowerOut_ - valueOut_;
3568 }
3569 return 0;
3570 } else {
3571 return -1;
3572 }
3573}
3574/* Expands out all possible combinations for a knapsack
3575 If buildObj NULL then just computes space needed - returns number elements
3576 On entry numberOutput is maximum allowed, on exit it is number needed or
3577 -1 (as will be number elements) if maximum exceeded. numberOutput will have at
3578 least space to return values which reconstruct input.
3579 Rows returned will be original rows but no entries will be returned for
3580 any rows all of whose entries are in knapsack. So up to user to allow for this.
3581 If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
3582 in expanded knapsack. Values in buildRow and buildElement;
3583*/
3584int
3585ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
3586 double * buildObj, CoinBigIndex * buildStart,
3587 int * buildRow, double * buildElement, int reConstruct) const
3588{
3589 int iRow;
3590 int iColumn;
3591 // Get column copy
3592 CoinPackedMatrix * columnCopy = matrix();
3593 // Get a row copy in standard format
3594 CoinPackedMatrix matrixByRow;
3595 matrixByRow.reverseOrderedCopyOf(*columnCopy);
3596 const double * elementByRow = matrixByRow.getElements();
3597 const int * column = matrixByRow.getIndices();
3598 const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3599 const int * rowLength = matrixByRow.getVectorLengths();
3600 CoinBigIndex j;
3601 int * whichColumn = new int [numberColumns_];
3602 int * whichRow = new int [numberRows_];
3603 int numJ = 0;
3604 // Get what other columns can compensate for
3605 double * lo = new double [numberRows_];
3606 double * high = new double [numberRows_];
3607 {
3608 // Use to get tight column bounds
3609 ClpSimplex tempModel(*this);
3610 tempModel.tightenPrimalBounds(0.0, 0, true);
3611 // Now another model without knapsacks
3612 int nCol = 0;
3613 for (iRow = 0; iRow < numberRows_; iRow++) {
3614 whichRow[iRow] = iRow;
3615 }
3616 for (iColumn = 0; iColumn < numberColumns_; iColumn++)
3617 whichColumn[iColumn] = -1;
3618 for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
3619 int iColumn = column[j];
3620 if (columnUpper_[iColumn] > columnLower_[iColumn]) {
3621 whichColumn[iColumn] = 0;
3622 } else {
3623 assert (!columnLower_[iColumn]); // fix later
3624 }
3625 }
3626 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3627 if (whichColumn[iColumn] < 0)
3628 whichColumn[nCol++] = iColumn;
3629 }
3630 ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
3631 // Row copy
3632 CoinPackedMatrix matrixByRow;
3633 matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
3634 const double * elementByRow = matrixByRow.getElements();
3635 const int * column = matrixByRow.getIndices();
3636 const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3637 const int * rowLength = matrixByRow.getVectorLengths();
3638 const double * columnLower = tempModel2.getColLower();
3639 const double * columnUpper = tempModel2.getColUpper();
3640 for (iRow = 0; iRow < numberRows_; iRow++) {
3641 lo[iRow] = -COIN_DBL_MAX;
3642 high[iRow] = COIN_DBL_MAX;
3643 if (rowLower_[iRow] > -1.0e20 || rowUpper_[iRow] < 1.0e20) {
3644
3645 // possible row
3646 int infiniteUpper = 0;
3647 int infiniteLower = 0;
3648 double maximumUp = 0.0;
3649 double maximumDown = 0.0;
3650 CoinBigIndex rStart = rowStart[iRow];
3651 CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3652 CoinBigIndex j;
3653 // Compute possible lower and upper ranges
3654
3655 for (j = rStart; j < rEnd; ++j) {
3656 double value = elementByRow[j];
3657 iColumn = column[j];
3658 if (value > 0.0) {
3659 if (columnUpper[iColumn] >= 1.0e20) {
3660 ++infiniteUpper;
3661 } else {
3662 maximumUp += columnUpper[iColumn] * value;
3663 }
3664 if (columnLower[iColumn] <= -1.0e20) {
3665 ++infiniteLower;
3666 } else {
3667 maximumDown += columnLower[iColumn] * value;
3668 }
3669 } else if (value < 0.0) {
3670 if (columnUpper[iColumn] >= 1.0e20) {
3671 ++infiniteLower;
3672 } else {
3673 maximumDown += columnUpper[iColumn] * value;
3674 }
3675 if (columnLower[iColumn] <= -1.0e20) {
3676 ++infiniteUpper;
3677 } else {
3678 maximumUp += columnLower[iColumn] * value;
3679 }
3680 }
3681 }
3682 // Build in a margin of error
3683 maximumUp += 1.0e-8 * fabs(maximumUp) + 1.0e-7;
3684 maximumDown -= 1.0e-8 * fabs(maximumDown) + 1.0e-7;
3685 // we want to save effective rhs
3686 double up = (infiniteUpper) ? COIN_DBL_MAX : maximumUp;
3687 double down = (infiniteLower) ? -COIN_DBL_MAX : maximumDown;
3688 if (up == COIN_DBL_MAX || rowLower_[iRow] == -COIN_DBL_MAX) {
3689 // However low we go it doesn't matter
3690 lo[iRow] = -COIN_DBL_MAX;
3691 } else {
3692 // If we go below this then can not be feasible
3693 lo[iRow] = rowLower_[iRow] - up;
3694 }
3695 if (down == -COIN_DBL_MAX || rowUpper_[iRow] == COIN_DBL_MAX) {
3696 // However high we go it doesn't matter
3697 high[iRow] = COIN_DBL_MAX;
3698 } else {
3699 // If we go above this then can not be feasible
3700 high[iRow] = rowUpper_[iRow] - down;
3701 }
3702 }
3703 }
3704 }
3705 numJ = 0;
3706 for (iColumn = 0; iColumn < numberColumns_; iColumn++)
3707 whichColumn[iColumn] = -1;
3708 int * markRow = new int [numberRows_];
3709 for (iRow = 0; iRow < numberRows_; iRow++)
3710 markRow[iRow] = 1;
3711 for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
3712 int iColumn = column[j];
3713 if (columnUpper_[iColumn] > columnLower_[iColumn]) {
3714 whichColumn[iColumn] = numJ;
3715 numJ++;
3716 }
3717 }
3718 /* mark rows
3719 -n in knapsack and n other variables
3720 1 no entries
3721 n+1000 not involved in knapsack but n entries
3722 0 only in knapsack
3723 */
3724 for (iRow = 0; iRow < numberRows_; iRow++) {
3725 int type = 1;
3726 for (j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
3727 int iColumn = column[j];
3728 if (whichColumn[iColumn] >= 0) {
3729 if (type == 1) {
3730 type = 0;
3731 } else if (type > 0) {
3732 assert (type > 1000);
3733 type = -(type - 1000);
3734 }
3735 } else if (type == 1) {
3736 type = 1001;
3737 } else if (type < 0) {
3738 type --;
3739 } else if (type == 0) {
3740 type = -1;
3741 } else {
3742 assert (type > 1000);
3743 type++;
3744 }
3745 }
3746 markRow[iRow] = type;
3747 if (type < 0 && type > -30 && false)
3748 printf("markrow on row %d is %d\n", iRow, markRow[iRow]);
3749 }
3750 int * bound = new int [numberColumns_+1];
3751 int * stack = new int [numberColumns_+1];
3752 int * flip = new int [numberColumns_+1];
3753 double * offset = new double[numberColumns_+1];
3754 double * size = new double [numberColumns_+1];
3755 double * rhsOffset = new double[numberRows_];
3756 int * build = new int[numberColumns_];
3757 int maxNumber = numberOutput;
3758 numJ = 0;
3759 double minSize = rowLower_[knapsackRow];
3760 double maxSize = rowUpper_[knapsackRow];
3761 double knapsackOffset = 0.0;
3762 for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
3763 int iColumn = column[j];
3764 double lowerColumn = columnLower_[iColumn];
3765 double upperColumn = columnUpper_[iColumn];
3766 if (lowerColumn == upperColumn)
3767 continue;
3768 double gap = upperColumn - lowerColumn;
3769 if (gap > 1.0e8)
3770 gap = 1.0e8;
3771 assert (fabs(floor(gap + 0.5) - gap) < 1.0e-5);
3772 whichColumn[numJ] = iColumn;
3773 bound[numJ] = static_cast<int> (gap);
3774 if (elementByRow[j] > 0.0) {
3775 flip[numJ] = 1;
3776 offset[numJ] = lowerColumn;
3777 size[numJ++] = elementByRow[j];
3778 } else {
3779 flip[numJ] = -1;
3780 offset[numJ] = upperColumn;
3781 size[numJ++] = -elementByRow[j];
3782 lowerColumn = upperColumn;
3783 }
3784 knapsackOffset += elementByRow[j] * lowerColumn;
3785 }
3786 int jRow;
3787 for (iRow = 0; iRow < numberRows_; iRow++)
3788 whichRow[iRow] = iRow;
3789 ClpSimplex smallModel(this, numberRows_, whichRow, numJ, whichColumn, true, true, true);
3790 // modify rhs to allow for nonzero lower bounds
3791 //double * rowLower = smallModel.rowLower();
3792 //double * rowUpper = smallModel.rowUpper();
3793 //const double * columnLower = smallModel.columnLower();
3794 //const double * columnUpper = smallModel.columnUpper();
3795 const CoinPackedMatrix * matrix = smallModel.matrix();
3796 const double * element = matrix->getElements();
3797 const int * row = matrix->getIndices();
3798 const CoinBigIndex * columnStart = matrix->getVectorStarts();
3799 const int * columnLength = matrix->getVectorLengths();
3800 const double * objective = smallModel.objective();
3801 //double objectiveOffset=0.0;
3802 // would use for fixed?
3803 CoinZeroN(rhsOffset, numberRows_);
3804 double * rowActivity = smallModel.primalRowSolution();
3805 CoinZeroN(rowActivity, numberRows_);
3806 maxSize -= knapsackOffset;
3807 minSize -= knapsackOffset;
3808 // now generate
3809 int i;
3810 int iStack = numJ;
3811 for (i = 0; i < numJ; i++) {
3812 stack[i] = 0;
3813 }
3814 double tooMuch = 10.0 * maxSize + 10000;
3815 stack[numJ] = 1;
3816 size[numJ] = tooMuch;
3817 bound[numJ] = 0;
3818 double sum = tooMuch;
3819 // allow for all zero being OK
3820 stack[numJ-1] = -1;
3821 sum -= size[numJ-1];
3822 numberOutput = 0;
3823 int nelCreate = 0;
3824 /* typeRun is - 0 for initial sizes
3825 1 for build
3826 2 for reconstruct
3827 */
3828 int typeRun = buildObj ? 1 : 0;
3829 if (reConstruct >= 0) {
3830 assert (buildRow && buildElement);
3831 typeRun = 2;
3832 }
3833 if (typeRun == 1)
3834 buildStart[0] = 0;
3835 while (iStack >= 0) {
3836 if (sum >= minSize && sum <= maxSize) {
3837 double checkSize = 0.0;
3838 bool good = true;
3839 int nRow = 0;
3840 double obj = 0.0;
3841 CoinZeroN(rowActivity, numberRows_);
3842 for (iColumn = 0; iColumn < numJ; iColumn++) {
3843 int iValue = stack[iColumn];
3844 if (iValue > bound[iColumn]) {
3845 good = false;
3846 break;
3847 } else {
3848 double realValue = offset[iColumn] + flip[iColumn] * iValue;
3849 if (realValue) {
3850 obj += objective[iColumn] * realValue;
3851 for (CoinBigIndex j = columnStart[iColumn];
3852 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3853 double value = element[j] * realValue;
3854 int kRow = row[j];
3855 if (rowActivity[kRow]) {
3856 rowActivity[kRow] += value;
3857 if (!rowActivity[kRow])
3858 rowActivity[kRow] = 1.0e-100;
3859 } else {
3860 build[nRow++] = kRow;
3861 rowActivity[kRow] = value;
3862 }
3863 }
3864 }
3865 }
3866 }
3867 if (good) {
3868 for (jRow = 0; jRow < nRow; jRow++) {
3869 int kRow = build[jRow];
3870 double value = rowActivity[kRow];
3871 if (value > high[kRow] || value < lo[kRow]) {
3872 good = false;
3873 break;
3874 }
3875 }
3876 }
3877 if (good) {
3878 if (typeRun == 1) {
3879 buildObj[numberOutput] = obj;
3880 for (jRow = 0; jRow < nRow; jRow++) {
3881 int kRow = build[jRow];
3882 double value = rowActivity[kRow];
3883 if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
3884 buildElement[nelCreate] = value;
3885 buildRow[nelCreate++] = kRow;
3886 }
3887 }
3888 buildStart[numberOutput+1] = nelCreate;
3889 } else if (!typeRun) {
3890 for (jRow = 0; jRow < nRow; jRow++) {
3891 int kRow = build[jRow];
3892 double value = rowActivity[kRow];
3893 if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
3894 nelCreate++;
3895 }
3896 }
3897 }
3898 if (typeRun == 2 && reConstruct == numberOutput) {
3899 // build and exit
3900 nelCreate = 0;
3901 for (iColumn = 0; iColumn < numJ; iColumn++) {
3902 int iValue = stack[iColumn];
3903 double realValue = offset[iColumn] + flip[iColumn] * iValue;
3904 if (realValue) {
3905 buildRow[nelCreate] = whichColumn[iColumn];
3906 buildElement[nelCreate++] = realValue;
3907 }
3908 }
3909 numberOutput = 1;
3910 for (i = 0; i < numJ; i++) {
3911 bound[i] = 0;
3912 }
3913 break;
3914 }
3915 numberOutput++;
3916 if (numberOutput > maxNumber) {
3917 nelCreate = -numberOutput;
3918 numberOutput = -1;
3919 for (i = 0; i < numJ; i++) {
3920 bound[i] = 0;
3921 }
3922 break;
3923 } else if (typeRun == 1 && numberOutput == maxNumber) {
3924 // On second run
3925 for (i = 0; i < numJ; i++) {
3926 bound[i] = 0;
3927 }
3928 break;
3929 }
3930 for (int j = 0; j < numJ; j++) {
3931 checkSize += stack[j] * size[j];
3932 }
3933 assert (fabs(sum - checkSize) < 1.0e-3);
3934 }
3935 for (jRow = 0; jRow < nRow; jRow++) {
3936 int kRow = build[jRow];
3937 rowActivity[kRow] = 0.0;
3938 }
3939 }
3940 if (sum > maxSize || stack[iStack] > bound[iStack]) {
3941 sum -= size[iStack] * stack[iStack];
3942 stack[iStack--] = 0;
3943 if (iStack >= 0) {
3944 stack[iStack] ++;
3945 sum += size[iStack];
3946 }
3947 } else {
3948 // must be less
3949 // add to last possible
3950 iStack = numJ - 1;
3951 sum += size[iStack];
3952 stack[iStack]++;
3953 }
3954 }
3955 //printf("%d will be created\n",numberOutput);
3956 delete [] whichColumn;
3957 delete [] whichRow;
3958 delete [] bound;
3959 delete [] stack;
3960 delete [] flip;
3961 delete [] size;
3962 delete [] offset;
3963 delete [] rhsOffset;
3964 delete [] build;
3965 delete [] markRow;
3966 delete [] lo;
3967 delete [] high;
3968 return nelCreate;
3969}
3970// Quick try at cleaning up duals if postsolve gets wrong
3971void
3972ClpSimplexOther::cleanupAfterPostsolve()
3973{
3974 // First mark singleton equality rows
3975 char * mark = new char [ numberRows_];
3976 memset(mark, 0, numberRows_);
3977 const int * row = matrix_->getIndices();
3978 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3979 const int * columnLength = matrix_->getVectorLengths();
3980 const double * element = matrix_->getElements();
3981 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3982 for (CoinBigIndex j = columnStart[iColumn];
3983 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3984 int iRow = row[j];
3985 if (mark[iRow])
3986 mark[iRow] = 2;
3987 else
3988 mark[iRow] = 1;
3989 }
3990 }
3991 // for now just == rows
3992 for (int iRow = 0; iRow < numberRows_; iRow++) {
3993 if (rowUpper_[iRow] > rowLower_[iRow])
3994 mark[iRow] = 3;
3995 }
3996 double dualTolerance = dblParam_[ClpDualTolerance];
3997 double primalTolerance = dblParam_[ClpPrimalTolerance];
3998 int numberCleaned = 0;
3999 double maxmin = optimizationDirection_;
4000 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4001 double dualValue = reducedCost_[iColumn] * maxmin;
4002 double primalValue = columnActivity_[iColumn];
4003 double lower = columnLower_[iColumn];
4004 double upper = columnUpper_[iColumn];
4005 int way = 0;
4006 switch(getColumnStatus(iColumn)) {
4007
4008 case basic:
4009 // dual should be zero
4010 if (dualValue > dualTolerance) {
4011 way = -1;
4012 } else if (dualValue < -dualTolerance) {
4013 way = 1;
4014 }
4015 break;
4016 case ClpSimplex::isFixed:
4017 break;
4018 case atUpperBound:
4019 // dual should not be positive
4020 if (dualValue > dualTolerance) {
4021 way = -1;
4022 }
4023 break;
4024 case atLowerBound:
4025 // dual should not be negative
4026 if (dualValue < -dualTolerance) {
4027 way = 1;
4028 }
4029 break;
4030 case superBasic:
4031 case isFree:
4032 if (primalValue < upper - primalTolerance) {
4033 // dual should not be negative
4034 if (dualValue < -dualTolerance) {
4035 way = 1;
4036 }
4037 }
4038 if (primalValue > lower + primalTolerance) {
4039 // dual should not be positive
4040 if (dualValue > dualTolerance) {
4041 way = -1;
4042 }
4043 }
4044 break;
4045 }
4046 if (way) {
4047 // see if can find singleton row
4048 for (CoinBigIndex j = columnStart[iColumn];
4049 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4050 int iRow = row[j];
4051 if (mark[iRow] == 1) {
4052 double value = element[j];
4053 // dj - addDual*value == 0.0
4054 double addDual = dualValue / value;
4055 dual_[iRow] += addDual;
4056 reducedCost_[iColumn] = 0.0;
4057 numberCleaned++;
4058 break;
4059 }
4060 }
4061 }
4062 }
4063 delete [] mark;
4064#ifdef CLP_INVESTIGATE
4065 printf("cleanupAfterPostsolve cleaned up %d columns\n", numberCleaned);
4066#endif
4067 // Redo
4068 memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
4069 matrix_->transposeTimes(-1.0, dual_, reducedCost_);
4070 checkSolutionInternal();
4071}
4072// Returns gub version of model or NULL
4073ClpSimplex *
4074ClpSimplexOther::gubVersion(int * whichRows, int * whichColumns,
4075 int neededGub,
4076 int factorizationFrequency)
4077{
4078 // find gub
4079 int numberRows = this->numberRows();
4080 int numberColumns = this->numberColumns();
4081 int iRow, iColumn;
4082 int * columnIsGub = new int [numberColumns];
4083 const double * columnLower = this->columnLower();
4084 const double * columnUpper = this->columnUpper();
4085 int numberFixed=0;
4086 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4087 if (columnUpper[iColumn] == columnLower[iColumn]) {
4088 columnIsGub[iColumn]=-2;
4089 numberFixed++;
4090 } else if (columnLower[iColumn]>=0) {
4091 columnIsGub[iColumn]=-1;
4092 } else {
4093 columnIsGub[iColumn]=-3;
4094 }
4095 }
4096 CoinPackedMatrix * matrix = this->matrix();
4097 // get row copy
4098 CoinPackedMatrix rowCopy = *matrix;
4099 rowCopy.reverseOrdering();
4100 const int * column = rowCopy.getIndices();
4101 const int * rowLength = rowCopy.getVectorLengths();
4102 const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
4103 const double * element = rowCopy.getElements();
4104 int numberNonGub = 0;
4105 int numberEmpty = numberRows;
4106 int * rowIsGub = new int [numberRows];
4107 int smallestGubRow=-1;
4108 int count=numberColumns+1;
4109 double * rowLower = this->rowLower();
4110 double * rowUpper = this->rowUpper();
4111 // make sure we can get rid of upper bounds
4112 double * fixedRow = new double [numberRows];
4113 for (iRow = 0 ; iRow < numberRows ; iRow++) {
4114 double sumFixed=0.0;
4115 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
4116 int iColumn = column[j];
4117 double value = columnLower[iColumn];
4118 if (value)
4119 sumFixed += element[j] * value;
4120 }
4121 fixedRow[iRow]=rowUpper[iRow]-sumFixed;
4122 }
4123 for (iRow = numberRows - 1; iRow >= 0; iRow--) {
4124 bool gubRow = true;
4125 int numberInRow=0;
4126 double sumFixed=0.0;
4127 double gap = fixedRow[iRow]-1.0e-12;
4128 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
4129 int iColumn = column[j];
4130 if (columnIsGub[iColumn]!=-2) {
4131 if (element[j] != 1.0||columnIsGub[iColumn]==-3||
4132 columnUpper[iColumn]-columnLower[iColumn]<gap) {
4133 gubRow = false;
4134 break;
4135 } else {
4136 numberInRow++;
4137 if (columnIsGub[iColumn] >= 0) {
4138 gubRow = false;
4139 break;
4140 }
4141 }
4142 } else {
4143 sumFixed += columnLower[iColumn]*element[j];
4144 }
4145 }
4146 if (!gubRow) {
4147 whichRows[numberNonGub++] = iRow;
4148 rowIsGub[iRow] = -1;
4149 } else if (numberInRow) {
4150 if (numberInRow<count) {
4151 count = numberInRow;
4152 smallestGubRow=iRow;
4153 }
4154 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
4155 int iColumn = column[j];
4156 if (columnIsGub[iColumn]!=-2)
4157 columnIsGub[iColumn] = iRow;
4158 }
4159 rowIsGub[iRow] = 0;
4160 } else {
4161 // empty row!
4162 whichRows[--numberEmpty] = iRow;
4163 rowIsGub[iRow] = -2;
4164 if (sumFixed>rowUpper[iRow]+1.0e-4||
4165 sumFixed<rowLower[iRow]-1.0e-4) {
4166 fprintf(stderr,"******** No infeasible empty rows - please!\n");
4167 abort();
4168 }
4169 }
4170 }
4171 delete [] fixedRow;
4172 char message[100];
4173 int numberGub = numberEmpty - numberNonGub;
4174 if (numberGub >= neededGub) {
4175 sprintf(message,"%d gub rows", numberGub);
4176 handler_->message(CLP_GENERAL2, messages_)
4177 << message << CoinMessageEol;
4178 int numberNormal = 0;
4179 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4180 if (columnIsGub[iColumn] < 0 && columnIsGub[iColumn] !=-2) {
4181 whichColumns[numberNormal++] = iColumn;
4182 }
4183 }
4184 if (!numberNormal) {
4185 sprintf(message,"Putting back one gub row to make non-empty");
4186 handler_->message(CLP_GENERAL2, messages_)
4187 << message << CoinMessageEol;
4188 rowIsGub[smallestGubRow]=-1;
4189 whichRows[numberNonGub++] = smallestGubRow;
4190 for (int j = rowStart[smallestGubRow];
4191 j < rowStart[smallestGubRow] + rowLength[smallestGubRow]; j++) {
4192 int iColumn = column[j];
4193 if (columnIsGub[iColumn]>=0) {
4194 columnIsGub[iColumn]=-4;
4195 whichColumns[numberNormal++] = iColumn;
4196 }
4197 }
4198 }
4199 std::sort(whichRows,whichRows+numberNonGub);
4200 std::sort(whichColumns,whichColumns+numberNormal);
4201 double * lower = CoinCopyOfArray(this->rowLower(),numberRows);
4202 double * upper = CoinCopyOfArray(this->rowUpper(),numberRows);
4203 // leave empty rows at end
4204 numberEmpty = numberRows-numberEmpty;
4205 const int * row = matrix->getIndices();
4206 const int * columnLength = matrix->getVectorLengths();
4207 const CoinBigIndex * columnStart = matrix->getVectorStarts();
4208 const double * elementByColumn = matrix->getElements();
4209 // Fixed at end
4210 int put2 = numberColumns-numberFixed;
4211 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4212 if (columnIsGub[iColumn] ==-2) {
4213 whichColumns[put2++] = iColumn;
4214 double value = columnLower[iColumn];
4215 for (int j = columnStart[iColumn];
4216 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4217 int iRow = row[j];
4218 if (lower[iRow]>-1.0e20)
4219 lower[iRow] -= value*element[j];
4220 if (upper[iRow]<1.0e20)
4221 upper[iRow] -= value*element[j];
4222 }
4223 }
4224 }
4225 int put = numberNormal;
4226 ClpSimplex * model2 =
4227 new ClpSimplex(this, numberNonGub, whichRows , numberNormal, whichColumns);
4228 // scale
4229 double * scaleArray = new double [numberRows];
4230 for (int i=0;i<numberRows;i++) {
4231 scaleArray[i]=1.0;
4232 if (rowIsGub[i]==-1) {
4233 double largest = 1.0e-30;
4234 double smallest = 1.0e30;
4235 for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
4236 int iColumn = column[j];
4237 if (columnIsGub[iColumn]!=-2) {
4238 double value =fabs(element[j]);
4239 largest = CoinMax(value,largest);
4240 smallest = CoinMin(value,smallest);
4241 }
4242 }
4243 double scale = CoinMax(0.001,1.0/sqrt(largest*smallest));
4244 scaleArray[i]=scale;
4245 if (lower[i]>-1.0e30)
4246 lower[i] *= scale;
4247 if (upper[i]<1.0e30)
4248 upper[i] *= scale;
4249 }
4250 }
4251 // scale partial matrix
4252 {
4253 CoinPackedMatrix * matrix = model2->matrix();
4254 const int * row = matrix->getIndices();
4255 const int * columnLength = matrix->getVectorLengths();
4256 const CoinBigIndex * columnStart = matrix->getVectorStarts();
4257 double * element = matrix->getMutableElements();
4258 for (int i=0;i<numberNormal;i++) {
4259 for (int j = columnStart[i];
4260 j < columnStart[i] + columnLength[i]; j++) {
4261 int iRow = row[j];
4262 iRow = whichRows[iRow];
4263 double scaleBy = scaleArray[iRow];
4264 element[j] *= scaleBy;
4265 }
4266 }
4267 }
4268 // adjust rhs
4269 double * rowLower = model2->rowLower();
4270 double * rowUpper = model2->rowUpper();
4271 for (int i=0;i<numberNonGub;i++) {
4272 int iRow = whichRows[i];
4273 rowLower[i] = lower[iRow];
4274 rowUpper[i] = upper[iRow];
4275 }
4276 int numberGubColumns = numberColumns - put - numberFixed;
4277 CoinBigIndex numberElements=0;
4278 int * temp1 = new int [numberRows+1];
4279 // get counts
4280 memset(temp1,0,numberRows*sizeof(int));
4281 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4282 int iGub = columnIsGub[iColumn];
4283 if (iGub>=0) {
4284 numberElements += columnLength[iColumn]-1;
4285 temp1[iGub]++;
4286 }
4287 }
4288 /* Optional but means can eventually simplify coding
4289 could even add in fixed slacks to deal with
4290 singularities - but should not be necessary */
4291 int numberSlacks=0;
4292 for (int i = 0; i < numberRows; i++) {
4293 if (rowIsGub[i]>=0) {
4294 if (lower[i]<upper[i]) {
4295 numberSlacks++;
4296 temp1[i]++;
4297 }
4298 }
4299 }
4300 int * gubStart = new int [numberGub+1];
4301 numberGub=0;
4302 gubStart[0]=0;
4303 for (int i = 0; i < numberRows; i++) {
4304 if (rowIsGub[i]>=0) {
4305 rowIsGub[i]=numberGub;
4306 gubStart[numberGub+1]=gubStart[numberGub]+temp1[i];
4307 temp1[numberGub]=0;
4308 lower[numberGub]=lower[i];
4309 upper[numberGub]=upper[i];
4310 whichRows[numberNonGub+numberGub]=i;
4311 numberGub++;
4312 }
4313 }
4314 int numberGubColumnsPlus = numberGubColumns + numberSlacks;
4315 double * lowerColumn2 = new double [numberGubColumnsPlus];
4316 CoinFillN(lowerColumn2, numberGubColumnsPlus, 0.0);
4317 double * upperColumn2 = new double [numberGubColumnsPlus];
4318 CoinFillN(upperColumn2, numberGubColumnsPlus, COIN_DBL_MAX);
4319 int * start2 = new int[numberGubColumnsPlus+1];
4320 int * row2 = new int[numberElements];
4321 double * element2 = new double[numberElements];
4322 double * cost2 = new double [numberGubColumnsPlus];
4323 CoinFillN(cost2, numberGubColumnsPlus, 0.0);
4324 const double * cost = this->objective();
4325 put = numberNormal;
4326 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4327 int iGub = columnIsGub[iColumn];
4328 if (iGub>=0) {
4329 // TEMP
4330 //this->setColUpper(iColumn,COIN_DBL_MAX);
4331 iGub = rowIsGub[iGub];
4332 assert (iGub>=0);
4333 int kPut = put+gubStart[iGub]+temp1[iGub];
4334 temp1[iGub]++;
4335 whichColumns[kPut]=iColumn;
4336 }
4337 }
4338 for (int i = 0; i < numberRows; i++) {
4339 if (rowIsGub[i]>=0) {
4340 int iGub = rowIsGub[i];
4341 if (lower[iGub]<upper[iGub]) {
4342 int kPut = put+gubStart[iGub]+temp1[iGub];
4343 temp1[iGub]++;
4344 whichColumns[kPut]=iGub+numberColumns;
4345 }
4346 }
4347 }
4348 //this->primal(1); // TEMP
4349 // redo rowIsGub to give lookup
4350 for (int i=0;i<numberRows;i++)
4351 rowIsGub[i]=-1;
4352 for (int i=0;i<numberNonGub;i++)
4353 rowIsGub[whichRows[i]]=i;
4354 start2[0]=0;
4355 numberElements = 0;
4356 for (int i=0;i<numberGubColumnsPlus;i++) {
4357 int iColumn = whichColumns[put++];
4358 if (iColumn<numberColumns) {
4359 cost2[i] = cost[iColumn];
4360 lowerColumn2[i] = columnLower[iColumn];
4361 upperColumn2[i] = columnUpper[iColumn];
4362 upperColumn2[i] = COIN_DBL_MAX;
4363 for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4364 int iRow = row[j];
4365 double scaleBy = scaleArray[iRow];
4366 iRow = rowIsGub[iRow];
4367 if (iRow >= 0) {
4368 row2[numberElements] = iRow;
4369 element2[numberElements++] = elementByColumn[j]*scaleBy;
4370 }
4371 }
4372 } else {
4373 // slack
4374 int iGub = iColumn-numberColumns;
4375 double slack = upper[iGub]-lower[iGub];
4376 assert (upper[iGub]<1.0e20);
4377 lower[iGub]=upper[iGub];
4378 cost2[i] = 0;
4379 lowerColumn2[i] = 0;
4380 upperColumn2[i] = slack;
4381 upperColumn2[i] = COIN_DBL_MAX;
4382 }
4383 start2[i+1] = numberElements;
4384 }
4385 // clean up bounds on variables
4386 for (int iSet=0;iSet<numberGub;iSet++) {
4387 double lowerValue=0.0;
4388 for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) {
4389 lowerValue += lowerColumn2[i];
4390 }
4391 assert (lowerValue<upper[iSet]+1.0e-6);
4392 double gap = CoinMax(0.0,upper[iSet]-lowerValue);
4393 for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) {
4394 if (upperColumn2[i]<1.0e30) {
4395 upperColumn2[i] = CoinMin(upperColumn2[i],
4396 lowerColumn2[i]+gap);
4397 }
4398 }
4399 }
4400 sprintf(message,"** Before adding matrix there are %d rows and %d columns",
4401 model2->numberRows(), model2->numberColumns());
4402 handler_->message(CLP_GENERAL2, messages_)
4403 << message << CoinMessageEol;
4404 delete [] scaleArray;
4405 delete [] temp1;
4406 model2->setFactorizationFrequency(factorizationFrequency);
4407 ClpDynamicMatrix * newMatrix =
4408 new ClpDynamicMatrix(model2, numberGub,
4409 numberGubColumnsPlus, gubStart,
4410 lower, upper,
4411 start2, row2, element2, cost2,
4412 lowerColumn2, upperColumn2);
4413 delete [] gubStart;
4414 delete [] lowerColumn2;
4415 delete [] upperColumn2;
4416 delete [] start2;
4417 delete [] row2;
4418 delete [] element2;
4419 delete [] cost2;
4420 delete [] lower;
4421 delete [] upper;
4422 model2->replaceMatrix(newMatrix,true);
4423#ifdef EVERY_ITERATION
4424 {
4425 ClpDynamicMatrix * gubMatrix =
4426 dynamic_cast< ClpDynamicMatrix*>(model2->clpMatrix());
4427 assert(gubMatrix);
4428 gubMatrix->writeMps("gub.mps");
4429 }
4430#endif
4431 delete [] columnIsGub;
4432 delete [] rowIsGub;
4433 newMatrix->switchOffCheck();
4434#ifdef EVERY_ITERATION
4435 newMatrix->setRefreshFrequency(1/*000*/);
4436#else
4437 newMatrix->setRefreshFrequency(1000);
4438#endif
4439 sprintf(message,
4440 "** While after adding matrix there are %d rows and %d columns",
4441 model2->numberRows(), model2->numberColumns());
4442 handler_->message(CLP_GENERAL2, messages_)
4443 << message << CoinMessageEol;
4444 model2->setSpecialOptions(4); // exactly to bound
4445 // Scaling off (done by hand)
4446 model2->scaling(0);
4447 return model2;
4448 } else {
4449 delete [] columnIsGub;
4450 delete [] rowIsGub;
4451 return NULL;
4452 }
4453}
4454// Sets basis from original
4455void
4456ClpSimplexOther::setGubBasis(ClpSimplex &original,const int * whichRows,
4457 const int * whichColumns)
4458{
4459 ClpDynamicMatrix * gubMatrix =
4460 dynamic_cast< ClpDynamicMatrix*>(clpMatrix());
4461 assert(gubMatrix);
4462 int numberGubColumns = gubMatrix->numberGubColumns();
4463 int numberNormal = gubMatrix->firstDynamic();
4464 //int lastOdd = gubMatrix->firstAvailable();
4465 //int numberTotalColumns = numberNormal + numberGubColumns;
4466 //assert (numberTotalColumns==numberColumns+numberSlacks);
4467 int numberRows = original.numberRows();
4468 int numberColumns = original.numberColumns();
4469 int * columnIsGub = new int [numberColumns];
4470 int numberNonGub = gubMatrix->numberStaticRows();
4471 //assert (firstOdd==numberNormal);
4472 double * solution = primalColumnSolution();
4473 double * originalSolution = original.primalColumnSolution();
4474 const double * upperSet = gubMatrix->upperSet();
4475 // Column copy of GUB part
4476 int numberSets = gubMatrix->numberSets();
4477 const int * startSet = gubMatrix->startSets();
4478 const CoinBigIndex * columnStart = gubMatrix->startColumn();
4479 const double * columnLower = gubMatrix->columnLower();
4480#ifdef TRY_IMPROVE
4481 const double * columnUpper = gubMatrix->columnUpper();
4482 const double * lowerSet = gubMatrix->lowerSet();
4483 const double * element = gubMatrix->element();
4484 const int * row = gubMatrix->row();
4485 bool allPositive=true;
4486 double * rowActivity = new double[numberNonGub];
4487 memset(rowActivity, 0, numberNonGub*sizeof(double));
4488 {
4489 // Non gub contribution
4490 const double * element = matrix_->getElements();
4491 const int * row = matrix_->getIndices();
4492 const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4493 const int * columnLength = matrix_->getVectorLengths();
4494 for (int i=0;i<numberNormal;i++) {
4495 int iColumn = whichColumns[i];
4496 double value = originalSolution[iColumn];
4497 if (value) {
4498 for (CoinBigIndex j = columnStart[i];
4499 j < columnStart[i] + columnLength[i]; j++) {
4500 int iRow = row[j];
4501 rowActivity[iRow] += value*element[j];
4502 }
4503 }
4504 }
4505 }
4506 double * newSolution = new double [numberGubColumns];
4507 int * slacks = new int [numberSets];
4508 for (int i=0;i<numberSets;i++) {
4509 double sum=0.0;
4510 int iSlack=-1;
4511 for (int j=startSet[i];j<startSet[i+1];j++) {
4512 gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
4513 int iColumn = whichColumns[j+numberNormal];
4514 if (iColumn<numberColumns) {
4515 columnIsGub[iColumn] = whichRows[numberNonGub+i];
4516 double value = originalSolution[iColumn];
4517 sum += value;
4518 newSolution[j]=value;
4519 for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
4520 int iRow = row[k];
4521 rowActivity[iRow] += value*element[k];
4522 if (element[k] < 0.0)
4523 allPositive=false;
4524 }
4525 if (columnStart[j]==columnStart[j+1])
4526 iSlack=j;
4527 } else {
4528 newSolution[j]=0.0;
4529 iSlack=j;
4530 allPositive=false; // for now
4531 }
4532 }
4533 slacks[i]=iSlack;
4534 if (sum>upperSet[i]+1.0e-8) {
4535 double gap = sum-upperSet[i];
4536 if (iSlack>=0) {
4537 double value=newSolution[iSlack];
4538 if (value>0.0) {
4539 double down = CoinMin(gap,value);
4540 gap -= down;
4541 sum -= down;
4542 newSolution[iSlack] = value-down;
4543 }
4544 }
4545 if (gap>1.0e-8) {
4546 for (int j=startSet[i];j<startSet[i+1];j++) {
4547 int iColumn = whichColumns[j+numberNormal];
4548 if (newSolution[j]>0.0&&iColumn<numberColumns) {
4549 double value = newSolution[j];
4550 double down = CoinMin(gap,value);
4551 gap -= down;
4552 sum -= down;
4553 newSolution[iSlack] = value-down;
4554 for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
4555 int iRow = row[k];
4556 rowActivity[iRow] -= down*element[k];
4557 }
4558 }
4559 }
4560 }
4561 assert (gap<1.0e-8);
4562 } else if (sum<lowerSet[i]-1.0e-8) {
4563 double gap = lowerSet[i]-sum;
4564 if (iSlack>=0) {
4565 double value=newSolution[iSlack];
4566 if (value<columnUpper[iSlack]) {
4567 double up = CoinMin(gap,columnUpper[iSlack]-value);
4568 gap -= up;
4569 sum += up;
4570 newSolution[iSlack] = value+up;
4571 }
4572 }
4573 if (gap>1.0e-8) {
4574 for (int j=startSet[i];j<startSet[i+1];j++) {
4575 int iColumn = whichColumns[j+numberNormal];
4576 if (newSolution[j]<columnUpper[j]&&iColumn<numberColumns) {
4577 double value = newSolution[j];
4578 double up = CoinMin(gap,columnUpper[j]-value);
4579 gap -= up;
4580 sum += up;
4581 newSolution[iSlack] = value+up;
4582 for (CoinBigIndex k = columnStart[j]; k < columnStart[j+1] ; k++) {
4583 int iRow = row[k];
4584 rowActivity[iRow] += up*element[k];
4585 }
4586 }
4587 }
4588 }
4589 assert (gap<1.0e-8);
4590 }
4591 if (fabs(sum-upperSet[i])>1.0e-7)
4592 printf("Sum for set %d is %g - lower %g, upper %g\n",i,
4593 sum,lowerSet[i],upperSet[i]);
4594 }
4595 if (allPositive) {
4596 // See if we can improve solution
4597 // first reduce if over
4598 double * gaps = new double [numberNonGub];
4599 double direction = optimizationDirection_;
4600 const double * cost = gubMatrix->cost();
4601 bool over=false;
4602 for (int i=0;i<numberNonGub;i++) {
4603 double activity = rowActivity[i];
4604 gaps[i]=0.0;
4605 if (activity>rowUpper_[i]+1.0e-6) {
4606 gaps[i]=activity-rowUpper_[i];
4607 over=true;
4608 }
4609 }
4610 double * weights = new double [numberGubColumns];
4611 int * which = new int [numberGubColumns];
4612 int * whichSet = new int [numberGubColumns];
4613 if (over) {
4614 int n=0;
4615 for (int i=0;i<numberSets;i++) {
4616 int iSlack = slacks[i];
4617 if (iSlack<0||newSolution[iSlack]>upperSet[i]-1.0e-8)
4618 continue;
4619 double slackCost = cost[iSlack]*direction;
4620 for (int j=startSet[i];j<startSet[i+1];j++) {
4621 whichSet[j]=i;
4622 double value = newSolution[j];
4623 double thisCost = cost[j]*direction;
4624 if (value>columnLower[j]&&j!=iSlack) {
4625 if(thisCost<slackCost) {
4626 double sum = 1.0e-30;
4627 for (CoinBigIndex k = columnStart[j];
4628 k < columnStart[j+1] ; k++) {
4629 int iRow = row[k];
4630 sum += gaps[iRow]*element[k];
4631 }
4632 which[n]=j;
4633 // big drop and small difference in cost better
4634 weights[n++]=(slackCost-thisCost)/sum;
4635 } else {
4636 // slack better anyway
4637 double move = value-columnLower[j];
4638 newSolution[iSlack]=CoinMin(upperSet[i],
4639 newSolution[iSlack]+move);
4640 newSolution[j]=columnLower[j];
4641 for (CoinBigIndex k = columnStart[j];
4642 k < columnStart[j+1] ; k++) {
4643 int iRow = row[k];
4644 rowActivity[iRow] -= move*element[k];
4645 }
4646 }
4647 }
4648 }
4649 }
4650 // sort
4651 CoinSort_2(weights,weights+n,which);
4652 for (int i=0;i<n;i++) {
4653 int j= which[i];
4654 int iSet = whichSet[j];
4655 int iSlack = slacks[iSet];
4656 assert (iSlack>=0);
4657 double move = 0.0;
4658 for (CoinBigIndex k = columnStart[j];
4659 k < columnStart[j+1] ; k++) {
4660 int iRow = row[k];
4661 if(rowActivity[iRow]-rowUpper_[iRow]>move*element[k]) {
4662 move = (rowActivity[iRow]-rowUpper_[iRow])/element[k];
4663 }
4664 }
4665 move=CoinMin(move,newSolution[j]-columnLower[j]);
4666 if (move) {
4667 newSolution[j] -= move;
4668 newSolution[iSlack] += move;
4669 for (CoinBigIndex k = columnStart[j];
4670 k < columnStart[j+1] ; k++) {
4671 int iRow = row[k];
4672 rowActivity[iRow] -= move*element[k];
4673 }
4674 }
4675 }
4676 }
4677 delete [] whichSet;
4678 delete [] which;
4679 delete [] weights;
4680 delete [] gaps;
4681 // redo original status!
4682 for (int i=0;i<numberSets;i++) {
4683 int numberBasic=0;
4684 int numberNewBasic=0;
4685 int j1=-1;
4686 int j2=-1;
4687 for (int j=startSet[i];j<startSet[i+1];j++) {
4688 if (newSolution[j]>columnLower[j]) {
4689 numberNewBasic++;
4690 j2=j;
4691 }
4692 int iOrig = whichColumns[j+numberNormal];
4693 if (iOrig<numberColumns) {
4694 if (original.getColumnStatus(iOrig)!=ClpSimplex::atLowerBound) {
4695 numberBasic++;
4696 j1=j;
4697 }
4698 } else {
4699 int iSet = iOrig - numberColumns;
4700 int iRow = whichRows[iSet+numberNonGub];
4701 if (original.getRowStatus(iRow)==ClpSimplex::basic) {
4702 numberBasic++;
4703 j1=j;
4704 abort();
4705 }
4706 }
4707 }
4708 if (numberBasic==1&&numberNewBasic==1&&
4709 j1!=j2) {
4710 int iOrig1=whichColumns[j1+numberNormal];
4711 int iOrig2=whichColumns[j2+numberNormal];
4712 ClpSimplex::Status status1 = original.getColumnStatus(iOrig1);
4713 ClpSimplex::Status status2 = original.getColumnStatus(iOrig2);
4714 originalSolution[iOrig1] = newSolution[j1];
4715 originalSolution[iOrig2] = newSolution[j2];
4716 original.setColumnStatus(iOrig1,status2);
4717 original.setColumnStatus(iOrig2,status1);
4718 }
4719 }
4720 }
4721 delete [] newSolution;
4722 delete [] slacks;
4723 delete [] rowActivity;
4724#else
4725 for (int i=0;i<numberSets;i++) {
4726 for (int j=startSet[i];j<startSet[i+1];j++) {
4727 gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
4728 int iColumn = whichColumns[j+numberNormal];
4729 if (iColumn<numberColumns) {
4730 columnIsGub[iColumn] = whichRows[numberNonGub+i];
4731 }
4732 }
4733 }
4734#endif
4735 int * numberKey = new int [numberRows];
4736 memset(numberKey,0,numberRows*sizeof(int));
4737 for (int i=0;i<numberGubColumns;i++) {
4738 int iOrig = whichColumns[i+numberNormal];
4739 if (iOrig<numberColumns) {
4740 if (original.getColumnStatus(iOrig)==ClpSimplex::basic) {
4741 int iRow = columnIsGub[iOrig];
4742 assert (iRow>=0);
4743 numberKey[iRow]++;
4744 }
4745 } else {
4746 // Set slack
4747 int iSet = iOrig - numberColumns;
4748 int iRow = whichRows[iSet+numberNonGub];
4749 if (original.getRowStatus(iRow)==ClpSimplex::basic)
4750 numberKey[iRow]++;
4751 }
4752 }
4753 /* Before going into cleanMatrix we need
4754 gub status set (inSmall just means basic and active)
4755 row status set
4756 */
4757 for (int i = 0; i < numberSets; i++) {
4758 gubMatrix->setStatus(i,ClpSimplex::isFixed);
4759 }
4760 for (int i = 0; i < numberGubColumns; i++) {
4761 int iOrig = whichColumns[i+numberNormal];
4762 if (iOrig<numberColumns) {
4763 ClpSimplex::Status status = original.getColumnStatus(iOrig);
4764 if (status==ClpSimplex::atUpperBound) {
4765 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atUpperBound);
4766 } else if (status==ClpSimplex::atLowerBound) {
4767 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound);
4768 } else if (status==ClpSimplex::basic) {
4769 int iRow = columnIsGub[iOrig];
4770 assert (iRow>=0);
4771 assert(numberKey[iRow]);
4772 if (numberKey[iRow]==1)
4773 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey);
4774 else
4775 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall);
4776 }
4777 } else {
4778 // slack
4779 int iSet = iOrig - numberColumns;
4780 int iRow = whichRows[iSet+numberNonGub];
4781 if (original.getRowStatus(iRow)==ClpSimplex::basic
4782#ifdef TRY_IMPROVE
4783 ||newSolution[i]>columnLower[i]+1.0e-8
4784#endif
4785 ) {
4786 assert(numberKey[iRow]);
4787 if (numberKey[iRow]==1)
4788 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey);
4789 else
4790 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall);
4791 } else {
4792 gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound);
4793 }
4794 }
4795 }
4796 // deal with sets without key
4797 for (int i = 0; i < numberSets; i++) {
4798 int iRow = whichRows[numberNonGub+i];
4799 if (!numberKey[iRow]) {
4800 double upper = upperSet[i]-1.0e-7;
4801 if (original.getRowStatus(iRow)==ClpSimplex::basic)
4802 gubMatrix->setStatus(i,ClpSimplex::basic);
4803 // If not at lb make key otherwise one with smallest number els
4804 double largest=0.0;
4805 int fewest=numberRows+1;
4806 int chosen=-1;
4807 for (int j=startSet[i];j<startSet[i+1];j++) {
4808 int length=columnStart[j+1]-columnStart[j];
4809 int iOrig = whichColumns[j+numberNormal];
4810 double value;
4811 if (iOrig<numberColumns) {
4812#ifdef TRY_IMPROVE
4813 value=newSolution[j]-columnLower[j];
4814#else
4815 value = originalSolution[iOrig]-columnLower[j];
4816#endif
4817 if (value>upper)
4818 gubMatrix->setStatus(i,ClpSimplex::atLowerBound);
4819 } else {
4820 // slack - take value as 0.0 as will win on length
4821 value=0.0;
4822 }
4823 if (value>largest+1.0e-8) {
4824 largest=value;
4825 fewest=length;
4826 chosen=j;
4827 } else if (fabs(value-largest)<=1.0e-8&&length<fewest) {
4828 largest=value;
4829 fewest=length;
4830 chosen=j;
4831 }
4832 }
4833 assert(chosen>=0);
4834 if (gubMatrix->getStatus(i)!=ClpSimplex::basic) {
4835 // set as key
4836 for (int j=startSet[i];j<startSet[i+1];j++) {
4837 if (j!=chosen)
4838 gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
4839 else
4840 gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::soloKey);
4841 }
4842 }
4843 }
4844 }
4845 for (int i = 0; i < numberNormal; i++) {
4846 int iOrig = whichColumns[i];
4847 setColumnStatus(i,original.getColumnStatus(iOrig));
4848 solution[i]=originalSolution[iOrig];
4849 }
4850 for (int i = 0; i < numberNonGub; i++) {
4851 int iOrig = whichRows[i];
4852 setRowStatus(i,original.getRowStatus(iOrig));
4853 }
4854 // Fill in current matrix
4855 gubMatrix->initialProblem();
4856 delete [] numberKey;
4857 delete [] columnIsGub;
4858}
4859// Restores basis to original
4860void
4861ClpSimplexOther::getGubBasis(ClpSimplex &original,const int * whichRows,
4862 const int * whichColumns) const
4863{
4864 ClpDynamicMatrix * gubMatrix =
4865 dynamic_cast< ClpDynamicMatrix*>(clpMatrix());
4866 assert(gubMatrix);
4867 int numberGubColumns = gubMatrix->numberGubColumns();
4868 int numberNormal = gubMatrix->firstDynamic();
4869 //int lastOdd = gubMatrix->firstAvailable();
4870 //int numberRows = original.numberRows();
4871 int numberColumns = original.numberColumns();
4872 int numberNonGub = gubMatrix->numberStaticRows();
4873 //assert (firstOdd==numberNormal);
4874 double * solution = primalColumnSolution();
4875 double * originalSolution = original.primalColumnSolution();
4876 int numberSets = gubMatrix->numberSets();
4877 const double * cost = original.objective();
4878 int lastOdd = gubMatrix->firstAvailable();
4879 //assert (numberTotalColumns==numberColumns+numberSlacks);
4880 int numberRows = original.numberRows();
4881 //int numberStaticRows = gubMatrix->numberStaticRows();
4882 const int * startSet = gubMatrix->startSets();
4883 unsigned char * status = original.statusArray();
4884 unsigned char * rowStatus = status+numberColumns;
4885 //assert (firstOdd==numberNormal);
4886 for (int i=0;i<numberSets;i++) {
4887 int iRow = whichRows[i+numberNonGub];
4888 original.setRowStatus(iRow,ClpSimplex::atLowerBound);
4889 }
4890 const int * id = gubMatrix->id();
4891 const double * columnLower = gubMatrix->columnLower();
4892 const double * columnUpper = gubMatrix->columnUpper();
4893 for (int i = 0; i < numberGubColumns; i++) {
4894 int iOrig = whichColumns[i+numberNormal];
4895 if (iOrig<numberColumns) {
4896 if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) {
4897 originalSolution[iOrig] = columnUpper[i];
4898 status[iOrig] = 2;
4899 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound && columnLower) {
4900 originalSolution[iOrig] = columnLower[i];
4901 status[iOrig] = 3;
4902 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) {
4903 int iSet = gubMatrix->whichSet(i);
4904 originalSolution[iOrig] = gubMatrix->keyValue(iSet);
4905 status[iOrig] = 1;
4906 } else {
4907 originalSolution[iOrig] = 0.0;
4908 status[iOrig] = 4;
4909 }
4910 } else {
4911 // slack
4912 int iSet = iOrig - numberColumns;
4913 int iRow = whichRows[iSet+numberNonGub];
4914 if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) {
4915 original.setRowStatus(iRow,ClpSimplex::atLowerBound);
4916 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound) {
4917 original.setRowStatus(iRow,ClpSimplex::atUpperBound);
4918 } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) {
4919 original.setRowStatus(iRow,ClpSimplex::basic);
4920 }
4921 }
4922 }
4923 for (int i = 0; i < numberNormal; i++) {
4924 int iOrig = whichColumns[i];
4925 ClpSimplex::Status thisStatus = getStatus(i);
4926 if (thisStatus == ClpSimplex::basic)
4927 status[iOrig] = 1;
4928 else if (thisStatus == ClpSimplex::atLowerBound)
4929 status[iOrig] = 3;
4930 else if (thisStatus == ClpSimplex::atUpperBound)
4931 status[iOrig] = 2;
4932 else if (thisStatus == ClpSimplex::isFixed)
4933 status[iOrig] = 5;
4934 else
4935 abort();
4936 originalSolution[iOrig] = solution[i];
4937 }
4938 for (int i = numberNormal; i < lastOdd; i++) {
4939 int iOrig = whichColumns[id[i-numberNormal] + numberNormal];
4940 if (iOrig<numberColumns) {
4941 ClpSimplex::Status thisStatus = getStatus(i);
4942 if (thisStatus == ClpSimplex::basic)
4943 status[iOrig] = 1;
4944 else if (thisStatus == ClpSimplex::atLowerBound)
4945 status[iOrig] = 3;
4946 else if (thisStatus == ClpSimplex::atUpperBound)
4947 status[iOrig] = 2;
4948 else if (thisStatus == ClpSimplex::isFixed)
4949 status[iOrig] = 5;
4950 else
4951 abort();
4952 originalSolution[iOrig] = solution[i];
4953 } else {
4954 // slack (basic probably)
4955 int iSet = iOrig - numberColumns;
4956 int iRow = whichRows[iSet+numberNonGub];
4957 ClpSimplex::Status thisStatus = getStatus(i);
4958 if (thisStatus == ClpSimplex::atLowerBound)
4959 thisStatus = ClpSimplex::atUpperBound;
4960 else if (thisStatus == ClpSimplex::atUpperBound)
4961 thisStatus = ClpSimplex::atLowerBound;
4962 original.setRowStatus(iRow,thisStatus);
4963 }
4964 }
4965 for (int i = 0; i < numberNonGub; i++) {
4966 int iOrig = whichRows[i];
4967 ClpSimplex::Status thisStatus = getRowStatus(i);
4968 if (thisStatus == ClpSimplex::basic)
4969 rowStatus[iOrig] = 1;
4970 else if (thisStatus == ClpSimplex::atLowerBound)
4971 rowStatus[iOrig] = 3;
4972 else if (thisStatus == ClpSimplex::atUpperBound)
4973 rowStatus[iOrig] = 2;
4974 else if (thisStatus == ClpSimplex::isFixed)
4975 rowStatus[iOrig] = 5;
4976 else
4977 abort();
4978 }
4979 int * numberKey = new int [numberRows];
4980 memset(numberKey,0,numberRows*sizeof(int));
4981 for (int i=0;i<numberSets;i++) {
4982 int iRow = whichRows[i+numberNonGub];
4983 for (int j=startSet[i];j<startSet[i+1];j++) {
4984 int iOrig = whichColumns[j+numberNormal];
4985 if (iOrig<numberColumns) {
4986 if (original.getColumnStatus(iOrig)==ClpSimplex::basic) {
4987 numberKey[iRow]++;
4988 }
4989 } else {
4990 // slack
4991 if (original.getRowStatus(iRow)==ClpSimplex::basic)
4992 numberKey[iRow]++;
4993 }
4994 }
4995 }
4996 for (int i=0;i<numberSets;i++) {
4997 int iRow = whichRows[i+numberNonGub];
4998 if (!numberKey[iRow]) {
4999 original.setRowStatus(iRow,ClpSimplex::basic);
5000 }
5001 }
5002 delete [] numberKey;
5003 double objValue = 0.0;
5004 for (int i = 0; i < numberColumns; i++)
5005 objValue += cost[i] * originalSolution[i];
5006 //printf("objective value is %g\n", objValue);
5007}
5008