source: trunk/Clp/src/ClpSimplexOther.cpp @ 1525

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 149.1 KB
Line 
1/* $Id: ClpSimplexOther.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinPragma.hpp"
6
7#include <math.h>
8
9#include "CoinHelperFunctions.hpp"
10#include "ClpSimplexOther.hpp"
11#include "ClpSimplexDual.hpp"
12#include "ClpSimplexPrimal.hpp"
13#include "ClpEventHandler.hpp"
14#include "ClpHelperFunctions.hpp"
15#include "ClpFactorization.hpp"
16#include "ClpDualRowDantzig.hpp"
17#include "CoinPackedMatrix.hpp"
18#include "CoinIndexedVector.hpp"
19#include "CoinBuild.hpp"
20#include "CoinMpsIO.hpp"
21#include "ClpMessage.hpp"
22#include <cfloat>
23#include <cassert>
24#include <string>
25#include <stdio.h>
26#include <iostream>
27/* Dual ranging.
28   This computes increase/decrease in cost for each given variable and corresponding
29   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
30   and numberColumns.. for artificials/slacks.
31   For non-basic variables the sequence number will be that of the non-basic variables.
32
33   Up to user to provide correct length arrays.
34
35*/
36void ClpSimplexOther::dualRanging(int numberCheck, const int * which,
37                                  double * costIncreased, int * sequenceIncreased,
38                                  double * costDecreased, int * sequenceDecreased,
39                                  double * valueIncrease, double * valueDecrease)
40{
41     rowArray_[1]->clear();
42     columnArray_[1]->clear();
43     // long enough for rows+columns
44     assert(rowArray_[3]->capacity() >= numberRows_ + numberColumns_);
45     rowArray_[3]->clear();
46     int * backPivot = rowArray_[3]->getIndices();
47     int i;
48     for ( i = 0; i < numberRows_ + numberColumns_; i++) {
49          backPivot[i] = -1;
50     }
51     for (i = 0; i < numberRows_; i++) {
52          int iSequence = pivotVariable_[i];
53          backPivot[iSequence] = i;
54     }
55     // dualTolerance may be zero if from CBC.  In fact use that fact
56     bool inCBC = !dualTolerance_;
57     if (inCBC)
58          assert (integerType_);
59     dualTolerance_ = dblParam_[ClpDualTolerance];
60     double * arrayX = rowArray_[0]->denseVector();
61     for ( i = 0; i < numberCheck; i++) {
62          rowArray_[0]->clear();
63          //rowArray_[0]->checkClear();
64          //rowArray_[1]->checkClear();
65          //columnArray_[1]->checkClear();
66          columnArray_[0]->clear();
67          //columnArray_[0]->checkClear();
68          int iSequence = which[i];
69          if (iSequence < 0) {
70               costIncreased[i] = 0.0;
71               sequenceIncreased[i] = -1;
72               costDecreased[i] = 0.0;
73               sequenceDecreased[i] = -1;
74               continue;
75          }
76          double costIncrease = COIN_DBL_MAX;
77          double costDecrease = COIN_DBL_MAX;
78          int sequenceIncrease = -1;
79          int sequenceDecrease = -1;
80          if (valueIncrease) {
81               assert (valueDecrease);
82               valueIncrease[i] = iSequence < numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_];
83               valueDecrease[i] = valueIncrease[i];
84          }
85
86          switch(getStatus(iSequence)) {
87
88          case basic: {
89               // non-trvial
90               // Get pivot row
91               int iRow = backPivot[iSequence];
92               assert (iRow >= 0);
93               double plusOne = 1.0;
94               rowArray_[0]->createPacked(1, &iRow, &plusOne);
95               factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
96               // put row of tableau in rowArray[0] and columnArray[0]
97               matrix_->transposeTimes(this, -1.0,
98                                       rowArray_[0], columnArray_[1], columnArray_[0]);
99               double alphaIncrease;
100               double alphaDecrease;
101               // do ratio test up and down
102               checkDualRatios(rowArray_[0], columnArray_[0], costIncrease, sequenceIncrease, alphaIncrease,
103                               costDecrease, sequenceDecrease, alphaDecrease);
104               if (!inCBC) {
105                    if (valueIncrease) {
106                         if (sequenceIncrease >= 0)
107                              valueIncrease[i] = primalRanging1(sequenceIncrease, iSequence);
108                         if (sequenceDecrease >= 0)
109                              valueDecrease[i] = primalRanging1(sequenceDecrease, iSequence);
110                    }
111               } else {
112                    int number = rowArray_[0]->getNumElements();
113                    double scale2 = 0.0;
114                    int j;
115                    for (j = 0; j < number; j++) {
116                         scale2 += arrayX[j] * arrayX[j];
117                    }
118                    scale2 = 1.0 / sqrt(scale2);
119                    valueIncrease[i] = scale2;
120                    if (sequenceIncrease >= 0) {
121                         double djValue = dj_[sequenceIncrease];
122                         if (fabs(djValue) > 10.0 * dualTolerance_) {
123                              // we are going to use for cutoff so be exact
124                              costIncrease = fabs(djValue / alphaIncrease);
125                              /* Not sure this is good idea as I don't think correct e.g.
126                                 suppose a continuous variable has dj slightly greater. */
127                              if(false && sequenceIncrease < numberColumns_ && integerType_[sequenceIncrease]) {
128                                   // can improve
129                                   double movement = (columnScale_ == NULL) ? 1.0 :
130                                                     rhsScale_ * inverseColumnScale_[sequenceIncrease];
131                                   costIncrease = CoinMax(fabs(djValue * movement), costIncrease);
132                              }
133                         } else {
134                              costIncrease = 0.0;
135                         }
136                    }
137                    if (sequenceDecrease >= 0) {
138                         double djValue = dj_[sequenceDecrease];
139                         if (fabs(djValue) > 10.0 * dualTolerance_) {
140                              // we are going to use for cutoff so be exact
141                              costDecrease = fabs(djValue / alphaDecrease);
142                              if(sequenceDecrease < numberColumns_ && integerType_[sequenceDecrease]) {
143                                   // can improve
144                                   double movement = (columnScale_ == NULL) ? 1.0 :
145                                                     rhsScale_ * inverseColumnScale_[sequenceDecrease];
146                                   costDecrease = CoinMax(fabs(djValue * movement), costDecrease);
147                              }
148                         } else {
149                              costDecrease = 0.0;
150                         }
151                    }
152                    costIncrease *= scale2;
153                    costDecrease *= scale2;
154               }
155          }
156          break;
157          case isFixed:
158               break;
159          case isFree:
160          case superBasic:
161               costIncrease = 0.0;
162               costDecrease = 0.0;
163               sequenceIncrease = iSequence;
164               sequenceDecrease = iSequence;
165               break;
166          case atUpperBound:
167               costIncrease = CoinMax(0.0, -dj_[iSequence]);
168               sequenceIncrease = iSequence;
169               if (valueIncrease)
170                    valueIncrease[i] = primalRanging1(iSequence, iSequence);
171               break;
172          case atLowerBound:
173               costDecrease = CoinMax(0.0, dj_[iSequence]);
174               sequenceDecrease = iSequence;
175               if (valueIncrease)
176                    valueDecrease[i] = primalRanging1(iSequence, iSequence);
177               break;
178          }
179          double scaleFactor;
180          if (rowScale_) {
181               if (iSequence < numberColumns_)
182                    scaleFactor = 1.0 / (objectiveScale_ * columnScale_[iSequence]);
183               else
184                    scaleFactor = rowScale_[iSequence-numberColumns_] / objectiveScale_;
185          } else {
186               scaleFactor = 1.0 / objectiveScale_;
187          }
188          if (costIncrease < 1.0e30)
189               costIncrease *= scaleFactor;
190          if (costDecrease < 1.0e30)
191               costDecrease *= scaleFactor;
192          if (optimizationDirection_ == 1.0) {
193               costIncreased[i] = costIncrease;
194               sequenceIncreased[i] = sequenceIncrease;
195               costDecreased[i] = costDecrease;
196               sequenceDecreased[i] = sequenceDecrease;
197          } else if (optimizationDirection_ == -1.0) {
198               costIncreased[i] = costDecrease;
199               sequenceIncreased[i] = sequenceDecrease;
200               costDecreased[i] = costIncrease;
201               sequenceDecreased[i] = sequenceIncrease;
202               if (valueIncrease) {
203                    double temp = valueIncrease[i];
204                    valueIncrease[i] = valueDecrease[i];
205                    valueDecrease[i] = temp;
206               }
207          } else if (optimizationDirection_ == 0.0) {
208               // !!!!!! ???
209               costIncreased[i] = COIN_DBL_MAX;
210               sequenceIncreased[i] = -1;
211               costDecreased[i] = COIN_DBL_MAX;
212               sequenceDecreased[i] = -1;
213          } else {
214               abort();
215          }
216     }
217     rowArray_[0]->clear();
218     //rowArray_[1]->clear();
219     //columnArray_[1]->clear();
220     columnArray_[0]->clear();
221     //rowArray_[3]->clear();
222     if (!optimizationDirection_)
223          printf("*** ????? Ranging with zero optimization costs\n");
224}
225/*
226   Row array has row part of pivot row
227   Column array has column part.
228   This is used in dual ranging
229*/
230void
231ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
232                                 CoinIndexedVector * columnArray,
233                                 double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
234                                 double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
235{
236     double acceptablePivot = 1.0e-9;
237     double * work;
238     int number;
239     int * which;
240     int iSection;
241
242     double thetaDown = 1.0e31;
243     double thetaUp = 1.0e31;
244     int sequenceDown = -1;
245     int sequenceUp = -1;
246     double alphaDown = 0.0;
247     double alphaUp = 0.0;
248
249     int addSequence;
250
251     for (iSection = 0; iSection < 2; iSection++) {
252
253          int i;
254          if (!iSection) {
255               work = rowArray->denseVector();
256               number = rowArray->getNumElements();
257               which = rowArray->getIndices();
258               addSequence = numberColumns_;
259          } else {
260               work = columnArray->denseVector();
261               number = columnArray->getNumElements();
262               which = columnArray->getIndices();
263               addSequence = 0;
264          }
265
266          for (i = 0; i < number; i++) {
267               int iSequence = which[i];
268               int iSequence2 = iSequence + addSequence;
269               double alpha = work[i];
270               if (fabs(alpha) < acceptablePivot)
271                    continue;
272               double oldValue = dj_[iSequence2];
273
274               switch(getStatus(iSequence2)) {
275
276               case basic:
277                    break;
278               case ClpSimplex::isFixed:
279                    break;
280               case isFree:
281               case superBasic:
282                    // treat dj as if zero
283                    thetaDown = 0.0;
284                    thetaUp = 0.0;
285                    sequenceDown = iSequence2;
286                    sequenceUp = iSequence2;
287                    break;
288               case atUpperBound:
289                    if (alpha > 0.0) {
290                         // test up
291                         if (oldValue + thetaUp * alpha > dualTolerance_) {
292                              thetaUp = (dualTolerance_ - oldValue) / alpha;
293                              sequenceUp = iSequence2;
294                              alphaUp = alpha;
295                         }
296                    } else {
297                         // test down
298                         if (oldValue - thetaDown * alpha > dualTolerance_) {
299                              thetaDown = -(dualTolerance_ - oldValue) / alpha;
300                              sequenceDown = iSequence2;
301                              alphaDown = alpha;
302                         }
303                    }
304                    break;
305               case atLowerBound:
306                    if (alpha < 0.0) {
307                         // test up
308                         if (oldValue + thetaUp * alpha < - dualTolerance_) {
309                              thetaUp = -(dualTolerance_ + oldValue) / alpha;
310                              sequenceUp = iSequence2;
311                              alphaUp = alpha;
312                         }
313                    } else {
314                         // test down
315                         if (oldValue - thetaDown * alpha < -dualTolerance_) {
316                              thetaDown = (dualTolerance_ + oldValue) / alpha;
317                              sequenceDown = iSequence2;
318                              alphaDown = alpha;
319                         }
320                    }
321                    break;
322               }
323          }
324     }
325     if (sequenceUp >= 0) {
326          costIncrease = thetaUp;
327          sequenceIncrease = sequenceUp;
328          alphaIncrease = alphaUp;
329     }
330     if (sequenceDown >= 0) {
331          costDecrease = thetaDown;
332          sequenceDecrease = sequenceDown;
333          alphaDecrease = alphaDown;
334     }
335}
336/** Primal ranging.
337    This computes increase/decrease in value for each given variable and corresponding
338    sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
339    and numberColumns.. for artificials/slacks.
340    For basic variables the sequence number will be that of the basic variables.
341
342    Up to user to provide correct length arrays.
343
344    When here - guaranteed optimal
345*/
346void
347ClpSimplexOther::primalRanging(int numberCheck, const int * which,
348                               double * valueIncreased, int * sequenceIncreased,
349                               double * valueDecreased, int * sequenceDecreased)
350{
351     rowArray_[0]->clear();
352     rowArray_[1]->clear();
353     lowerIn_ = -COIN_DBL_MAX;
354     upperIn_ = COIN_DBL_MAX;
355     valueIn_ = 0.0;
356     for ( int i = 0; i < numberCheck; i++) {
357          int iSequence = which[i];
358          double valueIncrease = COIN_DBL_MAX;
359          double valueDecrease = COIN_DBL_MAX;
360          int sequenceIncrease = -1;
361          int sequenceDecrease = -1;
362
363          switch(getStatus(iSequence)) {
364
365          case basic:
366          case isFree:
367          case superBasic:
368               // Easy
369               valueDecrease = CoinMax(0.0, upper_[iSequence] - solution_[iSequence]);
370               valueIncrease = CoinMax(0.0, solution_[iSequence] - lower_[iSequence]);
371               sequenceDecrease = iSequence;
372               sequenceIncrease = iSequence;
373               break;
374          case isFixed:
375          case atUpperBound:
376          case atLowerBound: {
377               // Non trivial
378               // Other bound is ignored
379               unpackPacked(rowArray_[1], iSequence);
380               factorization_->updateColumn(rowArray_[2], rowArray_[1]);
381               // Get extra rows
382               matrix_->extendUpdated(this, rowArray_[1], 0);
383               // do ratio test
384               checkPrimalRatios(rowArray_[1], 1);
385               if (pivotRow_ >= 0) {
386                    valueIncrease = theta_;
387                    sequenceIncrease = pivotVariable_[pivotRow_];
388               }
389               checkPrimalRatios(rowArray_[1], -1);
390               if (pivotRow_ >= 0) {
391                    valueDecrease = theta_;
392                    sequenceDecrease = pivotVariable_[pivotRow_];
393               }
394               rowArray_[1]->clear();
395          }
396          break;
397          }
398          double scaleFactor;
399          if (rowScale_) {
400               if (iSequence < numberColumns_)
401                    scaleFactor = columnScale_[iSequence] / rhsScale_;
402               else
403                    scaleFactor = 1.0 / (rowScale_[iSequence-numberColumns_] * rhsScale_);
404          } else {
405               scaleFactor = 1.0 / rhsScale_;
406          }
407          if (valueIncrease < 1.0e30)
408               valueIncrease *= scaleFactor;
409          else
410               valueIncrease = COIN_DBL_MAX;
411          if (valueDecrease < 1.0e30)
412               valueDecrease *= scaleFactor;
413          else
414               valueDecrease = COIN_DBL_MAX;
415          valueIncreased[i] = valueIncrease;
416          sequenceIncreased[i] = sequenceIncrease;
417          valueDecreased[i] = valueDecrease;
418          sequenceDecreased[i] = sequenceDecrease;
419     }
420}
421// Returns new value of whichOther when whichIn enters basis
422double
423ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
424{
425     rowArray_[0]->clear();
426     rowArray_[1]->clear();
427     int iSequence = whichIn;
428     double newValue = solution_[whichOther];
429     double alphaOther = 0.0;
430     Status status = getStatus(iSequence);
431     assert (status == atLowerBound || status == atUpperBound);
432     int wayIn = (status == atLowerBound) ? 1 : -1;
433
434     switch(getStatus(iSequence)) {
435
436     case basic:
437     case isFree:
438     case superBasic:
439          assert (whichIn == whichOther);
440          // Easy
441          newValue = wayIn > 0 ? upper_[iSequence] : lower_[iSequence];
442          break;
443     case isFixed:
444     case atUpperBound:
445     case atLowerBound:
446          // Non trivial
447     {
448          // Other bound is ignored
449          unpackPacked(rowArray_[1], iSequence);
450          factorization_->updateColumn(rowArray_[2], rowArray_[1]);
451          // Get extra rows
452          matrix_->extendUpdated(this, rowArray_[1], 0);
453          // do ratio test
454          double acceptablePivot = 1.0e-7;
455          double * work = rowArray_[1]->denseVector();
456          int number = rowArray_[1]->getNumElements();
457          int * which = rowArray_[1]->getIndices();
458
459          // we may need to swap sign
460          double way = wayIn;
461          double theta = 1.0e30;
462          for (int iIndex = 0; iIndex < number; iIndex++) {
463
464               int iRow = which[iIndex];
465               double alpha = work[iIndex] * way;
466               int iPivot = pivotVariable_[iRow];
467               if (iPivot == whichOther) {
468                    alphaOther = alpha;
469                    continue;
470               }
471               double oldValue = solution_[iPivot];
472               if (fabs(alpha) > acceptablePivot) {
473                    if (alpha > 0.0) {
474                         // basic variable going towards lower bound
475                         double bound = lower_[iPivot];
476                         oldValue -= bound;
477                         if (oldValue - theta * alpha < 0.0) {
478                              theta = CoinMax(0.0, oldValue / alpha);
479                         }
480                    } else {
481                         // basic variable going towards upper bound
482                         double bound = upper_[iPivot];
483                         oldValue = oldValue - bound;
484                         if (oldValue - theta * alpha > 0.0) {
485                              theta = CoinMax(0.0, oldValue / alpha);
486                         }
487                    }
488               }
489          }
490          if (whichIn != whichOther) {
491               if (theta < 1.0e30)
492                    newValue -= theta * alphaOther;
493               else
494                    newValue = alphaOther > 0.0 ? -1.0e30 : 1.0e30;
495          } else {
496               newValue += theta * wayIn;
497          }
498     }
499     rowArray_[1]->clear();
500     break;
501     }
502     double scaleFactor;
503     if (rowScale_) {
504          if (whichOther < numberColumns_)
505               scaleFactor = columnScale_[whichOther] / rhsScale_;
506          else
507               scaleFactor = 1.0 / (rowScale_[whichOther-numberColumns_] * rhsScale_);
508     } else {
509          scaleFactor = 1.0 / rhsScale_;
510     }
511     if (newValue < 1.0e29)
512          if (newValue > -1.0e29)
513               newValue *= scaleFactor;
514          else
515               newValue = -COIN_DBL_MAX;
516     else
517          newValue = COIN_DBL_MAX;
518     return newValue;
519}
520/*
521   Row array has pivot column
522   This is used in primal ranging
523*/
524void
525ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
526                                   int direction)
527{
528     // sequence stays as row number until end
529     pivotRow_ = -1;
530     double acceptablePivot = 1.0e-7;
531     double * work = rowArray->denseVector();
532     int number = rowArray->getNumElements();
533     int * which = rowArray->getIndices();
534
535     // we need to swap sign if going down
536     double way = direction;
537     theta_ = 1.0e30;
538     for (int iIndex = 0; iIndex < number; iIndex++) {
539
540          int iRow = which[iIndex];
541          double alpha = work[iIndex] * way;
542          int iPivot = pivotVariable_[iRow];
543          double oldValue = solution_[iPivot];
544          if (fabs(alpha) > acceptablePivot) {
545               if (alpha > 0.0) {
546                    // basic variable going towards lower bound
547                    double bound = lower_[iPivot];
548                    oldValue -= bound;
549                    if (oldValue - theta_ * alpha < 0.0) {
550                         pivotRow_ = iRow;
551                         theta_ = CoinMax(0.0, oldValue / alpha);
552                    }
553               } else {
554                    // basic variable going towards upper bound
555                    double bound = upper_[iPivot];
556                    oldValue = oldValue - bound;
557                    if (oldValue - theta_ * alpha > 0.0) {
558                         pivotRow_ = iRow;
559                         theta_ = CoinMax(0.0, oldValue / alpha);
560                    }
561               }
562          }
563     }
564}
565/* Write the basis in MPS format to the specified file.
566   If writeValues true writes values of structurals
567   (and adds VALUES to end of NAME card)
568
569   Row and column names may be null.
570   formatType is
571   <ul>
572   <li> 0 - normal
573   <li> 1 - extra accuracy
574   <li> 2 - IEEE hex (later)
575   </ul>
576
577   Returns non-zero on I/O error
578
579   This is based on code contributed by Thorsten Koch
580*/
581int
582ClpSimplexOther::writeBasis(const char *filename,
583                            bool writeValues,
584                            int formatType) const
585{
586     formatType = CoinMax(0, formatType);
587     formatType = CoinMin(2, formatType);
588     if (!writeValues)
589          formatType = 0;
590     // See if INTEL if IEEE
591     if (formatType == 2) {
592          // test intel here and add 1 if not intel
593          double value = 1.0;
594          char x[8];
595          memcpy(x, &value, 8);
596          if (x[0] == 63) {
597               formatType ++; // not intel
598          } else {
599               assert (x[0] == 0);
600          }
601     }
602
603     char number[20];
604     FILE * fp = fopen(filename, "w");
605     if (!fp)
606          return -1;
607
608     // NAME card
609
610     if (strcmp(strParam_[ClpProbName].c_str(), "") == 0) {
611          fprintf(fp, "NAME          BLANK      ");
612     } else {
613          fprintf(fp, "NAME          %s       ", strParam_[ClpProbName].c_str());
614     }
615     if (formatType >= 2)
616          fprintf(fp, "FREEIEEE");
617     else if (writeValues)
618          fprintf(fp, "VALUES");
619     // finish off name
620     fprintf(fp, "\n");
621     int iRow = 0;
622     for(int iColumn = 0; iColumn < numberColumns_; iColumn++) {
623          bool printit = false;
624          if( getColumnStatus(iColumn) == ClpSimplex::basic) {
625               printit = true;
626               // Find non basic row
627               for(; iRow < numberRows_; iRow++) {
628                    if (getRowStatus(iRow) != ClpSimplex::basic)
629                         break;
630               }
631               if (lengthNames_) {
632                    if (iRow != numberRows_) {
633                         fprintf(fp, " %s %-8s       %s",
634                                 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
635                                 columnNames_[iColumn].c_str(),
636                                 rowNames_[iRow].c_str());
637                         iRow++;
638                    } else {
639                         // Allow for too many basics!
640                         fprintf(fp, " BS %-8s       ",
641                                 columnNames_[iColumn].c_str());
642                         // Dummy row name if values
643                         if (writeValues)
644                              fprintf(fp, "      _dummy_");
645                    }
646               } else {
647                    // no names
648                    if (iRow != numberRows_) {
649                         fprintf(fp, " %s C%7.7d     R%7.7d",
650                                 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
651                                 iColumn, iRow);
652                         iRow++;
653                    } else {
654                         // Allow for too many basics!
655                         fprintf(fp, " BS C%7.7d", iColumn);
656                         // Dummy row name if values
657                         if (writeValues)
658                              fprintf(fp, "      _dummy_");
659                    }
660               }
661          } else  {
662               if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
663                    printit = true;
664                    if (lengthNames_)
665                         fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
666                    else
667                         fprintf(fp, " UL C%7.7d", iColumn);
668                    // Dummy row name if values
669                    if (writeValues)
670                         fprintf(fp, "      _dummy_");
671               }
672          }
673          if (printit && writeValues) {
674               // add value
675               CoinConvertDouble(0, formatType, columnActivity_[iColumn], number);
676               fprintf(fp, "     %s", number);
677          }
678          if (printit)
679               fprintf(fp, "\n");
680     }
681     fprintf(fp, "ENDATA\n");
682     fclose(fp);
683     return 0;
684}
685// Read a basis from the given filename
686int
687ClpSimplexOther::readBasis(const char *fileName)
688{
689     int status = 0;
690     bool canOpen = false;
691     if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
692          // stdin
693          canOpen = true;
694     } else {
695          FILE *fp = fopen(fileName, "r");
696          if (fp) {
697               // can open - lets go for it
698               fclose(fp);
699               canOpen = true;
700          } else {
701               handler_->message(CLP_UNABLE_OPEN, messages_)
702                         << fileName << CoinMessageEol;
703               return -1;
704          }
705     }
706     CoinMpsIO m;
707     m.passInMessageHandler(handler_);
708     *m.messagesPointer() = coinMessages();
709     bool savePrefix = m.messageHandler()->prefix();
710     m.messageHandler()->setPrefix(handler_->prefix());
711     status = m.readBasis(fileName, "", columnActivity_, status_ + numberColumns_,
712                          status_,
713                          columnNames_, numberColumns_,
714                          rowNames_, numberRows_);
715     m.messageHandler()->setPrefix(savePrefix);
716     if (status >= 0) {
717          if (!status) {
718               // set values
719               int iColumn, iRow;
720               for (iRow = 0; iRow < numberRows_; iRow++) {
721                    if (getRowStatus(iRow) == atLowerBound)
722                         rowActivity_[iRow] = rowLower_[iRow];
723                    else if (getRowStatus(iRow) == atUpperBound)
724                         rowActivity_[iRow] = rowUpper_[iRow];
725               }
726               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
727                    if (getColumnStatus(iColumn) == atLowerBound)
728                         columnActivity_[iColumn] = columnLower_[iColumn];
729                    else if (getColumnStatus(iColumn) == atUpperBound)
730                         columnActivity_[iColumn] = columnUpper_[iColumn];
731               }
732          } else {
733               memset(rowActivity_, 0, numberRows_ * sizeof(double));
734               matrix_->times(-1.0, columnActivity_, rowActivity_);
735          }
736     } else {
737          // errors
738          handler_->message(CLP_IMPORT_ERRORS, messages_)
739                    << status << fileName << CoinMessageEol;
740     }
741     return status;
742}
743/* Creates dual of a problem if looks plausible
744   (defaults will always create model)
745   fractionRowRanges is fraction of rows allowed to have ranges
746   fractionColumnRanges is fraction of columns allowed to have ranges
747*/
748ClpSimplex *
749ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const
750{
751     const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
752     bool changed = false;
753     int numberChanged = 0;
754     int iColumn;
755     // check if we need to change bounds to rows
756     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
757          if (columnUpper_[iColumn] < 1.0e20 &&
758                    columnLower_[iColumn] > -1.0e20) {
759               changed = true;
760               numberChanged++;
761          }
762     }
763     int iRow;
764     int numberExtraRows = 0;
765     if (numberChanged <= fractionColumnRanges * numberColumns_) {
766          for (iRow = 0; iRow < numberRows_; iRow++) {
767               if (rowLower_[iRow] > -1.0e20 &&
768                         rowUpper_[iRow] < 1.0e20) {
769                    if (rowUpper_[iRow] != rowLower_[iRow])
770                         numberExtraRows++;
771               }
772          }
773          if (numberExtraRows > fractionRowRanges * numberRows_)
774               return NULL;
775     } else {
776          return NULL;
777     }
778     if (changed) {
779          ClpSimplex * model3 = new ClpSimplex(*model2);
780          CoinBuild build;
781          double one = 1.0;
782          int numberColumns = model3->numberColumns();
783          const double * columnLower = model3->columnLower();
784          const double * columnUpper = model3->columnUpper();
785          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
786               if (columnUpper[iColumn] < 1.0e20 &&
787                         columnLower[iColumn] > -1.0e20) {
788                    if (fabs(columnLower[iColumn]) < fabs(columnUpper[iColumn])) {
789                         double value = columnUpper[iColumn];
790                         model3->setColumnUpper(iColumn, COIN_DBL_MAX);
791                         build.addRow(1, &iColumn, &one, -COIN_DBL_MAX, value);
792                    } else {
793                         double value = columnLower[iColumn];
794                         model3->setColumnLower(iColumn, -COIN_DBL_MAX);
795                         build.addRow(1, &iColumn, &one, value, COIN_DBL_MAX);
796                    }
797               }
798          }
799          model3->addRows(build);
800          model2 = model3;
801     }
802     int numberColumns = model2->numberColumns();
803     const double * columnLower = model2->columnLower();
804     const double * columnUpper = model2->columnUpper();
805     int numberRows = model2->numberRows();
806     double * rowLower = CoinCopyOfArray(model2->rowLower(), numberRows);
807     double * rowUpper = CoinCopyOfArray(model2->rowUpper(), numberRows);
808
809     const double * objective = model2->objective();
810     CoinPackedMatrix * matrix = model2->matrix();
811     // get transpose
812     CoinPackedMatrix rowCopy = *matrix;
813     const int * row = matrix->getIndices();
814     const int * columnLength = matrix->getVectorLengths();
815     const CoinBigIndex * columnStart = matrix->getVectorStarts();
816     const double * elementByColumn = matrix->getElements();
817     double objOffset = 0.0;
818     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
819          double offset = 0.0;
820          double objValue = optimizationDirection_ * objective[iColumn];
821          if (columnUpper[iColumn] > 1.0e20) {
822               if (columnLower[iColumn] > -1.0e20)
823                    offset = columnLower[iColumn];
824          } else if (columnLower[iColumn] < -1.0e20) {
825               offset = columnUpper[iColumn];
826          } else {
827               // taken care of before
828               abort();
829          }
830          if (offset) {
831               objOffset += offset * objValue;
832               for (CoinBigIndex j = columnStart[iColumn];
833                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
834                    int iRow = row[j];
835                    if (rowLower[iRow] > -1.0e20)
836                         rowLower[iRow] -= offset * elementByColumn[j];
837                    if (rowUpper[iRow] < 1.0e20)
838                         rowUpper[iRow] -= offset * elementByColumn[j];
839               }
840          }
841     }
842     int * which = new int[numberRows+numberExtraRows];
843     rowCopy.reverseOrdering();
844     rowCopy.transpose();
845     double * fromRowsLower = new double[numberRows+numberExtraRows];
846     double * fromRowsUpper = new double[numberRows+numberExtraRows];
847     double * newObjective = new double[numberRows+numberExtraRows];
848     double * fromColumnsLower = new double[numberColumns];
849     double * fromColumnsUpper = new double[numberColumns];
850     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
851          double objValue = optimizationDirection_ * objective[iColumn];
852          // Offset is already in
853          if (columnUpper[iColumn] > 1.0e20) {
854               if (columnLower[iColumn] > -1.0e20) {
855                    fromColumnsLower[iColumn] = -COIN_DBL_MAX;
856                    fromColumnsUpper[iColumn] = objValue;
857               } else {
858                    // free
859                    fromColumnsLower[iColumn] = objValue;
860                    fromColumnsUpper[iColumn] = objValue;
861               }
862          } else if (columnLower[iColumn] < -1.0e20) {
863               fromColumnsLower[iColumn] = objValue;
864               fromColumnsUpper[iColumn] = COIN_DBL_MAX;
865          } else {
866               abort();
867          }
868     }
869     int kRow = 0;
870     int kExtraRow = numberRows;
871     for (iRow = 0; iRow < numberRows; iRow++) {
872          if (rowLower[iRow] < -1.0e20) {
873               assert (rowUpper[iRow] < 1.0e20);
874               newObjective[kRow] = -rowUpper[iRow];
875               fromRowsLower[kRow] = -COIN_DBL_MAX;
876               fromRowsUpper[kRow] = 0.0;
877               which[kRow] = iRow;
878               kRow++;
879          } else if (rowUpper[iRow] > 1.0e20) {
880               newObjective[kRow] = -rowLower[iRow];
881               fromRowsLower[kRow] = 0.0;
882               fromRowsUpper[kRow] = COIN_DBL_MAX;
883               which[kRow] = iRow;
884               kRow++;
885          } else {
886               if (rowUpper[iRow] == rowLower[iRow]) {
887                    newObjective[kRow] = -rowLower[iRow];
888                    fromRowsLower[kRow] = -COIN_DBL_MAX;;
889                    fromRowsUpper[kRow] = COIN_DBL_MAX;
890                    which[kRow] = iRow;
891                    kRow++;
892               } else {
893                    // range
894                    newObjective[kRow] = -rowUpper[iRow];
895                    fromRowsLower[kRow] = -COIN_DBL_MAX;
896                    fromRowsUpper[kRow] = 0.0;
897                    which[kRow] = iRow;
898                    kRow++;
899                    newObjective[kExtraRow] = -rowLower[iRow];
900                    fromRowsLower[kExtraRow] = 0.0;
901                    fromRowsUpper[kExtraRow] = COIN_DBL_MAX;
902                    which[kExtraRow] = iRow;
903                    kExtraRow++;
904               }
905          }
906     }
907     if (numberExtraRows) {
908          CoinPackedMatrix newCopy;
909          newCopy.setExtraGap(0.0);
910          newCopy.setExtraMajor(0.0);
911          newCopy.submatrixOfWithDuplicates(rowCopy, kExtraRow, which);
912          rowCopy = newCopy;
913     }
914     ClpSimplex * modelDual = new ClpSimplex();
915     modelDual->loadProblem(rowCopy, fromRowsLower, fromRowsUpper, newObjective,
916                            fromColumnsLower, fromColumnsUpper);
917     modelDual->setObjectiveOffset(objOffset);
918     modelDual->setDualBound(model2->dualBound());
919     modelDual->setInfeasibilityCost(model2->infeasibilityCost());
920     modelDual->setDualTolerance(model2->dualTolerance());
921     modelDual->setPrimalTolerance(model2->primalTolerance());
922     modelDual->setPerturbation(model2->perturbation());
923     modelDual->setSpecialOptions(model2->specialOptions());
924     modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
925     delete [] fromRowsLower;
926     delete [] fromRowsUpper;
927     delete [] fromColumnsLower;
928     delete [] fromColumnsUpper;
929     delete [] newObjective;
930     delete [] which;
931     delete [] rowLower;
932     delete [] rowUpper;
933     if (changed)
934          delete model2;
935     modelDual->createStatus();
936     return modelDual;
937}
938// Restores solution from dualized problem
939int
940ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem)
941{
942     int returnCode = 0;;
943     createStatus();
944     // Number of rows in dual problem was original number of columns
945     assert (numberColumns_ == dualProblem->numberRows());
946     // If slack on d-row basic then column at bound otherwise column basic
947     // If d-column basic then rhs tight
948     int numberBasic = 0;
949     int iRow, iColumn = 0;
950     // Get number of extra rows from ranges
951     int numberExtraRows = 0;
952     for (iRow = 0; iRow < numberRows_; iRow++) {
953          if (rowLower_[iRow] > -1.0e20 &&
954                    rowUpper_[iRow] < 1.0e20) {
955               if (rowUpper_[iRow] != rowLower_[iRow])
956                    numberExtraRows++;
957          }
958     }
959     const double * objective = this->objective();
960     const double * dualDual = dualProblem->dualRowSolution();
961     const double * dualDj = dualProblem->dualColumnSolution();
962     const double * dualSol = dualProblem->primalColumnSolution();
963     const double * dualActs = dualProblem->primalRowSolution();
964#if 0
965     ClpSimplex thisCopy = *this;
966     thisCopy.dual(); // for testing
967     const double * primalDual = thisCopy.dualRowSolution();
968     const double * primalDj = thisCopy.dualColumnSolution();
969     const double * primalSol = thisCopy.primalColumnSolution();
970     const double * primalActs = thisCopy.primalRowSolution();
971     char ss[] = {'F', 'B', 'U', 'L', 'S', 'F'};
972     printf ("Dual problem row info %d rows\n", dualProblem->numberRows());
973     for (iRow = 0; iRow < dualProblem->numberRows(); iRow++)
974          printf("%d at %c primal %g dual %g\n",
975                 iRow, ss[dualProblem->getRowStatus(iRow)],
976                 dualActs[iRow], dualDual[iRow]);
977     printf ("Dual problem column info %d columns\n", dualProblem->numberColumns());
978     for (iColumn = 0; iColumn < dualProblem->numberColumns(); iColumn++)
979          printf("%d at %c primal %g dual %g\n",
980                 iColumn, ss[dualProblem->getColumnStatus(iColumn)],
981                 dualSol[iColumn], dualDj[iColumn]);
982     printf ("Primal problem row info %d rows\n", thisCopy.numberRows());
983     for (iRow = 0; iRow < thisCopy.numberRows(); iRow++)
984          printf("%d at %c primal %g dual %g\n",
985                 iRow, ss[thisCopy.getRowStatus(iRow)],
986                 primalActs[iRow], primalDual[iRow]);
987     printf ("Primal problem column info %d columns\n", thisCopy.numberColumns());
988     for (iColumn = 0; iColumn < thisCopy.numberColumns(); iColumn++)
989          printf("%d at %c primal %g dual %g\n",
990                 iColumn, ss[thisCopy.getColumnStatus(iColumn)],
991                 primalSol[iColumn], primalDj[iColumn]);
992#endif
993     // position at bound information
994     int jColumn = numberRows_;
995     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
996          double objValue = optimizationDirection_ * objective[iColumn];
997          Status status = dualProblem->getRowStatus(iColumn);
998          double otherValue = COIN_DBL_MAX;
999          if (columnUpper_[iColumn] < 1.0e20 &&
1000                    columnLower_[iColumn] > -1.0e20) {
1001               if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) {
1002                    otherValue = columnUpper_[iColumn] + dualDj[jColumn];
1003               } else {
1004                    otherValue = columnLower_[iColumn] + dualDj[jColumn];
1005               }
1006               jColumn++;
1007          }
1008          if (status == basic) {
1009               // column is at bound
1010               if (otherValue == COIN_DBL_MAX) {
1011                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1012                    if (columnUpper_[iColumn] > 1.0e20) {
1013                         if (columnLower_[iColumn] > -1.0e20) {
1014                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1015                                   setColumnStatus(iColumn, atLowerBound);
1016                              else
1017                                   setColumnStatus(iColumn, isFixed);
1018                              columnActivity_[iColumn] = columnLower_[iColumn];
1019                         } else {
1020                              // free
1021                              setColumnStatus(iColumn, isFree);
1022                              columnActivity_[iColumn] = 0.0;
1023                         }
1024                    } else {
1025                         setColumnStatus(iColumn, atUpperBound);
1026                         columnActivity_[iColumn] = columnUpper_[iColumn];
1027                    }
1028               } else {
1029                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1030                    //printf("other dual sol %g\n",otherValue);
1031                    if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1032                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1033                              setColumnStatus(iColumn, atLowerBound);
1034                         else
1035                              setColumnStatus(iColumn, isFixed);
1036                         columnActivity_[iColumn] = columnLower_[iColumn];
1037                    } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1038                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1039                              setColumnStatus(iColumn, atUpperBound);
1040                         else
1041                              setColumnStatus(iColumn, isFixed);
1042                         columnActivity_[iColumn] = columnUpper_[iColumn];
1043                    } else {
1044                         abort();
1045                    }
1046               }
1047          } else {
1048               if (otherValue == COIN_DBL_MAX) {
1049                    // column basic
1050                    setColumnStatus(iColumn, basic);
1051                    numberBasic++;
1052                    if (columnLower_[iColumn] > -1.0e20) {
1053                         columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn];
1054                    } else if (columnUpper_[iColumn] < 1.0e20) {
1055                         columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn];
1056                    } else {
1057                         columnActivity_[iColumn] = -dualDual[iColumn];
1058                    }
1059                    reducedCost_[iColumn] = 0.0;
1060               } else {
1061                    // may be at other bound
1062                    //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1063                    if (dualProblem->getColumnStatus(jColumn - 1) != basic) {
1064                         // column basic
1065                         setColumnStatus(iColumn, basic);
1066                         numberBasic++;
1067                         //printf("Col %d otherV %g dualDual %g\n",iColumn,
1068                         // otherValue,dualDual[iColumn]);
1069                         columnActivity_[iColumn] = -dualDual[iColumn];
1070                         columnActivity_[iColumn] = otherValue;
1071                         reducedCost_[iColumn] = 0.0;
1072                    } else {
1073                         reducedCost_[iColumn] = objValue - dualActs[iColumn];
1074                         if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1075                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1076                                   setColumnStatus(iColumn, atLowerBound);
1077                              else
1078                                   setColumnStatus(iColumn, isFixed);
1079                              columnActivity_[iColumn] = columnLower_[iColumn];
1080                         } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1081                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1082                                   setColumnStatus(iColumn, atUpperBound);
1083                              else
1084                                   setColumnStatus(iColumn, isFixed);
1085                              columnActivity_[iColumn] = columnUpper_[iColumn];
1086                         } else {
1087                              abort();
1088                         }
1089                    }
1090               }
1091          }
1092     }
1093     // now rows
1094     int kExtraRow = jColumn;
1095     int numberRanges = 0;
1096     for (iRow = 0; iRow < numberRows_; iRow++) {
1097          Status status = dualProblem->getColumnStatus(iRow);
1098          if (status == basic) {
1099               // row is at bound
1100               dual_[iRow] = dualSol[iRow];;
1101          } else {
1102               // row basic
1103               setRowStatus(iRow, basic);
1104               numberBasic++;
1105               dual_[iRow] = 0.0;
1106          }
1107          if (rowLower_[iRow] < -1.0e20) {
1108               if (status == basic) {
1109                    rowActivity_[iRow] = rowUpper_[iRow];
1110                    setRowStatus(iRow, atUpperBound);
1111               } else {
1112                    assert (dualDj[iRow] < 1.0e-5);
1113                    rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
1114               }
1115          } else if (rowUpper_[iRow] > 1.0e20) {
1116               if (status == basic) {
1117                    rowActivity_[iRow] = rowLower_[iRow];
1118                    setRowStatus(iRow, atLowerBound);
1119               } else {
1120                    rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
1121                    assert (dualDj[iRow] > -1.0e-5);
1122               }
1123          } else {
1124               if (rowUpper_[iRow] == rowLower_[iRow]) {
1125                    rowActivity_[iRow] = rowLower_[iRow];
1126                    if (status == basic) {
1127                         setRowStatus(iRow, isFixed);
1128                    }
1129               } else {
1130                    // range
1131                    numberRanges++;
1132                    Status statusL = dualProblem->getColumnStatus(kExtraRow);
1133                    //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
1134                    //     iRow,status,kExtraRow,statusL, dualSol[iRow],
1135                    //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
1136                    if (status == basic) {
1137                         assert (statusL != basic);
1138                         rowActivity_[iRow] = rowUpper_[iRow];
1139                         setRowStatus(iRow, atUpperBound);
1140                    } else if (statusL == basic) {
1141                         numberBasic--; // already counted
1142                         rowActivity_[iRow] = rowLower_[iRow];
1143                         setRowStatus(iRow, atLowerBound);
1144                         dual_[iRow] = dualSol[kExtraRow];;
1145                    } else {
1146                         rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
1147                         assert (dualDj[iRow] < 1.0e-5);
1148                         // row basic
1149                         //setRowStatus(iRow,basic);
1150                         //numberBasic++;
1151                         dual_[iRow] = 0.0;
1152                    }
1153                    kExtraRow++;
1154               }
1155          }
1156     }
1157     if (numberBasic != numberRows_) {
1158          printf("Bad basis - ranges - coding needed\n");
1159          assert (numberRanges);
1160          abort();
1161     }
1162     if (optimizationDirection_ < 0.0) {
1163          for (iRow = 0; iRow < numberRows_; iRow++) {
1164               dual_[iRow] = -dual_[iRow];
1165          }
1166     }
1167     // redo row activities
1168     memset(rowActivity_, 0, numberRows_ * sizeof(double));
1169     matrix_->times(1.0, columnActivity_, rowActivity_);
1170     // redo reduced costs
1171     memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
1172     matrix_->transposeTimes(-1.0, dual_, reducedCost_);
1173     checkSolutionInternal();
1174     if (sumDualInfeasibilities_ > 1.0e-5 || sumPrimalInfeasibilities_ > 1.0e-5) {
1175          returnCode = 1;
1176#ifdef CLP_INVESTIGATE
1177          printf("There are %d dual infeasibilities summing to %g ",
1178                 numberDualInfeasibilities_, sumDualInfeasibilities_);
1179          printf("and %d primal infeasibilities summing to %g\n",
1180                 numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
1181#endif
1182     }
1183     // Below will go to ..DEBUG later
1184#if 1 //ndef NDEBUG
1185     // Check if correct
1186     double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
1187     double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
1188     double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
1189     double * dual = CoinCopyOfArray(dual_, numberRows_);
1190     this->dual(); //primal();
1191     CoinRelFltEq eq(1.0e-5);
1192     for (iRow = 0; iRow < numberRows_; iRow++) {
1193          assert(eq(dual[iRow], dual_[iRow]));
1194     }
1195     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1196          assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
1197     }
1198     for (iRow = 0; iRow < numberRows_; iRow++) {
1199          assert(eq(rowActivity[iRow], rowActivity_[iRow]));
1200     }
1201     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1202          assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
1203     }
1204     delete [] columnActivity;
1205     delete [] rowActivity;
1206     delete [] reducedCost;
1207     delete [] dual;
1208#endif
1209     return returnCode;
1210}
1211/* Does very cursory presolve.
1212   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
1213*/
1214ClpSimplex *
1215ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
1216                        int & nBound, bool moreBounds, bool tightenBounds)
1217{
1218     //#define CHECK_STATUS
1219#ifdef CHECK_STATUS
1220     {
1221          int n = 0;
1222          int i;
1223          for (i = 0; i < numberColumns_; i++)
1224               if (getColumnStatus(i) == ClpSimplex::basic)
1225                    n++;
1226          for (i = 0; i < numberRows_; i++)
1227               if (getRowStatus(i) == ClpSimplex::basic)
1228                    n++;
1229          assert (n == numberRows_);
1230     }
1231#endif
1232
1233     const double * element = matrix_->getElements();
1234     const int * row = matrix_->getIndices();
1235     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1236     const int * columnLength = matrix_->getVectorLengths();
1237
1238     CoinZeroN(rhs, numberRows_);
1239     int iColumn;
1240     int iRow;
1241     CoinZeroN(whichRow, numberRows_);
1242     int * backColumn = whichColumn + numberColumns_;
1243     int numberRows2 = 0;
1244     int numberColumns2 = 0;
1245     double offset = 0.0;
1246     const double * objective = this->objective();
1247     double * solution = columnActivity_;
1248     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1249          double lower = columnLower_[iColumn];
1250          double upper = columnUpper_[iColumn];
1251          if (upper > lower || getColumnStatus(iColumn) == ClpSimplex::basic) {
1252               backColumn[iColumn] = numberColumns2;
1253               whichColumn[numberColumns2++] = iColumn;
1254               for (CoinBigIndex j = columnStart[iColumn];
1255                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1256                    int iRow = row[j];
1257                    int n = whichRow[iRow];
1258                    if (n == 0 && element[j])
1259                         whichRow[iRow] = -iColumn - 1;
1260                    else if (n < 0)
1261                         whichRow[iRow] = 2;
1262               }
1263          } else {
1264               // fixed
1265               backColumn[iColumn] = -1;
1266               solution[iColumn] = upper;
1267               if (upper) {
1268                    offset += objective[iColumn] * upper;
1269                    for (CoinBigIndex j = columnStart[iColumn];
1270                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1271                         int iRow = row[j];
1272                         double value = element[j];
1273                         rhs[iRow] += upper * value;
1274                    }
1275               }
1276          }
1277     }
1278     int returnCode = 0;
1279     double tolerance = primalTolerance();
1280     nBound = 2 * numberRows_;
1281     for (iRow = 0; iRow < numberRows_; iRow++) {
1282          int n = whichRow[iRow];
1283          if (n > 0) {
1284               whichRow[numberRows2++] = iRow;
1285          } else if (n < 0) {
1286               //whichRow[numberRows2++]=iRow;
1287               //continue;
1288               // Can only do in certain circumstances as we don't know current value
1289               if (rowLower_[iRow] == rowUpper_[iRow] || getRowStatus(iRow) == ClpSimplex::basic) {
1290                    // save row and column for bound
1291                    whichRow[--nBound] = iRow;
1292                    whichRow[nBound+numberRows_] = -n - 1;
1293               } else if (moreBounds) {
1294                    // save row and column for bound
1295                    whichRow[--nBound] = iRow;
1296                    whichRow[nBound+numberRows_] = -n - 1;
1297               } else {
1298                    whichRow[numberRows2++] = iRow;
1299               }
1300          } else {
1301               // empty
1302               double rhsValue = rhs[iRow];
1303               if (rhsValue < rowLower_[iRow] - tolerance || rhsValue > rowUpper_[iRow] + tolerance) {
1304                    returnCode = 1; // infeasible
1305               }
1306          }
1307     }
1308     ClpSimplex * small = NULL;
1309     if (!returnCode) {
1310          small = new ClpSimplex(this, numberRows2, whichRow,
1311                                 numberColumns2, whichColumn, true, false);
1312#if 0
1313          ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
1314          if (rowCopy) {
1315               assert(!small->rowCopy());
1316               small->setNewRowCopy(new ClpPackedMatrix(*rowCopy, numberRows2, whichRow,
1317                                    numberColumns2, whichColumn));
1318          }
1319#endif
1320          // Set some stuff
1321          small->setDualBound(dualBound_);
1322          small->setInfeasibilityCost(infeasibilityCost_);
1323          small->setSpecialOptions(specialOptions_);
1324          small->setPerturbation(perturbation_);
1325          small->defaultFactorizationFrequency();
1326          small->setAlphaAccuracy(alphaAccuracy_);
1327          // If no rows left then no tightening!
1328          if (!numberRows2 || !numberColumns2)
1329               tightenBounds = false;
1330
1331          int numberElements = getNumElements();
1332          int numberElements2 = small->getNumElements();
1333          small->setObjectiveOffset(objectiveOffset() - offset);
1334          handler_->message(CLP_CRUNCH_STATS, messages_)
1335                    << numberRows2 << -(numberRows_ - numberRows2)
1336                    << numberColumns2 << -(numberColumns_ - numberColumns2)
1337                    << numberElements2 << -(numberElements - numberElements2)
1338                    << CoinMessageEol;
1339          // And set objective value to match
1340          small->setObjectiveValue(this->objectiveValue());
1341          double * rowLower2 = small->rowLower();
1342          double * rowUpper2 = small->rowUpper();
1343          int jRow;
1344          for (jRow = 0; jRow < numberRows2; jRow++) {
1345               iRow = whichRow[jRow];
1346               if (rowLower2[jRow] > -1.0e20)
1347                    rowLower2[jRow] -= rhs[iRow];
1348               if (rowUpper2[jRow] < 1.0e20)
1349                    rowUpper2[jRow] -= rhs[iRow];
1350          }
1351          // and bounds
1352          double * columnLower2 = small->columnLower();
1353          double * columnUpper2 = small->columnUpper();
1354          const char * integerInformation = integerType_;
1355          for (jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1356               iRow = whichRow[jRow];
1357               iColumn = whichRow[jRow+numberRows_];
1358               double lowerRow = rowLower_[iRow];
1359               if (lowerRow > -1.0e20)
1360                    lowerRow -= rhs[iRow];
1361               double upperRow = rowUpper_[iRow];
1362               if (upperRow < 1.0e20)
1363                    upperRow -= rhs[iRow];
1364               int jColumn = backColumn[iColumn];
1365               double lower = columnLower2[jColumn];
1366               double upper = columnUpper2[jColumn];
1367               double value = 0.0;
1368               for (CoinBigIndex j = columnStart[iColumn];
1369                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1370                    if (iRow == row[j]) {
1371                         value = element[j];
1372                         break;
1373                    }
1374               }
1375               assert (value);
1376               // convert rowLower and Upper to implied bounds on column
1377               double newLower = -COIN_DBL_MAX;
1378               double newUpper = COIN_DBL_MAX;
1379               if (value > 0.0) {
1380                    if (lowerRow > -1.0e20)
1381                         newLower = lowerRow / value;
1382                    if (upperRow < 1.0e20)
1383                         newUpper = upperRow / value;
1384               } else {
1385                    if (upperRow < 1.0e20)
1386                         newLower = upperRow / value;
1387                    if (lowerRow > -1.0e20)
1388                         newUpper = lowerRow / value;
1389               }
1390               if (integerInformation && integerInformation[iColumn]) {
1391                    if (newLower - floor(newLower) < 10.0 * tolerance)
1392                         newLower = floor(newLower);
1393                    else
1394                         newLower = ceil(newLower);
1395                    if (ceil(newUpper) - newUpper < 10.0 * tolerance)
1396                         newUpper = ceil(newUpper);
1397                    else
1398                         newUpper = floor(newUpper);
1399               }
1400               newLower = CoinMax(lower, newLower);
1401               newUpper = CoinMin(upper, newUpper);
1402               if (newLower > newUpper + tolerance) {
1403                    //printf("XXYY inf on bound\n");
1404                    returnCode = 1;
1405               }
1406               columnLower2[jColumn] = newLower;
1407               columnUpper2[jColumn] = CoinMax(newLower, newUpper);
1408               if (getRowStatus(iRow) != ClpSimplex::basic) {
1409                    if (getColumnStatus(iColumn) == ClpSimplex::basic) {
1410                         if (columnLower2[jColumn] == columnUpper2[jColumn]) {
1411                              // can only get here if will be fixed
1412                              small->setColumnStatus(jColumn, ClpSimplex::isFixed);
1413                         } else {
1414                              // solution is valid
1415                              if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) <
1416                                        fabs(columnActivity_[iColumn] - columnUpper2[jColumn]))
1417                                   small->setColumnStatus(jColumn, ClpSimplex::atLowerBound);
1418                              else
1419                                   small->setColumnStatus(jColumn, ClpSimplex::atUpperBound);
1420                         }
1421                    } else {
1422                         //printf("what now neither basic\n");
1423                    }
1424               }
1425          }
1426          if (returnCode) {
1427               delete small;
1428               small = NULL;
1429          } else if (tightenBounds && integerInformation) {
1430               // See if we can tighten any bounds
1431               // use rhs for upper and small duals for lower
1432               double * up = rhs;
1433               double * lo = small->dualRowSolution();
1434               const double * element = small->clpMatrix()->getElements();
1435               const int * row = small->clpMatrix()->getIndices();
1436               const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1437               //const int * columnLength = small->clpMatrix()->getVectorLengths();
1438               CoinZeroN(lo, numberRows2);
1439               CoinZeroN(up, numberRows2);
1440               for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1441                    double upper = columnUpper2[iColumn];
1442                    double lower = columnLower2[iColumn];
1443                    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1444                    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1445                         int iRow = row[j];
1446                         double value = element[j];
1447                         if (value > 0.0) {
1448                              if (upper < 1.0e20)
1449                                   up[iRow] += upper * value;
1450                              else
1451                                   up[iRow] = COIN_DBL_MAX;
1452                              if (lower > -1.0e20)
1453                                   lo[iRow] += lower * value;
1454                              else
1455                                   lo[iRow] = -COIN_DBL_MAX;
1456                         } else {
1457                              if (upper < 1.0e20)
1458                                   lo[iRow] += upper * value;
1459                              else
1460                                   lo[iRow] = -COIN_DBL_MAX;
1461                              if (lower > -1.0e20)
1462                                   up[iRow] += lower * value;
1463                              else
1464                                   up[iRow] = COIN_DBL_MAX;
1465                         }
1466                    }
1467               }
1468               double * rowLower2 = small->rowLower();
1469               double * rowUpper2 = small->rowUpper();
1470               bool feasible = true;
1471               // make safer
1472               for (int iRow = 0; iRow < numberRows2; iRow++) {
1473                    double lower = lo[iRow];
1474                    if (lower > rowUpper2[iRow] + tolerance) {
1475                         feasible = false;
1476                         break;
1477                    } else {
1478                         lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance;
1479                    }
1480                    double upper = up[iRow];
1481                    if (upper < rowLower2[iRow] - tolerance) {
1482                         feasible = false;
1483                         break;
1484                    } else {
1485                         up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
1486                    }
1487               }
1488               if (!feasible) {
1489                    delete small;
1490                    small = NULL;
1491               } else {
1492                    // and tighten
1493                    for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1494                         if (integerInformation[whichColumn[iColumn]]) {
1495                              double upper = columnUpper2[iColumn];
1496                              double lower = columnLower2[iColumn];
1497                              double newUpper = upper;
1498                              double newLower = lower;
1499                              double difference = upper - lower;
1500                              if (lower > -1000.0 && upper < 1000.0) {
1501                                   for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1502                                        int iRow = row[j];
1503                                        double value = element[j];
1504                                        if (value > 0.0) {
1505                                             double upWithOut = up[iRow] - value * difference;
1506                                             if (upWithOut < 0.0) {
1507                                                  newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1508                                             }
1509                                             double lowWithOut = lo[iRow] + value * difference;
1510                                             if (lowWithOut > 0.0) {
1511                                                  newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1512                                             }
1513                                        } else {
1514                                             double upWithOut = up[iRow] + value * difference;
1515                                             if (upWithOut < 0.0) {
1516                                                  newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1517                                             }
1518                                             double lowWithOut = lo[iRow] - value * difference;
1519                                             if (lowWithOut > 0.0) {
1520                                                  newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1521                                             }
1522                                        }
1523                                   }
1524                                   if (newLower > lower || newUpper < upper) {
1525                                        if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1526                                             newUpper = floor(newUpper);
1527                                        else
1528                                             newUpper = floor(newUpper + 0.5);
1529                                        if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1530                                             newLower = ceil(newLower);
1531                                        else
1532                                             newLower = ceil(newLower - 0.5);
1533                                        // change may be too small - check
1534                                        if (newLower > lower || newUpper < upper) {
1535                                             if (newUpper >= newLower) {
1536                                                  // Could also tighten in this
1537                                                  //printf("%d bounds %g %g tightened to %g %g\n",
1538                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1539                                                  //     newLower,newUpper);
1540#if 1
1541                                                  columnUpper2[iColumn] = newUpper;
1542                                                  columnLower2[iColumn] = newLower;
1543                                                  columnUpper_[whichColumn[iColumn]] = newUpper;
1544                                                  columnLower_[whichColumn[iColumn]] = newLower;
1545#endif
1546                                                  // and adjust bounds on rows
1547                                                  newUpper -= upper;
1548                                                  newLower -= lower;
1549                                                  for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1550                                                       int iRow = row[j];
1551                                                       double value = element[j];
1552                                                       if (value > 0.0) {
1553                                                            up[iRow] += newUpper * value;
1554                                                            lo[iRow] += newLower * value;
1555                                                       } else {
1556                                                            lo[iRow] += newUpper * value;
1557                                                            up[iRow] += newLower * value;
1558                                                       }
1559                                                  }
1560                                             } else {
1561                                                  // infeasible
1562                                                  //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1563                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1564                                                  //     newLower,newUpper);
1565#if 1
1566                                                  delete small;
1567                                                  small = NULL;
1568                                                  break;
1569#endif
1570                                             }
1571                                        }
1572                                   }
1573                              }
1574                         }
1575                    }
1576               }
1577          }
1578     }
1579#if 0
1580     if (small) {
1581          static int which = 0;
1582          which++;
1583          char xxxx[20];
1584          sprintf(xxxx, "bad%d.mps", which);
1585          small->writeMps(xxxx, 0, 1);
1586          sprintf(xxxx, "largebad%d.mps", which);
1587          writeMps(xxxx, 0, 1);
1588          printf("bad%d %x old size %d %d new %d %d\n", which, small,
1589                 numberRows_, numberColumns_, small->numberRows(), small->numberColumns());
1590#if 0
1591          for (int i = 0; i < numberColumns_; i++)
1592               printf("Bound %d %g %g\n", i, columnLower_[i], columnUpper_[i]);
1593          for (int i = 0; i < numberRows_; i++)
1594               printf("Row bound %d %g %g\n", i, rowLower_[i], rowUpper_[i]);
1595#endif
1596     }
1597#endif
1598#ifdef CHECK_STATUS
1599     {
1600          int n = 0;
1601          int i;
1602          for (i = 0; i < small->numberColumns(); i++)
1603               if (small->getColumnStatus(i) == ClpSimplex::basic)
1604                    n++;
1605          for (i = 0; i < small->numberRows(); i++)
1606               if (small->getRowStatus(i) == ClpSimplex::basic)
1607                    n++;
1608          assert (n == small->numberRows());
1609     }
1610#endif
1611     return small;
1612}
1613/* After very cursory presolve.
1614   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
1615*/
1616void
1617ClpSimplexOther::afterCrunch(const ClpSimplex & small,
1618                             const int * whichRow,
1619                             const int * whichColumn, int nBound)
1620{
1621#ifndef NDEBUG
1622     for (int i = 0; i < small.numberRows(); i++)
1623          assert (whichRow[i] >= 0 && whichRow[i] < numberRows_);
1624     for (int i = 0; i < small.numberColumns(); i++)
1625          assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_);
1626#endif
1627     getbackSolution(small, whichRow, whichColumn);
1628     // and deal with status for bounds
1629     const double * element = matrix_->getElements();
1630     const int * row = matrix_->getIndices();
1631     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1632     const int * columnLength = matrix_->getVectorLengths();
1633     double tolerance = primalTolerance();
1634     double djTolerance = dualTolerance();
1635     for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1636          int iRow = whichRow[jRow];
1637          int iColumn = whichRow[jRow+numberRows_];
1638          if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1639               double lower = columnLower_[iColumn];
1640               double upper = columnUpper_[iColumn];
1641               double value = columnActivity_[iColumn];
1642               double djValue = reducedCost_[iColumn];
1643               dual_[iRow] = 0.0;
1644               if (upper > lower) {
1645                    if (value < lower + tolerance && djValue > -djTolerance) {
1646                         setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1647                         setRowStatus(iRow, ClpSimplex::basic);
1648                    } else if (value > upper - tolerance && djValue < djTolerance) {
1649                         setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1650                         setRowStatus(iRow, ClpSimplex::basic);
1651                    } else {
1652                         // has to be basic
1653                         setColumnStatus(iColumn, ClpSimplex::basic);
1654                         reducedCost_[iColumn] = 0.0;
1655                         double value = 0.0;
1656                         for (CoinBigIndex j = columnStart[iColumn];
1657                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1658                              if (iRow == row[j]) {
1659                                   value = element[j];
1660                                   break;
1661                              }
1662                         }
1663                         dual_[iRow] = djValue / value;
1664                         if (rowUpper_[iRow] > rowLower_[iRow]) {
1665                              if (fabs(rowActivity_[iRow] - rowLower_[iRow]) <
1666                                        fabs(rowActivity_[iRow] - rowUpper_[iRow]))
1667                                   setRowStatus(iRow, ClpSimplex::atLowerBound);
1668                              else
1669                                   setRowStatus(iRow, ClpSimplex::atUpperBound);
1670                         } else {
1671                              setRowStatus(iRow, ClpSimplex::isFixed);
1672                         }
1673                    }
1674               } else {
1675                    // row can always be basic
1676                    setRowStatus(iRow, ClpSimplex::basic);
1677               }
1678          } else {
1679               // row can always be basic
1680               setRowStatus(iRow, ClpSimplex::basic);
1681          }
1682     }
1683     //#ifndef NDEBUG
1684#if 0
1685     if  (small.status() == 0) {
1686          int n = 0;
1687          int i;
1688          for (i = 0; i < numberColumns; i++)
1689               if (getColumnStatus(i) == ClpSimplex::basic)
1690                    n++;
1691          for (i = 0; i < numberRows; i++)
1692               if (getRowStatus(i) == ClpSimplex::basic)
1693                    n++;
1694          assert (n == numberRows);
1695     }
1696#endif
1697}
1698/* Tightens integer bounds - returns number tightened or -1 if infeasible
1699 */
1700int
1701ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1702{
1703     // See if we can tighten any bounds
1704     // use rhs for upper and small duals for lower
1705     double * up = rhsSpace;
1706     double * lo = dual_;
1707     const double * element = matrix_->getElements();
1708     const int * row = matrix_->getIndices();
1709     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1710     const int * columnLength = matrix_->getVectorLengths();
1711     CoinZeroN(lo, numberRows_);
1712     CoinZeroN(up, numberRows_);
1713     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1714          double upper = columnUpper_[iColumn];
1715          double lower = columnLower_[iColumn];
1716          //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1717          for (CoinBigIndex j = columnStart[iColumn];
1718                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1719               int iRow = row[j];
1720               double value = element[j];
1721               if (value > 0.0) {
1722                    if (upper < 1.0e20)
1723                         up[iRow] += upper * value;
1724                    else
1725                         up[iRow] = COIN_DBL_MAX;
1726                    if (lower > -1.0e20)
1727                         lo[iRow] += lower * value;
1728                    else
1729                         lo[iRow] = -COIN_DBL_MAX;
1730               } else {
1731                    if (upper < 1.0e20)
1732                         lo[iRow] += upper * value;
1733                    else
1734                         lo[iRow] = -COIN_DBL_MAX;
1735                    if (lower > -1.0e20)
1736                         up[iRow] += lower * value;
1737                    else
1738                         up[iRow] = COIN_DBL_MAX;
1739               }
1740          }
1741     }
1742     bool feasible = true;
1743     // make safer
1744     double tolerance = primalTolerance();
1745     for (int iRow = 0; iRow < numberRows_; iRow++) {
1746          double lower = lo[iRow];
1747          if (lower > rowUpper_[iRow] + tolerance) {
1748               feasible = false;
1749               break;
1750          } else {
1751               lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance;
1752          }
1753          double upper = up[iRow];
1754          if (upper < rowLower_[iRow] - tolerance) {
1755               feasible = false;
1756               break;
1757          } else {
1758               up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance;
1759          }
1760     }
1761     int numberTightened = 0;
1762     if (!feasible) {
1763          return -1;
1764     } else if (integerType_) {
1765          // and tighten
1766          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1767               if (integerType_[iColumn]) {
1768                    double upper = columnUpper_[iColumn];
1769                    double lower = columnLower_[iColumn];
1770                    double newUpper = upper;
1771                    double newLower = lower;
1772                    double difference = upper - lower;
1773                    if (lower > -1000.0 && upper < 1000.0) {
1774                         for (CoinBigIndex j = columnStart[iColumn];
1775                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1776                              int iRow = row[j];
1777                              double value = element[j];
1778                              if (value > 0.0) {
1779                                   double upWithOut = up[iRow] - value * difference;
1780                                   if (upWithOut < 0.0) {
1781                                        newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1782                                   }
1783                                   double lowWithOut = lo[iRow] + value * difference;
1784                                   if (lowWithOut > 0.0) {
1785                                        newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1786                                   }
1787                              } else {
1788                                   double upWithOut = up[iRow] + value * difference;
1789                                   if (upWithOut < 0.0) {
1790                                        newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1791                                   }
1792                                   double lowWithOut = lo[iRow] - value * difference;
1793                                   if (lowWithOut > 0.0) {
1794                                        newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1795                                   }
1796                              }
1797                         }
1798                         if (newLower > lower || newUpper < upper) {
1799                              if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1800                                   newUpper = floor(newUpper);
1801                              else
1802                                   newUpper = floor(newUpper + 0.5);
1803                              if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1804                                   newLower = ceil(newLower);
1805                              else
1806                                   newLower = ceil(newLower - 0.5);
1807                              // change may be too small - check
1808                              if (newLower > lower || newUpper < upper) {
1809                                   if (newUpper >= newLower) {
1810                                        numberTightened++;
1811                                        //printf("%d bounds %g %g tightened to %g %g\n",
1812                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1813                                        //     newLower,newUpper);
1814                                        columnUpper_[iColumn] = newUpper;
1815                                        columnLower_[iColumn] = newLower;
1816                                        // and adjust bounds on rows
1817                                        newUpper -= upper;
1818                                        newLower -= lower;
1819                                        for (CoinBigIndex j = columnStart[iColumn];
1820                                                  j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1821                                             int iRow = row[j];
1822                                             double value = element[j];
1823                                             if (value > 0.0) {
1824                                                  up[iRow] += newUpper * value;
1825                                                  lo[iRow] += newLower * value;
1826                                             } else {
1827                                                  lo[iRow] += newUpper * value;
1828                                                  up[iRow] += newLower * value;
1829                                             }
1830                                        }
1831                                   } else {
1832                                        // infeasible
1833                                        //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1834                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1835                                        //     newLower,newUpper);
1836                                        return -1;
1837                                   }
1838                              }
1839                         }
1840                    }
1841               }
1842          }
1843     }
1844     return numberTightened;
1845}
1846/* Parametrics
1847   This is an initial slow version.
1848   The code uses current bounds + theta * change (if change array not NULL)
1849   and similarly for objective.
1850   It starts at startingTheta and returns ending theta in endingTheta.
1851   If reportIncrement 0.0 it will report on any movement
1852   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1853   If it can not reach input endingTheta return code will be 1 for infeasible,
1854   2 for unbounded, if error on ranges -1,  otherwise 0.
1855   Normal report is just theta and objective but
1856   if event handler exists it may do more
1857   On exit endingTheta is maximum reached (can be used for next startingTheta)
1858*/
1859int
1860ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement,
1861                             const double * changeLowerBound, const double * changeUpperBound,
1862                             const double * changeLowerRhs, const double * changeUpperRhs,
1863                             const double * changeObjective)
1864{
1865     bool needToDoSomething = true;
1866     bool canTryQuick = (reportIncrement) ? true : false;
1867     // Save copy of model
1868     ClpSimplex copyModel = *this;
1869     int savePerturbation = perturbation_;
1870     perturbation_ = 102; // switch off
1871     while (needToDoSomething) {
1872          needToDoSomething = false;
1873          algorithm_ = -1;
1874
1875          // save data
1876          ClpDataSave data = saveData();
1877          int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
1878          int iRow, iColumn;
1879          double * chgUpper = NULL;
1880          double * chgLower = NULL;
1881          double * chgObjective = NULL;
1882
1883          // Dantzig (as will not be used) (out later)
1884          ClpDualRowPivot * savePivot = dualRowPivot_;
1885          dualRowPivot_ = new ClpDualRowDantzig();
1886
1887          if (!returnCode) {
1888               // Find theta when bounds will cross over and create arrays
1889               int numberTotal = numberRows_ + numberColumns_;
1890               chgLower = new double[numberTotal];
1891               memset(chgLower, 0, numberTotal * sizeof(double));
1892               chgUpper = new double[numberTotal];
1893               memset(chgUpper, 0, numberTotal * sizeof(double));
1894               chgObjective = new double[numberTotal];
1895               memset(chgObjective, 0, numberTotal * sizeof(double));
1896               assert (!rowScale_);
1897               double maxTheta = 1.0e50;
1898               if (changeLowerRhs || changeUpperRhs) {
1899                    for (iRow = 0; iRow < numberRows_; iRow++) {
1900                         double lower = rowLower_[iRow];
1901                         double upper = rowUpper_[iRow];
1902                         if (lower > upper) {
1903                              maxTheta = -1.0;
1904                              break;
1905                         }
1906                         double changeLower = (changeLowerRhs) ? changeLowerRhs[iRow] : 0.0;
1907                         double changeUpper = (changeUpperRhs) ? changeUpperRhs[iRow] : 0.0;
1908                         if (lower > -1.0e20 && upper < 1.0e20) {
1909                              if (lower + maxTheta * changeLower > upper + maxTheta * changeUpper) {
1910                                   maxTheta = (upper - lower) / (changeLower - changeUpper);
1911                              }
1912                         }
1913                         if (lower > -1.0e20) {
1914                              lower_[numberColumns_+iRow] += startingTheta * changeLower;
1915                              chgLower[numberColumns_+iRow] = changeLower;
1916                         }
1917                         if (upper < 1.0e20) {
1918                              upper_[numberColumns_+iRow] += startingTheta * changeUpper;
1919                              chgUpper[numberColumns_+iRow] = changeUpper;
1920                         }
1921                    }
1922               }
1923               if (maxTheta > 0.0) {
1924                    if (changeLowerBound || changeUpperBound) {
1925                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1926                              double lower = columnLower_[iColumn];
1927                              double upper = columnUpper_[iColumn];
1928                              if (lower > upper) {
1929                                   maxTheta = -1.0;
1930                                   break;
1931                              }
1932                              double changeLower = (changeLowerBound) ? changeLowerBound[iColumn] : 0.0;
1933                              double changeUpper = (changeUpperBound) ? changeUpperBound[iColumn] : 0.0;
1934                              if (lower > -1.0e20 && upper < 1.0e20) {
1935                                   if (lower + maxTheta * changeLower > upper + maxTheta * changeUpper) {
1936                                        maxTheta = (upper - lower) / (changeLower - changeUpper);
1937                                   }
1938                              }
1939                              if (lower > -1.0e20) {
1940                                   lower_[iColumn] += startingTheta * changeLower;
1941                                   chgLower[iColumn] = changeLower;
1942                              }
1943                              if (upper < 1.0e20) {
1944                                   upper_[iColumn] += startingTheta * changeUpper;
1945                                   chgUpper[iColumn] = changeUpper;
1946                              }
1947                         }
1948                    }
1949                    if (maxTheta == 1.0e50)
1950                         maxTheta = COIN_DBL_MAX;
1951               }
1952               if (maxTheta < 0.0) {
1953                    // bad ranges or initial
1954                    returnCode = -1;
1955               }
1956               endingTheta = CoinMin(endingTheta, maxTheta);
1957               if (endingTheta < startingTheta) {
1958                    // bad initial
1959                    returnCode = -2;
1960               }
1961          }
1962          double saveEndingTheta = endingTheta;
1963          if (!returnCode) {
1964               if (changeObjective) {
1965                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1966                         chgObjective[iColumn] = changeObjective[iColumn];
1967                         cost_[iColumn] += startingTheta * changeObjective[iColumn];
1968                    }
1969               }
1970               double * saveDuals = NULL;
1971               reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
1972               assert (!problemStatus_);
1973               // Now do parametrics
1974               printf("at starting theta of %g objective value is %g\n", startingTheta,
1975                      objectiveValue());
1976               while (!returnCode) {
1977                    assert (reportIncrement);
1978                    returnCode = parametricsLoop(startingTheta, endingTheta, reportIncrement,
1979                                                 chgLower, chgUpper, chgObjective, data,
1980                                                 canTryQuick);
1981                    if (!returnCode) {
1982                         //double change = endingTheta-startingTheta;
1983                         startingTheta = endingTheta;
1984                         endingTheta = saveEndingTheta;
1985                         //for (int i=0;i<numberTotal;i++) {
1986                         //lower_[i] += change*chgLower[i];
1987                         //upper_[i] += change*chgUpper[i];
1988                         //cost_[i] += change*chgObjective[i];
1989                         //}
1990                         printf("at theta of %g objective value is %g\n", startingTheta,
1991                                objectiveValue());
1992                         if (startingTheta >= endingTheta)
1993                              break;
1994                    } else if (returnCode == -1) {
1995                         // trouble - do external solve
1996                         needToDoSomething = true;
1997                    } else {
1998                         abort();
1999                    }
2000               }
2001          }
2002          reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
2003
2004          delete dualRowPivot_;
2005          dualRowPivot_ = savePivot;
2006          // Restore any saved stuff
2007          restoreData(data);
2008          if (needToDoSomething) {
2009               double saveStartingTheta = startingTheta; // known to be feasible
2010               int cleanedUp = 1;
2011               while (cleanedUp) {
2012                    // tweak
2013                    if (cleanedUp == 1) {
2014                         if (!reportIncrement)
2015                              startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
2016                         else
2017                              startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
2018                    } else {
2019                         // restoring to go slowly
2020                         startingTheta = saveStartingTheta;
2021                    }
2022                    // only works if not scaled
2023                    int i;
2024                    const double * obj1 = objective();
2025                    double * obj2 = copyModel.objective();
2026                    const double * lower1 = columnLower_;
2027                    double * lower2 = copyModel.columnLower();
2028                    const double * upper1 = columnUpper_;
2029                    double * upper2 = copyModel.columnUpper();
2030                    for (i = 0; i < numberColumns_; i++) {
2031                         obj2[i] = obj1[i] + startingTheta * chgObjective[i];
2032                         lower2[i] = lower1[i] + startingTheta * chgLower[i];
2033                         upper2[i] = upper1[i] + startingTheta * chgUpper[i];
2034                    }
2035                    lower1 = rowLower_;
2036                    lower2 = copyModel.rowLower();
2037                    upper1 = rowUpper_;
2038                    upper2 = copyModel.rowUpper();
2039                    for (i = 0; i < numberRows_; i++) {
2040                         lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
2041                         upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
2042                    }
2043                    copyModel.dual();
2044                    if (copyModel.problemStatus()) {
2045                         printf("Can not get to theta of %g\n", startingTheta);
2046                         canTryQuick = false; // do slowly to get exact amount
2047                         // back to last known good
2048                         if (cleanedUp == 1)
2049                              cleanedUp = 2;
2050                         else
2051                              abort();
2052                    } else {
2053                         // and move stuff back
2054                         int numberTotal = numberRows_ + numberColumns_;
2055                         CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
2056                         CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
2057                         CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
2058                         cleanedUp = 0;
2059                    }
2060               }
2061          }
2062          delete [] chgLower;
2063          delete [] chgUpper;
2064          delete [] chgObjective;
2065     }
2066     perturbation_ = savePerturbation;
2067     return problemStatus_;
2068}
2069int
2070ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta, double reportIncrement,
2071                                 const double * changeLower, const double * changeUpper,
2072                                 const double * changeObjective, ClpDataSave & data,
2073                                 bool canTryQuick)
2074{
2075     // stuff is already at starting
2076     // For this crude version just try and go to end
2077     double change = 0.0;
2078     if (reportIncrement && canTryQuick) {
2079          endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
2080          change = endingTheta - startingTheta;
2081     }
2082     int numberTotal = numberRows_ + numberColumns_;
2083     int i;
2084     for ( i = 0; i < numberTotal; i++) {
2085          lower_[i] += change * changeLower[i];
2086          upper_[i] += change * changeUpper[i];
2087          switch(getStatus(i)) {
2088
2089          case basic:
2090          case isFree:
2091          case superBasic:
2092               break;
2093          case isFixed:
2094          case atUpperBound:
2095               solution_[i] = upper_[i];
2096               break;
2097          case atLowerBound:
2098               solution_[i] = lower_[i];
2099               break;
2100          }
2101          cost_[i] += change * changeObjective[i];
2102     }
2103     problemStatus_ = -1;
2104
2105     // This says whether to restore things etc
2106     // startup will have factorized so can skip
2107     int factorType = 0;
2108     // Start check for cycles
2109     progress_.startCheck();
2110     // Say change made on first iteration
2111     changeMade_ = 1;
2112     /*
2113       Status of problem:
2114       0 - optimal
2115       1 - infeasible
2116       2 - unbounded
2117       -1 - iterating
2118       -2 - factorization wanted
2119       -3 - redo checking without factorization
2120       -4 - looks infeasible
2121     */
2122     while (problemStatus_ < 0) {
2123          int iRow, iColumn;
2124          // clear
2125          for (iRow = 0; iRow < 4; iRow++) {
2126               rowArray_[iRow]->clear();
2127          }
2128
2129          for (iColumn = 0; iColumn < 2; iColumn++) {
2130               columnArray_[iColumn]->clear();
2131          }
2132
2133          // give matrix (and model costs and bounds a chance to be
2134          // refreshed (normally null)
2135          matrix_->refresh(this);
2136          // may factorize, checks if problem finished
2137          statusOfProblemInParametrics(factorType, data);
2138          // Say good factorization
2139          factorType = 1;
2140          if (data.sparseThreshold_) {
2141               // use default at present
2142               factorization_->sparseThreshold(0);
2143               factorization_->goSparse();
2144          }
2145
2146          // exit if victory declared
2147          if (problemStatus_ >= 0)
2148               break;
2149
2150          // test for maximum iterations
2151          if (hitMaximumIterations()) {
2152               problemStatus_ = 3;
2153               break;
2154          }
2155          // Check event
2156          {
2157               int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2158               if (status >= 0) {
2159                    problemStatus_ = 5;
2160                    secondaryStatus_ = ClpEventHandler::endOfFactorization;
2161                    break;
2162               }
2163          }
2164          // Do iterations
2165          if (canTryQuick) {
2166               double * saveDuals = NULL;
2167               reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
2168          } else {
2169               whileIterating(startingTheta,  endingTheta, reportIncrement,
2170                              changeLower, changeUpper,
2171                              changeObjective);
2172          }
2173     }
2174     if (!problemStatus_) {
2175          theta_ = change + startingTheta;
2176          eventHandler_->event(ClpEventHandler::theta);
2177          return 0;
2178     } else if (problemStatus_ == 10) {
2179          return -1;
2180     } else {
2181          return problemStatus_;
2182     }
2183}
2184/* Checks if finished.  Updates status */
2185void
2186ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
2187{
2188     if (type == 2) {
2189          // trouble - go to recovery
2190          problemStatus_ = 10;
2191          return;
2192     }
2193     if (problemStatus_ > -3 || factorization_->pivots()) {
2194          // factorize
2195          // later on we will need to recover from singularities
2196          // also we could skip if first time
2197          if (type) {
2198               // is factorization okay?
2199               if (internalFactorize(1)) {
2200                    // trouble - go to recovery
2201                    problemStatus_ = 10;
2202                    return;
2203               }
2204          }
2205          if (problemStatus_ != -4 || factorization_->pivots() > 10)
2206               problemStatus_ = -3;
2207     }
2208     // at this stage status is -3 or -4 if looks infeasible
2209     // get primal and dual solutions
2210     gutsOfSolution(NULL, NULL);
2211     double realDualInfeasibilities = sumDualInfeasibilities_;
2212     // If bad accuracy treat as singular
2213     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
2214          // trouble - go to recovery
2215          problemStatus_ = 10;
2216          return;
2217     } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
2218          // Can reduce tolerance
2219          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
2220          factorization_->pivotTolerance(newTolerance);
2221     }
2222     // Check if looping
2223     int loop;
2224     if (type != 2)
2225          loop = progress_.looping();
2226     else
2227          loop = -1;
2228     if (loop >= 0) {
2229          problemStatus_ = loop; //exit if in loop
2230          if (!problemStatus_) {
2231               // declaring victory
2232               numberPrimalInfeasibilities_ = 0;
2233               sumPrimalInfeasibilities_ = 0.0;
2234          } else {
2235               problemStatus_ = 10; // instead - try other algorithm
2236          }
2237          return;
2238     } else if (loop < -1) {
2239          // something may have changed
2240          gutsOfSolution(NULL, NULL);
2241     }
2242     progressFlag_ = 0; //reset progress flag
2243     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
2244          handler_->message(CLP_SIMPLEX_STATUS, messages_)
2245                    << numberIterations_ << objectiveValue();
2246          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
2247                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
2248          handler_->printing(sumDualInfeasibilities_ > 0.0)
2249                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
2250          handler_->printing(numberDualInfeasibilitiesWithoutFree_
2251                             < numberDualInfeasibilities_)
2252                    << numberDualInfeasibilitiesWithoutFree_;
2253          handler_->message() << CoinMessageEol;
2254     }
2255     /* If we are primal feasible and any dual infeasibilities are on
2256        free variables then it is better to go to primal */
2257     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
2258               numberDualInfeasibilities_) {
2259          problemStatus_ = 10;
2260          return;
2261     }
2262
2263     // check optimal
2264     // give code benefit of doubt
2265     if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
2266               sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
2267          // say optimal (with these bounds etc)
2268          numberDualInfeasibilities_ = 0;
2269          sumDualInfeasibilities_ = 0.0;
2270          numberPrimalInfeasibilities_ = 0;
2271          sumPrimalInfeasibilities_ = 0.0;
2272     }
2273     if (dualFeasible() || problemStatus_ == -4) {
2274          progress_.modifyObjective(objectiveValue_
2275                                    - sumDualInfeasibilities_ * dualBound_);
2276     }
2277     if (numberPrimalInfeasibilities_) {
2278          if (problemStatus_ == -4 || problemStatus_ == -5) {
2279               problemStatus_ = 1; // infeasible
2280          }
2281     } else if (numberDualInfeasibilities_) {
2282          // clean up
2283          problemStatus_ = 10;
2284     } else {
2285          problemStatus_ = 0;
2286     }
2287     lastGoodIteration_ = numberIterations_;
2288     if (problemStatus_ < 0) {
2289          sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
2290          if (sumDualInfeasibilities_)
2291               numberDualInfeasibilities_ = 1;
2292     }
2293     // Allow matrices to be sorted etc
2294     int fake = -999; // signal sort
2295     matrix_->correctSequence(this, fake, fake);
2296}
2297/* This has the flow between re-factorizations
2298   Reasons to come out:
2299   -1 iterations etc
2300   -2 inaccuracy
2301   -3 slight inaccuracy (and done iterations)
2302   +0 looks optimal (might be unbounded - but we will investigate)
2303   +1 looks infeasible
2304   +3 max iterations
2305*/
2306int
2307ClpSimplexOther::whileIterating(double startingTheta, double & endingTheta, double /*reportIncrement*/,
2308                                const double * changeLower, const double * changeUpper,
2309                                const double * changeObjective)
2310{
2311     {
2312          int i;
2313          for (i = 0; i < 4; i++) {
2314               rowArray_[i]->clear();
2315          }
2316          for (i = 0; i < 2; i++) {
2317               columnArray_[i]->clear();
2318          }
2319     }
2320     // if can't trust much and long way from optimal then relax
2321     if (largestPrimalError_ > 10.0)
2322          factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
2323     else
2324          factorization_->relaxAccuracyCheck(1.0);
2325     // status stays at -1 while iterating, >=0 finished, -2 to invert
2326     // status -3 to go to top without an invert
2327     int returnCode = -1;
2328     double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
2329     double useTheta = startingTheta;
2330     double * primalChange = new double[numberRows_];
2331     double * dualChange = new double[numberColumns_];
2332     int numberTotal = numberColumns_ + numberRows_;
2333     int iSequence;
2334     // See if bounds
2335     int type = 0;
2336     for (iSequence = 0; iSequence < numberTotal; iSequence++) {
2337          if (changeLower[iSequence] || changeUpper[iSequence]) {
2338               type = 1;
2339               break;
2340          }
2341     }
2342     // See if objective
2343     for (iSequence = 0; iSequence < numberTotal; iSequence++) {
2344          if (changeObjective[iSequence]) {
2345               type |= 2;
2346               break;
2347          }
2348     }
2349     assert (type);
2350     while (problemStatus_ == -1) {
2351          double increaseTheta = CoinMin(endingTheta - useTheta, 1.0e50);
2352
2353          // Get theta for bounds - we know can't crossover
2354          int pivotType = nextTheta(type, increaseTheta, primalChange, dualChange,
2355                                    changeLower, changeUpper, changeObjective);
2356          if (pivotType)
2357               abort();
2358          // choose row to go out
2359          // dualRow will go to virtual row pivot choice algorithm
2360          reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
2361          if (pivotRow_ >= 0) {
2362               // we found a pivot row
2363               if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
2364                    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
2365                              << pivotRow_
2366                              << CoinMessageEol;
2367               }
2368               // check accuracy of weights
2369               dualRowPivot_->checkAccuracy();
2370               // Get good size for pivot
2371               // Allow first few iterations to take tiny
2372               double acceptablePivot = 1.0e-9;
2373               if (numberIterations_ > 100)
2374                    acceptablePivot = 1.0e-8;
2375               if (factorization_->pivots() > 10 ||
2376                         (factorization_->pivots() && saveSumDual))
2377                    acceptablePivot = 1.0e-5; // if we have iterated be more strict
2378               else if (factorization_->pivots() > 5)
2379                    acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
2380               else if (factorization_->pivots())
2381                    acceptablePivot = 1.0e-8; // relax
2382               double bestPossiblePivot = 1.0;
2383               // get sign for finding row of tableau
2384               // normal iteration
2385               // create as packed
2386               double direction = directionOut_;
2387               rowArray_[0]->createPacked(1, &pivotRow_, &direction);
2388               factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
2389               // put row of tableau in rowArray[0] and columnArray[0]
2390               matrix_->transposeTimes(this, -1.0,
2391                                       rowArray_[0], rowArray_[3], columnArray_[0]);
2392               // do ratio test for normal iteration
2393               bestPossiblePivot = reinterpret_cast<ClpSimplexDual *> ( this)->dualColumn(rowArray_[0],
2394                                   columnArray_[0], columnArray_[1],
2395                                   rowArray_[3], acceptablePivot, NULL);
2396               if (sequenceIn_ >= 0) {
2397                    // normal iteration
2398                    // update the incoming column
2399                    double btranAlpha = -alpha_ * directionOut_; // for check
2400                    unpackPacked(rowArray_[1]);
2401                    // moved into updateWeights factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2402                    // and update dual weights (can do in parallel - with extra array)
2403                    alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
2404                                                          rowArray_[2],
2405                                                          rowArray_[3],
2406                                                          rowArray_[1]);
2407                    // see if update stable
2408#ifdef CLP_DEBUG
2409                    if ((handler_->logLevel() & 32))
2410                         printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
2411#endif
2412                    double checkValue = 1.0e-7;
2413                    // if can't trust much and long way from optimal then relax
2414                    if (largestPrimalError_ > 10.0)
2415                         checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
2416                    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
2417                              fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
2418                         handler_->message(CLP_DUAL_CHECK, messages_)
2419                                   << btranAlpha
2420                                   << alpha_
2421                                   << CoinMessageEol;
2422                         if (factorization_->pivots()) {
2423                              dualRowPivot_->unrollWeights();
2424                              problemStatus_ = -2; // factorize now
2425                              rowArray_[0]->clear();
2426                              rowArray_[1]->clear();
2427                              columnArray_[0]->clear();
2428                              returnCode = -2;
2429                              break;
2430                         } else {
2431                              // take on more relaxed criterion
2432                              double test;
2433                              if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
2434                                   test = 1.0e-1 * fabs(alpha_);
2435                              else
2436                                   test = 1.0e-4 * (1.0 + fabs(alpha_));
2437                              if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
2438                                        fabs(btranAlpha - alpha_) > test) {
2439                                   dualRowPivot_->unrollWeights();
2440                                   // need to reject something
2441                                   char x = isColumn(sequenceOut_) ? 'C' : 'R';
2442                                   handler_->message(CLP_SIMPLEX_FLAG, messages_)
2443                                             << x << sequenceWithin(sequenceOut_)
2444                                             << CoinMessageEol;
2445                                   setFlagged(sequenceOut_);
2446                                   progress_.clearBadTimes();
2447                                   lastBadIteration_ = numberIterations_; // say be more cautious
2448                                   rowArray_[0]->clear();
2449                                   rowArray_[1]->clear();
2450                                   columnArray_[0]->clear();
2451                                   if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
2452                                        //printf("I think should declare infeasible\n");
2453                                        problemStatus_ = 1;
2454                                        returnCode = 1;
2455                                        break;
2456                                   }
2457                                   continue;
2458                              }
2459                         }
2460                    }
2461                    // update duals BEFORE replaceColumn so can do updateColumn
2462                    double objectiveChange = 0.0;
2463                    // do duals first as variables may flip bounds
2464                    // rowArray_[0] and columnArray_[0] may have flips
2465                    // so use rowArray_[3] for work array from here on
2466                    int nswapped = 0;
2467                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
2468                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
2469                    nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
2470                               rowArray_[2], theta_,
2471                               objectiveChange, false);
2472
2473                    // which will change basic solution
2474                    if (nswapped) {
2475                         factorization_->updateColumn(rowArray_[3], rowArray_[2]);
2476                         dualRowPivot_->updatePrimalSolution(rowArray_[2],
2477                                                             1.0, objectiveChange);
2478                         // recompute dualOut_
2479                         valueOut_ = solution_[sequenceOut_];
2480                         if (directionOut_ < 0) {
2481                              dualOut_ = valueOut_ - upperOut_;
2482                         } else {
2483                              dualOut_ = lowerOut_ - valueOut_;
2484                         }
2485                    }
2486                    // amount primal will move
2487                    double movement = -dualOut_ * directionOut_ / alpha_;
2488                    // so objective should increase by fabs(dj)*movement
2489                    // but we already have objective change - so check will be good
2490                    if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
2491#ifdef CLP_DEBUG
2492                         if (handler_->logLevel() & 32)
2493                              printf("movement %g, swap change %g, rest %g  * %g\n",
2494                                     objectiveChange + fabs(movement * dualIn_),
2495                                     objectiveChange, movement, dualIn_);
2496#endif
2497                         if(factorization_->pivots()) {
2498                              // going backwards - factorize
2499                              dualRowPivot_->unrollWeights();
2500                              problemStatus_ = -2; // factorize now
2501                              returnCode = -2;
2502                              break;
2503                         }
2504                    }
2505                    CoinAssert(fabs(dualOut_) < 1.0e50);
2506                    // if stable replace in basis
2507                    int updateStatus = factorization_->replaceColumn(this,
2508                                       rowArray_[2],
2509                                       rowArray_[1],
2510                                       pivotRow_,
2511                                       alpha_);
2512                    // if no pivots, bad update but reasonable alpha - take and invert
2513                    if (updateStatus == 2 &&
2514                              !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
2515                         updateStatus = 4;
2516                    if (updateStatus == 1 || updateStatus == 4) {
2517                         // slight error
2518                         if (factorization_->pivots() > 5 || updateStatus == 4) {
2519                              problemStatus_ = -2; // factorize now
2520                              returnCode = -3;
2521                         }
2522                    } else if (updateStatus == 2) {
2523                         // major error
2524                         dualRowPivot_->unrollWeights();
2525                         // later we may need to unwind more e.g. fake bounds
2526                         if (factorization_->pivots()) {
2527                              problemStatus_ = -2; // factorize now
2528                              returnCode = -2;
2529                              break;
2530                         } else {
2531                              // need to reject something
2532                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
2533                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
2534                                        << x << sequenceWithin(sequenceOut_)
2535                                        << CoinMessageEol;
2536                              setFlagged(sequenceOut_);
2537                              progress_.clearBadTimes();
2538                              lastBadIteration_ = numberIterations_; // say be more cautious
2539                              rowArray_[0]->clear();
2540                              rowArray_[1]->clear();
2541                              columnArray_[0]->clear();
2542                              // make sure dual feasible
2543                              // look at all rows and columns
2544                              double objectiveChange = 0.0;
2545                              reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
2546                                        0.0, objectiveChange, true);
2547                              continue;
2548                         }
2549                    } else if (updateStatus == 3) {
2550                         // out of memory
2551                         // increase space if not many iterations
2552                         if (factorization_->pivots() <
2553                                   0.5 * factorization_->maximumPivots() &&
2554                                   factorization_->pivots() < 200)
2555                              factorization_->areaFactor(
2556                                   factorization_->areaFactor() * 1.1);
2557                         problemStatus_ = -2; // factorize now
2558                    } else if (updateStatus == 5) {
2559                         problemStatus_ = -2; // factorize now
2560                    }
2561                    // update primal solution
2562                    if (theta_ < 0.0) {
2563#ifdef CLP_DEBUG
2564                         if (handler_->logLevel() & 32)
2565                              printf("negative theta %g\n", theta_);
2566#endif
2567                         theta_ = 0.0;
2568                    }
2569                    // do actual flips
2570                    reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
2571                    //rowArray_[1]->expand();
2572                    dualRowPivot_->updatePrimalSolution(rowArray_[1],
2573                                                        movement,
2574                                                        objectiveChange);
2575                    // modify dualout
2576                    dualOut_ /= alpha_;
2577                    dualOut_ *= -directionOut_;
2578                    //setStatus(sequenceIn_,basic);
2579                    dj_[sequenceIn_] = 0.0;
2580                    double oldValue = valueIn_;
2581                    if (directionIn_ == -1) {
2582                         // as if from upper bound
2583                         valueIn_ = upperIn_ + dualOut_;
2584                    } else {
2585                         // as if from lower bound
2586                         valueIn_ = lowerIn_ + dualOut_;
2587                    }
2588                    objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
2589                    // outgoing
2590                    // set dj to zero unless values pass
2591                    if (directionOut_ > 0) {
2592                         valueOut_ = lowerOut_;
2593                         dj_[sequenceOut_] = theta_;
2594                    } else {
2595                         valueOut_ = upperOut_;
2596                         dj_[sequenceOut_] = -theta_;
2597                    }
2598                    solution_[sequenceOut_] = valueOut_;
2599                    int whatNext = housekeeping(objectiveChange);
2600                    // and set bounds correctly
2601                    reinterpret_cast<ClpSimplexDual *> ( this)->originalBound(sequenceIn_);
2602                    reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
2603                    if (whatNext == 1) {
2604                         problemStatus_ = -2; // refactorize
2605                    } else if (whatNext == 2) {
2606                         // maximum iterations or equivalent
2607                         problemStatus_ = 3;
2608                         returnCode = 3;
2609                         break;
2610                    }
2611                    // Check event
2612                    {
2613                         int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2614                         if (status >= 0) {
2615                              problemStatus_ = 5;
2616                              secondaryStatus_ = ClpEventHandler::endOfIteration;
2617                              returnCode = 4;
2618                              break;
2619                         }
2620                    }
2621               } else {
2622                    // no incoming column is valid
2623                    pivotRow_ = -1;
2624#ifdef CLP_DEBUG
2625                    if (handler_->logLevel() & 32)
2626                         printf("** no column pivot\n");
2627#endif
2628                    if (factorization_->pivots() < 5) {
2629                         // If not in branch and bound etc save ray
2630                         if ((specialOptions_&(1024 | 4096)) == 0) {
2631                              // create ray anyway
2632                              delete [] ray_;
2633                              ray_ = new double [ numberRows_];
2634                              rowArray_[0]->expand(); // in case packed
2635                              ClpDisjointCopyN(rowArray_[0]->denseVector(), numberRows_, ray_);
2636                         }
2637                         // If we have just factorized and infeasibility reasonable say infeas
2638                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
2639                              if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3
2640                                        || (specialOptions_ & 64) == 0) {
2641                                   // say infeasible
2642                                   problemStatus_ = 1;
2643                                   // unless primal feasible!!!!
2644                                   //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
2645                                   //     numberDualInfeasibilities_,sumDualInfeasibilities_);
2646                                   if (numberDualInfeasibilities_)
2647                                        problemStatus_ = 10;
2648                                   rowArray_[0]->clear();
2649                                   columnArray_[0]->clear();
2650                                   returnCode = 1;
2651                                   break;
2652                              }
2653                         }
2654                         // If special option set - put off as long as possible
2655                         if ((specialOptions_ & 64) == 0) {
2656                              problemStatus_ = -4; //say looks infeasible
2657                         } else {
2658                              // flag
2659                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
2660                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
2661                                        << x << sequenceWithin(sequenceOut_)
2662                                        << CoinMessageEol;
2663                              setFlagged(sequenceOut_);
2664                              if (!factorization_->pivots()) {
2665                                   rowArray_[0]->clear();
2666                                   columnArray_[0]->clear();
2667                                   continue;
2668                              }
2669                         }
2670                    }
2671                    rowArray_[0]->clear();
2672                    columnArray_[0]->clear();
2673                    returnCode = 1;
2674                    break;
2675               }
2676          } else {
2677               // no pivot row
2678#ifdef CLP_DEBUG
2679               if (handler_->logLevel() & 32)
2680                    printf("** no row pivot\n");
2681#endif
2682               int numberPivots = factorization_->pivots();
2683               bool specialCase;
2684               int useNumberFake;
2685               returnCode = 0;
2686               if (numberPivots < 20 &&
2687                         (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
2688                         && dualBound_ > 1.0e8) {
2689                    specialCase = true;
2690                    // as dual bound high - should be okay
2691                    useNumberFake = 0;
2692               } else {
2693                    specialCase = false;
2694                    useNumberFake = numberFake_;
2695               }
2696               if (!numberPivots || specialCase) {
2697                    // may have crept through - so may be optimal
2698                    // check any flagged variables
2699                    int iRow;
2700                    for (iRow = 0; iRow < numberRows_; iRow++) {
2701                         int iPivot = pivotVariable_[iRow];
2702                         if (flagged(iPivot))
2703                              break;
2704                    }
2705                    if (iRow < numberRows_ && numberPivots) {
2706                         // try factorization
2707                         returnCode = -2;
2708                    }
2709
2710                    if (useNumberFake || numberDualInfeasibilities_) {
2711                         // may be dual infeasible
2712                         problemStatus_ = -5;
2713                    } else {
2714                         if (iRow < numberRows_) {
2715                              problemStatus_ = -5;
2716                         } else {
2717                              if (numberPivots) {
2718                                   // objective may be wrong
2719                                   objectiveValue_ = innerProduct(cost_,
2720                                                                  numberColumns_ + numberRows_,
2721                                                                  solution_);
2722                                   objectiveValue_ += objective_->nonlinearOffset();
2723                                   objectiveValue_ /= (objectiveScale_ * rhsScale_);
2724                                   if ((specialOptions_ & 16384) == 0) {
2725                                        // and dual_ may be wrong (i.e. for fixed or basic)
2726                                        CoinIndexedVector * arrayVector = rowArray_[1];
2727                                        arrayVector->clear();
2728                                        int iRow;
2729                                        double * array = arrayVector->denseVector();
2730                                        /* Use dual_ instead of array
2731                                           Even though dual_ is only numberRows_ long this is
2732                                           okay as gets permuted to longer rowArray_[2]
2733                                        */
2734                                        arrayVector->setDenseVector(dual_);
2735                                        int * index = arrayVector->getIndices();
2736                                        int number = 0;
2737                                        for (iRow = 0; iRow < numberRows_; iRow++) {
2738                                             int iPivot = pivotVariable_[iRow];
2739                                             double value = cost_[iPivot];
2740                                             dual_[iRow] = value;
2741                                             if (value) {
2742                                                  index[number++] = iRow;
2743                                             }
2744                                        }
2745                                        arrayVector->setNumElements(number);
2746                                        // Extended duals before "updateTranspose"
2747                                        matrix_->dualExpanded(this, arrayVector, NULL, 0);
2748                                        // Btran basic costs
2749                                        rowArray_[2]->clear();
2750                                        factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
2751                                        // and return vector
2752                                        arrayVector->setDenseVector(array);
2753                                   }
2754                              }
2755                              problemStatus_ = 0;
2756                              sumPrimalInfeasibilities_ = 0.0;
2757                              if ((specialOptions_&(1024 + 16384)) != 0) {
2758                                   CoinIndexedVector * arrayVector = rowArray_[1];
2759                                   arrayVector->clear();
2760                                   double * rhs = arrayVector->denseVector();
2761                                   times(1.0, solution_, rhs);
2762                                   bool bad2 = false;
2763                                   int i;
2764                                   for ( i = 0; i < numberRows_; i++) {
2765                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
2766                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
2767                                             bad2 = true;
2768                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
2769                                        }
2770                                        rhs[i] = 0.0;
2771                                   }
2772                                   for ( i = 0; i < numberColumns_; i++) {
2773                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
2774                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
2775                                             bad2 = true;
2776                                        }
2777                                   }
2778                                   if (bad2) {
2779                                        problemStatus_ = -3;
2780                                        returnCode = -2;
2781                                        // Force to re-factorize early next time
2782                                        int numberPivots = factorization_->pivots();
2783                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2784                                   }
2785                              }
2786                         }
2787                    }
2788               } else {
2789                    problemStatus_ = -3;
2790                    returnCode = -2;
2791                    // Force to re-factorize early next time
2792                    int numberPivots = factorization_->pivots();
2793                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2794               }
2795               break;
2796          }
2797     }
2798     delete [] primalChange;
2799     delete [] dualChange;
2800     return returnCode;
2801}
2802// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
2803int
2804ClpSimplexOther::nextTheta(int type, double maxTheta, double * primalChange, double * /*dualChange*/,
2805                           const double * changeLower, const double * changeUpper,
2806                           const double * /*changeObjective*/)
2807{
2808     int numberTotal = numberColumns_ + numberRows_;
2809     int iSequence;
2810     int iRow;
2811     theta_ = maxTheta;
2812     bool toLower = false;
2813     if ((type & 1) != 0) {
2814          // get change
2815          for (iSequence = 0; iSequence < numberTotal; iSequence++) {
2816               primalChange[iSequence] = 0.0;
2817               switch(getStatus(iSequence)) {
2818
2819               case basic:
2820               case isFree:
2821               case superBasic:
2822                    break;
2823               case isFixed:
2824               case atUpperBound:
2825                    primalChange[iSequence] = changeUpper[iSequence];
2826                    break;
2827               case atLowerBound:
2828                    primalChange[iSequence] = changeLower[iSequence];
2829                    break;
2830               }
2831          }
2832          // use array
2833          double * array = rowArray_[1]->denseVector();
2834          times(1.0, primalChange, array);
2835          int * index = rowArray_[1]->getIndices();
2836          int number = 0;
2837          for (iRow = 0; iRow < numberRows_; iRow++) {
2838               double value = array[iRow];
2839               if (value) {
2840                    array[iRow] = value;
2841                    index[number++] = iRow;
2842               }
2843          }
2844          // ftran it
2845          rowArray_[1]->setNumElements(number);
2846          factorization_->updateColumn(rowArray_[0], rowArray_[1]);
2847          number = rowArray_[1]->getNumElements();
2848          pivotRow_ = -1;
2849          for (iRow = 0; iRow < number; iRow++) {
2850               int iPivot = index[iRow];
2851               iSequence = pivotVariable_[iPivot];
2852               // solution value will be sol - theta*alpha
2853               // bounds will be bounds + change *theta
2854               double currentSolution = solution_[iSequence];
2855               double currentLower = lower_[iSequence];
2856               double currentUpper = upper_[iSequence];
2857               double alpha = array[iPivot];
2858               assert (currentSolution >= currentLower - primalTolerance_);
2859               assert (currentSolution <= currentUpper + primalTolerance_);
2860               double thetaCoefficient;
2861               double hitsLower = COIN_DBL_MAX;
2862               thetaCoefficient = changeLower[iSequence] + alpha;
2863               if (fabs(thetaCoefficient) > 1.0e-8)
2864                    hitsLower = (currentSolution - currentLower) / thetaCoefficient;
2865               if (hitsLower < 0.0) {
2866                    // does not hit - but should we check further
2867                    hitsLower = COIN_DBL_MAX;
2868               }
2869               double hitsUpper = COIN_DBL_MAX;
2870               thetaCoefficient = changeUpper[iSequence] + alpha;
2871               if (fabs(thetaCoefficient) > 1.0e-8)
2872                    hitsUpper = (currentSolution - currentUpper) / thetaCoefficient;
2873               if (hitsUpper < 0.0) {
2874                    // does not hit - but should we check further
2875                    hitsUpper = COIN_DBL_MAX;
2876               }
2877               if (CoinMin(hitsLower, hitsUpper) < theta_) {
2878                    theta_ = CoinMin(hitsLower, hitsUpper);
2879                    toLower = hitsLower < hitsUpper;
2880                    pivotRow_ = iPivot;
2881               }
2882          }
2883     }
2884     if ((type & 2) != 0) {
2885          abort();
2886     }
2887     if (pivotRow_ >= 0) {
2888          sequenceOut_ = pivotVariable_[pivotRow_];
2889          valueOut_ = solution_[sequenceOut_];
2890          lowerOut_ = lower_[sequenceOut_];
2891          upperOut_ = upper_[sequenceOut_];
2892          if (!toLower) {
2893               directionOut_ = -1;
2894               dualOut_ = valueOut_ - upperOut_;
2895          } else if (valueOut_ < lowerOut_) {
2896               directionOut_ = 1;
2897               dualOut_ = lowerOut_ - valueOut_;
2898          }
2899          return 0;
2900     } else {
2901          return -1;
2902     }
2903}
2904/* Expands out all possible combinations for a knapsack
2905   If buildObj NULL then just computes space needed - returns number elements
2906   On entry numberOutput is maximum allowed, on exit it is number needed or
2907   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
2908   least space to return values which reconstruct input.
2909   Rows returned will be original rows but no entries will be returned for
2910   any rows all of whose entries are in knapsack.  So up to user to allow for this.
2911   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
2912   in expanded knapsack.  Values in buildRow and buildElement;
2913*/
2914int
2915ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
2916                                double * buildObj, CoinBigIndex * buildStart,
2917                                int * buildRow, double * buildElement, int reConstruct) const
2918{
2919     int iRow;
2920     int iColumn;
2921     // Get column copy
2922     CoinPackedMatrix * columnCopy = matrix();
2923     // Get a row copy in standard format
2924     CoinPackedMatrix matrixByRow;
2925     matrixByRow.reverseOrderedCopyOf(*columnCopy);
2926     const double * elementByRow = matrixByRow.getElements();
2927     const int * column = matrixByRow.getIndices();
2928     const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2929     const int * rowLength = matrixByRow.getVectorLengths();
2930     CoinBigIndex j;
2931     int * whichColumn = new int [numberColumns_];
2932     int * whichRow = new int [numberRows_];
2933     int numJ = 0;
2934     // Get what other columns can compensate for
2935     double * lo = new double [numberRows_];
2936     double * high = new double [numberRows_];
2937     {
2938          // Use to get tight column bounds
2939          ClpSimplex tempModel(*this);
2940          tempModel.tightenPrimalBounds(0.0, 0, true);
2941          // Now another model without knapsacks
2942          int nCol = 0;
2943          for (iRow = 0; iRow < numberRows_; iRow++) {
2944               whichRow[iRow] = iRow;
2945          }
2946          for (iColumn = 0; iColumn < numberColumns_; iColumn++)
2947               whichColumn[iColumn] = -1;
2948          for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
2949               int iColumn = column[j];
2950               if (columnUpper_[iColumn] > columnLower_[iColumn]) {
2951                    whichColumn[iColumn] = 0;
2952               } else {
2953                    assert (!columnLower_[iColumn]); // fix later
2954               }
2955          }
2956          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2957               if (whichColumn[iColumn] < 0)
2958                    whichColumn[nCol++] = iColumn;
2959          }
2960          ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
2961          // Row copy
2962          CoinPackedMatrix matrixByRow;
2963          matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
2964          const double * elementByRow = matrixByRow.getElements();
2965          const int * column = matrixByRow.getIndices();
2966          const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2967          const int * rowLength = matrixByRow.getVectorLengths();
2968          const double * columnLower = tempModel2.getColLower();
2969          const double * columnUpper = tempModel2.getColUpper();
2970          for (iRow = 0; iRow < numberRows_; iRow++) {
2971               lo[iRow] = -COIN_DBL_MAX;
2972               high[iRow] = COIN_DBL_MAX;
2973               if (rowLower_[iRow] > -1.0e20 || rowUpper_[iRow] < 1.0e20) {
2974
2975                    // possible row
2976                    int infiniteUpper = 0;
2977                    int infiniteLower = 0;
2978                    double maximumUp = 0.0;
2979                    double maximumDown = 0.0;
2980                    CoinBigIndex rStart = rowStart[iRow];
2981                    CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
2982                    CoinBigIndex j;
2983                    // Compute possible lower and upper ranges
2984
2985                    for (j = rStart; j < rEnd; ++j) {
2986                         double value = elementByRow[j];
2987                         iColumn = column[j];
2988                         if (value > 0.0) {
2989                              if (columnUpper[iColumn] >= 1.0e20) {
2990                                   ++infiniteUpper;
2991                              } else {
2992                                   maximumUp += columnUpper[iColumn] * value;
2993                              }
2994                              if (columnLower[iColumn] <= -1.0e20) {
2995                                   ++infiniteLower;
2996                              } else {
2997                                   maximumDown += columnLower[iColumn] * value;
2998                              }
2999                         } else if (value < 0.0) {
3000                              if (columnUpper[iColumn] >= 1.0e20) {
3001                                   ++infiniteLower;
3002                              } else {
3003                                   maximumDown += columnUpper[iColumn] * value;
3004                              }
3005                              if (columnLower[iColumn] <= -1.0e20) {
3006                                   ++infiniteUpper;
3007                              } else {
3008                                   maximumUp += columnLower[iColumn] * value;
3009                              }
3010                         }
3011                    }
3012                    // Build in a margin of error
3013                    maximumUp += 1.0e-8 * fabs(maximumUp) + 1.0e-7;
3014                    maximumDown -= 1.0e-8 * fabs(maximumDown) + 1.0e-7;
3015                    // we want to save effective rhs
3016                    double up = (infiniteUpper) ? COIN_DBL_MAX : maximumUp;
3017                    double down = (infiniteLower) ? -COIN_DBL_MAX : maximumDown;
3018                    if (up == COIN_DBL_MAX || rowLower_[iRow] == -COIN_DBL_MAX) {
3019                         // However low we go it doesn't matter
3020                         lo[iRow] = -COIN_DBL_MAX;
3021                    } else {
3022                         // If we go below this then can not be feasible
3023                         lo[iRow] = rowLower_[iRow] - up;
3024                    }
3025                    if (down == -COIN_DBL_MAX || rowUpper_[iRow] == COIN_DBL_MAX) {
3026                         // However high we go it doesn't matter
3027                         high[iRow] = COIN_DBL_MAX;
3028                    } else {
3029                         // If we go above this then can not be feasible
3030                         high[iRow] = rowUpper_[iRow] - down;
3031                    }
3032               }
3033          }
3034     }
3035     numJ = 0;
3036     for (iColumn = 0; iColumn < numberColumns_; iColumn++)
3037          whichColumn[iColumn] = -1;
3038     int * markRow = new int [numberRows_];
3039     for (iRow = 0; iRow < numberRows_; iRow++)
3040          markRow[iRow] = 1;
3041     for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
3042          int iColumn = column[j];
3043          if (columnUpper_[iColumn] > columnLower_[iColumn]) {
3044               whichColumn[iColumn] = numJ;
3045               numJ++;
3046          }
3047     }
3048     /* mark rows
3049        -n in knapsack and n other variables
3050        1 no entries
3051        n+1000 not involved in knapsack but n entries
3052        0 only in knapsack
3053     */
3054     for (iRow = 0; iRow < numberRows_; iRow++) {
3055          int type = 1;
3056          for (j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
3057               int iColumn = column[j];
3058               if (whichColumn[iColumn] >= 0) {
3059                    if (type == 1) {
3060                         type = 0;
3061                    } else if (type > 0) {
3062                         assert (type > 1000);
3063                         type = -(type - 1000);
3064                    }
3065               } else if (type == 1) {
3066                    type = 1001;
3067               } else if (type < 0) {
3068                    type --;
3069               } else if (type == 0) {
3070                    type = -1;
3071               } else {
3072                    assert (type > 1000);
3073                    type++;
3074               }
3075          }
3076          markRow[iRow] = type;
3077          if (type < 0 && type > -30 && false)
3078               printf("markrow on row %d is %d\n", iRow, markRow[iRow]);
3079     }
3080     int * bound = new int [numberColumns_+1];
3081     int * stack = new int [numberColumns_+1];
3082     int * flip = new int [numberColumns_+1];
3083     double * offset = new double[numberColumns_+1];
3084     double * size = new double [numberColumns_+1];
3085     double * rhsOffset = new double[numberRows_];
3086     int * build = new int[numberColumns_];
3087     int maxNumber = numberOutput;
3088     numJ = 0;
3089     double minSize = rowLower_[knapsackRow];
3090     double maxSize = rowUpper_[knapsackRow];
3091     double knapsackOffset = 0.0;
3092     for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
3093          int iColumn = column[j];
3094          double lowerColumn = columnLower_[iColumn];
3095          double upperColumn = columnUpper_[iColumn];
3096          if (lowerColumn == upperColumn)
3097               continue;
3098          double gap = upperColumn - lowerColumn;
3099          if (gap > 1.0e8)
3100               gap = 1.0e8;
3101          assert (fabs(floor(gap + 0.5) - gap) < 1.0e-5);
3102          whichColumn[numJ] = iColumn;
3103          bound[numJ] = static_cast<int> (gap);
3104          if (elementByRow[j] > 0.0) {
3105               flip[numJ] = 1;
3106               offset[numJ] = lowerColumn;
3107               size[numJ++] = elementByRow[j];
3108          } else {
3109               flip[numJ] = -1;
3110               offset[numJ] = upperColumn;
3111               size[numJ++] = -elementByRow[j];
3112               lowerColumn = upperColumn;
3113          }
3114          knapsackOffset += elementByRow[j] * lowerColumn;
3115     }
3116     int jRow;
3117     for (iRow = 0; iRow < numberRows_; iRow++)
3118          whichRow[iRow] = iRow;
3119     ClpSimplex smallModel(this, numberRows_, whichRow, numJ, whichColumn, true, true, true);
3120     // modify rhs to allow for nonzero lower bounds
3121     //double * rowLower = smallModel.rowLower();
3122     //double * rowUpper = smallModel.rowUpper();
3123     //const double * columnLower = smallModel.columnLower();
3124     //const double * columnUpper = smallModel.columnUpper();
3125     const CoinPackedMatrix * matrix = smallModel.matrix();
3126     const double * element = matrix->getElements();
3127     const int * row = matrix->getIndices();
3128     const CoinBigIndex * columnStart = matrix->getVectorStarts();
3129     const int * columnLength = matrix->getVectorLengths();
3130     const double * objective = smallModel.objective();
3131     //double objectiveOffset=0.0;
3132     // would use for fixed?
3133     CoinZeroN(rhsOffset, numberRows_);
3134     double * rowActivity = smallModel.primalRowSolution();
3135     CoinZeroN(rowActivity, numberRows_);
3136     maxSize -= knapsackOffset;
3137     minSize -= knapsackOffset;
3138     // now generate
3139     int i;
3140     int iStack = numJ;
3141     for (i = 0; i < numJ; i++) {
3142          stack[i] = 0;
3143     }
3144     double tooMuch = 10.0 * maxSize + 10000;
3145     stack[numJ] = 1;
3146     size[numJ] = tooMuch;
3147     bound[numJ] = 0;
3148     double sum = tooMuch;
3149     // allow for all zero being OK
3150     stack[numJ-1] = -1;
3151     sum -= size[numJ-1];
3152     numberOutput = 0;
3153     int nelCreate = 0;
3154     /* typeRun is - 0 for initial sizes
3155                     1 for build
3156          2 for reconstruct
3157     */
3158     int typeRun = buildObj ? 1 : 0;
3159     if (reConstruct >= 0) {
3160          assert (buildRow && buildElement);
3161          typeRun = 2;
3162     }
3163     if (typeRun == 1)
3164          buildStart[0] = 0;
3165     while (iStack >= 0) {
3166          if (sum >= minSize && sum <= maxSize) {
3167               double checkSize = 0.0;
3168               bool good = true;
3169               int nRow = 0;
3170               double obj = 0.0;
3171               CoinZeroN(rowActivity, numberRows_);
3172               for (iColumn = 0; iColumn < numJ; iColumn++) {
3173                    int iValue = stack[iColumn];
3174                    if (iValue > bound[iColumn]) {
3175                         good = false;
3176                         break;
3177                    } else {
3178                         double realValue = offset[iColumn] + flip[iColumn] * iValue;
3179                         if (realValue) {
3180                              obj += objective[iColumn] * realValue;
3181                              for (CoinBigIndex j = columnStart[iColumn];
3182                                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3183                                   double value = element[j] * realValue;
3184                                   int kRow = row[j];
3185                                   if (rowActivity[kRow]) {
3186                                        rowActivity[kRow] += value;
3187                                        if (!rowActivity[kRow])
3188                                             rowActivity[kRow] = 1.0e-100;
3189                                   } else {
3190                                        build[nRow++] = kRow;
3191                                        rowActivity[kRow] = value;
3192                                   }
3193                              }
3194                         }
3195                    }
3196               }
3197               if (good) {
3198                    for (jRow = 0; jRow < nRow; jRow++) {
3199                         int kRow = build[jRow];
3200                         double value = rowActivity[kRow];
3201                         if (value > high[kRow] || value < lo[kRow]) {
3202                              good = false;
3203                              break;
3204                         }
3205                    }
3206               }
3207               if (good) {
3208                    if (typeRun == 1) {
3209                         buildObj[numberOutput] = obj;
3210                         for (jRow = 0; jRow < nRow; jRow++) {
3211                              int kRow = build[jRow];
3212                              double value = rowActivity[kRow];
3213                              if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
3214                                   buildElement[nelCreate] = value;
3215                                   buildRow[nelCreate++] = kRow;
3216                              }
3217                         }
3218                         buildStart[numberOutput+1] = nelCreate;
3219                    } else if (!typeRun) {
3220                         for (jRow = 0; jRow < nRow; jRow++) {
3221                              int kRow = build[jRow];
3222                              double value = rowActivity[kRow];
3223                              if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
3224                                   nelCreate++;
3225                              }
3226                         }
3227                    }
3228                    if (typeRun == 2 && reConstruct == numberOutput) {
3229                         // build and exit
3230                         nelCreate = 0;
3231                         for (iColumn = 0; iColumn < numJ; iColumn++) {
3232                              int iValue = stack[iColumn];
3233                              double realValue = offset[iColumn] + flip[iColumn] * iValue;
3234                              if (realValue) {
3235                                   buildRow[nelCreate] = whichColumn[iColumn];
3236                                   buildElement[nelCreate++] = realValue;
3237                              }
3238                         }
3239                         numberOutput = 1;
3240                         for (i = 0; i < numJ; i++) {
3241                              bound[i] = 0;
3242                         }
3243                         break;
3244                    }
3245                    numberOutput++;
3246                    if (numberOutput > maxNumber) {
3247                         nelCreate = -numberOutput;
3248                         numberOutput = -1;
3249                         for (i = 0; i < numJ; i++) {
3250                              bound[i] = 0;
3251                         }
3252                         break;
3253                    } else if (typeRun == 1 && numberOutput == maxNumber) {
3254                         // On second run
3255                         for (i = 0; i < numJ; i++) {
3256                              bound[i] = 0;
3257                         }
3258                         break;
3259                    }
3260                    for (int j = 0; j < numJ; j++) {
3261                         checkSize += stack[j] * size[j];
3262                    }
3263                    assert (fabs(sum - checkSize) < 1.0e-3);
3264               }
3265               for (jRow = 0; jRow < nRow; jRow++) {
3266                    int kRow = build[jRow];
3267                    rowActivity[kRow] = 0.0;
3268               }
3269          }
3270          if (sum > maxSize || stack[iStack] > bound[iStack]) {
3271               sum -= size[iStack] * stack[iStack];
3272               stack[iStack--] = 0;
3273               if (iStack >= 0) {
3274                    stack[iStack] ++;
3275                    sum += size[iStack];
3276               }
3277          } else {
3278               // must be less
3279               // add to last possible
3280               iStack = numJ - 1;
3281               sum += size[iStack];
3282               stack[iStack]++;
3283          }
3284     }
3285     //printf("%d will be created\n",numberOutput);
3286     delete [] whichColumn;
3287     delete [] whichRow;
3288     delete [] bound;
3289     delete [] stack;
3290     delete [] flip;
3291     delete [] size;
3292     delete [] offset;
3293     delete [] rhsOffset;
3294     delete [] build;
3295     delete [] markRow;
3296     delete [] lo;
3297     delete [] high;
3298     return nelCreate;
3299}
3300// Quick try at cleaning up duals if postsolve gets wrong
3301void
3302ClpSimplexOther::cleanupAfterPostsolve()
3303{
3304     // First mark singleton equality rows
3305     char * mark = new char [ numberRows_];
3306     memset(mark, 0, numberRows_);
3307     const int * row = matrix_->getIndices();
3308     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3309     const int * columnLength = matrix_->getVectorLengths();
3310     const double * element = matrix_->getElements();
3311     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3312          for (CoinBigIndex j = columnStart[iColumn];
3313                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3314               int iRow = row[j];
3315               if (mark[iRow])
3316                    mark[iRow] = 2;
3317               else
3318                    mark[iRow] = 1;
3319          }
3320     }
3321     // for now just == rows
3322     for (int iRow = 0; iRow < numberRows_; iRow++) {
3323          if (rowUpper_[iRow] > rowLower_[iRow])
3324               mark[iRow] = 3;
3325     }
3326     double dualTolerance = dblParam_[ClpDualTolerance];
3327     double primalTolerance = dblParam_[ClpPrimalTolerance];
3328     int numberCleaned = 0;
3329     double maxmin = optimizationDirection_;
3330     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3331          double dualValue = reducedCost_[iColumn] * maxmin;
3332          double primalValue = columnActivity_[iColumn];
3333          double lower = columnLower_[iColumn];
3334          double upper = columnUpper_[iColumn];
3335          int way = 0;
3336          switch(getColumnStatus(iColumn)) {
3337
3338          case basic:
3339               // dual should be zero
3340               if (dualValue > dualTolerance) {
3341                    way = -1;
3342               } else if (dualValue < -dualTolerance) {
3343                    way = 1;
3344               }
3345               break;
3346          case ClpSimplex::isFixed:
3347               break;
3348          case atUpperBound:
3349               // dual should not be positive
3350               if (dualValue > dualTolerance) {
3351                    way = -1;
3352               }
3353               break;
3354          case atLowerBound:
3355               // dual should not be negative
3356               if (dualValue < -dualTolerance) {
3357                    way = 1;
3358               }
3359               break;
3360          case superBasic:
3361          case isFree:
3362               if (primalValue < upper - primalTolerance) {
3363                    // dual should not be negative
3364                    if (dualValue < -dualTolerance) {
3365                         way = 1;
3366                    }
3367               }
3368               if (primalValue > lower + primalTolerance) {
3369                    // dual should not be positive
3370                    if (dualValue > dualTolerance) {
3371                         way = -1;
3372                    }
3373               }
3374               break;
3375          }
3376          if (way) {
3377               // see if can find singleton row
3378               for (CoinBigIndex j = columnStart[iColumn];
3379                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3380                    int iRow = row[j];
3381                    if (mark[iRow] == 1) {
3382                         double value = element[j];
3383                         // dj - addDual*value == 0.0
3384                         double addDual = dualValue / value;
3385                         dual_[iRow] += addDual;
3386                         reducedCost_[iColumn] = 0.0;
3387                         numberCleaned++;
3388                         break;
3389                    }
3390               }
3391          }
3392     }
3393     delete [] mark;
3394#ifdef CLP_INVESTIGATE
3395     printf("cleanupAfterPostsolve cleaned up %d columns\n", numberCleaned);
3396#endif
3397     // Redo
3398     memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
3399     matrix_->transposeTimes(-1.0, dual_, reducedCost_);
3400     checkSolutionInternal();
3401}
Note: See TracBrowser for help on using the repository browser.