Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpSimplexOther.cpp

    r2384 r2385  
    4747
    4848*/
    49 void ClpSimplexOther::dualRanging(int numberCheck, const int * which,
    50                                   double * costIncreased, int * sequenceIncreased,
    51                                   double * costDecreased, int * sequenceDecreased,
    52                                   double * valueIncrease, double * valueDecrease)
     49void ClpSimplexOther::dualRanging(int numberCheck, const int *which,
     50  double *costIncreased, int *sequenceIncreased,
     51  double *costDecreased, int *sequenceDecreased,
     52  double *valueIncrease, double *valueDecrease)
    5353{
    54      rowArray_[1]->clear();
     54  rowArray_[1]->clear();
    5555#ifdef LONG_REGION_2
    56      rowArray_[2]->clear();
     56  rowArray_[2]->clear();
    5757#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           }
     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    }
    100100
    101           switch(getStatus(iSequence)) {
     101    switch (getStatus(iSequence)) {
    102102
    103           case basic: {
    104                // non-trvial
    105                // Get pivot row
    106                int iRow = backPivot[iSequence];
    107                assert (iRow >= 0);
     103    case basic: {
     104      // non-trvial
     105      // Get pivot row
     106      int iRow = backPivot[iSequence];
     107      assert(iRow >= 0);
    108108#ifndef COIN_FAC_NEW
    109                double plusOne = 1.0;
    110                rowArray_[0]->createPacked(1, &iRow, &plusOne);
     109      double plusOne = 1.0;
     110      rowArray_[0]->createPacked(1, &iRow, &plusOne);
    111111#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]);
     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]);
    124124#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();
     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();
    141141#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++) {
     142        const int *index = rowArray_[0]->getIndices();
     143#endif
     144        double scale2 = 0.0;
     145        int j;
     146        for (j = 0; j < number; j++) {
    147147#ifndef COIN_FAC_NEW
    148                          scale2 += arrayX[j] * arrayX[j];
     148          scale2 += arrayX[j] * arrayX[j];
    149149#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.
     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.
    162162                                 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 :
    166                                                      rhsScale_ * inverseColumnScale_[sequenceIncrease];
    167                                    costIncrease = CoinMax(fabs(djValue * movement), costIncrease);
    168                               }
    169                          } else {
    170                               costIncrease = 0.0;
    171                          }
    172                     }
    173                     if (sequenceDecrease >= 0) {
    174                          double djValue = dj_[sequenceDecrease];
    175                          if (fabs(djValue) > 10.0 * dualTolerance_) {
    176                               // we are going to use for cutoff so be exact
    177                               costDecrease = fabs(djValue / alphaDecrease);
    178                               if(sequenceDecrease < numberColumns_ && integerType_[sequenceDecrease]) {
    179                                    // can improve
    180                                    double movement = (columnScale_ == NULL) ? 1.0 :
    181                                                      rhsScale_ * inverseColumnScale_[sequenceDecrease];
    182                                    costDecrease = CoinMax(fabs(djValue * movement), costDecrease);
    183                               }
    184                          } else {
    185                               costDecrease = 0.0;
    186                          }
    187                     }
    188                     costIncrease *= scale2;
    189                     costDecrease *= scale2;
    190                }
    191           }
    192           break;
    193           case isFixed:
    194                break;
    195           case isFree:
    196           case superBasic:
    197                costIncrease = 0.0;
    198                costDecrease = 0.0;
    199                sequenceIncrease = iSequence;
    200                sequenceDecrease = iSequence;
    201                break;
    202           case atUpperBound:
    203                costIncrease = CoinMax(0.0, -dj_[iSequence]);
    204                sequenceIncrease = iSequence;
    205                if (valueIncrease)
    206                     valueIncrease[i] = primalRanging1(iSequence, iSequence);
    207                break;
    208           case atLowerBound:
    209                costDecrease = CoinMax(0.0, dj_[iSequence]);
    210                sequenceDecrease = iSequence;
    211                if (valueIncrease)
    212                     valueDecrease[i] = primalRanging1(iSequence, iSequence);
    213                break;
    214           }
    215           double scaleFactor;
    216           if (rowScale_) {
    217                if (iSequence < numberColumns_)
    218                     scaleFactor = 1.0 / (objectiveScale_ * columnScale_[iSequence]);
    219                else
    220                     scaleFactor = rowScale_[iSequence-numberColumns_] / objectiveScale_;
     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            }
    221168          } else {
    222                scaleFactor = 1.0 / objectiveScale_;
    223           }
    224           if (costIncrease < 1.0e30)
    225                costIncrease *= scaleFactor;
    226           if (costDecrease < 1.0e30)
    227                costDecrease *= scaleFactor;
    228           if (optimizationDirection_ == 1.0) {
    229                costIncreased[i] = costIncrease;
    230                sequenceIncreased[i] = sequenceIncrease;
    231                costDecreased[i] = costDecrease;
    232                sequenceDecreased[i] = sequenceDecrease;
    233           } else if (optimizationDirection_ == -1.0) {
    234                costIncreased[i] = costDecrease;
    235                sequenceIncreased[i] = sequenceDecrease;
    236                costDecreased[i] = costIncrease;
    237                sequenceDecreased[i] = sequenceIncrease;
    238                if (valueIncrease) {
    239                     double temp = valueIncrease[i];
    240                     valueIncrease[i] = valueDecrease[i];
    241                     valueDecrease[i] = temp;
    242                }
    243           } else if (optimizationDirection_ == 0.0) {
    244                // !!!!!! ???
    245                costIncreased[i] = COIN_DBL_MAX;
    246                sequenceIncreased[i] = -1;
    247                costDecreased[i] = COIN_DBL_MAX;
    248                sequenceDecreased[i] = -1;
     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            }
    249182          } else {
    250                abort();
    251           }
    252      }
    253      rowArray_[0]->clear();
    254      //rowArray_[1]->clear();
    255      //columnArray_[1]->clear();
    256      columnArray_[0]->clear();
    257      delete [] backPivot;
    258      if (!optimizationDirection_)
    259           printf("*** ????? Ranging with zero optimization costs\n");
     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");
    260257}
    261258/*
     
    264261   This is used in dual ranging
    265262*/
    266 void
    267 ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
    268                                  CoinIndexedVector * columnArray,
    269                                  double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
    270                                  double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
     263void ClpSimplexOther::checkDualRatios(CoinIndexedVector *rowArray,
     264  CoinIndexedVector *columnArray,
     265  double &costIncrease, int &sequenceIncrease, double &alphaIncrease,
     266  double &costDecrease, int &sequenceDecrease, double &alphaDecrease)
    271267{
    272      double acceptablePivot = 1.0e-9;
    273      double * work;
    274      int number;
    275      int * which;
    276      int iSection;
     268  double acceptablePivot = 1.0e-9;
     269  double *work;
     270  int number;
     271  int *which;
     272  int iSection;
    277273
    278      double thetaDown = 1.0e31;
    279      double thetaUp = 1.0e31;
    280      int sequenceDown = -1;
    281      int sequenceUp = -1;
    282      double alphaDown = 0.0;
    283      double alphaUp = 0.0;
     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;
    284280
    285      int addSequence;
     281  int addSequence;
    286282
    287      for (iSection = 0; iSection < 2; iSection++) {
     283  for (iSection = 0; iSection < 2; iSection++) {
    288284
    289           int i;
    290           if (!iSection) {
    291                work = rowArray->denseVector();
    292                number = rowArray->getNumElements();
    293                which = rowArray->getIndices();
    294                addSequence = numberColumns_;
    295           } else {
    296                work = columnArray->denseVector();
    297                number = columnArray->getNumElements();
    298                which = columnArray->getIndices();
    299                addSequence = 0;
    300           }
     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    }
    301297
    302           for (i = 0; i < number; i++) {
    303                int iSequence = which[i];
    304                int iSequence2 = iSequence + addSequence;
     298    for (i = 0; i < number; i++) {
     299      int iSequence = which[i];
     300      int iSequence2 = iSequence + addSequence;
    305301#ifndef COIN_FAC_NEW
    306                double alpha = work[i];
     302      double alpha = work[i];
    307303#else
    308                double alpha = !addSequence ? work[i] : work[iSequence];
    309 #endif
    310                if (fabs(alpha) < acceptablePivot)
    311                     continue;
    312                double oldValue = dj_[iSequence2];
     304      double alpha = !addSequence ? work[i] : work[iSequence];
     305#endif
     306      if (fabs(alpha) < acceptablePivot)
     307        continue;
     308      double oldValue = dj_[iSequence2];
    313309
    314                switch(getStatus(iSequence2)) {
     310      switch (getStatus(iSequence2)) {
    315311
    316                case basic:
    317                     break;
    318                case ClpSimplex::isFixed:
    319                     break;
    320                case isFree:
    321                case superBasic:
    322                     // treat dj as if zero
    323                     thetaDown = 0.0;
    324                     thetaUp = 0.0;
    325                     sequenceDown = iSequence2;
    326                     sequenceUp = iSequence2;
    327                     break;
    328                case atUpperBound:
    329                     if (alpha > 0.0) {
    330                          // test up
    331                          if (oldValue + thetaUp * alpha > dualTolerance_) {
    332                               thetaUp = (dualTolerance_ - oldValue) / alpha;
    333                               sequenceUp = iSequence2;
    334                               alphaUp = alpha;
    335                          }
    336                     } else {
    337                          // test down
    338                          if (oldValue - thetaDown * alpha > dualTolerance_) {
    339                               thetaDown = -(dualTolerance_ - oldValue) / alpha;
    340                               sequenceDown = iSequence2;
    341                               alphaDown = alpha;
    342                          }
    343                     }
    344                     break;
    345                case atLowerBound:
    346                     if (alpha < 0.0) {
    347                          // test up
    348                          if (oldValue + thetaUp * alpha < - dualTolerance_) {
    349                               thetaUp = -(dualTolerance_ + oldValue) / alpha;
    350                               sequenceUp = iSequence2;
    351                               alphaUp = alpha;
    352                          }
    353                     } else {
    354                          // test down
    355                          if (oldValue - thetaDown * alpha < -dualTolerance_) {
    356                               thetaDown = (dualTolerance_ + oldValue) / alpha;
    357                               sequenceDown = iSequence2;
    358                               alphaDown = alpha;
    359                          }
    360                     }
    361                     break;
    362                }
    363           }
    364      }
    365      if (sequenceUp >= 0) {
    366           costIncrease = thetaUp;
    367           sequenceIncrease = sequenceUp;
    368           alphaIncrease = alphaUp;
    369      }
    370      if (sequenceDown >= 0) {
    371           costDecrease = thetaDown;
    372           sequenceDecrease = sequenceDown;
    373           alphaDecrease = alphaDown;
    374      }
     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  }
    375371}
    376372/** Primal ranging.
     
    384380    When here - guaranteed optimal
    385381*/
    386 void
    387 ClpSimplexOther::primalRanging(int numberCheck, const int * which,
    388                                double * valueIncreased, int * sequenceIncreased,
    389                                double * valueDecreased, int * sequenceDecreased)
     382void ClpSimplexOther::primalRanging(int numberCheck, const int *which,
     383  double *valueIncreased, int *sequenceIncreased,
     384  double *valueDecreased, int *sequenceDecreased)
    390385{
    391      rowArray_[0]->clear();
    392      rowArray_[1]->clear();
    393      lowerIn_ = -COIN_DBL_MAX;
    394      upperIn_ = COIN_DBL_MAX;
    395      valueIn_ = 0.0;
    396      for ( int i = 0; i < numberCheck; i++) {
    397           int iSequence = which[i];
    398           double valueIncrease = COIN_DBL_MAX;
    399           double valueDecrease = COIN_DBL_MAX;
    400           int sequenceIncrease = -1;
    401           int sequenceDecrease = -1;
     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;
    402397
    403           switch(getStatus(iSequence)) {
     398    switch (getStatus(iSequence)) {
    404399
    405           case basic:
    406           case isFree:
    407           case superBasic:
    408                // Easy
    409                valueDecrease = CoinMax(0.0, upper_[iSequence] - solution_[iSequence]);
    410                valueIncrease = CoinMax(0.0, solution_[iSequence] - lower_[iSequence]);
    411                sequenceDecrease = iSequence;
    412                sequenceIncrease = iSequence;
    413                break;
    414           case isFixed:
    415           case atUpperBound:
    416           case atLowerBound: {
    417                // Non trivial
    418                // Other bound is ignored
     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
    419414#ifndef COIN_FAC_NEW
    420                unpackPacked(rowArray_[1], iSequence);
     415      unpackPacked(rowArray_[1], iSequence);
    421416#else
    422                unpack(rowArray_[1], iSequence);
    423 #endif
    424                factorization_->updateColumn(rowArray_[2], rowArray_[1]);
    425                // Get extra rows
    426                matrix_->extendUpdated(this, rowArray_[1], 0);
    427                // do ratio test
    428                checkPrimalRatios(rowArray_[1], 1);
    429                if (pivotRow_ >= 0) {
    430                     valueIncrease = theta_;
    431                     sequenceIncrease = pivotVariable_[pivotRow_];
    432                }
    433                checkPrimalRatios(rowArray_[1], -1);
    434                if (pivotRow_ >= 0) {
    435                     valueDecrease = theta_;
    436                     sequenceDecrease = pivotVariable_[pivotRow_];
    437                }
    438                rowArray_[1]->clear();
    439           }
    440           break;
    441           }
    442           double scaleFactor;
    443           if (rowScale_) {
    444                if (iSequence < numberColumns_)
    445                     scaleFactor = columnScale_[iSequence] / rhsScale_;
    446                else
    447                     scaleFactor = 1.0 / (rowScale_[iSequence-numberColumns_] * rhsScale_);
    448           } else {
    449                scaleFactor = 1.0 / rhsScale_;
    450           }
    451           if (valueIncrease < 1.0e30)
    452                valueIncrease *= scaleFactor;
    453           else
    454                valueIncrease = COIN_DBL_MAX;
    455           if (valueDecrease < 1.0e30)
    456                valueDecrease *= scaleFactor;
    457           else
    458                valueDecrease = COIN_DBL_MAX;
    459           valueIncreased[i] = valueIncrease;
    460           sequenceIncreased[i] = sequenceIncrease;
    461           valueDecreased[i] = valueDecrease;
    462           sequenceDecreased[i] = sequenceDecrease;
    463      }
     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  }
    464458}
    465459// Returns new value of whichOther when whichIn enters basis
     
    467461ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
    468462{
    469      rowArray_[0]->clear();
    470      rowArray_[1]->clear();
    471      int iSequence = whichIn;
    472      double newValue = solution_[whichOther];
    473      double alphaOther = 0.0;
    474      Status status = getStatus(iSequence);
    475      assert (status == atLowerBound || status == atUpperBound);
    476      int wayIn = (status == atLowerBound) ? 1 : -1;
     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;
    477471
    478      switch(getStatus(iSequence)) {
     472  switch (getStatus(iSequence)) {
    479473
    480      case basic:
    481      case isFree:
    482      case superBasic:
    483           assert (whichIn == whichOther);
    484           // Easy
    485           newValue = wayIn > 0 ? upper_[iSequence] : lower_[iSequence];
    486           break;
    487      case isFixed:
    488      case atUpperBound:
    489      case atLowerBound:
    490           // Non trivial
    491      {
    492           // Other bound is ignored
     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
    493487#ifndef COIN_FAC_NEW
    494           unpackPacked(rowArray_[1], iSequence);
     488      unpackPacked(rowArray_[1], iSequence);
    495489#else
    496           unpack(rowArray_[1], iSequence);
    497 #endif
    498           factorization_->updateColumn(rowArray_[2], rowArray_[1]);
    499           // Get extra rows
    500           matrix_->extendUpdated(this, rowArray_[1], 0);
    501           // do ratio test
    502           double acceptablePivot = 1.0e-7;
    503           double * work = rowArray_[1]->denseVector();
    504           int number = rowArray_[1]->getNumElements();
    505           int * which = rowArray_[1]->getIndices();
     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();
    506500
    507           // we may need to swap sign
    508           double way = wayIn;
    509           double theta = 1.0e30;
    510           for (int iIndex = 0; iIndex < number; iIndex++) {
     501      // we may need to swap sign
     502      double way = wayIn;
     503      double theta = 1.0e30;
     504      for (int iIndex = 0; iIndex < number; iIndex++) {
    511505
    512                int iRow = which[iIndex];
     506        int iRow = which[iIndex];
    513507#ifndef COIN_FAC_NEW
    514                double alpha = work[iIndex] * way;
     508        double alpha = work[iIndex] * way;
    515509#else
    516                double alpha = work[iRow] * way;
    517 #endif
    518                int iPivot = pivotVariable_[iRow];
    519                if (iPivot == whichOther) {
    520                     alphaOther = alpha;
    521                     continue;
    522                }
    523                double oldValue = solution_[iPivot];
    524                if (fabs(alpha) > acceptablePivot) {
    525                     if (alpha > 0.0) {
    526                          // basic variable going towards lower bound
    527                          double bound = lower_[iPivot];
    528                          oldValue -= bound;
    529                          if (oldValue - theta * alpha < 0.0) {
    530                               theta = CoinMax(0.0, oldValue / alpha);
    531                          }
    532                     } else {
    533                          // basic variable going towards upper bound
    534                          double bound = upper_[iPivot];
    535                          oldValue = oldValue - bound;
    536                          if (oldValue - theta * alpha > 0.0) {
    537                               theta = CoinMax(0.0, oldValue / alpha);
    538                          }
    539                     }
    540                }
    541           }
    542           if (whichIn != whichOther) {
    543                if (theta < 1.0e30)
    544                     newValue -= theta * alphaOther;
    545                else
    546                     newValue = alphaOther > 0.0 ? -1.0e30 : 1.0e30;
     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            }
    547526          } else {
    548                newValue += theta * wayIn;
    549           }
    550      }
    551      rowArray_[1]->clear();
    552      break;
    553      }
    554      double scaleFactor;
    555      if (rowScale_) {
    556           if (whichOther < numberColumns_)
    557                scaleFactor = columnScale_[whichOther] / rhsScale_;
    558           else
    559                scaleFactor = 1.0 / (rowScale_[whichOther-numberColumns_] * rhsScale_);
    560      } else {
    561           scaleFactor = 1.0 / rhsScale_;
    562      }
    563      if (newValue < 1.0e29)
    564           if (newValue > -1.0e29)
    565                newValue *= scaleFactor;
    566           else
    567                newValue = -COIN_DBL_MAX;
    568      else
    569           newValue = COIN_DBL_MAX;
    570      return newValue;
     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;
    571565}
    572566/*
     
    574568   This is used in primal ranging
    575569*/
    576 void
    577 ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
    578                                    int direction)
     570void ClpSimplexOther::checkPrimalRatios(CoinIndexedVector *rowArray,
     571  int direction)
    579572{
    580      // sequence stays as row number until end
    581      pivotRow_ = -1;
    582      double acceptablePivot = 1.0e-7;
    583      double * work = rowArray->denseVector();
    584      int number = rowArray->getNumElements();
    585      int * which = rowArray->getIndices();
     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();
    586579
    587      // we need to swap sign if going down
    588      double way = direction;
    589      theta_ = 1.0e30;
    590      for (int iIndex = 0; iIndex < number; iIndex++) {
     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++) {
    591584
    592           int iRow = which[iIndex];
     585    int iRow = which[iIndex];
    593586#ifndef COIN_FAC_NEW
    594           double alpha = work[iIndex] * way;
     587    double alpha = work[iIndex] * way;
    595588#else
    596           double alpha = work[iRow] * way;
    597 #endif
    598           int iPivot = pivotVariable_[iRow];
    599           double oldValue = solution_[iPivot];
    600           if (fabs(alpha) > acceptablePivot) {
    601                if (alpha > 0.0) {
    602                     // basic variable going towards lower bound
    603                     double bound = lower_[iPivot];
    604                     oldValue -= bound;
    605                     if (oldValue - theta_ * alpha < 0.0) {
    606                          pivotRow_ = iRow;
    607                          theta_ = CoinMax(0.0, oldValue / alpha);
    608                     }
    609                } else {
    610                     // basic variable going towards upper bound
    611                     double bound = upper_[iPivot];
    612                     oldValue = oldValue - bound;
    613                     if (oldValue - theta_ * alpha > 0.0) {
    614                          pivotRow_ = iRow;
    615                          theta_ = CoinMax(0.0, oldValue / alpha);
    616                     }
    617                }
    618           }
    619      }
     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  }
    620613}
    621614/* Write the basis in MPS format to the specified file.
     
    635628   This is based on code contributed by Thorsten Koch
    636629*/
    637 int
    638 ClpSimplexOther::writeBasis(const char *filename,
    639                             bool writeValues,
    640                             int formatType) const
     630int ClpSimplexOther::writeBasis(const char *filename,
     631  bool writeValues,
     632  int formatType) const
    641633{
    642      formatType = CoinMax(0, formatType);
    643      formatType = CoinMin(2, formatType);
    644      if (!writeValues)
    645           formatType = 0;
    646      // See if INTEL if IEEE
    647      if (formatType == 2) {
    648           // test intel here and add 1 if not intel
    649           double value = 1.0;
    650           char x[8];
    651           memcpy(x, &value, 8);
    652           if (x[0] == 63) {
    653                formatType ++; // not intel
    654           } else {
    655                assert (x[0] == 0);
    656           }
    657      }
     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  }
    658650
    659      char number[20];
    660      FILE * fp = fopen(filename, "w");
    661      if (!fp)
    662           return -1;
     651  char number[20];
     652  FILE *fp = fopen(filename, "w");
     653  if (!fp)
     654    return -1;
    663655
    664      // NAME card
     656  // NAME card
    665657
    666      // Set locale so won't get , instead of .
    667      char * saveLocale = strdup(setlocale(LC_ALL,NULL));
    668      setlocale(LC_ALL,"C");
    669      if (strcmp(strParam_[ClpProbName].c_str(), "") == 0) {
    670           fprintf(fp, "NAME          BLANK      ");
    671      } else {
    672           fprintf(fp, "NAME          %s       ", strParam_[ClpProbName].c_str());
    673      }
    674      if (formatType >= 2)
    675           fprintf(fp, "FREEIEEE");
    676      else if (writeValues)
    677           fprintf(fp, "VALUES");
    678      // finish off name
    679      fprintf(fp, "\n");
    680      int iRow = 0;
    681      for(int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    682           bool printit = false;
    683           if( getColumnStatus(iColumn) == ClpSimplex::basic) {
    684                printit = true;
    685                // Find non basic row
    686                for(; iRow < numberRows_; iRow++) {
    687                     if (getRowStatus(iRow) != ClpSimplex::basic)
    688                          break;
    689                }
    690                if (lengthNames_) {
    691                     if (iRow != numberRows_) {
    692                          fprintf(fp, " %s %-8s       %s",
    693                                  getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
    694                                  columnNames_[iColumn].c_str(),
    695                                  rowNames_[iRow].c_str());
    696                          iRow++;
    697                     } else {
    698                          // Allow for too many basics!
    699                          fprintf(fp, " BS %-8s       ",
    700                                  columnNames_[iColumn].c_str());
    701                          // Dummy row name if values
    702                          if (writeValues)
    703                               fprintf(fp, "      _dummy_");
    704                     }
    705                } else {
    706                     // no names
    707                     if (iRow != numberRows_) {
    708                          fprintf(fp, " %s C%7.7d     R%7.7d",
    709                                  getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
    710                                  iColumn, iRow);
    711                          iRow++;
    712                     } else {
    713                          // Allow for too many basics!
    714                          fprintf(fp, " BS C%7.7d", iColumn);
    715                          // Dummy row name if values
    716                          if (writeValues)
    717                               fprintf(fp, "      _dummy_");
    718                     }
    719                }
    720           } else  {
    721                if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
    722                     printit = true;
    723                     if (lengthNames_)
    724                          fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
    725                     else
    726                          fprintf(fp, " UL C%7.7d", iColumn);
    727                     // Dummy row name if values
    728                     if (writeValues)
    729                          fprintf(fp, "      _dummy_");
    730                } else if( (getColumnStatus(iColumn) == ClpSimplex::superBasic||
    731                            getColumnStatus(iColumn) == ClpSimplex::isFree)&&
    732                           writeValues) {
    733                     printit = true;
    734                     if (lengthNames_)
    735                          fprintf(fp, " BS %s", columnNames_[iColumn].c_str());
    736                     else
    737                          fprintf(fp, " BS C%7.7d", iColumn);
    738                     // Dummy row name if values
    739                     if (writeValues)
    740                          fprintf(fp, "      _dummy_");
    741                }
    742           }
    743           if (printit && writeValues) {
    744                // add value
    745                CoinConvertDouble(0, formatType, columnActivity_[iColumn], number);
    746                fprintf(fp, "     %s", number);
    747           }
    748           if (printit)
    749                fprintf(fp, "\n");
    750      }
    751      fprintf(fp, "ENDATA\n");
    752      fclose(fp);
    753      setlocale(LC_ALL,saveLocale);
    754      free(saveLocale);
    755      return 0;
     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;
    756746}
    757747// Read a basis from the given filename
    758 int
    759 ClpSimplexOther::readBasis(const char *fileName)
     748int ClpSimplexOther::readBasis(const char *fileName)
    760749{
    761      int status = 0;
    762      if (strcmp(fileName, "-") != 0 && strcmp(fileName, "stdin") != 0) {
    763           FILE *fp = fopen(fileName, "r");
    764           if (fp) {
    765                // can open - lets go for it
    766                fclose(fp);
    767           } else {
    768                handler_->message(CLP_UNABLE_OPEN, messages_)
    769                          << fileName << CoinMessageEol;
    770                return -1;
    771           }
    772      }
    773      CoinMpsIO m;
    774      m.passInMessageHandler(handler_);
    775      *m.messagesPointer() = coinMessages();
    776      bool savePrefix = m.messageHandler()->prefix();
    777      m.messageHandler()->setPrefix(handler_->prefix());
    778      status = m.readBasis(fileName, "", columnActivity_, status_ + numberColumns_,
    779                           status_,
    780                           columnNames_, numberColumns_,
    781                           rowNames_, numberRows_);
    782      m.messageHandler()->setPrefix(savePrefix);
    783      if (status >= 0) {
    784           if (!status) {
    785                // set values
    786                int iColumn, iRow;
    787                for (iRow = 0; iRow < numberRows_; iRow++) {
    788                     if (getRowStatus(iRow) == atLowerBound)
    789                          rowActivity_[iRow] = rowLower_[iRow];
    790                     else if (getRowStatus(iRow) == atUpperBound)
    791                          rowActivity_[iRow] = rowUpper_[iRow];
    792                }
    793                for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    794                     if (getColumnStatus(iColumn) == atLowerBound)
    795                          columnActivity_[iColumn] = columnLower_[iColumn];
    796                     else if (getColumnStatus(iColumn) == atUpperBound)
    797                          columnActivity_[iColumn] = columnUpper_[iColumn];
    798                }
    799           } else {
    800                memset(rowActivity_, 0, numberRows_ * sizeof(double));
    801                matrix_->times(-1.0, columnActivity_, rowActivity_);
    802           }
    803      } else {
    804           // errors
    805           handler_->message(CLP_IMPORT_ERRORS, messages_)
    806                     << status << fileName << CoinMessageEol;
    807      }
    808      return status;
     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;
    809798}
    810799/* Creates dual of a problem if looks plausible
     
    816805ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const
    817806{
    818      const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
    819      bool changed = false;
    820      int numberChanged = 0;
    821      int numberFreeColumnsInPrimal=0;
    822      int iColumn;
    823      // check if we need to change bounds to rows
    824      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    825        if (columnUpper_[iColumn] < 1.0e20) {
    826          if (columnLower_[iColumn] > -1.0e20) {
    827                changed = true;
    828                numberChanged++;
    829           }
    830        } else if (columnLower_[iColumn] < -1.0e20) {
    831          numberFreeColumnsInPrimal++;
    832        }
    833      }
    834      int iRow;
    835      int numberExtraRows = 0;
    836      int numberFreeColumnsInDual=0;
    837      if (numberChanged <= fractionColumnRanges * numberColumns_) {
    838           for (iRow = 0; iRow < numberRows_; iRow++) {
    839                if (rowLower_[iRow] > -1.0e20 &&
    840                          rowUpper_[iRow] < 1.0e20) {
    841                     if (rowUpper_[iRow] != rowLower_[iRow])
    842                          numberExtraRows++;
    843                     else
    844                       numberFreeColumnsInDual++;
    845                }
    846           }
    847           if (numberExtraRows > fractionRowRanges * numberRows_)
    848                return NULL;
    849      } else {
    850           return NULL;
    851      }
    852      printf("would have %d free columns in primal, %d in dual\n",
    853             numberFreeColumnsInPrimal,numberFreeColumnsInDual);
    854      if (4*(numberFreeColumnsInDual-numberFreeColumnsInPrimal)>
    855          numberColumns_&&fractionRowRanges<1.0)
    856        return NULL; //dangerous (well anyway in dual)
    857      if (changed) {
    858           ClpSimplex * model3 = new ClpSimplex(*model2);
    859           CoinBuild build;
    860           double one = 1.0;
    861           int numberColumns = model3->numberColumns();
    862           const double * columnLower = model3->columnLower();
    863           const double * columnUpper = model3->columnUpper();
    864           for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    865                if (columnUpper[iColumn] < 1.0e20 &&
    866                          columnLower[iColumn] > -1.0e20) {
    867                     if (fabs(columnLower[iColumn]) < fabs(columnUpper[iColumn])) {
    868                          double value = columnUpper[iColumn];
    869                          model3->setColumnUpper(iColumn, COIN_DBL_MAX);
    870                          build.addRow(1, &iColumn, &one, -COIN_DBL_MAX, value);
    871                     } else {
    872                          double value = columnLower[iColumn];
    873                          model3->setColumnLower(iColumn, -COIN_DBL_MAX);
    874                          build.addRow(1, &iColumn, &one, value, COIN_DBL_MAX);
    875                     }
    876                }
    877           }
    878           model3->addRows(build);
    879           model2 = model3;
    880      }
    881      int numberColumns = model2->numberColumns();
    882      const double * columnLower = model2->columnLower();
    883      const double * columnUpper = model2->columnUpper();
    884      int numberRows = model2->numberRows();
    885      double * rowLower = CoinCopyOfArray(model2->rowLower(), numberRows);
    886      double * rowUpper = CoinCopyOfArray(model2->rowUpper(), numberRows);
     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);
    887873
    888      const double * objective = model2->objective();
    889      CoinPackedMatrix * matrix = model2->matrix();
    890      // get transpose
    891      CoinPackedMatrix rowCopy = *matrix;
    892      const int * row = matrix->getIndices();
    893      const int * columnLength = matrix->getVectorLengths();
    894      const CoinBigIndex * columnStart = matrix->getVectorStarts();
    895      const double * elementByColumn = matrix->getElements();
    896      double objOffset = 0.0;
    897      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    898           double offset = 0.0;
    899           double objValue = optimizationDirection_ * objective[iColumn];
    900           if (columnUpper[iColumn] > 1.0e20) {
    901                if (columnLower[iColumn] > -1.0e20)
    902                     offset = columnLower[iColumn];
    903           } else if (columnLower[iColumn] < -1.0e20) {
    904                offset = columnUpper[iColumn];
    905           } else {
    906                // taken care of before
    907                abort();
    908           }
    909           if (offset) {
    910                objOffset += offset * objValue;
    911                for (CoinBigIndex j = columnStart[iColumn];
    912                          j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    913                     int iRow = row[j];
    914                     if (rowLower[iRow] > -1.0e20)
    915                          rowLower[iRow] -= offset * elementByColumn[j];
    916                     if (rowUpper[iRow] < 1.0e20)
    917                          rowUpper[iRow] -= offset * elementByColumn[j];
    918                }
    919           }
    920      }
    921      int * which = new int[numberRows+numberExtraRows];
    922      rowCopy.reverseOrdering();
    923      rowCopy.transpose();
    924      double * fromRowsLower = new double[numberRows+numberExtraRows];
    925      double * fromRowsUpper = new double[numberRows+numberExtraRows];
    926      double * newObjective = new double[numberRows+numberExtraRows];
    927      double * fromColumnsLower = new double[numberColumns];
    928      double * fromColumnsUpper = new double[numberColumns];
    929      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    930           double objValue = optimizationDirection_ * objective[iColumn];
    931           // Offset is already in
    932           if (columnUpper[iColumn] > 1.0e20) {
    933                if (columnLower[iColumn] > -1.0e20) {
    934                     fromColumnsLower[iColumn] = -COIN_DBL_MAX;
    935                     fromColumnsUpper[iColumn] = objValue;
    936                } else {
    937                     // free
    938                     fromColumnsLower[iColumn] = objValue;
    939                     fromColumnsUpper[iColumn] = objValue;
    940                }
    941           } else if (columnLower[iColumn] < -1.0e20) {
    942                fromColumnsLower[iColumn] = objValue;
    943                fromColumnsUpper[iColumn] = COIN_DBL_MAX;
    944           } else {
    945                abort();
    946           }
    947      }
    948      int kRow = 0;
    949      int kExtraRow = numberRows;
    950      for (iRow = 0; iRow < numberRows; iRow++) {
    951           if (rowLower[iRow] < -1.0e20) {
    952                assert (rowUpper[iRow] < 1.0e20);
    953                newObjective[kRow] = -rowUpper[iRow];
    954                fromRowsLower[kRow] = -COIN_DBL_MAX;
    955                fromRowsUpper[kRow] = 0.0;
    956                which[kRow] = iRow;
    957                kRow++;
    958           } else if (rowUpper[iRow] > 1.0e20) {
    959                newObjective[kRow] = -rowLower[iRow];
    960                fromRowsLower[kRow] = 0.0;
    961                fromRowsUpper[kRow] = COIN_DBL_MAX;
    962                which[kRow] = iRow;
    963                kRow++;
    964           } else {
    965                if (rowUpper[iRow] == rowLower[iRow]) {
    966                     newObjective[kRow] = -rowLower[iRow];
    967                     fromRowsLower[kRow] = -COIN_DBL_MAX;;
    968                     fromRowsUpper[kRow] = COIN_DBL_MAX;
    969                     which[kRow] = iRow;
    970                     kRow++;
    971                } else {
    972                     // range
    973                     newObjective[kRow] = -rowUpper[iRow];
    974                     fromRowsLower[kRow] = -COIN_DBL_MAX;
    975                     fromRowsUpper[kRow] = 0.0;
    976                     which[kRow] = iRow;
    977                     kRow++;
    978                     newObjective[kExtraRow] = -rowLower[iRow];
    979                     fromRowsLower[kExtraRow] = 0.0;
    980                     fromRowsUpper[kExtraRow] = COIN_DBL_MAX;
    981                     which[kExtraRow] = iRow;
    982                     kExtraRow++;
    983                }
    984           }
    985      }
    986      if (numberExtraRows) {
    987           CoinPackedMatrix newCopy;
    988           newCopy.setExtraGap(0.0);
    989           newCopy.setExtraMajor(0.0);
    990           newCopy.submatrixOfWithDuplicates(rowCopy, kExtraRow, which);
    991           rowCopy = newCopy;
    992      }
    993      ClpSimplex * modelDual = new ClpSimplex();
    994      modelDual->passInEventHandler(eventHandler_);
    995      modelDual->loadProblem(rowCopy, fromRowsLower, fromRowsUpper, newObjective,
    996                             fromColumnsLower, fromColumnsUpper);
    997      modelDual->setObjectiveOffset(objOffset);
    998      modelDual->setDualBound(model2->dualBound());
    999      modelDual->setInfeasibilityCost(model2->infeasibilityCost());
    1000      modelDual->setDualTolerance(model2->dualTolerance());
    1001      modelDual->setPrimalTolerance(model2->primalTolerance());
    1002      modelDual->setPerturbation(model2->perturbation());
    1003      modelDual->setSpecialOptions(model2->specialOptions());
    1004      modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
    1005      modelDual->setMaximumIterations(model2->maximumIterations());
    1006      modelDual->setFactorizationFrequency(model2->factorizationFrequency());
    1007      modelDual->setLogLevel(model2->logLevel());
    1008      delete [] fromRowsLower;
    1009      delete [] fromRowsUpper;
    1010      delete [] fromColumnsLower;
    1011      delete [] fromColumnsUpper;
    1012      delete [] newObjective;
    1013      delete [] which;
    1014      delete [] rowLower;
    1015      delete [] rowUpper;
    1016      if (changed)
    1017           delete model2;
    1018      modelDual->createStatus();
    1019      return modelDual;
     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;
    10201007}
    10211008// Restores solution from dualized problem
    1022 int
    1023 ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem,
    1024                                  bool checkAccuracy)
     1009int ClpSimplexOther::restoreFromDual(const ClpSimplex *dualProblem,
     1010  bool checkAccuracy)
    10251011{
    1026      int returnCode = 0;;
    1027      createStatus();
    1028      // Number of rows in dual problem was original number of columns
    1029      assert (numberColumns_ == dualProblem->numberRows());
    1030      // If slack on d-row basic then column at bound otherwise column basic
    1031      // If d-column basic then rhs tight
    1032      int numberBasic = 0;
    1033      int iRow, iColumn = 0;
    1034      // Get number of extra rows from ranges
    1035      int numberExtraRows = 0;
    1036      for (iRow = 0; iRow < numberRows_; iRow++) {
    1037           if (rowLower_[iRow] > -1.0e20 &&
    1038                     rowUpper_[iRow] < 1.0e20) {
    1039                if (rowUpper_[iRow] != rowLower_[iRow])
    1040                     numberExtraRows++;
    1041           }
    1042      }
    1043      const double * objective = this->objective();
    1044      const double * dualDual = dualProblem->dualRowSolution();
    1045      const double * dualDj = dualProblem->dualColumnSolution();
    1046      const double * dualSol = dualProblem->primalColumnSolution();
    1047      const double * dualActs = dualProblem->primalRowSolution();
     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();
    10481034#if 0
    10491035     ClpSimplex thisCopy = *this;
     
    10751061                 primalSol[iColumn], primalDj[iColumn]);
    10761062#endif
    1077      // position at bound information
    1078      int jColumn = numberRows_;
    1079      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1080           double objValue = optimizationDirection_ * objective[iColumn];
    1081           Status status = dualProblem->getRowStatus(iColumn);
    1082           double otherValue = COIN_DBL_MAX;
    1083           if (columnUpper_[iColumn] < 1.0e20 &&
    1084                     columnLower_[iColumn] > -1.0e20) {
    1085                if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) {
    1086                     otherValue = columnUpper_[iColumn] + dualDj[jColumn];
    1087                } else {
    1088                     otherValue = columnLower_[iColumn] + dualDj[jColumn];
    1089                }
    1090                jColumn++;
    1091           }
    1092           if (status == basic) {
    1093                // column is at bound
    1094                if (otherValue == COIN_DBL_MAX) {
    1095                     reducedCost_[iColumn] = objValue - dualActs[iColumn];
    1096                     if (columnUpper_[iColumn] > 1.0e20) {
    1097                          if (columnLower_[iColumn] > -1.0e20) {
    1098                               if (columnUpper_[iColumn] > columnLower_[iColumn])
    1099                                    setColumnStatus(iColumn, atLowerBound);
    1100                               else
    1101                                    setColumnStatus(iColumn, isFixed);
    1102                               columnActivity_[iColumn] = columnLower_[iColumn];
    1103                          } else {
    1104                               // free
    1105                               setColumnStatus(iColumn, isFree);
    1106                               columnActivity_[iColumn] = 0.0;
    1107                          }
    1108                     } else {
    1109                          setColumnStatus(iColumn, atUpperBound);
    1110                          columnActivity_[iColumn] = columnUpper_[iColumn];
    1111                     }
    1112                } else {
    1113                     reducedCost_[iColumn] = objValue - dualActs[iColumn];
    1114                     //printf("other dual sol %g\n",otherValue);
    1115                     if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
    1116                          if (columnUpper_[iColumn] > columnLower_[iColumn])
    1117                               setColumnStatus(iColumn, atLowerBound);
    1118                          else
    1119                               setColumnStatus(iColumn, isFixed);
    1120                          columnActivity_[iColumn] = columnLower_[iColumn];
    1121                     } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
    1122                          if (columnUpper_[iColumn] > columnLower_[iColumn])
    1123                               setColumnStatus(iColumn, atUpperBound);
    1124                          else
    1125                               setColumnStatus(iColumn, isFixed);
    1126                          columnActivity_[iColumn] = columnUpper_[iColumn];
    1127                     } else {
    1128                          setColumnStatus(iColumn, superBasic);
    1129                          columnActivity_[iColumn] = otherValue;
    1130                     }
    1131                }
     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];
    11321088          } else {
    1133                if (otherValue == COIN_DBL_MAX) {
    1134                     // column basic
    1135                     setColumnStatus(iColumn, basic);
    1136                     numberBasic++;
    1137                     if (columnLower_[iColumn] > -1.0e20) {
    1138                          columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn];
    1139                     } else if (columnUpper_[iColumn] < 1.0e20) {
    1140                          columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn];
    1141                     } else {
    1142                          columnActivity_[iColumn] = -dualDual[iColumn];
    1143                     }
    1144                     reducedCost_[iColumn] = 0.0;
    1145                } else {
    1146                     // may be at other bound
    1147                     //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
    1148                     if (dualProblem->getColumnStatus(jColumn - 1) != basic) {
    1149                          // column basic
    1150                          setColumnStatus(iColumn, basic);
    1151                          numberBasic++;
    1152                          //printf("Col %d otherV %g dualDual %g\n",iColumn,
    1153                          // otherValue,dualDual[iColumn]);
    1154                          columnActivity_[iColumn] = -dualDual[iColumn];
    1155                          columnActivity_[iColumn] = otherValue;
    1156                          reducedCost_[iColumn] = 0.0;
    1157                     } else {
    1158                          reducedCost_[iColumn] = objValue - dualActs[iColumn];
    1159                          if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
    1160                               if (columnUpper_[iColumn] > columnLower_[iColumn])
    1161                                    setColumnStatus(iColumn, atLowerBound);
    1162                               else
    1163                                    setColumnStatus(iColumn, isFixed);
    1164                               columnActivity_[iColumn] = columnLower_[iColumn];
    1165                          } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
    1166                               if (columnUpper_[iColumn] > columnLower_[iColumn])
    1167                                    setColumnStatus(iColumn, atUpperBound);
    1168                               else
    1169                                    setColumnStatus(iColumn, isFixed);
    1170                               columnActivity_[iColumn] = columnUpper_[iColumn];
    1171                          } else {
    1172                               setColumnStatus(iColumn, superBasic);
    1173                               columnActivity_[iColumn] = otherValue;
    1174                          }
    1175                     }
    1176                }
    1177           }
    1178      }
    1179      // now rows
    1180      int kExtraRow = jColumn;
    1181      int numberRanges = 0;
    1182      for (iRow = 0; iRow < numberRows_; iRow++) {
    1183           Status status = dualProblem->getColumnStatus(iRow);
    1184           if (status == basic) {
    1185                // row is at bound
    1186                dual_[iRow] = dualSol[iRow];;
     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];
    11871156          } else {
    1188                // row basic
    1189                setRowStatus(iRow, basic);
    1190                numberBasic++;
    1191                dual_[iRow] = 0.0;
    1192           }
    1193           if (rowLower_[iRow] < -1.0e20) {
    1194                if (status == basic) {
    1195                     rowActivity_[iRow] = rowUpper_[iRow];
    1196                     setRowStatus(iRow, atUpperBound);
    1197                } else {
    1198                     // might be stopped assert (dualDj[iRow] < 1.0e-5);
    1199                     rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
    1200                }
    1201           } else if (rowUpper_[iRow] > 1.0e20) {
    1202                if (status == basic) {
    1203                     rowActivity_[iRow] = rowLower_[iRow];
    1204                     setRowStatus(iRow, atLowerBound);
    1205                } else {
    1206                     rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
    1207                     // might be stopped assert (dualDj[iRow] > -1.0e-5);
    1208                }
    1209           } else {
    1210                if (rowUpper_[iRow] == rowLower_[iRow]) {
    1211                     rowActivity_[iRow] = rowLower_[iRow];
    1212                     if (status == basic) {
    1213                          setRowStatus(iRow, isFixed);
    1214                     }
    1215                } else {
    1216                     // range
    1217                     numberRanges++;
    1218                     Status statusL = dualProblem->getColumnStatus(kExtraRow);
    1219                     //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
    1220                     //     iRow,status,kExtraRow,statusL, dualSol[iRow],
    1221                     //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
    1222                     if (status == basic) {
    1223                          // might be stopped assert (statusL != basic);
    1224                          rowActivity_[iRow] = rowUpper_[iRow];
    1225                          setRowStatus(iRow, atUpperBound);
    1226                     } else if (statusL == basic) {
    1227                          numberBasic--; // already counted
    1228                          rowActivity_[iRow] = rowLower_[iRow];
    1229                          setRowStatus(iRow, atLowerBound);
    1230                          dual_[iRow] = dualSol[kExtraRow];;
    1231                     } else {
    1232                          rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
    1233                          // might be stopped assert (dualDj[iRow] < 1.0e-5);
    1234                          // row basic
    1235                          //setRowStatus(iRow,basic);
    1236                          //numberBasic++;
    1237                          dual_[iRow] = 0.0;
    1238                     }
    1239                     kExtraRow++;
    1240                }
    1241           }
    1242      }
    1243      if (numberBasic != numberRows_ && 0) {
    1244           printf("Bad basis - ranges - coding needed\n");
    1245           assert (numberRanges);
    1246           abort();
    1247      }
    1248      if (optimizationDirection_ < 0.0) {
    1249           for (iRow = 0; iRow < numberRows_; iRow++) {
    1250                dual_[iRow] = -dual_[iRow];
    1251           }
    1252      }
    1253      // redo row activities
    1254      memset(rowActivity_, 0, numberRows_ * sizeof(double));
    1255      matrix_->times(1.0, columnActivity_, rowActivity_);
    1256      // redo reduced costs
    1257      memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
    1258      matrix_->transposeTimes(-1.0, dual_, reducedCost_);
    1259      checkSolutionInternal();
    1260      if (sumDualInfeasibilities_ > 1.0e-5 || sumPrimalInfeasibilities_ > 1.0e-5) {
    1261           returnCode = 1;
     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;
    12621249#ifdef CLP_INVESTIGATE
    1263           printf("There are %d dual infeasibilities summing to %g ",
    1264                  numberDualInfeasibilities_, sumDualInfeasibilities_);
    1265           printf("and %d primal infeasibilities summing to %g\n",
    1266                  numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
    1267 #endif
    1268      }
    1269      // Below will go to ..DEBUG later
     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
    12701257#if 1 //ndef NDEBUG
    1271      if (checkAccuracy) {
    1272        // Check if correct
    1273        double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
    1274        double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
    1275        double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
    1276        double * dual = CoinCopyOfArray(dual_, numberRows_);
    1277        this->dual(); //primal();
    1278        CoinRelFltEq eq(1.0e-5);
    1279        for (iRow = 0; iRow < numberRows_; iRow++) {
    1280         assert(eq(dual[iRow], dual_[iRow]));
    1281        }
    1282        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1283         assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
    1284        }
    1285        for (iRow = 0; iRow < numberRows_; iRow++) {
    1286         assert(eq(rowActivity[iRow], rowActivity_[iRow]));
    1287        }
    1288        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1289         assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
    1290        }
    1291        delete [] columnActivity;
    1292        delete [] rowActivity;
    1293        delete [] reducedCost;
    1294        delete [] dual;
    1295      }
    1296 #endif
    1297      return returnCode;
     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;
    12981285}
    12991286/* Sets solution in dualized problem
    13001287   non-zero return code indicates minor problems
    13011288*/
    1302 int
    1303 ClpSimplexOther::setInDual(ClpSimplex * dualProblem)
     1289int ClpSimplexOther::setInDual(ClpSimplex *dualProblem)
    13041290{
    13051291  // Number of rows in dual problem was original number of columns
    1306   assert (numberColumns_ == dualProblem->numberRows());
     1292  assert(numberColumns_ == dualProblem->numberRows());
    13071293  // out If slack on d-row basic then column at bound otherwise column basic
    13081294  // out If d-column basic then rhs tight
     
    13101296  // if column basic then slack on d-row at bound
    13111297  // if rhs non-basic then d-column basic
    1312   // if rhs basic then d-column ? 
     1298  // if rhs basic then d-column ?
    13131299  int numberBasic = 0;
    13141300  int iRow, iColumn = 0;
     
    13171303  //double * dualDual = dualProblem->dualRowSolution();
    13181304  //double * dualDj = dualProblem->dualColumnSolution();
    1319   double * dualSol = dualProblem->primalColumnSolution();
     1305  double *dualSol = dualProblem->primalColumnSolution();
    13201306  //double * dualActs = dualProblem->primalRowSolution();
    1321   const double * lower = dualProblem->columnLower();
    1322   const double * upper = dualProblem->columnUpper();
     1307  const double *lower = dualProblem->columnLower();
     1308  const double *upper = dualProblem->columnUpper();
    13231309  // position at bound information
    13241310  int jColumn = numberRows_;
     
    13271313    Status statusD = dualProblem->getRowStatus(iColumn);
    13281314    Status statusDJ = dualProblem->getColumnStatus(jColumn);
    1329     if (status==atLowerBound||
    1330         status==isFixed||
    1331         status==atUpperBound) {
    1332       dualProblem->setRowStatus(iColumn,basic);
     1315    if (status == atLowerBound || status == isFixed || status == atUpperBound) {
     1316      dualProblem->setRowStatus(iColumn, basic);
    13331317      numberBasic++;
    1334       if (columnUpper_[iColumn] < 1.0e20 &&
    1335           columnLower_[iColumn] > -1.0e20) {
    1336         bool mainLower =(fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn]));
    1337         // fix this
    1338         if (mainLower) {
    1339           if (status==atUpperBound) {
    1340             dualProblem->setColumnStatus(jColumn,atUpperBound);
    1341           } else {
    1342             dualProblem->setColumnStatus(jColumn,atUpperBound);
    1343           }
    1344         } else {
    1345           if (status==atUpperBound) {
    1346             dualProblem->setColumnStatus(jColumn,atLowerBound);
    1347           } else {
    1348             dualProblem->setColumnStatus(jColumn,atLowerBound);
    1349           }
    1350         }
    1351         assert(statusDJ == dualProblem->getColumnStatus(jColumn));
    1352         jColumn++;
    1353       }
    1354     } else if (status==isFree) {
    1355       dualProblem->setRowStatus(iColumn,basic);
     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);
    13561339      numberBasic++;
    13571340    } else {
    1358       assert (status==basic);
     1341      assert(status == basic);
    13591342      //numberBasic++;
    13601343    }
     
    13681351      // dual variable is at bound
    13691352      if (!lower[iRow]) {
    1370         dualProblem->setColumnStatus(iRow,atLowerBound);
     1353        dualProblem->setColumnStatus(iRow, atLowerBound);
    13711354      } else if (!upper[iRow]) {
    1372         dualProblem->setColumnStatus(iRow,atUpperBound);
     1355        dualProblem->setColumnStatus(iRow, atUpperBound);
    13731356      } else {
    1374         dualProblem->setColumnStatus(iRow,isFree);
    1375         dualSol[iRow]=0.0;
     1357        dualProblem->setColumnStatus(iRow, isFree);
     1358        dualSol[iRow] = 0.0;
    13761359      }
    13771360    } else {
    13781361      // dual variable is basic
    1379       dualProblem->setColumnStatus(iRow,basic);
     1362      dualProblem->setColumnStatus(iRow, basic);
    13801363      numberBasic++;
    13811364    }
    1382     if (rowLower_[iRow] < -1.0e20 &&
    1383         rowUpper_[iRow] > 1.0e20) {
     1365    if (rowLower_[iRow] < -1.0e20 && rowUpper_[iRow] > 1.0e20) {
    13841366      if (rowUpper_[iRow] != rowLower_[iRow]) {
    1385         printf("can't handle ranges yet\n");
    1386         abort();
     1367        printf("can't handle ranges yet\n");
     1368        abort();
    13871369      }
    13881370    }
     
    13991381*/
    14001382ClpSimplex *
    1401 ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
    1402                         int & nBound, bool moreBounds, bool tightenBounds)
     1383ClpSimplexOther::crunch(double *rhs, int *whichRow, int *whichColumn,
     1384  int &nBound, bool moreBounds, bool tightenBounds)
    14031385{
    1404      //#define CHECK_STATUS
     1386  //#define CHECK_STATUS
    14051387#ifdef CHECK_STATUS
    1406      {
    1407           int n = 0;
    1408           int i;
    1409           for (i = 0; i < numberColumns_; i++)
    1410                if (getColumnStatus(i) == ClpSimplex::basic)
    1411                     n++;
    1412           for (i = 0; i < numberRows_; i++)
    1413                if (getRowStatus(i) == ClpSimplex::basic)
    1414                     n++;
    1415           assert (n == numberRows_);
    1416      }
     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  }
    14171399#endif
    14181400
    1419      const double * element = matrix_->getElements();
    1420      const int * row = matrix_->getIndices();
    1421      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1422      const int * columnLength = matrix_->getVectorLengths();
     1401  const double *element = matrix_->getElements();
     1402  const int *row = matrix_->getIndices();
     1403  const CoinBigIndex *columnStart = matrix_->getVectorStarts();
     1404  const int *columnLength = matrix_->getVectorLengths();
    14231405
    1424      CoinZeroN(rhs, numberRows_);
    1425      int iColumn;
    1426      int iRow;
    1427      CoinZeroN(whichRow, numberRows_);
    1428      int * backColumn = whichColumn + numberColumns_;
    1429      int numberRows2 = 0;
    1430      int numberColumns2 = 0;
    1431      double offset = 0.0;
    1432      const double * objective = this->objective();
    1433      double * solution = columnActivity_;
    1434      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1435           double lower = columnLower_[iColumn];
    1436           double upper = columnUpper_[iColumn];
    1437           if (upper > lower || getColumnStatus(iColumn) == ClpSimplex::basic) {
    1438                backColumn[iColumn] = numberColumns2;
    1439                whichColumn[numberColumns2++] = iColumn;
    1440                for (CoinBigIndex j = columnStart[iColumn];
    1441                          j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1442                     int iRow = row[j];
    1443                     int n = whichRow[iRow];
    1444                     if (n == 0 && element[j])
    1445                          whichRow[iRow] = -iColumn - 1;
    1446                     else if (n < 0)
    1447                          whichRow[iRow] = 2;
    1448                }
    1449           } else {
    1450                // fixed
    1451                backColumn[iColumn] = -1;
    1452                solution[iColumn] = upper;
    1453                if (upper) {
    1454                     offset += objective[iColumn] * upper;
    1455                     for (CoinBigIndex j = columnStart[iColumn];
    1456                               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1457                          int iRow = row[j];
    1458                          double value = element[j];
    1459                          rhs[iRow] += upper * value;
    1460                     }
    1461                }
    1462           }
    1463      }
    1464      int returnCode = 0;
    1465      double tolerance = primalTolerance();
    1466      nBound = 2 * numberRows_;
    1467      for (iRow = 0; iRow < numberRows_; iRow++) {
    1468           int n = whichRow[iRow];
    1469           if (n > 0) {
    1470                whichRow[numberRows2++] = iRow;
    1471           } else if (n < 0) {
    1472                //whichRow[numberRows2++]=iRow;
    1473                //continue;
    1474                // Can only do in certain circumstances as we don't know current value
    1475                if (rowLower_[iRow] == rowUpper_[iRow] || getRowStatus(iRow) == ClpSimplex::basic) {
    1476                     // save row and column for bound
    1477                     whichRow[--nBound] = iRow;
    1478                     whichRow[nBound+numberRows_] = -n - 1;
    1479                } else if (moreBounds) {
    1480                     // save row and column for bound
    1481                     whichRow[--nBound] = iRow;
    1482                     whichRow[nBound+numberRows_] = -n - 1;
    1483                } else {
    1484                     whichRow[numberRows2++] = iRow;
    1485                }
    1486           } else {
    1487                // empty
    1488                double rhsValue = rhs[iRow];
    1489                if (rhsValue < rowLower_[iRow] - tolerance || rhsValue > rowUpper_[iRow] + tolerance) {
    1490                     returnCode = 1; // infeasible
    1491                }
    1492           }
    1493      }
    1494      ClpSimplex * small = NULL;
    1495      if (!returnCode) {
    1496        //printf("CRUNCH from (%d,%d) to (%d,%d)\n",
    1497        //     numberRows_,numberColumns_,numberRows2,numberColumns2);
    1498           small = new ClpSimplex(this, numberRows2, whichRow,
    1499                                  numberColumns2, whichColumn, true, false);
     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);
    15001482#if 0
    15011483          ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
     
    15061488          }
    15071489#endif
    1508           // Set some stuff
    1509           small->setDualBound(dualBound_);
    1510           small->setInfeasibilityCost(infeasibilityCost_);
    1511           small->setSpecialOptions(specialOptions_);
    1512           small->setPerturbation(perturbation_);
    1513           small->defaultFactorizationFrequency();
    1514           small->setAlphaAccuracy(alphaAccuracy_);
    1515           // If no rows left then no tightening!
    1516           if (!numberRows2 || !numberColumns2)
    1517                tightenBounds = false;
     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;
    15181500
    1519           CoinBigIndex numberElements = getNumElements();
    1520           CoinBigIndex numberElements2 = small->getNumElements();
    1521           small->setObjectiveOffset(objectiveOffset() - offset);
    1522           handler_->message(CLP_CRUNCH_STATS, messages_)
    1523                     << numberRows2 << -(numberRows_ - numberRows2)
    1524                     << numberColumns2 << -(numberColumns_ - numberColumns2)
    1525                     << numberElements2 << -(numberElements - numberElements2)
    1526                     << CoinMessageEol;
    1527           // And set objective value to match
    1528           small->setObjectiveValue(this->objectiveValue());
    1529           double * rowLower2 = small->rowLower();
    1530           double * rowUpper2 = small->rowUpper();
    1531           int jRow;
    1532           for (jRow = 0; jRow < numberRows2; jRow++) {
    1533                iRow = whichRow[jRow];
    1534                if (rowLower2[jRow] > -1.0e20)
    1535                     rowLower2[jRow] -= rhs[iRow];
    1536                if (rowUpper2[jRow] < 1.0e20)
    1537                     rowUpper2[jRow] -= rhs[iRow];
    1538           }
    1539           // and bounds
    1540           double * columnLower2 = small->columnLower();
    1541           double * columnUpper2 = small->columnUpper();
    1542           const char * integerInformation = integerType_;
    1543           for (jRow = nBound; jRow < 2 * numberRows_; jRow++) {
    1544                iRow = whichRow[jRow];
    1545                iColumn = whichRow[jRow+numberRows_];
    1546                double lowerRow = rowLower_[iRow];
    1547                if (lowerRow > -1.0e20)
    1548                     lowerRow -= rhs[iRow];
    1549                double upperRow = rowUpper_[iRow];
    1550                if (upperRow < 1.0e20)
    1551                     upperRow -= rhs[iRow];
    1552                int jColumn = backColumn[iColumn];
    1553                double lower = columnLower2[jColumn];
    1554                double upper = columnUpper2[jColumn];
    1555                double value = 0.0;
    1556                for (CoinBigIndex j = columnStart[iColumn];
    1557                          j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1558                     if (iRow == row[j]) {
    1559                          value = element[j];
    1560                          break;
     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                      }
    15611735                    }
    1562                }
    1563                assert (value);
    1564                // convert rowLower and Upper to implied bounds on column
    1565                double newLower = -COIN_DBL_MAX;
    1566                double newUpper = COIN_DBL_MAX;
    1567                if (value > 0.0) {
    1568                     if (lowerRow > -1.0e20)
    1569                          newLower = lowerRow / value;
    1570                     if (upperRow < 1.0e20)
    1571                          newUpper = upperRow / value;
    1572                } else {
    1573                     if (upperRow < 1.0e20)
    1574                          newLower = upperRow / value;
    1575                     if (lowerRow > -1.0e20)
    1576                          newUpper = lowerRow / value;
    1577                }
    1578                if (integerInformation && integerInformation[iColumn]) {
    1579                     if (newLower - floor(newLower) < 10.0 * tolerance)
    1580                          newLower = floor(newLower);
    1581                     else
    1582                          newLower = ceil(newLower);
    1583                     if (ceil(newUpper) - newUpper < 10.0 * tolerance)
    1584                          newUpper = ceil(newUpper);
    1585                     else
    1586                          newUpper = floor(newUpper);
    1587                }
    1588                newLower = CoinMax(lower, newLower);
    1589                newUpper = CoinMin(upper, newUpper);
    1590                if (newLower > newUpper + tolerance) {
    1591                     //printf("XXYY inf on bound\n");
    1592                     returnCode = 1;
    1593                }
    1594                columnLower2[jColumn] = newLower;
    1595                columnUpper2[jColumn] = CoinMax(newLower, newUpper);
    1596                if (getRowStatus(iRow) != ClpSimplex::basic) {
    1597                     if (getColumnStatus(iColumn) == ClpSimplex::basic) {
    1598                          if (columnLower2[jColumn] == columnUpper2[jColumn]) {
    1599                               // can only get here if will be fixed
    1600                               small->setColumnStatus(jColumn, ClpSimplex::isFixed);
    1601                          } else {
    1602                               // solution is valid
    1603                               if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) <
    1604                                         fabs(columnActivity_[iColumn] - columnUpper2[jColumn]))
    1605                                    small->setColumnStatus(jColumn, ClpSimplex::atLowerBound);
    1606                               else
    1607                                    small->setColumnStatus(jColumn, ClpSimplex::atUpperBound);
    1608                          }
    1609                     } else {
    1610                          //printf("what now neither basic\n");
    1611                     }
    1612                }
    1613           }
    1614           if (returnCode) {
    1615                delete small;
    1616                small = NULL;
    1617           } else if (tightenBounds && integerInformation) {
    1618                // See if we can tighten any bounds
    1619                // use rhs for upper and small duals for lower
    1620                double * up = rhs;
    1621                double * lo = small->dualRowSolution();
    1622                const double * element = small->clpMatrix()->getElements();
    1623                const int * row = small->clpMatrix()->getIndices();
    1624                const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
    1625                //const int * columnLength = small->clpMatrix()->getVectorLengths();
    1626                CoinZeroN(lo, numberRows2);
    1627                CoinZeroN(up, numberRows2);
    1628                for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
    1629                     double upper = columnUpper2[iColumn];
    1630                     double lower = columnLower2[iColumn];
    1631                     //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
    1632                     for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
    1633                          int iRow = row[j];
    1634                          double value = element[j];
    1635                          if (value > 0.0) {
    1636                               if (upper < 1.0e20)
    1637                                    up[iRow] += upper * value;
    1638                               else
    1639                                    up[iRow] = COIN_DBL_MAX;
    1640                               if (lower > -1.0e20)
    1641                                    lo[iRow] += lower * value;
    1642                               else
    1643                                    lo[iRow] = -COIN_DBL_MAX;
    1644                          } else {
    1645                               if (upper < 1.0e20)
    1646                                    lo[iRow] += upper * value;
    1647                               else
    1648                                    lo[iRow] = -COIN_DBL_MAX;
    1649                               if (lower > -1.0e20)
    1650                                    up[iRow] += lower * value;
    1651                               else
    1652                                    up[iRow] = COIN_DBL_MAX;
    1653                          }
    1654                     }
    1655                }
    1656                double * rowLower2 = small->rowLower();
    1657                double * rowUpper2 = small->rowUpper();
    1658                bool feasible = true;
    1659                // make safer
    1660                for (int iRow = 0; iRow < numberRows2; iRow++) {
    1661                     double lower = lo[iRow];
    1662                     if (lower > rowUpper2[iRow] + tolerance) {
    1663                          feasible = false;
    1664                          break;
    1665                     } else {
    1666                          lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance;
    1667                     }
    1668                     double upper = up[iRow];
    1669                     if (upper < rowLower2[iRow] - tolerance) {
    1670                          feasible = false;
    1671                          break;
    1672                     } else {
    1673                          up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
    1674                     }
    1675                     // tighten row bounds
    1676                     if (lower>-1.0e10)
    1677                       rowLower2[iRow] = CoinMax(rowLower2[iRow],
    1678                                                 lower - 1.0e-6*(1.0+fabs(lower)));
    1679                     if (upper<1.0e10)
    1680                       rowUpper2[iRow] = CoinMin(rowUpper2[iRow],
    1681                                                 upper + 1.0e-6*(1.0+fabs(upper)));
    1682                }
    1683                if (!feasible) {
     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
    16841742                    delete small;
    16851743                    small = NULL;
    1686                } else {
    1687                     // and tighten
    1688                     for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
    1689                          if (integerInformation[whichColumn[iColumn]]) {
    1690                               double upper = columnUpper2[iColumn];
    1691                               double lower = columnLower2[iColumn];
    1692                               double newUpper = upper;
    1693                               double newLower = lower;
    1694                               double difference = upper - lower;
    1695                               if (lower > -1000.0 && upper < 1000.0) {
    1696                                    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
    1697                                         int iRow = row[j];
    1698                                         double value = element[j];
    1699                                         if (value > 0.0) {
    1700                                              double upWithOut = up[iRow] - value * difference;
    1701                                              if (upWithOut < 0.0) {
    1702                                                   newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
    1703                                              }
    1704                                              double lowWithOut = lo[iRow] + value * difference;
    1705                                              if (lowWithOut > 0.0) {
    1706                                                   newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
    1707                                              }
    1708                                         } else {
    1709                                              double upWithOut = up[iRow] + value * difference;
    1710                                              if (upWithOut < 0.0) {
    1711                                                   newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
    1712                                              }
    1713                                              double lowWithOut = lo[iRow] - value * difference;
    1714                                              if (lowWithOut > 0.0) {
    1715                                                   newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
    1716                                              }
    1717                                         }
    1718                                    }
    1719                                    if (newLower > lower || newUpper < upper) {
    1720                                         if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
    1721                                              newUpper = floor(newUpper);
    1722                                         else
    1723                                              newUpper = floor(newUpper + 0.5);
    1724                                         if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
    1725                                              newLower = ceil(newLower);
    1726                                         else
    1727                                              newLower = ceil(newLower - 0.5);
    1728                                         // change may be too small - check
    1729                                         if (newLower > lower || newUpper < upper) {
    1730                                              if (newUpper >= newLower) {
    1731                                                   // Could also tighten in this
    1732                                                   //printf("%d bounds %g %g tightened to %g %g\n",
    1733                                                   //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
    1734                                                   //     newLower,newUpper);
    1735 #if 1
    1736                                                   columnUpper2[iColumn] = newUpper;
    1737                                                   columnLower2[iColumn] = newLower;
    1738                                                   columnUpper_[whichColumn[iColumn]] = newUpper;
    1739                                                   columnLower_[whichColumn[iColumn]] = newLower;
    1740 #endif
    1741                                                   // and adjust bounds on rows
    1742                                                   newUpper -= upper;
    1743                                                   newLower -= lower;
    1744                                                   for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
    1745                                                        int iRow = row[j];
    1746                                                        double value = element[j];
    1747                                                        if (value > 0.0) {
    1748                                                             up[iRow] += newUpper * value;
    1749                                                             lo[iRow] += newLower * value;
    1750                                                        } else {
    1751                                                             lo[iRow] += newUpper * value;
    1752                                                             up[iRow] += newLower * value;
    1753                                                        }
    1754                                                   }
    1755                                              } else {
    1756                                                   // infeasible
    1757                                                   //printf("%d bounds infeasible %g %g tightened to %g %g\n",
    1758                                                   //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
    1759                                                   //     newLower,newUpper);
    1760 #if 1
    1761                                                   delete small;
    1762                                                   small = NULL;
    1763                                                   break;
    1764 #endif
    1765                                              }
    1766                                         }
    1767                                    }
    1768                               }
    1769                          }
    1770                     }
    1771                }
    1772           }
    1773      }
     1744                    break;
     1745#endif
     1746                  }
     1747                }
     1748              }
     1749            }
     1750          }
     1751        }
     1752      }
     1753    }
     1754  }
    17741755#if 0
    17751756     if (small) {
     
    17921773#endif
    17931774#ifdef CHECK_STATUS
    1794      {
    1795           int n = 0;
    1796           int i;
    1797           for (i = 0; i < small->numberColumns(); i++)
    1798                if (small->getColumnStatus(i) == ClpSimplex::basic)
    1799                     n++;
    1800           for (i = 0; i < small->numberRows(); i++)
    1801                if (small->getRowStatus(i) == ClpSimplex::basic)
    1802                     n++;
    1803           assert (n == small->numberRows());
    1804      }
    1805 #endif
    1806      return small;
     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;
    18071788}
    18081789/* After very cursory presolve.
    18091790   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
    18101791*/
    1811 void
    1812 ClpSimplexOther::afterCrunch(const ClpSimplex & small,
    1813                              const int * whichRow,
    1814                              const int * whichColumn, int nBound)
     1792void ClpSimplexOther::afterCrunch(const ClpSimplex &small,
     1793  const int *whichRow,
     1794  const int *whichColumn, int nBound)
    18151795{
    18161796#ifndef NDEBUG
    1817      for (int i = 0; i < small.numberRows(); i++)
    1818           assert (whichRow[i] >= 0 && whichRow[i] < numberRows_);
    1819      for (int i = 0; i < small.numberColumns(); i++)
    1820           assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_);
    1821 #endif
    1822      getbackSolution(small, whichRow, whichColumn);
    1823      // and deal with status for bounds
    1824      const double * element = matrix_->getElements();
    1825      const int * row = matrix_->getIndices();
    1826      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1827      const int * columnLength = matrix_->getVectorLengths();
    1828      double tolerance = primalTolerance();
    1829      double djTolerance = dualTolerance();
    1830      for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) {
    1831           int iRow = whichRow[jRow];
    1832           int iColumn = whichRow[jRow+numberRows_];
    1833           if (getColumnStatus(iColumn) != ClpSimplex::basic) {
    1834                double lower = columnLower_[iColumn];
    1835                double upper = columnUpper_[iColumn];
    1836                double value = columnActivity_[iColumn];
    1837                double djValue = reducedCost_[iColumn];
    1838                dual_[iRow] = 0.0;
    1839                if (upper > lower) {
    1840                     if (value < lower + tolerance && djValue > -djTolerance) {
    1841                          setColumnStatus(iColumn, ClpSimplex::atLowerBound);
    1842                          setRowStatus(iRow, ClpSimplex::basic);
    1843                     } else if (value > upper - tolerance && djValue < djTolerance) {
    1844                          setColumnStatus(iColumn, ClpSimplex::atUpperBound);
    1845                          setRowStatus(iRow, ClpSimplex::basic);
    1846                     } else {
    1847                          // has to be basic
    1848                          setColumnStatus(iColumn, ClpSimplex::basic);
    1849                          reducedCost_[iColumn] = 0.0;
    1850                          double value = 0.0;
    1851                          for (CoinBigIndex j = columnStart[iColumn];
    1852                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1853                               if (iRow == row[j]) {
    1854                                    value = element[j];
    1855                                    break;
    1856                               }
    1857                          }
    1858                          dual_[iRow] = djValue / value;
    1859                          if (rowUpper_[iRow] > rowLower_[iRow]) {
    1860                               if (fabs(rowActivity_[iRow] - rowLower_[iRow]) <
    1861                                         fabs(rowActivity_[iRow] - rowUpper_[iRow]))
    1862                                    setRowStatus(iRow, ClpSimplex::atLowerBound);
    1863                               else
    1864                                    setRowStatus(iRow, ClpSimplex::atUpperBound);
    1865                          } else {
    1866                               setRowStatus(iRow, ClpSimplex::isFixed);
    1867                          }
    1868                     }
    1869                } else {
    1870                     // row can always be basic
    1871                     setRowStatus(iRow, ClpSimplex::basic);
    1872                }
     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);
    18731844          } else {
    1874                // row can always be basic
    1875                setRowStatus(iRow, ClpSimplex::basic);
    1876           }
    1877      }
    1878      //#ifndef NDEBUG
     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
    18791858#if 0
    18801859     if  (small.status() == 0) {
     
    18931872/* Tightens integer bounds - returns number tightened or -1 if infeasible
    18941873 */
    1895 int
    1896 ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
     1874int ClpSimplexOther::tightenIntegerBounds(double *rhsSpace)
    18971875{
    1898      // See if we can tighten any bounds
    1899      // use rhs for upper and small duals for lower
    1900      double * up = rhsSpace;
    1901      double * lo = dual_;
    1902      const double * element = matrix_->getElements();
    1903      const int * row = matrix_->getIndices();
    1904      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    1905      const int * columnLength = matrix_->getVectorLengths();
    1906      CoinZeroN(lo, numberRows_);
    1907      CoinZeroN(up, numberRows_);
    1908      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1909           double upper = columnUpper_[iColumn];
    1910           double lower = columnLower_[iColumn];
    1911           //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
     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) {
    19121947          for (CoinBigIndex j = columnStart[iColumn];
    1913                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1914                int iRow = row[j];
    1915                double value = element[j];
    1916                if (value > 0.0) {
    1917                     if (upper < 1.0e20)
    1918                          up[iRow] += upper * value;
    1919                     else
    1920                          up[iRow] = COIN_DBL_MAX;
    1921                     if (lower > -1.0e20)
    1922                          lo[iRow] += lower * value;
    1923                     else
    1924                          lo[iRow] = -COIN_DBL_MAX;
    1925                } else {
    1926                     if (upper < 1.0e20)
    1927                          lo[iRow] += upper * value;
    1928                     else
    1929                          lo[iRow] = -COIN_DBL_MAX;
    1930                     if (lower > -1.0e20)
    1931                          up[iRow] += lower * value;
    1932                     else
    1933                          up[iRow] = COIN_DBL_MAX;
    1934                }
    1935           }
    1936      }
    1937      bool feasible = true;
    1938      // make safer
    1939      double tolerance = primalTolerance();
    1940      for (int iRow = 0; iRow < numberRows_; iRow++) {
    1941           double lower = lo[iRow];
    1942           if (lower > rowUpper_[iRow] + tolerance) {
    1943                feasible = false;
    1944                break;
    1945           } else {
    1946                lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance;
    1947           }
    1948           double upper = up[iRow];
    1949           if (upper < rowLower_[iRow] - tolerance) {
    1950                feasible = false;
    1951                break;
    1952           } else {
    1953                up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance;
    1954           }
    1955      }
    1956      int numberTightened = 0;
    1957      if (!feasible) {
    1958           return -1;
    1959      } else if (integerType_) {
    1960           // and tighten
    1961           for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1962                if (integerType_[iColumn]) {
    1963                     double upper = columnUpper_[iColumn];
    1964                     double lower = columnLower_[iColumn];
    1965                     double newUpper = upper;
    1966                     double newLower = lower;
    1967                     double difference = upper - lower;
    1968                     if (lower > -1000.0 && upper < 1000.0) {
    1969                          for (CoinBigIndex j = columnStart[iColumn];
    1970                                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    1971                               int iRow = row[j];
    1972                               double value = element[j];
    1973                               if (value > 0.0) {
    1974                                    double upWithOut = up[iRow] - value * difference;
    1975                                    if (upWithOut < 0.0) {
    1976                                         newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
    1977                                    }
    1978                                    double lowWithOut = lo[iRow] + value * difference;
    1979                                    if (lowWithOut > 0.0) {
    1980                                         newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
    1981                                    }
    1982                               } else {
    1983                                    double upWithOut = up[iRow] + value * difference;
    1984                                    if (upWithOut < 0.0) {
    1985                                         newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
    1986                                    }
    1987                                    double lowWithOut = lo[iRow] - value * difference;
    1988                                    if (lowWithOut > 0.0) {
    1989                                         newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
    1990                                    }
    1991                               }
    1992                          }
    1993                          if (newLower > lower || newUpper < upper) {
    1994                               if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
    1995                                    newUpper = floor(newUpper);
    1996                               else
    1997                                    newUpper = floor(newUpper + 0.5);
    1998                               if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
    1999                                    newLower = ceil(newLower);
    2000                               else
    2001                                    newLower = ceil(newLower - 0.5);
    2002                               // change may be too small - check
    2003                               if (newLower > lower || newUpper < upper) {
    2004                                    if (newUpper >= newLower) {
    2005                                         numberTightened++;
    2006                                         //printf("%d bounds %g %g tightened to %g %g\n",
    2007                                         //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
    2008                                         //     newLower,newUpper);
    2009                                         columnUpper_[iColumn] = newUpper;
    2010                                         columnLower_[iColumn] = newLower;
    2011                                         // and adjust bounds on rows
    2012                                         newUpper -= upper;
    2013                                         newLower -= lower;
    2014                                         for (CoinBigIndex j = columnStart[iColumn];
    2015                                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2016                                              int iRow = row[j];
    2017                                              double value = element[j];
    2018                                              if (value > 0.0) {
    2019                                                   up[iRow] += newUpper * value;
    2020                                                   lo[iRow] += newLower * value;
    2021                                              } else {
    2022                                                   lo[iRow] += newUpper * value;
    2023                                                   up[iRow] += newLower * value;
    2024                                              }
    2025                                         }
    2026                                    } else {
    2027                                         // infeasible
    2028                                         //printf("%d bounds infeasible %g %g tightened to %g %g\n",
    2029                                         //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
    2030                                         //     newLower,newUpper);
    2031                                         return -1;
    2032                                    }
    2033                               }
    2034                          }
    2035                     }
    2036                }
    2037           }
    2038      }
    2039      return numberTightened;
     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;
    20402018}
    20412019/* Parametrics
     
    20522030   On exit endingTheta is maximum reached (can be used for next startingTheta)
    20532031*/
    2054 int
    2055 ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement,
    2056                              const double * lowerChangeBound, const double * upperChangeBound,
    2057                              const double * lowerChangeRhs, const double * upperChangeRhs,
    2058                              const double * changeObjective)
     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)
    20592036{
    2060      bool needToDoSomething = true;
    2061      bool canTryQuick = (reportIncrement) ? true : false;
    2062      // Save copy of model
    2063      ClpSimplex copyModel = *this;
    2064      int savePerturbation = perturbation_;
    2065      perturbation_ = 102; // switch off
    2066      while (needToDoSomething) {
    2067           needToDoSomething = false;
    2068           algorithm_ = -1;
     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;
    20692046
    2070           // save data
    2071           ClpDataSave data = saveData();
    2072           // Dantzig
    2073           ClpDualRowPivot * savePivot = dualRowPivot_;
    2074           dualRowPivot_ = new ClpDualRowDantzig();
    2075           dualRowPivot_->setModel(this);
    2076           int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
    2077           int iRow, iColumn;
    2078           double * chgUpper = NULL;
    2079           double * chgLower = NULL;
    2080           double * chgObjective = NULL;
     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;
    20812058
     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);
    20822199
    2083           if (!returnCode) {
    2084                // Find theta when bounds will cross over and create arrays
    2085                int numberTotal = numberRows_ + numberColumns_;
    2086                chgLower = new double[numberTotal];
    2087                memset(chgLower, 0, numberTotal * sizeof(double));
    2088                chgUpper = new double[numberTotal];
    2089                memset(chgUpper, 0, numberTotal * sizeof(double));
    2090                chgObjective = new double[numberTotal];
    2091                memset(chgObjective, 0, numberTotal * sizeof(double));
    2092                assert (!rowScale_);
    2093                double maxTheta = 1.0e50;
    2094                if (lowerChangeRhs || upperChangeRhs) {
    2095                     for (iRow = 0; iRow < numberRows_; iRow++) {
    2096                          double lower = rowLower_[iRow];
    2097                          double upper = rowUpper_[iRow];
    2098                          if (lower > upper) {
    2099                               maxTheta = -1.0;
    2100                               break;
    2101                          }
    2102                          double lowerChange = (lowerChangeRhs) ? lowerChangeRhs[iRow] : 0.0;
    2103                          double upperChange = (upperChangeRhs) ? upperChangeRhs[iRow] : 0.0;
    2104                          if (lower > -1.0e20 && upper < 1.0e20) {
    2105                               if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
    2106                                    maxTheta = (upper - lower) / (lowerChange - upperChange);
    2107                               }
    2108                          }
    2109                          if (lower > -1.0e20) {
    2110                               lower_[numberColumns_+iRow] += startingTheta * lowerChange;
    2111                               chgLower[numberColumns_+iRow] = lowerChange;
    2112                          }
    2113                          if (upper < 1.0e20) {
    2114                               upper_[numberColumns_+iRow] += startingTheta * upperChange;
    2115                               chgUpper[numberColumns_+iRow] = upperChange;
    2116                          }
    2117                     }
    2118                }
    2119                if (maxTheta > 0.0) {
    2120                     if (lowerChangeBound || upperChangeBound) {
    2121                          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2122                               double lower = columnLower_[iColumn];
    2123                               double upper = columnUpper_[iColumn];
    2124                               if (lower > upper) {
    2125                                    maxTheta = -1.0;
    2126                                    break;
    2127                               }
    2128                               double lowerChange = (lowerChangeBound) ? lowerChangeBound[iColumn] : 0.0;
    2129                               double upperChange = (upperChangeBound) ? upperChangeBound[iColumn] : 0.0;
    2130                               if (lower > -1.0e20 && upper < 1.0e20) {
    2131                                    if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
    2132                                         maxTheta = (upper - lower) / (lowerChange - upperChange);
    2133                                    }
    2134                               }
    2135                               if (lower > -1.0e20) {
    2136                                    lower_[iColumn] += startingTheta * lowerChange;
    2137                                    chgLower[iColumn] = lowerChange;
    2138                               }
    2139                               if (upper < 1.0e20) {
    2140                                    upper_[iColumn] += startingTheta * upperChange;
    2141                                    chgUpper[iColumn] = upperChange;
    2142                               }
    2143                          }
    2144                     }
    2145                     if (maxTheta == 1.0e50)
    2146                          maxTheta = COIN_DBL_MAX;
    2147                }
    2148                if (maxTheta < 0.0) {
    2149                     // bad ranges or initial
    2150                     returnCode = -1;
    2151                }
    2152                if (maxTheta < endingTheta) {
    2153                  char line[100];
    2154                  sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n",
    2155                          endingTheta,maxTheta);
    2156                  handler_->message(CLP_GENERAL,messages_)
    2157                    << line << CoinMessageEol;
    2158                  endingTheta = maxTheta;
    2159                }
    2160                if (endingTheta < startingTheta) {
    2161                     // bad initial
    2162                     returnCode = -2;
    2163                }
    2164           }
    2165           double saveEndingTheta = endingTheta;
    2166           if (!returnCode) {
    2167                if (changeObjective) {
    2168                     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2169                          chgObjective[iColumn] = changeObjective[iColumn];
    2170                          cost_[iColumn] += startingTheta * changeObjective[iColumn];
    2171                     }
    2172                }
    2173                double * saveDuals = NULL;
    2174                reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
    2175                assert (!problemStatus_);
    2176                for (int i=0;i<numberRows_+numberColumns_;i++)
    2177                  setFakeBound(i, noFake);
    2178                // Now do parametrics
    2179                handler_->message(CLP_PARAMETRICS_STATS, messages_)
    2180                  << startingTheta << objectiveValue() << CoinMessageEol;
    2181                while (!returnCode) {
    2182                     //assert (reportIncrement);
    2183                  parametricsData paramData;
    2184                  paramData.startingTheta=startingTheta;
    2185                  paramData.endingTheta=endingTheta;
    2186                  paramData.maxTheta=COIN_DBL_MAX;
    2187                  paramData.lowerChange = chgLower;
    2188                  paramData.upperChange = chgUpper;
    2189                     returnCode = parametricsLoop(paramData, reportIncrement,
    2190                                                  chgLower, chgUpper, chgObjective, data,
    2191                                                  canTryQuick);
    2192                  startingTheta=paramData.startingTheta;
    2193                  endingTheta=paramData.endingTheta;
    2194                     if (!returnCode) {
    2195                          //double change = endingTheta-startingTheta;
    2196                          startingTheta = endingTheta;
    2197                          endingTheta = saveEndingTheta;
    2198                          //for (int i=0;i<numberTotal;i++) {
    2199                          //lower_[i] += change*chgLower[i];
    2200                          //upper_[i] += change*chgUpper[i];
    2201                          //cost_[i] += change*chgObjective[i];
    2202                          //}
    2203                          handler_->message(CLP_PARAMETRICS_STATS, messages_)
    2204                            << startingTheta << objectiveValue() << CoinMessageEol;
    2205                          if (startingTheta >= endingTheta)
    2206                               break;
    2207                     } else if (returnCode == -1) {
    2208                          // trouble - do external solve
    2209                          needToDoSomething = true;
    2210                     } else if (problemStatus_==1) {
    2211                       // can't move any further
    2212                       if (!canTryQuick) {
    2213                         handler_->message(CLP_PARAMETRICS_STATS, messages_)
    2214                           << endingTheta << objectiveValue() << CoinMessageEol;
    2215                         problemStatus_=0;
    2216                       }
    2217                     } else {
    2218                          abort();
    2219                     }
    2220                }
    2221           }
    2222           reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
    2223 
    2224           delete dualRowPivot_;
    2225           dualRowPivot_ = savePivot;
    2226           // Restore any saved stuff
    2227           restoreData(data);
    2228           if (needToDoSomething) {
    2229                double saveStartingTheta = startingTheta; // known to be feasible
    2230                int cleanedUp = 1;
    2231                while (cleanedUp) {
    2232                     // tweak
    2233                     if (cleanedUp == 1) {
    2234                          if (!reportIncrement)
    2235                               startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
    2236                          else
    2237                               startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
    2238                     } else {
    2239                          // restoring to go slowly
    2240                          startingTheta = saveStartingTheta;
    2241                     }
    2242                     // only works if not scaled
    2243                     int i;
    2244                     const double * obj1 = objective();
    2245                     double * obj2 = copyModel.objective();
    2246                     const double * lower1 = columnLower_;
    2247                     double * lower2 = copyModel.columnLower();
    2248                     const double * upper1 = columnUpper_;
    2249                     double * upper2 = copyModel.columnUpper();
    2250                     for (i = 0; i < numberColumns_; i++) {
    2251                          obj2[i] = obj1[i] + startingTheta * chgObjective[i];
    2252                          lower2[i] = lower1[i] + startingTheta * chgLower[i];
    2253                          upper2[i] = upper1[i] + startingTheta * chgUpper[i];
    2254                     }
    2255                     lower1 = rowLower_;
    2256                     lower2 = copyModel.rowLower();
    2257                     upper1 = rowUpper_;
    2258                     upper2 = copyModel.rowUpper();
    2259                     for (i = 0; i < numberRows_; i++) {
    2260                          lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
    2261                          upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
    2262                     }
    2263                     copyModel.dual();
    2264                     if (copyModel.problemStatus()) {
    2265                       char line[100];
    2266                       sprintf(line,"Can not get to theta of %g\n", startingTheta);
    2267                       handler_->message(CLP_GENERAL,messages_)
    2268                         << line << CoinMessageEol;
    2269                          canTryQuick = false; // do slowly to get exact amount
    2270                          // back to last known good
    2271                          if (cleanedUp == 1)
    2272                               cleanedUp = 2;
    2273                          else
    2274                               abort();
    2275                     } else {
    2276                          // and move stuff back
    2277                          int numberTotal = numberRows_ + numberColumns_;
    2278                          CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
    2279                          CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
    2280                          CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
    2281                          cleanedUp = 0;
    2282                     }
    2283                }
    2284           }
    2285           delete [] chgLower;
    2286           delete [] chgUpper;
    2287           delete [] chgObjective;
    2288      }
    2289      perturbation_ = savePerturbation;
    2290      char line[100];
    2291      sprintf(line,"Ending theta %g\n", endingTheta);
    2292      handler_->message(CLP_GENERAL,messages_)
    2293        << line << CoinMessageEol;
    2294      return problemStatus_;
     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_;
    22952271}
    22962272/* Version of parametrics which reads from file
    22972273   See CbcClpParam.cpp for details of format
    22982274   Returns -2 if unable to open file */
    2299 int
    2300 ClpSimplexOther::parametrics(const char * dataFile)
     2275int ClpSimplexOther::parametrics(const char *dataFile)
    23012276{
    2302   int returnCode=-2;
     2277  int returnCode = -2;
    23032278  FILE *fp = fopen(dataFile, "r");
    23042279  char line[200];
     
    23102285
    23112286  if (!fgets(line, 200, fp)) {
    2312     sprintf(line,"Empty parametrics file %s?",dataFile);
    2313     handler_->message(CLP_GENERAL,messages_)
     2287    sprintf(line, "Empty parametrics file %s?", dataFile);
     2288    handler_->message(CLP_GENERAL, messages_)
    23142289      << line << CoinMessageEol;
    23152290    fclose(fp);
    23162291    return -2;
    23172292  }
    2318   char * pos = line;
    2319   char * put = line;
     2293  char *pos = line;
     2294  char *put = line;
    23202295  while (*pos >= ' ' && *pos != '\n') {
    23212296    if (*pos != ' ' && *pos != '\t') {
    2322       *put = static_cast<char>(tolower(*pos));
     2297      *put = static_cast< char >(tolower(*pos));
    23232298      put++;
    23242299    }
     
    23272302  *put = '\0';
    23282303  pos = line;
    2329   double startTheta=0.0;
    2330   double endTheta=0.0;
    2331   double intervalTheta=COIN_DBL_MAX;
    2332   int detail=0;
     2304  double startTheta = 0.0;
     2305  double endTheta = 0.0;
     2306  double intervalTheta = COIN_DBL_MAX;
     2307  int detail = 0;
    23332308  bool good = true;
    23342309  while (good) {
    2335     good=false;
     2310    good = false;
    23362311    // check ROWS
    2337     char * comma = strchr(pos, ',');
     2312    char *comma = strchr(pos, ',');
    23382313    if (!comma)
    23392314      break;
    23402315    *comma = '\0';
    2341     if (strcmp(pos,"rows"))
     2316    if (strcmp(pos, "rows"))
    23422317      break;
    23432318    *comma = ',';
    2344     pos = comma+1;
     2319    pos = comma + 1;
    23452320    // check lower theta
    23462321    comma = strchr(pos, ',');
     
    23502325    startTheta = atof(pos);
    23512326    *comma = ',';
    2352     pos = comma+1;
     2327    pos = comma + 1;
    23532328    // check upper theta
    23542329    comma = strchr(pos, ',');
    2355     good=true;
     2330    good = true;
    23562331    if (comma)
    23572332      *comma = '\0';
     
    23592334    if (comma) {
    23602335      *comma = ',';
    2361       pos = comma+1;
     2336      pos = comma + 1;
    23622337      comma = strchr(pos, ',');
    23632338      if (comma)
    2364         *comma = '\0';
     2339        *comma = '\0';
    23652340      intervalTheta = atof(pos);
    23662341      if (comma) {
    2367         *comma = ',';
    2368         pos = comma+1;
    2369         comma = strchr(pos, ',');
    2370         if (comma)
    2371           *comma = '\0';
    2372         detail = atoi(pos);
    2373         if (comma)
    2374         *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 = ',';
    23752350      }
    23762351    }
     
    23782353  }
    23792354  if (good) {
    2380     if (startTheta<0.0||
    2381         startTheta>endTheta||
    2382         intervalTheta<0.0)
    2383       good=false;
    2384     if (detail<0||detail>1)
    2385       good=false;
    2386   }
    2387   if (intervalTheta>=endTheta)
    2388     intervalTheta=0.0;
     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;
    23892362  if (!good) {
    2390     sprintf(line,"Odd first line %s on file %s?",line,dataFile);
    2391     handler_->message(CLP_GENERAL,messages_)
     2363    sprintf(line, "Odd first line %s on file %s?", line, dataFile);
     2364    handler_->message(CLP_GENERAL, messages_)
    23922365      << line << CoinMessageEol;
    23932366    fclose(fp);
     
    23952368  }
    23962369  if (!fgets(line, 200, fp)) {
    2397     sprintf(line,"Not enough records on parametrics file %s?",dataFile);
    2398     handler_->message(CLP_GENERAL,messages_)
     2370    sprintf(line, "Not enough records on parametrics file %s?", dataFile);
     2371    handler_->message(CLP_GENERAL, messages_)
    23992372      << line << CoinMessageEol;
    24002373    fclose(fp);
    24012374    return -2;
    24022375  }
    2403   double * lowerRowMove = NULL;
    2404   double * upperRowMove = NULL;
    2405   double * lowerColumnMove = NULL;
    2406   double * upperColumnMove = NULL;
    2407   double * objectiveMove = NULL;
     2376  double *lowerRowMove = NULL;
     2377  double *upperRowMove = NULL;
     2378  double *lowerColumnMove = NULL;
     2379  double *upperColumnMove = NULL;
     2380  double *objectiveMove = NULL;
    24082381  char saveLine[200];
    2409   saveLine[0]='\0';
    2410   std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
    2411   int gotRow[] = { -1, -1, -1, -1, -1};
     2382  saveLine[0] = '\0';
     2383  std::string headingsRow[] = { "name", "number", "lower", "upper", "rhs" };
     2384  int gotRow[] = { -1, -1, -1, -1, -1 };
    24122385  int orderRow[5];
    24132386  assert(sizeof(gotRow) == sizeof(orderRow));
     
    24172390  while (*pos >= ' ' && *pos != '\n') {
    24182391    if (*pos != ' ' && *pos != '\t') {
    2419       *put = static_cast<char>(tolower(*pos));
     2392      *put = static_cast< char >(tolower(*pos));
    24202393      put++;
    24212394    }
     
    24262399  int i;
    24272400  good = true;
    2428   if (strncmp(line,"column",6)) {
     2401  if (strncmp(line, "column", 6)) {
    24292402    while (pos) {
    2430       char * comma = strchr(pos, ',');
     2403      char *comma = strchr(pos, ',');
    24312404      if (comma)
    2432         *comma = '\0';
    2433       for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
    2434         if (headingsRow[i] == pos) {
    2435           if (gotRow[i] < 0) {
    2436             orderRow[nAcross] = i;
    2437             gotRow[i] = nAcross++;
    2438           } else {
    2439             // duplicate
    2440             good = false;
    2441           }
    2442           break;
    2443         }
    2444       }
    2445       if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
    2446         good = false;
     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;
    24472420      if (comma) {
    2448         *comma = ',';
    2449         pos = comma + 1;
     2421        *comma = ',';
     2422        pos = comma + 1;
    24502423      } else {
    2451         break;
     2424        break;
    24522425      }
    24532426    }
     
    24582431    if (gotRow[0] >= 0 && !lengthNames())
    24592432      good = false;
    2460     if (gotRow[4]<0) {
     2433    if (gotRow[4] < 0) {
    24612434      if (gotRow[2] < 0 && gotRow[3] >= 0)
    2462         good = false;
     2435        good = false;
    24632436      else if (gotRow[3] < 0 && gotRow[2] >= 0)
    2464         good = false;
    2465     } else if (gotRow[2]>=0||gotRow[3]>=0) {
     2437        good = false;
     2438    } else if (gotRow[2] >= 0 || gotRow[3] >= 0) {
    24662439      good = false;
    24672440    }
    24682441    if (good) {
    2469       char ** rowNames = new char * [numberRows_];
     2442      char **rowNames = new char *[numberRows_];
    24702443      int iRow;
    24712444      for (iRow = 0; iRow < numberRows_; iRow++) {
    2472         rowNames[iRow] =
    2473           CoinStrdup(rowName(iRow).c_str());
    2474       }
    2475       lowerRowMove = new double [numberRows_];
    2476       memset(lowerRowMove,0,numberRows_*sizeof(double));
    2477       upperRowMove = new double [numberRows_];
    2478       memset(upperRowMove,0,numberRows_*sizeof(double));
     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));
    24792451      int nLine = 0;
    24802452      int nBadLine = 0;
    24812453      int nBadName = 0;
    24822454      while (fgets(line, 200, fp)) {
    2483         if (!strncmp(line, "ENDATA", 6)||
    2484             !strncmp(line, "COLUMN",6))
    2485           break;
    2486         nLine++;
    2487         iRow = -1;
    2488         double upper = 0.0;
    2489         double lower = 0.0;
    2490         char * pos = line;
    2491         char * put = line;
    2492         while (*pos >= ' ' && *pos != '\n') {
    2493           if (*pos != ' ' && *pos != '\t') {
    2494             *put = *pos;
    2495             put++;
    2496           }
    2497           pos++;
    2498         }
    2499         *put = '\0';
    2500         pos = line;
    2501         for (int i = 0; i < nAcross; i++) {
    2502           char * comma = strchr(pos, ',');
    2503           if (comma) {
    2504             *comma = '\0';
    2505           } else if (i < nAcross - 1) {
    2506             nBadLine++;
    2507             break;
    2508           }
    2509           switch (orderRow[i]) {
    2510             // name
    2511           case 0:
    2512             // For large problems this could be slow
    2513             for (iRow = 0; iRow < numberRows_; iRow++) {
    2514               if (!strcmp(rowNames[iRow], pos))
    2515                 break;
    2516             }
    2517             if (iRow == numberRows_)
    2518               iRow = -1;
    2519             break;
    2520             // number
    2521           case 1:
    2522             iRow = atoi(pos);
    2523             if (iRow < 0 || iRow >= numberRows_)
    2524               iRow = -1;
    2525             break;
    2526             // lower
    2527           case 2:
    2528             upper = atof(pos);
    2529             break;
    2530             // upper
    2531           case 3:
    2532             lower = atof(pos);
    2533             break;
    2534             // rhs
    2535           case 4:
    2536             lower = atof(pos);
    2537             upper = lower;
    2538             break;
    2539           }
    2540           if (comma) {
    2541             *comma = ',';
    2542             pos = comma + 1;
    2543           }
    2544         }
    2545         if (iRow >= 0) {
    2546           if (rowLower_[iRow]>-1.0e20)
    2547             lowerRowMove[iRow] = lower;
    2548           else
    2549             lowerRowMove[iRow]=0.0;
    2550           if (rowUpper_[iRow]<1.0e20)
    2551             upperRowMove[iRow] = upper;
    2552           else
    2553             upperRowMove[iRow] = lower;
    2554         } else {
    2555           nBadName++;
    2556           if(saveLine[0]=='\0')
    2557             strcpy(saveLine,line);
    2558         }
    2559       }
    2560       sprintf(line,"%d Row fields and %d records", nAcross, nLine);
    2561       handler_->message(CLP_GENERAL,messages_)
    2562         << line << CoinMessageEol;
     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;
    25632534      if (nBadName) {
    2564         sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
    2565         handler_->message(CLP_GENERAL,messages_)
    2566           << line << CoinMessageEol;
    2567         returnCode=-1;
    2568         good=false;
     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;
    25692540      }
    25702541      for (iRow = 0; iRow < numberRows_; iRow++) {
    2571         free(rowNames[iRow]);
    2572       }
    2573       delete [] rowNames;
     2542        free(rowNames[iRow]);
     2543      }
     2544      delete[] rowNames;
    25742545    } else {
    2575       sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
    2576       handler_->message(CLP_GENERAL,messages_)
    2577         << line << CoinMessageEol;
    2578       returnCode=-1;
    2579       good=false;
    2580     }
    2581   }
    2582   if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) {
     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))) {
    25832554    if (!fgets(line, 200, fp)) {
    2584       sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
    2585       handler_->message(CLP_GENERAL,messages_)
    2586         << line << CoinMessageEol;
     2555      sprintf(line, "Not enough records on parametrics file %s after COLUMNS?", dataFile);
     2556      handler_->message(CLP_GENERAL, messages_)
     2557        << line << CoinMessageEol;
    25872558      fclose(fp);
    25882559      return -2;
    25892560    }
    2590     std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
    2591     saveLine[0]='\0';
    2592     int gotColumn[] = { -1, -1, -1, -1, -1};
     2561    std::string headingsColumn[] = { "name", "number", "lower", "upper", "objective" };
     2562    saveLine[0] = '\0';
     2563    int gotColumn[] = { -1, -1, -1, -1, -1 };
    25932564    int orderColumn[5];
    25942565    assert(sizeof(gotColumn) == sizeof(orderColumn));
     
    25982569    while (*pos >= ' ' && *pos != '\n') {
    25992570      if (*pos != ' ' && *pos != '\t') {
    2600         *put = static_cast<char>(tolower(*pos));
    2601         put++;
     2571        *put = static_cast< char >(tolower(*pos));
     2572        put++;
    26022573      }
    26032574      pos++;
     
    26062577    pos = line;
    26072578    int i;
    2608     if (strncmp(line,"endata",6)&&good) {
     2579    if (strncmp(line, "endata", 6) && good) {
    26092580      while (pos) {
    2610         char * comma = strchr(pos, ',');
    2611         if (comma)
    2612           *comma = '\0';
    2613         for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
    2614           if (headingsColumn[i] == pos) {
    2615             if (gotColumn[i] < 0) {
    2616               orderColumn[nAcross] = i;
    2617               gotColumn[i] = nAcross++;
    2618             } else {
    2619               // duplicate
    2620               good = false;
    2621             }
    2622             break;
    2623           }
    2624         }
    2625         if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
    2626           good = false;
    2627         if (comma) {
    2628           *comma = ',';
    2629           pos = comma + 1;
    2630         } else {
    2631           break;
    2632         }
     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        }
    26332604      }
    26342605      if (gotColumn[0] < 0 && gotColumn[1] < 0)
    2635         good = false;
     2606        good = false;
    26362607      if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
    2637         good = false;
     2608        good = false;
    26382609      if (gotColumn[0] >= 0 && !lengthNames())
    2639         good = false;
     2610        good = false;
    26402611      if (good) {
    2641         char ** columnNames = new char * [numberColumns_];
    2642         int iColumn;
    2643         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2644           columnNames[iColumn] =
    2645             CoinStrdup(columnName(iColumn).c_str());
    2646         }
    2647         lowerColumnMove = new double [numberColumns_];
    2648         memset(lowerColumnMove,0,numberColumns_*sizeof(double));
    2649         upperColumnMove = new double [numberColumns_];
    2650         memset(upperColumnMove,0,numberColumns_*sizeof(double));
    2651         objectiveMove = new double [numberColumns_];
    2652         memset(objectiveMove,0,numberColumns_*sizeof(double));
    2653         int nLine = 0;
    2654         int nBadLine = 0;
    2655         int nBadName = 0;
    2656         while (fgets(line, 200, fp)) {
    2657           if (!strncmp(line, "ENDATA", 6))
    2658             break;
    2659           nLine++;
    2660           iColumn = -1;
    2661           double upper = 0.0;
    2662           double lower = 0.0;
    2663           double obj =0.0;
    2664           char * pos = line;
    2665           char * put = line;
    2666           while (*pos >= ' ' && *pos != '\n') {
    2667             if (*pos != ' ' && *pos != '\t') {
    2668               *put = *pos;
    2669               put++;
    2670             }
    2671             pos++;
    2672           }
    2673           *put = '\0';
    2674           pos = line;
    2675           for (int i = 0; i < nAcross; i++) {
    2676             char * comma = strchr(pos, ',');
    2677             if (comma) {
    2678               *comma = '\0';
    2679             } else if (i < nAcross - 1) {
    2680               nBadLine++;
    2681               break;
    2682             }
    2683             switch (orderColumn[i]) {
    2684               // name
    2685             case 0:
    2686               // For large problems this could be slow
    2687               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2688                 if (!strcmp(columnNames[iColumn], pos))
    2689                   break;
    2690               }
    2691               if (iColumn == numberColumns_)
    2692                 iColumn = -1;
    2693               break;
    2694               // number
    2695             case 1:
    2696               iColumn = atoi(pos);
    2697               if (iColumn < 0 || iColumn >= numberColumns_)
    2698                 iColumn = -1;
    2699               break;
    2700               // lower
    2701             case 2:
    2702               upper = atof(pos);
    2703               break;
    2704               // upper
    2705             case 3:
    2706               lower = atof(pos);
    2707               break;
    2708               // objective
    2709             case 4:
    2710               obj = atof(pos);
    2711               upper = lower;
    2712               break;
    2713             }
    2714             if (comma) {
    2715               *comma = ',';
    2716               pos = comma + 1;
    2717             }
    2718           }
    2719           if (iColumn >= 0) {
    2720             if (columnLower_[iColumn]>-1.0e20)
    2721               lowerColumnMove[iColumn] = lower;
    2722             else
    2723               lowerColumnMove[iColumn]=0.0;
    2724             if (columnUpper_[iColumn]<1.0e20)
    2725               upperColumnMove[iColumn] = upper;
    2726             else
    2727               upperColumnMove[iColumn] = lower;
    2728             objectiveMove[iColumn] = obj;
    2729           } else {
    2730             nBadName++;
    2731             if(saveLine[0]=='\0')
    2732               strcpy(saveLine,line);
    2733           }
    2734         }
    2735         sprintf(line,"%d Column fields and %d records", nAcross, nLine);
    2736         handler_->message(CLP_GENERAL,messages_)
    2737           << line << CoinMessageEol;
    2738         if (nBadName) {
    2739           sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
    2740           handler_->message(CLP_GENERAL,messages_)
    2741             << line << CoinMessageEol;
    2742           returnCode=-1;
    2743           good=false;
    2744         }
    2745         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2746           free(columnNames[iColumn]);
    2747         }
    2748         delete [] columnNames;
     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;
    27492719      } else {
    2750         sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
    2751         handler_->message(CLP_GENERAL,messages_)
    2752           << line << CoinMessageEol;
    2753         returnCode=-1;
    2754         good=false;
    2755       }
    2756     }
    2757   }
    2758   returnCode=-1;
     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;
    27592729  if (good) {
    27602730    // clean arrays
    27612731    if (lowerRowMove) {
    2762       bool empty=true;
    2763       for (int i=0;i<numberRows_;i++) {
    2764         if (lowerRowMove[i]) {
    2765           empty=false;
    2766         break;
    2767         }
     2732      bool empty = true;
     2733      for (int i = 0; i < numberRows_; i++) {
     2734        if (lowerRowMove[i]) {
     2735          empty = false;
     2736          break;
     2737        }
    27682738      }
    27692739      if (empty) {
    2770         delete [] lowerRowMove;
    2771         lowerRowMove=NULL;
     2740        delete[] lowerRowMove;
     2741        lowerRowMove = NULL;
    27722742      }
    27732743    }
    27742744    if (upperRowMove) {
    2775       bool empty=true;
    2776       for (int i=0;i<numberRows_;i++) {
    2777         if (upperRowMove[i]) {
    2778           empty=false;
    2779         break;
    2780         }
     2745      bool empty = true;
     2746      for (int i = 0; i < numberRows_; i++) {
     2747        if (upperRowMove[i]) {
     2748          empty = false;
     2749          break;
     2750        }
    27812751      }
    27822752      if (empty) {
    2783         delete [] upperRowMove;
    2784         upperRowMove=NULL;
     2753        delete[] upperRowMove;
     2754        upperRowMove = NULL;
    27852755      }
    27862756    }
    27872757    if (lowerColumnMove) {
    2788       bool empty=true;
    2789       for (int i=0;i<numberColumns_;i++) {
    2790         if (lowerColumnMove[i]) {
    2791           empty=false;
    2792         break;
    2793         }
     2758      bool empty = true;
     2759      for (int i = 0; i < numberColumns_; i++) {
     2760        if (lowerColumnMove[i]) {
     2761          empty = false;
     2762          break;
     2763        }
    27942764      }
    27952765      if (empty) {
    2796         delete [] lowerColumnMove;
    2797         lowerColumnMove=NULL;
     2766        delete[] lowerColumnMove;
     2767        lowerColumnMove = NULL;
    27982768      }
    27992769    }
    28002770    if (upperColumnMove) {
    2801       bool empty=true;
    2802       for (int i=0;i<numberColumns_;i++) {
    2803         if (upperColumnMove[i]) {
    2804           empty=false;
    2805         break;
    2806         }
     2771      bool empty = true;
     2772      for (int i = 0; i < numberColumns_; i++) {
     2773        if (upperColumnMove[i]) {
     2774          empty = false;
     2775          break;
     2776        }
    28072777      }
    28082778      if (empty) {
    2809         delete [] upperColumnMove;
    2810         upperColumnMove=NULL;
     2779        delete[] upperColumnMove;
     2780        upperColumnMove = NULL;
    28112781      }
    28122782    }
    28132783    if (objectiveMove) {
    2814       bool empty=true;
    2815       for (int i=0;i<numberColumns_;i++) {
    2816         if (objectiveMove[i]) {
    2817           empty=false;
    2818         break;
    2819         }
     2784      bool empty = true;
     2785      for (int i = 0; i < numberColumns_; i++) {
     2786        if (objectiveMove[i]) {
     2787          empty = false;
     2788          break;
     2789        }
    28202790      }
    28212791      if (empty) {
    2822         delete [] objectiveMove;
    2823         objectiveMove=NULL;
     2792        delete[] objectiveMove;
     2793        objectiveMove = NULL;
    28242794      }
    28252795    }
     
    28272797    scalingFlag_ = 0;
    28282798    int saveLogLevel = handler_->logLevel();
    2829     if (detail>0&&!intervalTheta)
     2799    if (detail > 0 && !intervalTheta)
    28302800      handler_->setLogLevel(3);
    28312801    else
    28322802      handler_->setLogLevel(1);
    2833     returnCode = parametrics(startTheta,endTheta,intervalTheta,
    2834                              lowerColumnMove,upperColumnMove,
    2835                              lowerRowMove,upperRowMove,
    2836                              objectiveMove);
     2803    returnCode = parametrics(startTheta, endTheta, intervalTheta,
     2804      lowerColumnMove, upperColumnMove,
     2805      lowerRowMove, upperRowMove,
     2806      objectiveMove);
    28372807    scalingFlag_ = saveScaling;
    28382808    handler_->setLogLevel(saveLogLevel);
    28392809  }
    2840   delete [] lowerRowMove;
    2841   delete [] upperRowMove;
    2842   delete [] lowerColumnMove;
    2843   delete [] upperColumnMove;
    2844   delete [] objectiveMove;
     2810  delete[] lowerRowMove;
     2811  delete[] upperRowMove;
     2812  delete[] lowerColumnMove;
     2813  delete[] upperColumnMove;
     2814  delete[] objectiveMove;
    28452815  fclose(fp);
    28462816  return returnCode;
    28472817}
    2848 int
    2849 ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement,
    2850                                  const double * lowerChange, const double * upperChange,
    2851                                  const double * changeObjective, ClpDataSave & data,
    2852                                  bool canTryQuick)
     2818int ClpSimplexOther::parametricsLoop(parametricsData &paramData, double reportIncrement,
     2819  const double *lowerChange, const double *upperChange,
     2820  const double *changeObjective, ClpDataSave &data,
     2821  bool canTryQuick)
    28532822{
    28542823  double startingTheta = paramData.startingTheta;
    2855   double & endingTheta = paramData.endingTheta;
    2856      // stuff is already at starting
    2857      // For this crude version just try and go to end
    2858      double change = 0.0;
    2859      if (reportIncrement && canTryQuick) {
    2860           endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
    2861           change = endingTheta - startingTheta;
    2862      }
    2863      int numberTotal = numberRows_ + numberColumns_;
    2864      int i;
    2865      for ( i = 0; i < numberTotal; i++) {
    2866           lower_[i] += change * lowerChange[i];
    2867           upper_[i] += change * upperChange[i];
    2868           switch(getStatus(i)) {
     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)) {
    28692838
    2870           case basic:
    2871           case isFree:
    2872           case superBasic:
    2873                break;
    2874           case isFixed:
    2875           case atUpperBound:
    2876                solution_[i] = upper_[i];
    2877                break;
    2878           case atLowerBound:
    2879                solution_[i] = lower_[i];
    2880                break;
    2881           }
    2882           cost_[i] += change * changeObjective[i];
    2883      }
    2884      problemStatus_ = -1;
     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;
    28852854
    2886      // This says whether to restore things etc
    2887      // startup will have factorized so can skip
    2888      int factorType = 0;
    2889      // Start check for cycles
    2890      progress_.startCheck();
    2891      // Say change made on first iteration
    2892      changeMade_ = 1;
    2893      /*
     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  /*
    28942863       Status of problem:
    28952864       0 - optimal
     
    29012870       -4 - looks infeasible
    29022871     */
    2903      while (problemStatus_ < 0) {
    2904           int iRow, iColumn;
    2905           // clear
    2906           for (iRow = 0; iRow < 4; iRow++) {
    2907                rowArray_[iRow]->clear();
    2908           }
     2872  while (problemStatus_ < 0) {
     2873    int iRow, iColumn;
     2874    // clear
     2875    for (iRow = 0; iRow < 4; iRow++) {
     2876      rowArray_[iRow]->clear();
     2877    }
    29092878
    2910           for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
    2911                columnArray_[iColumn]->clear();
    2912           }
     2879    for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
     2880      columnArray_[iColumn]->clear();
     2881    }
    29132882
    2914           // give matrix (and model costs and bounds a chance to be
    2915           // refreshed (normally null)
    2916           matrix_->refresh(this);
    2917           // may factorize, checks if problem finished
    2918           statusOfProblemInParametrics(factorType, data);
    2919           // Say good factorization
    2920           factorType = 1;
    2921           if (data.sparseThreshold_) {
    2922                // use default at present
    2923                factorization_->sparseThreshold(0);
    2924                factorization_->goSparse();
    2925           }
     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    }
    29262895
    2927           // exit if victory declared
    2928           if (problemStatus_ >= 0 &&
    2929               (canTryQuick || startingTheta>=endingTheta-1.0e-7) )
    2930                break;
     2896    // exit if victory declared
     2897    if (problemStatus_ >= 0 && (canTryQuick || startingTheta >= endingTheta - 1.0e-7))
     2898      break;
    29312899
    2932           // test for maximum iterations
    2933           if (hitMaximumIterations()) {
    2934                problemStatus_ = 3;
    2935                break;
    2936           }
    2937           // Check event
    2938           {
    2939                int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
    2940                if (status >= 0) {
    2941                     problemStatus_ = 5;
    2942                     secondaryStatus_ = ClpEventHandler::endOfFactorization;
    2943                     break;
    2944                }
    2945           }
    2946           // Do iterations
    2947           problemStatus_=-1;
    2948           if (canTryQuick) {
    2949                double * saveDuals = NULL;
    2950                reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
    2951           } else {
    2952                whileIterating(paramData, reportIncrement,
    2953                               changeObjective);
    2954                startingTheta = endingTheta;
    2955           }
    2956      }
    2957      if (!problemStatus_) {
    2958           theta_ = change + startingTheta;
    2959           eventHandler_->event(ClpEventHandler::theta);
    2960           return 0;
    2961      } else if (problemStatus_ == 10) {
    2962           return -1;
    2963      } else {
    2964           return problemStatus_;
    2965      }
     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  }
    29662934}
    29672935/* Parametrics
     
    29732941   On exit endingTheta is maximum reached (can be used for next startingTheta)
    29742942*/
    2975 int
    2976 ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,
    2977                              const double * lowerChangeBound, const double * upperChangeBound,
    2978                              const double * lowerChangeRhs, const double * upperChangeRhs)
     2943int ClpSimplexOther::parametrics(double startingTheta, double &endingTheta,
     2944  const double *lowerChangeBound, const double *upperChangeBound,
     2945  const double *lowerChangeRhs, const double *upperChangeRhs)
    29792946{
    29802947  int savePerturbation = perturbation_;
     
    29852952  int numberDense = factorization_->numberDense();
    29862953  int length = numberRows_ + numberDense + maximumPivots;
    2987   assert (!rowArray_[4]);
    2988   rowArray_[4]=new CoinIndexedVector(length);
    2989   assert (!rowArray_[5]);
    2990   rowArray_[5]=new CoinIndexedVector(length);
     2954  assert(!rowArray_[4]);
     2955  rowArray_[4] = new CoinIndexedVector(length);
     2956  assert(!rowArray_[5]);
     2957  rowArray_[5] = new CoinIndexedVector(length);
    29912958
    29922959  // save data
    29932960  ClpDataSave data = saveData();
    29942961  int numberTotal = numberRows_ + numberColumns_;
    2995   int ratio = static_cast<int>((2*sizeof(int))/sizeof(double));
    2996   assert (ratio==1||ratio==2);
     2962  int ratio = static_cast< int >((2 * sizeof(int)) / sizeof(double));
     2963  assert(ratio == 1 || ratio == 2);
    29972964  // allow for unscaled - even if not needed
    2998   int lengthArrays = 4*numberTotal+(3*numberTotal+2)*ratio+2*numberRows_+1;
    2999   int unscaledChangesOffset=lengthArrays; // need extra copy of change
     2965  int lengthArrays = 4 * numberTotal + (3 * numberTotal + 2) * ratio + 2 * numberRows_ + 1;
     2966  int unscaledChangesOffset = lengthArrays; // need extra copy of change
    30002967  lengthArrays += numberTotal;
    3001  
     2968
    30022969  /*
    30032970    Save information and modify
    30042971  */
    3005   double * saveLower = new double [lengthArrays];
    3006   double * saveUpper = new double [lengthArrays];
    3007   double * lowerCopy = saveLower+2*numberTotal;
    3008   double * upperCopy = saveUpper+2*numberTotal;
    3009   double * lowerChange = saveLower+numberTotal;
    3010   double * upperChange = saveUpper+numberTotal;
    3011   double * lowerGap = saveLower+4*numberTotal;
    3012   double * lowerCoefficient = lowerGap+numberRows_;
    3013   double * upperGap = saveUpper+4*numberTotal;
    3014   double * upperCoefficient = upperGap+numberRows_;
    3015   int * lowerList = (reinterpret_cast<int *>(saveLower+4*numberTotal+2*numberRows_))+2;
    3016   int * upperList = (reinterpret_cast<int *>(saveUpper+4*numberTotal+2*numberRows_))+2;
    3017   int * lowerActive = lowerList+numberTotal+1;
    3018   int * upperActive = upperList+numberTotal+1;
     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;
    30192986  // To mark as odd
    3020   char * markDone = reinterpret_cast<char *>(lowerActive+numberTotal);
     2987  char *markDone = reinterpret_cast< char * >(lowerActive + numberTotal);
    30212988  //memset(markDone,0,numberTotal);
    3022   int * backwardBasic = upperActive+numberTotal;
     2989  int *backwardBasic = upperActive + numberTotal;
    30232990  parametricsData paramData;
    30242991  paramData.lowerChange = lowerChange;
    3025   paramData.lowerList=lowerList;
     2992  paramData.lowerList = lowerList;
    30262993  paramData.upperChange = upperChange;
    3027   paramData.upperList=upperList;
    3028   paramData.markDone=markDone;
    3029   paramData.backwardBasic=backwardBasic;
     2994  paramData.upperList = upperList;
     2995  paramData.markDone = markDone;
     2996  paramData.backwardBasic = backwardBasic;
    30302997  paramData.lowerActive = lowerActive;
    30312998  paramData.lowerGap = lowerGap;
     
    30343001  paramData.upperGap = upperGap;
    30353002  paramData.upperCoefficient = upperCoefficient;
    3036   paramData.unscaledChangesOffset = unscaledChangesOffset-numberTotal;
    3037   paramData.firstIteration=true;
     3003  paramData.unscaledChangesOffset = unscaledChangesOffset - numberTotal;
     3004  paramData.firstIteration = true;
    30383005  // Find theta when bounds will cross over and create arrays
    30393006  memset(lowerChange, 0, numberTotal * sizeof(double));
    30403007  memset(upperChange, 0, numberTotal * sizeof(double));
    30413008  if (lowerChangeBound)
    3042     memcpy(lowerChange,lowerChangeBound,numberColumns_*sizeof(double));
     3009    memcpy(lowerChange, lowerChangeBound, numberColumns_ * sizeof(double));
    30433010  if (upperChangeBound)
    3044     memcpy(upperChange,upperChangeBound,numberColumns_*sizeof(double));
     3011    memcpy(upperChange, upperChangeBound, numberColumns_ * sizeof(double));
    30453012  if (lowerChangeRhs)
    3046     memcpy(lowerChange+numberColumns_,
    3047            lowerChangeRhs,numberRows_*sizeof(double));
     3013    memcpy(lowerChange + numberColumns_,
     3014      lowerChangeRhs, numberRows_ * sizeof(double));
    30483015  if (upperChangeRhs)
    3049     memcpy(upperChange+numberColumns_,
    3050            upperChangeRhs,numberRows_*sizeof(double));
     3016    memcpy(upperChange + numberColumns_,
     3017      upperChangeRhs, numberRows_ * sizeof(double));
    30513018  // clean
    30523019  for (int iRow = 0; iRow < numberRows_; iRow++) {
    30533020    double lower = rowLower_[iRow];
    30543021    double upper = rowUpper_[iRow];
    3055     if (lower<-1.0e30)
    3056       lowerChange[numberColumns_+iRow]=0.0;
    3057     if (upper>1.0e30)
    3058       upperChange[numberColumns_+iRow]=0.0;
     3022    if (lower < -1.0e30)
     3023      lowerChange[numberColumns_ + iRow] = 0.0;
     3024    if (upper > 1.0e30)
     3025      upperChange[numberColumns_ + iRow] = 0.0;
    30593026  }
    30603027  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    30613028    double lower = columnLower_[iColumn];
    30623029    double upper = columnUpper_[iColumn];
    3063     if (lower<-1.0e30)
    3064       lowerChange[iColumn]=0.0;
    3065     if (upper>1.0e30)
    3066       upperChange[iColumn]=0.0;
     3030    if (lower < -1.0e30)
     3031      lowerChange[iColumn] = 0.0;
     3032    if (upper > 1.0e30)
     3033      upperChange[iColumn] = 0.0;
    30673034  }
    30683035  // save unscaled version of changes
    3069   memcpy(saveLower+unscaledChangesOffset,lowerChange,numberTotal*sizeof(double));
    3070   memcpy(saveUpper+unscaledChangesOffset,upperChange,numberTotal*sizeof(double));
    3071   int nLowerChange=0;
    3072   int nUpperChange=0;
    3073   for (int i=0;i<numberColumns_;i++) {
    3074     if (lowerChange[i]) { 
    3075       lowerList[nLowerChange++]=i;
    3076     }
    3077     if (upperChange[i]) { 
    3078       upperList[nUpperChange++]=i;
    3079     }
    3080   }
    3081   lowerList[-2]=nLowerChange;
    3082   upperList[-2]=nUpperChange;
    3083   for (int i=numberColumns_;i<numberTotal;i++) {
    3084     if (lowerChange[i]) { 
    3085       lowerList[nLowerChange++]=i;
    3086     }
    3087     if (upperChange[i]) { 
    3088       upperList[nUpperChange++]=i;
    3089     }
    3090   }
    3091   lowerList[-1]=nLowerChange;
    3092   upperList[-1]=nUpperChange;
    3093   memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double));
    3094   memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double));
    3095   memcpy(lowerCopy+numberColumns_,
    3096          rowLower_,numberRows_*sizeof(double));
    3097   memcpy(upperCopy+numberColumns_,
    3098          rowUpper_,numberRows_*sizeof(double));
     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));
    30993066  {
    31003067    //  extra for unscaled
    3101     double * unscaledCopy;
    3102     unscaledCopy = lowerCopy + numberTotal; 
    3103     memcpy(unscaledCopy,columnLower_,numberColumns_*sizeof(double));
    3104     memcpy(unscaledCopy+numberColumns_,
    3105            rowLower_,numberRows_*sizeof(double));
    3106     unscaledCopy = upperCopy + numberTotal; 
    3107     memcpy(unscaledCopy,columnUpper_,numberColumns_*sizeof(double));
    3108     memcpy(unscaledCopy+numberColumns_,
    3109            rowUpper_,numberRows_*sizeof(double));
    3110   }
    3111   int returnCode=0;
    3112   paramData.startingTheta=startingTheta;
    3113   paramData.endingTheta=endingTheta;
    3114   paramData.maxTheta=endingTheta;
     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;
    31153082  computeRhsEtc(paramData);
    3116   bool swapped=false;
    3117   // Dantzig 
     3083  bool swapped = false;
     3084  // Dantzig
    31183085#define ALL_DANTZIG
    31193086#ifdef ALL_DANTZIG
    3120   ClpDualRowPivot * savePivot = dualRowPivot_;
     3087  ClpDualRowPivot *savePivot = dualRowPivot_;
    31213088  dualRowPivot_ = new ClpDualRowDantzig();
    31223089  dualRowPivot_->setModel(this);
    31233090#else
    3124   ClpDualRowPivot * savePivot = NULL;
     3091  ClpDualRowPivot *savePivot = NULL;
    31253092#endif
    31263093  if (!returnCode) {
    3127     assert (objective_->type()==1);
     3094    assert(objective_->type() == 1);
    31283095    objective_->setType(2); // in case matrix empty
    3129     returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
     3096    returnCode = reinterpret_cast< ClpSimplexDual * >(this)->startupSolve(0, NULL, 0);
    31303097    objective_->setType(1);
    31313098    if (!returnCode) {
    3132       double saveDualBound=dualBound_;
    3133       dualBound_=CoinMax(dualBound_,1.0e15);
    3134       swapped=true;
    3135       double * temp;
    3136       memcpy(saveLower,lower_,numberTotal*sizeof(double));
    3137       temp=saveLower;
    3138       saveLower=lower_;
    3139       lower_=temp;
     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;
    31403107      //columnLowerWork_ = lower_;
    31413108      //rowLowerWork_ = lower_ + numberColumns_;
    3142       memcpy(saveUpper,upper_,numberTotal*sizeof(double));
    3143       temp=saveUpper;
    3144       saveUpper=upper_;
    3145       upper_=temp;
     3109      memcpy(saveUpper, upper_, numberTotal * sizeof(double));
     3110      temp = saveUpper;
     3111      saveUpper = upper_;