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

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

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