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

formatting

File:
1 edited

Legend:

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

    r2383 r2385  
    1111#include <math.h>
    1212
    13 #if SLIM_CLP==2
     13#if SLIM_CLP == 2
    1414#define SLIM_NOIO
    1515#endif
     
    4343//#############################################################################
    4444
    45 ClpSimplex::ClpSimplex (bool emptyMessages) :
    46 
    47      ClpModel(emptyMessages),
    48      bestPossibleImprovement_(0.0),
    49      zeroTolerance_(1.0e-13),
    50      columnPrimalSequence_(-2),
    51      rowPrimalSequence_(-2),
    52      bestObjectiveValue_(-COIN_DBL_MAX),
    53      moreSpecialOptions_(2),
    54      baseIteration_(0),
    55      vectorMode_(0),
    56      primalToleranceToGetOptimal_(-1.0),
    57      largeValue_(1.0e15),
    58      largestPrimalError_(0.0),
    59      largestDualError_(0.0),
    60      alphaAccuracy_(-1.0),
    61      dualBound_(1.0e10),
    62      alpha_(0.0),
    63      theta_(0.0),
    64      lowerIn_(0.0),
    65      valueIn_(0.0),
    66      upperIn_(-COIN_DBL_MAX),
    67      dualIn_(0.0),
    68      lowerOut_(-1),
    69      valueOut_(-1),
    70      upperOut_(-1),
    71      dualOut_(-1),
    72      dualTolerance_(1.0e-7),
    73      primalTolerance_(1.0e-7),
    74      sumDualInfeasibilities_(0.0),
    75      sumPrimalInfeasibilities_(0.0),
    76      infeasibilityCost_(1.0e10),
    77      sumOfRelaxedDualInfeasibilities_(0.0),
    78      sumOfRelaxedPrimalInfeasibilities_(0.0),
    79      acceptablePivot_(1.0e-8),
    80      lower_(NULL),
    81      rowLowerWork_(NULL),
    82      columnLowerWork_(NULL),
    83      upper_(NULL),
    84      rowUpperWork_(NULL),
    85      columnUpperWork_(NULL),
    86      cost_(NULL),
    87      rowObjectiveWork_(NULL),
    88      objectiveWork_(NULL),
    89      sequenceIn_(-1),
    90      directionIn_(-1),
    91      sequenceOut_(-1),
    92      directionOut_(-1),
    93      pivotRow_(-1),
    94      lastGoodIteration_(-100),
    95      dj_(NULL),
    96      rowReducedCost_(NULL),
    97      reducedCostWork_(NULL),
    98      solution_(NULL),
    99      rowActivityWork_(NULL),
    100      columnActivityWork_(NULL),
    101      numberDualInfeasibilities_(0),
    102      numberDualInfeasibilitiesWithoutFree_(0),
    103      numberPrimalInfeasibilities_(100),
    104      numberRefinements_(0),
    105      pivotVariable_(NULL),
    106      factorization_(NULL),
    107      savedSolution_(NULL),
    108      numberTimesOptimal_(0),
    109      disasterArea_(NULL),
    110      changeMade_(1),
    111      algorithm_(0),
    112      forceFactorization_(-1),
    113      perturbation_(100),
    114      nonLinearCost_(NULL),
    115      lastBadIteration_(-999999),
    116      lastFlaggedIteration_(-999999),
    117      numberFake_(0),
    118      numberChanged_(0),
    119      progressFlag_(0),
    120      firstFree_(-1),
    121      numberExtraRows_(0),
    122      maximumBasic_(0),
    123      dontFactorizePivots_(0),
    124      incomingInfeasibility_(1.0),
    125      allowedInfeasibility_(10.0),
    126      automaticScale_(0),
    127      maximumPerturbationSize_(0),
    128      perturbationArray_(NULL),
    129      baseModel_(NULL)
     45ClpSimplex::ClpSimplex(bool emptyMessages)
     46  :
     47
     48  ClpModel(emptyMessages)
     49  , bestPossibleImprovement_(0.0)
     50  , zeroTolerance_(1.0e-13)
     51  , columnPrimalSequence_(-2)
     52  , rowPrimalSequence_(-2)
     53  , bestObjectiveValue_(-COIN_DBL_MAX)
     54  , moreSpecialOptions_(2)
     55  , baseIteration_(0)
     56  , vectorMode_(0)
     57  , primalToleranceToGetOptimal_(-1.0)
     58  , largeValue_(1.0e15)
     59  , largestPrimalError_(0.0)
     60  , largestDualError_(0.0)
     61  , alphaAccuracy_(-1.0)
     62  , dualBound_(1.0e10)
     63  , alpha_(0.0)
     64  , theta_(0.0)
     65  , lowerIn_(0.0)
     66  , valueIn_(0.0)
     67  , upperIn_(-COIN_DBL_MAX)
     68  , dualIn_(0.0)
     69  , lowerOut_(-1)
     70  , valueOut_(-1)
     71  , upperOut_(-1)
     72  , dualOut_(-1)
     73  , dualTolerance_(1.0e-7)
     74  , primalTolerance_(1.0e-7)
     75  , sumDualInfeasibilities_(0.0)
     76  , sumPrimalInfeasibilities_(0.0)
     77  , infeasibilityCost_(1.0e10)
     78  , sumOfRelaxedDualInfeasibilities_(0.0)
     79  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     80  , acceptablePivot_(1.0e-8)
     81  , lower_(NULL)
     82  , rowLowerWork_(NULL)
     83  , columnLowerWork_(NULL)
     84  , upper_(NULL)
     85  , rowUpperWork_(NULL)
     86  , columnUpperWork_(NULL)
     87  , cost_(NULL)
     88  , rowObjectiveWork_(NULL)
     89  , objectiveWork_(NULL)
     90  , sequenceIn_(-1)
     91  , directionIn_(-1)
     92  , sequenceOut_(-1)
     93  , directionOut_(-1)
     94  , pivotRow_(-1)
     95  , lastGoodIteration_(-100)
     96  , dj_(NULL)
     97  , rowReducedCost_(NULL)
     98  , reducedCostWork_(NULL)
     99  , solution_(NULL)
     100  , rowActivityWork_(NULL)
     101  , columnActivityWork_(NULL)
     102  , numberDualInfeasibilities_(0)
     103  , numberDualInfeasibilitiesWithoutFree_(0)
     104  , numberPrimalInfeasibilities_(100)
     105  , numberRefinements_(0)
     106  , pivotVariable_(NULL)
     107  , factorization_(NULL)
     108  , savedSolution_(NULL)
     109  , numberTimesOptimal_(0)
     110  , disasterArea_(NULL)
     111  , changeMade_(1)
     112  , algorithm_(0)
     113  , forceFactorization_(-1)
     114  , perturbation_(100)
     115  , nonLinearCost_(NULL)
     116  , lastBadIteration_(-999999)
     117  , lastFlaggedIteration_(-999999)
     118  , numberFake_(0)
     119  , numberChanged_(0)
     120  , progressFlag_(0)
     121  , firstFree_(-1)
     122  , numberExtraRows_(0)
     123  , maximumBasic_(0)
     124  , dontFactorizePivots_(0)
     125  , incomingInfeasibility_(1.0)
     126  , allowedInfeasibility_(10.0)
     127  , automaticScale_(0)
     128  , maximumPerturbationSize_(0)
     129  , perturbationArray_(NULL)
     130  , baseModel_(NULL)
    130131#ifdef ABC_INHERIT
    131      ,abcSimplex_(NULL),
    132      abcState_(0)
    133 #endif
    134 {
    135      int i;
    136      for (i = 0; i < 6; i++) {
    137           rowArray_[i] = NULL;
    138           columnArray_[i] = NULL;
    139      }
    140      for (i = 0; i < 4; i++) {
    141           spareIntArray_[i] = 0;
    142           spareDoubleArray_[i] = 0.0;
    143      }
    144      saveStatus_ = NULL;
    145      // get an empty factorization so we can set tolerances etc
    146      getEmptyFactorization();
    147      // Say sparse
    148      factorization_->sparseThreshold(1);
    149      // say Steepest pricing
    150      dualRowPivot_ = new ClpDualRowSteepest();
    151      // say Steepest pricing
    152      primalColumnPivot_ = new ClpPrimalColumnSteepest();
    153      solveType_ = 1; // say simplex based life form
    154 
     132  , abcSimplex_(NULL)
     133  , abcState_(0)
     134#endif
     135{
     136  int i;
     137  for (i = 0; i < 6; i++) {
     138    rowArray_[i] = NULL;
     139    columnArray_[i] = NULL;
     140  }
     141  for (i = 0; i < 4; i++) {
     142    spareIntArray_[i] = 0;
     143    spareDoubleArray_[i] = 0.0;
     144  }
     145  saveStatus_ = NULL;
     146  // get an empty factorization so we can set tolerances etc
     147  getEmptyFactorization();
     148  // Say sparse
     149  factorization_->sparseThreshold(1);
     150  // say Steepest pricing
     151  dualRowPivot_ = new ClpDualRowSteepest();
     152  // say Steepest pricing
     153  primalColumnPivot_ = new ClpPrimalColumnSteepest();
     154  solveType_ = 1; // say simplex based life form
    155155}
    156156
    157157// Subproblem constructor
    158 ClpSimplex::ClpSimplex ( const ClpModel * rhs,
    159                          int numberRows, const int * whichRow,
    160                          int numberColumns, const int * whichColumn,
    161                          bool dropNames, bool dropIntegers, bool fixOthers)
    162      : ClpModel(rhs, numberRows, whichRow,
    163                 numberColumns, whichColumn, dropNames, dropIntegers),
    164      bestPossibleImprovement_(0.0),
    165      zeroTolerance_(1.0e-13),
    166      columnPrimalSequence_(-2),
    167      rowPrimalSequence_(-2),
    168      bestObjectiveValue_(-COIN_DBL_MAX),
    169      moreSpecialOptions_(2),
    170      baseIteration_(0),
    171      vectorMode_(0),
    172      primalToleranceToGetOptimal_(-1.0),
    173      largeValue_(1.0e15),
    174      largestPrimalError_(0.0),
    175      largestDualError_(0.0),
    176      alphaAccuracy_(-1.0),
    177      dualBound_(1.0e10),
    178      alpha_(0.0),
    179      theta_(0.0),
    180      lowerIn_(0.0),
    181      valueIn_(0.0),
    182      upperIn_(-COIN_DBL_MAX),
    183      dualIn_(0.0),
    184      lowerOut_(-1),
    185      valueOut_(-1),
    186      upperOut_(-1),
    187      dualOut_(-1),
    188      dualTolerance_(1.0e-7),
    189      primalTolerance_(1.0e-7),
    190      sumDualInfeasibilities_(0.0),
    191      sumPrimalInfeasibilities_(0.0),
    192      infeasibilityCost_(1.0e10),
    193      sumOfRelaxedDualInfeasibilities_(0.0),
    194      sumOfRelaxedPrimalInfeasibilities_(0.0),
    195      acceptablePivot_(1.0e-8),
    196      lower_(NULL),
    197      rowLowerWork_(NULL),
    198      columnLowerWork_(NULL),
    199      upper_(NULL),
    200      rowUpperWork_(NULL),
    201      columnUpperWork_(NULL),
    202      cost_(NULL),
    203      rowObjectiveWork_(NULL),
    204      objectiveWork_(NULL),
    205      sequenceIn_(-1),
    206      directionIn_(-1),
    207      sequenceOut_(-1),
    208      directionOut_(-1),
    209      pivotRow_(-1),
    210      lastGoodIteration_(-100),
    211      dj_(NULL),
    212      rowReducedCost_(NULL),
    213      reducedCostWork_(NULL),
    214      solution_(NULL),
    215      rowActivityWork_(NULL),
    216      columnActivityWork_(NULL),
    217      numberDualInfeasibilities_(0),
    218      numberDualInfeasibilitiesWithoutFree_(0),
    219      numberPrimalInfeasibilities_(100),
    220      numberRefinements_(0),
    221      pivotVariable_(NULL),
    222      factorization_(NULL),
    223      savedSolution_(NULL),
    224      numberTimesOptimal_(0),
    225      disasterArea_(NULL),
    226      changeMade_(1),
    227      algorithm_(0),
    228      forceFactorization_(-1),
    229      perturbation_(100),
    230      nonLinearCost_(NULL),
    231      lastBadIteration_(-999999),
    232      lastFlaggedIteration_(-999999),
    233      numberFake_(0),
    234      numberChanged_(0),
    235      progressFlag_(0),
    236      firstFree_(-1),
    237      numberExtraRows_(0),
    238      maximumBasic_(0),
    239      dontFactorizePivots_(0),
    240      incomingInfeasibility_(1.0),
    241      allowedInfeasibility_(10.0),
    242      automaticScale_(0),
    243      maximumPerturbationSize_(0),
    244      perturbationArray_(NULL),
    245     baseModel_(NULL)
     158ClpSimplex::ClpSimplex(const ClpModel *rhs,
     159  int numberRows, const int *whichRow,
     160  int numberColumns, const int *whichColumn,
     161  bool dropNames, bool dropIntegers, bool fixOthers)
     162  : ClpModel(rhs, numberRows, whichRow,
     163      numberColumns, whichColumn, dropNames, dropIntegers)
     164  , bestPossibleImprovement_(0.0)
     165  , zeroTolerance_(1.0e-13)
     166  , columnPrimalSequence_(-2)
     167  , rowPrimalSequence_(-2)
     168  , bestObjectiveValue_(-COIN_DBL_MAX)
     169  , moreSpecialOptions_(2)
     170  , baseIteration_(0)
     171  , vectorMode_(0)
     172  , primalToleranceToGetOptimal_(-1.0)
     173  , largeValue_(1.0e15)
     174  , largestPrimalError_(0.0)
     175  , largestDualError_(0.0)
     176  , alphaAccuracy_(-1.0)
     177  , dualBound_(1.0e10)
     178  , alpha_(0.0)
     179  , theta_(0.0)
     180  , lowerIn_(0.0)
     181  , valueIn_(0.0)
     182  , upperIn_(-COIN_DBL_MAX)
     183  , dualIn_(0.0)
     184  , lowerOut_(-1)
     185  , valueOut_(-1)
     186  , upperOut_(-1)
     187  , dualOut_(-1)
     188  , dualTolerance_(1.0e-7)
     189  , primalTolerance_(1.0e-7)
     190  , sumDualInfeasibilities_(0.0)
     191  , sumPrimalInfeasibilities_(0.0)
     192  , infeasibilityCost_(1.0e10)
     193  , sumOfRelaxedDualInfeasibilities_(0.0)
     194  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     195  , acceptablePivot_(1.0e-8)
     196  , lower_(NULL)
     197  , rowLowerWork_(NULL)
     198  , columnLowerWork_(NULL)
     199  , upper_(NULL)
     200  , rowUpperWork_(NULL)
     201  , columnUpperWork_(NULL)
     202  , cost_(NULL)
     203  , rowObjectiveWork_(NULL)
     204  , objectiveWork_(NULL)
     205  , sequenceIn_(-1)
     206  , directionIn_(-1)
     207  , sequenceOut_(-1)
     208  , directionOut_(-1)
     209  , pivotRow_(-1)
     210  , lastGoodIteration_(-100)
     211  , dj_(NULL)
     212  , rowReducedCost_(NULL)
     213  , reducedCostWork_(NULL)
     214  , solution_(NULL)
     215  , rowActivityWork_(NULL)
     216  , columnActivityWork_(NULL)
     217  , numberDualInfeasibilities_(0)
     218  , numberDualInfeasibilitiesWithoutFree_(0)
     219  , numberPrimalInfeasibilities_(100)
     220  , numberRefinements_(0)
     221  , pivotVariable_(NULL)
     222  , factorization_(NULL)
     223  , savedSolution_(NULL)
     224  , numberTimesOptimal_(0)
     225  , disasterArea_(NULL)
     226  , changeMade_(1)
     227  , algorithm_(0)
     228  , forceFactorization_(-1)
     229  , perturbation_(100)
     230  , nonLinearCost_(NULL)
     231  , lastBadIteration_(-999999)
     232  , lastFlaggedIteration_(-999999)
     233  , numberFake_(0)
     234  , numberChanged_(0)
     235  , progressFlag_(0)
     236  , firstFree_(-1)
     237  , numberExtraRows_(0)
     238  , maximumBasic_(0)
     239  , dontFactorizePivots_(0)
     240  , incomingInfeasibility_(1.0)
     241  , allowedInfeasibility_(10.0)
     242  , automaticScale_(0)
     243  , maximumPerturbationSize_(0)
     244  , perturbationArray_(NULL)
     245  , baseModel_(NULL)
    246246#ifdef ABC_INHERIT
    247      ,abcSimplex_(NULL),
    248     abcState_(0)
    249 #endif
    250 {
    251      int i;
    252      for (i = 0; i < 6; i++) {
    253           rowArray_[i] = NULL;
    254           columnArray_[i] = NULL;
    255      }
    256      for (i = 0; i < 4; i++) {
    257           spareIntArray_[i] = 0;
    258           spareDoubleArray_[i] = 0.0;
    259      }
    260      saveStatus_ = NULL;
    261      // get an empty factorization so we can set tolerances etc
    262      getEmptyFactorization();
    263      // say Steepest pricing
    264      dualRowPivot_ = new ClpDualRowSteepest();
    265      // say Steepest pricing
    266      primalColumnPivot_ = new ClpPrimalColumnSteepest();
    267      solveType_ = 1; // say simplex based life form
    268      if (fixOthers) {
    269           int numberOtherColumns = rhs->numberColumns();
    270           int numberOtherRows = rhs->numberRows();
    271           double * solution = new double [numberOtherColumns];
    272           CoinZeroN(solution, numberOtherColumns);
    273           int i;
    274           for (i = 0; i < numberColumns; i++) {
    275                int iColumn = whichColumn[i];
    276                if (solution[iColumn])
    277                     fixOthers = false; // duplicates
    278                solution[iColumn] = 1.0;
    279           }
    280           if (fixOthers) {
    281                const double * otherSolution = rhs->primalColumnSolution();
    282                const double * objective = rhs->objective();
    283                double offset = 0.0;
    284                for (i = 0; i < numberOtherColumns; i++) {
    285                     if (solution[i]) {
    286                          solution[i] = 0.0; // in
    287                     } else {
    288                          solution[i] = otherSolution[i];
    289                          offset += objective[i] * otherSolution[i];
    290                     }
    291                }
    292                double * rhsModification = new double [numberOtherRows];
    293                CoinZeroN(rhsModification, numberOtherRows);
    294                rhs->matrix()->times(solution, rhsModification) ;
    295                for ( i = 0; i < numberRows; i++) {
    296                     int iRow = whichRow[i];
    297                     if (rowLower_[i] > -1.0e20)
    298                          rowLower_[i] -= rhsModification[iRow];
    299                     if (rowUpper_[i] < 1.0e20)
    300                          rowUpper_[i] -= rhsModification[iRow];
    301                }
    302                delete [] rhsModification;
    303                setObjectiveOffset(rhs->objectiveOffset() - offset);
    304                // And set objective value to match
    305                setObjectiveValue(rhs->objectiveValue());
    306           }
    307           delete [] solution;
    308      }
     247  , abcSimplex_(NULL)
     248  , abcState_(0)
     249#endif
     250{
     251  int i;
     252  for (i = 0; i < 6; i++) {
     253    rowArray_[i] = NULL;
     254    columnArray_[i] = NULL;
     255  }
     256  for (i = 0; i < 4; i++) {
     257    spareIntArray_[i] = 0;
     258    spareDoubleArray_[i] = 0.0;
     259  }
     260  saveStatus_ = NULL;
     261  // get an empty factorization so we can set tolerances etc
     262  getEmptyFactorization();
     263  // say Steepest pricing
     264  dualRowPivot_ = new ClpDualRowSteepest();
     265  // say Steepest pricing
     266  primalColumnPivot_ = new ClpPrimalColumnSteepest();
     267  solveType_ = 1; // say simplex based life form
     268  if (fixOthers) {
     269    int numberOtherColumns = rhs->numberColumns();
     270    int numberOtherRows = rhs->numberRows();
     271    double *solution = new double[numberOtherColumns];
     272    CoinZeroN(solution, numberOtherColumns);
     273    int i;
     274    for (i = 0; i < numberColumns; i++) {
     275      int iColumn = whichColumn[i];
     276      if (solution[iColumn])
     277        fixOthers = false; // duplicates
     278      solution[iColumn] = 1.0;
     279    }
     280    if (fixOthers) {
     281      const double *otherSolution = rhs->primalColumnSolution();
     282      const double *objective = rhs->objective();
     283      double offset = 0.0;
     284      for (i = 0; i < numberOtherColumns; i++) {
     285        if (solution[i]) {
     286          solution[i] = 0.0; // in
     287        } else {
     288          solution[i] = otherSolution[i];
     289          offset += objective[i] * otherSolution[i];
     290        }
     291      }
     292      double *rhsModification = new double[numberOtherRows];
     293      CoinZeroN(rhsModification, numberOtherRows);
     294      rhs->matrix()->times(solution, rhsModification);
     295      for (i = 0; i < numberRows; i++) {
     296        int iRow = whichRow[i];
     297        if (rowLower_[i] > -1.0e20)
     298          rowLower_[i] -= rhsModification[iRow];
     299        if (rowUpper_[i] < 1.0e20)
     300          rowUpper_[i] -= rhsModification[iRow];
     301      }
     302      delete[] rhsModification;
     303      setObjectiveOffset(rhs->objectiveOffset() - offset);
     304      // And set objective value to match
     305      setObjectiveValue(rhs->objectiveValue());
     306    }
     307    delete[] solution;
     308  }
    309309}
    310310// Subproblem constructor
    311 ClpSimplex::ClpSimplex ( const ClpSimplex * rhs,
    312                          int numberRows, const int * whichRow,
    313                          int numberColumns, const int * whichColumn,
    314                          bool dropNames, bool dropIntegers, bool fixOthers)
    315      : ClpModel(rhs, numberRows, whichRow,
    316                 numberColumns, whichColumn, dropNames, dropIntegers),
    317      bestPossibleImprovement_(0.0),
    318      zeroTolerance_(1.0e-13),
    319      columnPrimalSequence_(-2),
    320      rowPrimalSequence_(-2),
    321      bestObjectiveValue_(-COIN_DBL_MAX),
    322      moreSpecialOptions_(2),
    323      baseIteration_(0),
    324      vectorMode_(0),
    325      primalToleranceToGetOptimal_(-1.0),
    326      largeValue_(1.0e15),
    327      largestPrimalError_(0.0),
    328      largestDualError_(0.0),
    329      alphaAccuracy_(-1.0),
    330      dualBound_(1.0e10),
    331      alpha_(0.0),
    332      theta_(0.0),
    333      lowerIn_(0.0),
    334      valueIn_(0.0),
    335      upperIn_(-COIN_DBL_MAX),
    336      dualIn_(0.0),
    337      lowerOut_(-1),
    338      valueOut_(-1),
    339      upperOut_(-1),
    340      dualOut_(-1),
    341      dualTolerance_(rhs->dualTolerance_),
    342      primalTolerance_(rhs->primalTolerance_),
    343      sumDualInfeasibilities_(0.0),
    344      sumPrimalInfeasibilities_(0.0),
    345      infeasibilityCost_(1.0e10),
    346      sumOfRelaxedDualInfeasibilities_(0.0),
    347      sumOfRelaxedPrimalInfeasibilities_(0.0),
    348      acceptablePivot_(1.0e-8),
    349      lower_(NULL),
    350      rowLowerWork_(NULL),
    351      columnLowerWork_(NULL),
    352      upper_(NULL),
    353      rowUpperWork_(NULL),
    354      columnUpperWork_(NULL),
    355      cost_(NULL),
    356      rowObjectiveWork_(NULL),
    357      objectiveWork_(NULL),
    358      sequenceIn_(-1),
    359      directionIn_(-1),
    360      sequenceOut_(-1),
    361      directionOut_(-1),
    362      pivotRow_(-1),
    363      lastGoodIteration_(-100),
    364      dj_(NULL),
    365      rowReducedCost_(NULL),
    366      reducedCostWork_(NULL),
    367      solution_(NULL),
    368      rowActivityWork_(NULL),
    369      columnActivityWork_(NULL),
    370      numberDualInfeasibilities_(0),
    371      numberDualInfeasibilitiesWithoutFree_(0),
    372      numberPrimalInfeasibilities_(100),
    373      numberRefinements_(0),
    374      pivotVariable_(NULL),
    375      factorization_(NULL),
    376      savedSolution_(NULL),
    377      numberTimesOptimal_(0),
    378      disasterArea_(NULL),
    379      changeMade_(1),
    380      algorithm_(0),
    381      forceFactorization_(-1),
    382      perturbation_(100),
    383      nonLinearCost_(NULL),
    384      lastBadIteration_(-999999),
    385      lastFlaggedIteration_(-999999),
    386      numberFake_(0),
    387      numberChanged_(0),
    388      progressFlag_(0),
    389      firstFree_(-1),
    390      numberExtraRows_(0),
    391      maximumBasic_(0),
    392      dontFactorizePivots_(0),
    393      incomingInfeasibility_(1.0),
    394      allowedInfeasibility_(10.0),
    395      automaticScale_(0),
    396      maximumPerturbationSize_(0),
    397      perturbationArray_(NULL),
    398     baseModel_(NULL)
     311ClpSimplex::ClpSimplex(const ClpSimplex *rhs,
     312  int numberRows, const int *whichRow,
     313  int numberColumns, const int *whichColumn,
     314  bool dropNames, bool dropIntegers, bool fixOthers)
     315  : ClpModel(rhs, numberRows, whichRow,
     316      numberColumns, whichColumn, dropNames, dropIntegers)
     317  , bestPossibleImprovement_(0.0)
     318  , zeroTolerance_(1.0e-13)
     319  , columnPrimalSequence_(-2)
     320  , rowPrimalSequence_(-2)
     321  , bestObjectiveValue_(-COIN_DBL_MAX)
     322  , moreSpecialOptions_(2)
     323  , baseIteration_(0)
     324  , vectorMode_(0)
     325  , primalToleranceToGetOptimal_(-1.0)
     326  , largeValue_(1.0e15)
     327  , largestPrimalError_(0.0)
     328  , largestDualError_(0.0)
     329  , alphaAccuracy_(-1.0)
     330  , dualBound_(1.0e10)
     331  , alpha_(0.0)
     332  , theta_(0.0)
     333  , lowerIn_(0.0)
     334  , valueIn_(0.0)
     335  , upperIn_(-COIN_DBL_MAX)
     336  , dualIn_(0.0)
     337  , lowerOut_(-1)
     338  , valueOut_(-1)
     339  , upperOut_(-1)
     340  , dualOut_(-1)
     341  , dualTolerance_(rhs->dualTolerance_)
     342  , primalTolerance_(rhs->primalTolerance_)
     343  , sumDualInfeasibilities_(0.0)
     344  , sumPrimalInfeasibilities_(0.0)
     345  , infeasibilityCost_(1.0e10)
     346  , sumOfRelaxedDualInfeasibilities_(0.0)
     347  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     348  , acceptablePivot_(1.0e-8)
     349  , lower_(NULL)
     350  , rowLowerWork_(NULL)
     351  , columnLowerWork_(NULL)
     352  , upper_(NULL)
     353  , rowUpperWork_(NULL)
     354  , columnUpperWork_(NULL)
     355  , cost_(NULL)
     356  , rowObjectiveWork_(NULL)
     357  , objectiveWork_(NULL)
     358  , sequenceIn_(-1)
     359  , directionIn_(-1)
     360  , sequenceOut_(-1)
     361  , directionOut_(-1)
     362  , pivotRow_(-1)
     363  , lastGoodIteration_(-100)
     364  , dj_(NULL)
     365  , rowReducedCost_(NULL)
     366  , reducedCostWork_(NULL)
     367  , solution_(NULL)
     368  , rowActivityWork_(NULL)
     369  , columnActivityWork_(NULL)
     370  , numberDualInfeasibilities_(0)
     371  , numberDualInfeasibilitiesWithoutFree_(0)
     372  , numberPrimalInfeasibilities_(100)
     373  , numberRefinements_(0)
     374  , pivotVariable_(NULL)
     375  , factorization_(NULL)
     376  , savedSolution_(NULL)
     377  , numberTimesOptimal_(0)
     378  , disasterArea_(NULL)
     379  , changeMade_(1)
     380  , algorithm_(0)
     381  , forceFactorization_(-1)
     382  , perturbation_(100)
     383  , nonLinearCost_(NULL)
     384  , lastBadIteration_(-999999)
     385  , lastFlaggedIteration_(-999999)
     386  , numberFake_(0)
     387  , numberChanged_(0)
     388  , progressFlag_(0)
     389  , firstFree_(-1)
     390  , numberExtraRows_(0)
     391  , maximumBasic_(0)
     392  , dontFactorizePivots_(0)
     393  , incomingInfeasibility_(1.0)
     394  , allowedInfeasibility_(10.0)
     395  , automaticScale_(0)
     396  , maximumPerturbationSize_(0)
     397  , perturbationArray_(NULL)
     398  , baseModel_(NULL)
    399399#ifdef ABC_INHERIT
    400      ,abcSimplex_(NULL),
    401      abcState_(rhs->abcState_)
    402 #endif
    403 {
    404      int i;
    405      for (i = 0; i < 6; i++) {
    406           rowArray_[i] = NULL;
    407           columnArray_[i] = NULL;
    408      }
    409      for (i = 0; i < 4; i++) {
    410           spareIntArray_[i] = 0;
    411           spareDoubleArray_[i] = 0.0;
    412      }
    413      saveStatus_ = NULL;
    414      factorization_ = new ClpFactorization(*rhs->factorization_, -numberRows_);
    415      //factorization_ = new ClpFactorization(*rhs->factorization_,
    416      //                         rhs->factorization_->goDenseThreshold());
    417      ClpPEDualRowSteepest * pivotDualPE =
    418           dynamic_cast< ClpPEDualRowSteepest*>(rhs->dualRowPivot_);
    419      if (pivotDualPE) {
    420        dualRowPivot_ = new ClpPEDualRowSteepest(pivotDualPE->psi());
    421      } else {
    422        ClpDualRowDantzig * pivot =
    423          dynamic_cast< ClpDualRowDantzig*>(rhs->dualRowPivot_);
    424        // say Steepest pricing
    425        if (!pivot)
    426          dualRowPivot_ = new ClpDualRowSteepest();
    427        else
    428          dualRowPivot_ = new ClpDualRowDantzig();
    429      }
    430      ClpPEPrimalColumnSteepest * pivotPrimalPE =
    431           dynamic_cast< ClpPEPrimalColumnSteepest*>(rhs->primalColumnPivot_);
    432      if (pivotPrimalPE) {
    433        primalColumnPivot_ = new ClpPEPrimalColumnSteepest(pivotPrimalPE->psi());
    434      } else {
    435        // say Steepest pricing
    436        primalColumnPivot_ = new ClpPrimalColumnSteepest();
    437      }
    438      solveType_ = 1; // say simplex based life form
    439      if (fixOthers) {
    440           int numberOtherColumns = rhs->numberColumns();
    441           int numberOtherRows = rhs->numberRows();
    442           double * solution = new double [numberOtherColumns];
    443           CoinZeroN(solution, numberOtherColumns);
    444           int i;
    445           for (i = 0; i < numberColumns; i++) {
    446                int iColumn = whichColumn[i];
    447                if (solution[iColumn])
    448                     fixOthers = false; // duplicates
    449                solution[iColumn] = 1.0;
    450           }
    451           if (fixOthers) {
    452                const double * otherSolution = rhs->primalColumnSolution();
    453                const double * objective = rhs->objective();
    454                double offset = 0.0;
    455                for (i = 0; i < numberOtherColumns; i++) {
    456                     if (solution[i]) {
    457                          solution[i] = 0.0; // in
    458                     } else {
    459                          solution[i] = otherSolution[i];
    460                          offset += objective[i] * otherSolution[i];
    461                     }
    462                }
    463                double * rhsModification = new double [numberOtherRows];
    464                CoinZeroN(rhsModification, numberOtherRows);
    465                rhs->matrix()->times(solution, rhsModification) ;
    466                for ( i = 0; i < numberRows; i++) {
    467                     int iRow = whichRow[i];
    468                     if (rowLower_[i] > -1.0e20)
    469                          rowLower_[i] -= rhsModification[iRow];
    470                     if (rowUpper_[i] < 1.0e20)
    471                          rowUpper_[i] -= rhsModification[iRow];
    472                }
    473                delete [] rhsModification;
    474                setObjectiveOffset(rhs->objectiveOffset() - offset);
    475                // And set objective value to match
    476                setObjectiveValue(rhs->objectiveValue());
    477           }
    478           delete [] solution;
    479      }
    480      if (rhs->maximumPerturbationSize_) {
    481           maximumPerturbationSize_ = 2 * numberColumns;
    482           perturbationArray_ = new double [maximumPerturbationSize_];
    483           for (i = 0; i < numberColumns; i++) {
    484                int iColumn = whichColumn[i];
    485                perturbationArray_[2*i] = rhs->perturbationArray_[2*iColumn];
    486                perturbationArray_[2*i+1] = rhs->perturbationArray_[2*iColumn+1];
    487           }
    488      }
     400  , abcSimplex_(NULL)
     401  , abcState_(rhs->abcState_)
     402#endif
     403{
     404  int i;
     405  for (i = 0; i < 6; i++) {
     406    rowArray_[i] = NULL;
     407    columnArray_[i] = NULL;
     408  }
     409  for (i = 0; i < 4; i++) {
     410    spareIntArray_[i] = 0;
     411    spareDoubleArray_[i] = 0.0;
     412  }
     413  saveStatus_ = NULL;
     414  factorization_ = new ClpFactorization(*rhs->factorization_, -numberRows_);
     415  //factorization_ = new ClpFactorization(*rhs->factorization_,
     416  //                            rhs->factorization_->goDenseThreshold());
     417  ClpPEDualRowSteepest *pivotDualPE = dynamic_cast< ClpPEDualRowSteepest * >(rhs->dualRowPivot_);
     418  if (pivotDualPE) {
     419    dualRowPivot_ = new ClpPEDualRowSteepest(pivotDualPE->psi());
     420  } else {
     421    ClpDualRowDantzig *pivot = dynamic_cast< ClpDualRowDantzig * >(rhs->dualRowPivot_);
     422    // say Steepest pricing
     423    if (!pivot)
     424      dualRowPivot_ = new ClpDualRowSteepest();
     425    else
     426      dualRowPivot_ = new ClpDualRowDantzig();
     427  }
     428  ClpPEPrimalColumnSteepest *pivotPrimalPE = dynamic_cast< ClpPEPrimalColumnSteepest * >(rhs->primalColumnPivot_);
     429  if (pivotPrimalPE) {
     430    primalColumnPivot_ = new ClpPEPrimalColumnSteepest(pivotPrimalPE->psi());
     431  } else {
     432    // say Steepest pricing
     433    primalColumnPivot_ = new ClpPrimalColumnSteepest();
     434  }
     435  solveType_ = 1; // say simplex based life form
     436  if (fixOthers) {
     437    int numberOtherColumns = rhs->numberColumns();
     438    int numberOtherRows = rhs->numberRows();
     439    double *solution = new double[numberOtherColumns];
     440    CoinZeroN(solution, numberOtherColumns);
     441    int i;
     442    for (i = 0; i < numberColumns; i++) {
     443      int iColumn = whichColumn[i];
     444      if (solution[iColumn])
     445        fixOthers = false; // duplicates
     446      solution[iColumn] = 1.0;
     447    }
     448    if (fixOthers) {
     449      const double *otherSolution = rhs->primalColumnSolution();
     450      const double *objective = rhs->objective();
     451      double offset = 0.0;
     452      for (i = 0; i < numberOtherColumns; i++) {
     453        if (solution[i]) {
     454          solution[i] = 0.0; // in
     455        } else {
     456          solution[i] = otherSolution[i];
     457          offset += objective[i] * otherSolution[i];
     458        }
     459      }
     460      double *rhsModification = new double[numberOtherRows];
     461      CoinZeroN(rhsModification, numberOtherRows);
     462      rhs->matrix()->times(solution, rhsModification);
     463      for (i = 0; i < numberRows; i++) {
     464        int iRow = whichRow[i];
     465        if (rowLower_[i] > -1.0e20)
     466          rowLower_[i] -= rhsModification[iRow];
     467        if (rowUpper_[i] < 1.0e20)
     468          rowUpper_[i] -= rhsModification[iRow];
     469      }
     470      delete[] rhsModification;
     471      setObjectiveOffset(rhs->objectiveOffset() - offset);
     472      // And set objective value to match
     473      setObjectiveValue(rhs->objectiveValue());
     474    }
     475    delete[] solution;
     476  }
     477  if (rhs->maximumPerturbationSize_) {
     478    maximumPerturbationSize_ = 2 * numberColumns;
     479    perturbationArray_ = new double[maximumPerturbationSize_];
     480    for (i = 0; i < numberColumns; i++) {
     481      int iColumn = whichColumn[i];
     482      perturbationArray_[2 * i] = rhs->perturbationArray_[2 * iColumn];
     483      perturbationArray_[2 * i + 1] = rhs->perturbationArray_[2 * iColumn + 1];
     484    }
     485  }
    489486}
    490487// Puts solution back small model
    491 void
    492 ClpSimplex::getbackSolution(const ClpSimplex & smallModel, const int * whichRow, const int * whichColumn)
    493 {
    494      setSumDualInfeasibilities(smallModel.sumDualInfeasibilities());
    495      setNumberDualInfeasibilities(smallModel.numberDualInfeasibilities());
    496      setSumPrimalInfeasibilities(smallModel.sumPrimalInfeasibilities());
    497      setNumberPrimalInfeasibilities(smallModel.numberPrimalInfeasibilities());
    498      setNumberIterations(smallModel.numberIterations());
    499      setProblemStatus(smallModel.status());
    500      setObjectiveValue(smallModel.objectiveValue());
    501      const double * solution2 = smallModel.primalColumnSolution();
    502      int i;
    503      int numberRows2 = smallModel.numberRows();
    504      int numberColumns2 = smallModel.numberColumns();
    505      const double * dj2 = smallModel.dualColumnSolution();
    506      for ( i = 0; i < numberColumns2; i++) {
    507           int iColumn = whichColumn[i];
    508           columnActivity_[iColumn] = solution2[i];
    509           reducedCost_[iColumn] = dj2[i];
    510           setStatus(iColumn, smallModel.getStatus(i));
    511      }
    512      const double * dual2 = smallModel.dualRowSolution();
    513      memset(dual_, 0, numberRows_ * sizeof(double));
    514      for (i = 0; i < numberRows2; i++) {
    515           int iRow = whichRow[i];
    516           setRowStatus(iRow, smallModel.getRowStatus(i));
    517           dual_[iRow] = dual2[i];
    518      }
    519      CoinZeroN(rowActivity_, numberRows_);
     488void ClpSimplex::getbackSolution(const ClpSimplex &smallModel, const int *whichRow, const int *whichColumn)
     489{
     490  setSumDualInfeasibilities(smallModel.sumDualInfeasibilities());
     491  setNumberDualInfeasibilities(smallModel.numberDualInfeasibilities());
     492  setSumPrimalInfeasibilities(smallModel.sumPrimalInfeasibilities());
     493  setNumberPrimalInfeasibilities(smallModel.numberPrimalInfeasibilities());
     494  setNumberIterations(smallModel.numberIterations());
     495  setProblemStatus(smallModel.status());
     496  setObjectiveValue(smallModel.objectiveValue());
     497  const double *solution2 = smallModel.primalColumnSolution();
     498  int i;
     499  int numberRows2 = smallModel.numberRows();
     500  int numberColumns2 = smallModel.numberColumns();
     501  const double *dj2 = smallModel.dualColumnSolution();
     502  for (i = 0; i < numberColumns2; i++) {
     503    int iColumn = whichColumn[i];
     504    columnActivity_[iColumn] = solution2[i];
     505    reducedCost_[iColumn] = dj2[i];
     506    setStatus(iColumn, smallModel.getStatus(i));
     507  }
     508  const double *dual2 = smallModel.dualRowSolution();
     509  memset(dual_, 0, numberRows_ * sizeof(double));
     510  for (i = 0; i < numberRows2; i++) {
     511    int iRow = whichRow[i];
     512    setRowStatus(iRow, smallModel.getRowStatus(i));
     513    dual_[iRow] = dual2[i];
     514  }
     515  CoinZeroN(rowActivity_, numberRows_);
    520516#if 0
    521517     if (!problemStatus_) {
     
    534530     }
    535531#endif
    536      matrix()->times(columnActivity_, rowActivity_) ;
     532  matrix()->times(columnActivity_, rowActivity_);
    537533}
    538534
    539535//-----------------------------------------------------------------------------
    540536
    541 ClpSimplex::~ClpSimplex ()
    542 {
    543      setPersistenceFlag(0);
    544      gutsOfDelete(0);
    545      delete nonLinearCost_;
     537ClpSimplex::~ClpSimplex()
     538{
     539  setPersistenceFlag(0);
     540  gutsOfDelete(0);
     541  delete nonLinearCost_;
    546542}
    547543//#############################################################################
    548 void ClpSimplex::setLargeValue( double value)
    549 {
    550      if (value > 0.0 && value < COIN_DBL_MAX)
    551           largeValue_ = value;
    552 }
    553 int
    554 ClpSimplex::gutsOfSolution ( double * givenDuals,
    555                              const double * givenPrimals,
    556                              bool valuesPass)
    557 {
    558 
    559 
    560      // if values pass, save values of basic variables
    561      double * save = NULL;
    562      double oldValue = 0.0;
    563      double oldLargestPrimalError=largestPrimalError_;
    564      double oldLargestDualError=largestDualError_;
    565      if (valuesPass) {
    566           assert(algorithm_ > 0); // only primal at present
    567           assert(nonLinearCost_);
    568           int iRow;
    569           checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    570           // get correct bounds on all variables
    571           nonLinearCost_->checkInfeasibilities(primalTolerance_);
    572           oldValue = nonLinearCost_->largestInfeasibility();
    573           save = new double[numberRows_];
    574           for (iRow = 0; iRow < numberRows_; iRow++) {
    575                int iPivot = pivotVariable_[iRow];
    576                save[iRow] = solution_[iPivot];
    577           }
    578      }
    579      // do work
    580      computePrimals(rowActivityWork_, columnActivityWork_);
    581      // If necessary - override results
    582      if (givenPrimals) {
    583           CoinMemcpyN(givenPrimals, numberColumns_, columnActivityWork_);
    584           memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
    585           times(-1.0, columnActivityWork_, rowActivityWork_);
    586      }
    587      double objectiveModification = 0.0;
    588      if (algorithm_ > 0 && nonLinearCost_ != NULL) {
     544void ClpSimplex::setLargeValue(double value)
     545{
     546  if (value > 0.0 && value < COIN_DBL_MAX)
     547    largeValue_ = value;
     548}
     549int ClpSimplex::gutsOfSolution(double *givenDuals,
     550  const double *givenPrimals,
     551  bool valuesPass)
     552{
     553
     554  // if values pass, save values of basic variables
     555  double *save = NULL;
     556  double oldValue = 0.0;
     557  double oldLargestPrimalError = largestPrimalError_;
     558  double oldLargestDualError = largestDualError_;
     559  if (valuesPass) {
     560    assert(algorithm_ > 0); // only primal at present
     561    assert(nonLinearCost_);
     562    int iRow;
     563    checkPrimalSolution(rowActivityWork_, columnActivityWork_);
     564    // get correct bounds on all variables
     565    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     566    oldValue = nonLinearCost_->largestInfeasibility();
     567    save = new double[numberRows_];
     568    for (iRow = 0; iRow < numberRows_; iRow++) {
     569      int iPivot = pivotVariable_[iRow];
     570      save[iRow] = solution_[iPivot];
     571    }
     572  }
     573  // do work
     574  computePrimals(rowActivityWork_, columnActivityWork_);
     575  // If necessary - override results
     576  if (givenPrimals) {
     577    CoinMemcpyN(givenPrimals, numberColumns_, columnActivityWork_);
     578    memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
     579    times(-1.0, columnActivityWork_, rowActivityWork_);
     580  }
     581  double objectiveModification = 0.0;
     582  if (algorithm_ > 0 && nonLinearCost_ != NULL) {
    589583#ifdef CLP_USER_DRIVEN
    590           eventHandler_->eventWithInfo(ClpEventHandler::goodFactorization,NULL);
    591 #endif
    592           // primal algorithm
    593           // get correct bounds on all variables
    594           // If  4 bit set - Force outgoing variables to exact bound (primal)
    595           if ((specialOptions_ & 4) == 0)
    596                nonLinearCost_->checkInfeasibilities(primalTolerance_);
    597           else
    598                nonLinearCost_->checkInfeasibilities(0.0);
    599           objectiveModification += nonLinearCost_->changeInCost();
    600           if (nonLinearCost_->numberInfeasibilities())
    601                if (handler_->detail(CLP_SIMPLEX_NONLINEAR, messages_) < 100) {
    602                     handler_->message(CLP_SIMPLEX_NONLINEAR, messages_)
    603                               << nonLinearCost_->changeInCost()
    604                               << nonLinearCost_->numberInfeasibilities()
    605                               << CoinMessageEol;
    606                }
    607      }
    608      if (valuesPass) {
    609           double badInfeasibility = nonLinearCost_->largestInfeasibility();
     584    eventHandler_->eventWithInfo(ClpEventHandler::goodFactorization, NULL);
     585#endif
     586    // primal algorithm
     587    // get correct bounds on all variables
     588    // If  4 bit set - Force outgoing variables to exact bound (primal)
     589    if ((specialOptions_ & 4) == 0)
     590      nonLinearCost_->checkInfeasibilities(primalTolerance_);
     591    else
     592      nonLinearCost_->checkInfeasibilities(0.0);
     593    objectiveModification += nonLinearCost_->changeInCost();
     594    if (nonLinearCost_->numberInfeasibilities())
     595      if (handler_->detail(CLP_SIMPLEX_NONLINEAR, messages_) < 100) {
     596        handler_->message(CLP_SIMPLEX_NONLINEAR, messages_)
     597          << nonLinearCost_->changeInCost()
     598          << nonLinearCost_->numberInfeasibilities()
     599          << CoinMessageEol;
     600      }
     601  }
     602  if (valuesPass) {
     603    double badInfeasibility = nonLinearCost_->largestInfeasibility();
    610604#ifdef CLP_DEBUG
    611           std::cout << "Largest given infeasibility " << oldValue
    612                     << " now " << nonLinearCost_->largestInfeasibility() << std::endl;
    613 #endif
    614           int numberOut = 0;
    615           // But may be very large rhs etc
    616           double useError = CoinMin(largestPrimalError_,
    617                                     1.0e5 / maximumAbsElement(solution_, numberRows_ + numberColumns_));
    618           if ((oldValue < incomingInfeasibility_ || badInfeasibility >
    619                     (CoinMax(10.0 * allowedInfeasibility_, 100.0 * oldValue)))
    620                     && (badInfeasibility > CoinMax(incomingInfeasibility_, allowedInfeasibility_) ||
    621                         useError > 1.0e-3)) {
    622                if (algorithm_>1) {
    623                  // nonlinear
    624                  //printf("Original largest infeas %g, now %g, primalError %g\n",
    625                  //     oldValue,nonLinearCost_->largestInfeasibility(),
    626                  //     largestPrimalError_);
    627                  //printf("going to all slack\n");
    628                  allSlackBasis(true);
    629                  CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
    630                  return 1;
    631                }
    632                //printf("Original largest infeas %g, now %g, primalError %g\n",
    633                //     oldValue,nonLinearCost_->largestInfeasibility(),
    634                //     largestPrimalError_);
    635                // throw out up to 1000 structurals
    636                int maxOut = (allowedInfeasibility_==10.0) ? 1000 : 100;
    637                int iRow;
    638                int * sort = new int[numberRows_];
    639                // first put back solution and store difference
    640                for (iRow = 0; iRow < numberRows_; iRow++) {
    641                     int iPivot = pivotVariable_[iRow];
    642                     double difference = fabs(solution_[iPivot] - save[iRow]);
    643                     solution_[iPivot] = save[iRow];
    644                     save[iRow] = difference;
    645                }
    646                int numberBasic = 0;
    647                for (iRow = 0; iRow < numberRows_; iRow++) {
    648                     int iPivot = pivotVariable_[iRow];
    649 
    650                     if (iPivot < numberColumns_) {
    651                          // column
    652                          double difference = save[iRow];
    653                          if (difference > 1.0e-4) {
    654                               sort[numberOut] = iRow;
    655                               save[numberOut++] = -difference;
    656                               if (getStatus(iPivot) == basic)
    657                                    numberBasic++;
    658                          }
    659                     }
    660                }
    661                if (!numberBasic) {
    662                     //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
     605    std::cout << "Largest given infeasibility " << oldValue
     606              << " now " << nonLinearCost_->largestInfeasibility() << std::endl;
     607#endif
     608    int numberOut = 0;
     609    // But may be very large rhs etc
     610    double useError = CoinMin(largestPrimalError_,
     611      1.0e5 / maximumAbsElement(solution_, numberRows_ + numberColumns_));
     612    if ((oldValue < incomingInfeasibility_ || badInfeasibility > (CoinMax(10.0 * allowedInfeasibility_, 100.0 * oldValue)))
     613      && (badInfeasibility > CoinMax(incomingInfeasibility_, allowedInfeasibility_) || useError > 1.0e-3)) {
     614      if (algorithm_ > 1) {
     615        // nonlinear
     616        //printf("Original largest infeas %g, now %g, primalError %g\n",
     617        //      oldValue,nonLinearCost_->largestInfeasibility(),
     618        //      largestPrimalError_);
     619        //printf("going to all slack\n");
     620        allSlackBasis(true);
     621        CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
     622        return 1;
     623      }
     624      //printf("Original largest infeas %g, now %g, primalError %g\n",
     625      //     oldValue,nonLinearCost_->largestInfeasibility(),
     626      //     largestPrimalError_);
     627      // throw out up to 1000 structurals
     628      int maxOut = (allowedInfeasibility_ == 10.0) ? 1000 : 100;
     629      int iRow;
     630      int *sort = new int[numberRows_];
     631      // first put back solution and store difference
     632      for (iRow = 0; iRow < numberRows_; iRow++) {
     633        int iPivot = pivotVariable_[iRow];
     634        double difference = fabs(solution_[iPivot] - save[iRow]);
     635        solution_[iPivot] = save[iRow];
     636        save[iRow] = difference;
     637      }
     638      int numberBasic = 0;
     639      for (iRow = 0; iRow < numberRows_; iRow++) {
     640        int iPivot = pivotVariable_[iRow];
     641
     642        if (iPivot < numberColumns_) {
     643          // column
     644          double difference = save[iRow];
     645          if (difference > 1.0e-4) {
     646            sort[numberOut] = iRow;
     647            save[numberOut++] = -difference;
     648            if (getStatus(iPivot) == basic)
     649              numberBasic++;
     650          }
     651        }
     652      }
     653      if (!numberBasic) {
     654        //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
    663655#if 0
    664656                    allSlackBasis(true);
    665657                    CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
    666658#else
    667                     // allow
    668                     numberOut = 0;
    669 #endif
    670                }
    671                CoinSort_2(save, save + numberOut, sort);
    672                numberOut = CoinMin(maxOut, numberOut);
    673                for (iRow = 0; iRow < numberOut; iRow++) {
    674                     int jRow = sort[iRow];
    675                     int iColumn = pivotVariable_[jRow];
    676                     setColumnStatus(iColumn, superBasic);
    677                     setRowStatus(jRow, basic);
    678                     pivotVariable_[jRow] = jRow + numberColumns_;
    679                     if (fabs(solution_[iColumn]) > 1.0e10) {
    680                          if (upper_[iColumn] < 0.0) {
    681                               solution_[iColumn] = upper_[iColumn];
    682                          } else if (lower_[iColumn] > 0.0) {
    683                               solution_[iColumn] = lower_[iColumn];
    684                          } else {
    685                               solution_[iColumn] = 0.0;
    686                          }
    687                     }
    688                }
    689                delete [] sort;
    690           }
    691           delete [] save;
    692           if (numberOut)
    693                return numberOut;
    694      }
    695      if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
    696           //printf("trying feas pump\n");
    697           const char * integerType = integerInformation();
    698           assert (integerType);
    699           assert (perturbationArray_);
    700           CoinZeroN(cost_, numberRows_ + numberColumns_);
    701           for (int i = 0; i < numberRows_ - numberRows_; i++) {
    702                int iSequence = pivotVariable_[i];
    703                if (iSequence < numberColumns_ && integerType[iSequence]) {
    704                     double lower = lower_[iSequence];
    705                     double upper = upper_[iSequence];
    706                     double value = solution_[iSequence];
    707                     if (value >= lower - primalTolerance_ &&
    708                               value <= upper + primalTolerance_) {
    709                          double sign;
    710                          if (value - lower < upper - value)
    711                               sign = 1.0;
    712                          else
    713                               sign = -1.0;
    714                          cost_[iSequence] = sign * perturbationArray_[iSequence];
    715                     }
    716                }
    717           }
    718      }
    719 #if CAN_HAVE_ZERO_OBJ>1
    720      if ((specialOptions_&16777216)==0) {
    721 #endif
    722      computeDuals(givenDuals);
    723      if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
    724           const char * integerType = integerInformation();
    725           // Need to do columns and rows to stay dual feasible
    726           for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
    727                if (integerType[iSequence] && getStatus(iSequence) != basic) {
    728                     double djValue = dj_[iSequence];
    729                     double change = 0.0;
    730                     if (getStatus(iSequence) == atLowerBound)
    731                          change = CoinMax(-djValue, 10.0 * perturbationArray_[iSequence]);
    732                     else if (getStatus(iSequence) == atUpperBound)
    733                          change = CoinMin(-djValue, -10.0 * perturbationArray_[iSequence]);
    734                     cost_[iSequence] = change;
    735                     dj_[iSequence] += change;
    736                }
    737           }
    738      }
    739 
    740      // now check solutions
    741      //checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    742      //checkDualSolution();
    743      checkBothSolutions();
    744      objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_);
    745 #if CAN_HAVE_ZERO_OBJ>1
    746      } else {
    747        checkPrimalSolution( rowActivityWork_, columnActivityWork_);
     659        // allow
     660        numberOut = 0;
     661#endif
     662      }
     663      CoinSort_2(save, save + numberOut, sort);
     664      numberOut = CoinMin(maxOut, numberOut);
     665      for (iRow = 0; iRow < numberOut; iRow++) {
     666        int jRow = sort[iRow];
     667        int iColumn = pivotVariable_[jRow];
     668        setColumnStatus(iColumn, superBasic);
     669        setRowStatus(jRow, basic);
     670        pivotVariable_[jRow] = jRow + numberColumns_;
     671        if (fabs(solution_[iColumn]) > 1.0e10) {
     672          if (upper_[iColumn] < 0.0) {
     673            solution_[iColumn] = upper_[iColumn];
     674          } else if (lower_[iColumn] > 0.0) {
     675            solution_[iColumn] = lower_[iColumn];
     676          } else {
     677            solution_[iColumn] = 0.0;
     678          }
     679        }
     680      }
     681      delete[] sort;
     682    }
     683    delete[] save;
     684    if (numberOut)
     685      return numberOut;
     686  }
     687  if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
     688    //printf("trying feas pump\n");
     689    const char *integerType = integerInformation();
     690    assert(integerType);
     691    assert(perturbationArray_);
     692    CoinZeroN(cost_, numberRows_ + numberColumns_);
     693    for (int i = 0; i < numberRows_ - numberRows_; i++) {
     694      int iSequence = pivotVariable_[i];
     695      if (iSequence < numberColumns_ && integerType[iSequence]) {
     696        double lower = lower_[iSequence];
     697        double upper = upper_[iSequence];
     698        double value = solution_[iSequence];
     699        if (value >= lower - primalTolerance_ && value <= upper + primalTolerance_) {
     700          double sign;
     701          if (value - lower < upper - value)
     702            sign = 1.0;
     703          else
     704            sign = -1.0;
     705          cost_[iSequence] = sign * perturbationArray_[iSequence];
     706        }
     707      }
     708    }
     709  }
     710#if CAN_HAVE_ZERO_OBJ > 1
     711  if ((specialOptions_ & 16777216) == 0) {
     712#endif
     713    computeDuals(givenDuals);
     714    if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
     715      const char *integerType = integerInformation();
     716      // Need to do columns and rows to stay dual feasible
     717      for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
     718        if (integerType[iSequence] && getStatus(iSequence) != basic) {
     719          double djValue = dj_[iSequence];
     720          double change = 0.0;
     721          if (getStatus(iSequence) == atLowerBound)
     722            change = CoinMax(-djValue, 10.0 * perturbationArray_[iSequence]);
     723          else if (getStatus(iSequence) == atUpperBound)
     724            change = CoinMin(-djValue, -10.0 * perturbationArray_[iSequence]);
     725          cost_[iSequence] = change;
     726          dj_[iSequence] += change;
     727        }
     728      }
     729    }
     730
     731    // now check solutions
     732    //checkPrimalSolution( rowActivityWork_, columnActivityWork_);
     733    //checkDualSolution();
     734    checkBothSolutions();
     735    objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_);
     736#if CAN_HAVE_ZERO_OBJ > 1
     737  } else {
     738    checkPrimalSolution(rowActivityWork_, columnActivityWork_);
    748739#ifndef COIN_REUSE_RANDOM
    749        memset(dj_,0,(numberRows_+numberColumns_)*sizeof(double));
     740    memset(dj_, 0, (numberRows_ + numberColumns_) * sizeof(double));
    750741#else
    751        for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    752          double value;
    753          switch (getStatus(iSequence)) {
    754          case atLowerBound:
    755            value=1.0e-9*(1.0+randomNumberGenerator_.randomDouble());
    756            break;
    757          case atUpperBound:
    758            value=-1.0e-9*(1.0+randomNumberGenerator_.randomDouble());
    759            break;
    760          default:
    761            value=0.0;
    762            break;
    763          }
    764        }
    765 #endif
    766        objectiveValue_=0.0;
    767      }
    768 #endif
    769      if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 ||
    770                                       largestDualError_ > 1.0e-2))
    771           handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
    772                     << largestPrimalError_
    773                     << largestDualError_
    774                     << CoinMessageEol;
    775      if (largestPrimalError_ > 1.0e-1 && numberRows_ > 100 && numberIterations_) {
    776           // Change factorization tolerance
    777           if (factorization_->zeroTolerance() > 1.0e-18)
    778                factorization_->zeroTolerance(1.0e-18);
    779      }
    780      int returnCode=0;
    781      bool notChanged=true;
    782      if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_<0 || factorization_->pivotTolerance()<0.9899999999) &&
    783          (oldLargestDualError||oldLargestPrimalError)) {
    784        double useOldDualError=oldLargestDualError;
    785        double useDualError=largestDualError_;
    786        if (algorithm_>0&&nonLinearCost_&&
    787            nonLinearCost_->sumInfeasibilities()) {
    788          double factor=CoinMax(1.0,CoinMin(1.0e3,infeasibilityCost_*1.0e-6));
    789          useOldDualError /= factor;
    790          useDualError /= factor;
    791        }
    792        if ((largestPrimalError_>1.0e3&&
    793            oldLargestPrimalError*1.0e2<largestPrimalError_)||
    794            (useDualError>1.0e3&&
    795             useOldDualError*1.0e2<useDualError)) {
    796          double pivotTolerance = factorization_->pivotTolerance();
    797          double factor=(largestPrimalError_>1.0e10||largestDualError_>1.0e10)
    798            ? 2.0 : 1.2;
    799          if (pivotTolerance<0.1)
    800            factorization_->pivotTolerance(0.1);
    801          else if (pivotTolerance<0.98999999)
    802            factorization_->pivotTolerance(CoinMin(0.99,pivotTolerance*factor));
    803          notChanged=pivotTolerance==factorization_->pivotTolerance();
     742    for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
     743      double value;
     744      switch (getStatus(iSequence)) {
     745      case atLowerBound:
     746        value = 1.0e-9 * (1.0 + randomNumberGenerator_.randomDouble());
     747        break;
     748      case atUpperBound:
     749        value = -1.0e-9 * (1.0 + randomNumberGenerator_.randomDouble());
     750        break;
     751      default:
     752        value = 0.0;
     753        break;
     754      }
     755    }
     756#endif
     757    objectiveValue_ = 0.0;
     758  }
     759#endif
     760  if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 || largestDualError_ > 1.0e-2))
     761    handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
     762      << largestPrimalError_
     763      << largestDualError_
     764      << CoinMessageEol;
     765  if (largestPrimalError_ > 1.0e-1 && numberRows_ > 100 && numberIterations_) {
     766    // Change factorization tolerance
     767    if (factorization_->zeroTolerance() > 1.0e-18)
     768      factorization_->zeroTolerance(1.0e-18);
     769  }
     770  int returnCode = 0;
     771  bool notChanged = true;
     772  if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_ < 0 || factorization_->pivotTolerance() < 0.9899999999) && (oldLargestDualError || oldLargestPrimalError)) {
     773    double useOldDualError = oldLargestDualError;
     774    double useDualError = largestDualError_;
     775    if (algorithm_ > 0 && nonLinearCost_ && nonLinearCost_->sumInfeasibilities()) {
     776      double factor = CoinMax(1.0, CoinMin(1.0e3, infeasibilityCost_ * 1.0e-6));
     777      useOldDualError /= factor;
     778      useDualError /= factor;
     779    }
     780    if ((largestPrimalError_ > 1.0e3 && oldLargestPrimalError * 1.0e2 < largestPrimalError_) || (useDualError > 1.0e3 && useOldDualError * 1.0e2 < useDualError)) {
     781      double pivotTolerance = factorization_->pivotTolerance();
     782      double factor = (largestPrimalError_ > 1.0e10 || largestDualError_ > 1.0e10)
     783        ? 2.0
     784        : 1.2;
     785      if (pivotTolerance < 0.1)
     786        factorization_->pivotTolerance(0.1);
     787      else if (pivotTolerance < 0.98999999)
     788        factorization_->pivotTolerance(CoinMin(0.99, pivotTolerance * factor));
     789      notChanged = pivotTolerance == factorization_->pivotTolerance();
    804790#ifdef CLP_USEFUL_PRINTOUT
    805          if (pivotTolerance<0.9899999) {
    806            printf("Changing pivot tolerance from %g to %g and backtracking\n",
    807                   pivotTolerance,factorization_->pivotTolerance());
    808          }
    809          printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n",
    810                 oldLargestPrimalError,largestPrimalError_,
    811                 oldLargestDualError,largestDualError_,
    812                 pivotTolerance);
    813 #endif
    814          if (pivotTolerance<0.9899999) {
    815            largestPrimalError_=0.0;
    816            largestDualError_=0.0;
    817            returnCode=-123456789;
    818          }
    819        }
    820      }
    821      if (progress_.iterationNumber_[0]>0&&
    822          progress_.iterationNumber_[CLP_PROGRESS-1]
    823          -progress_.iterationNumber_[0]<CLP_PROGRESS*3&&
    824          factorization_->pivotTolerance()<0.25&&notChanged) {
    825        double pivotTolerance = factorization_->pivotTolerance();
    826        factorization_->pivotTolerance(pivotTolerance*1.5);
     791      if (pivotTolerance < 0.9899999) {
     792        printf("Changing pivot tolerance from %g to %g and backtracking\n",
     793          pivotTolerance, factorization_->pivotTolerance());
     794      }
     795      printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n",
     796        oldLargestPrimalError, largestPrimalError_,
     797        oldLargestDualError, largestDualError_,
     798        pivotTolerance);
     799#endif
     800      if (pivotTolerance < 0.9899999) {
     801        largestPrimalError_ = 0.0;
     802        largestDualError_ = 0.0;
     803        returnCode = -123456789;
     804      }
     805    }
     806  }
     807  if (progress_.iterationNumber_[0] > 0 && progress_.iterationNumber_[CLP_PROGRESS - 1] - progress_.iterationNumber_[0] < CLP_PROGRESS * 3 && factorization_->pivotTolerance() < 0.25 && notChanged) {
     808    double pivotTolerance = factorization_->pivotTolerance();
     809    factorization_->pivotTolerance(pivotTolerance * 1.5);
    827810#ifdef CLP_USEFUL_PRINTOUT
    828        printf("Changing pivot tolerance from %g to %g - inverting too often\n",
    829               pivotTolerance,factorization_->pivotTolerance());
    830 #endif
    831      }
    832      // Switch off false values pass indicator
    833      if (!valuesPass && algorithm_ > 0)
    834           firstFree_ = -1;
    835      if (handler_->logLevel()==63)
    836      printf("end getsolution algorithm %d status %d npinf %d sum,relaxed %g,%g ndinf %d sum,relaxed %g,%g\n",
    837             algorithm_,problemStatus_,
    838             numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,sumOfRelaxedPrimalInfeasibilities_,
    839             numberDualInfeasibilities_,sumDualInfeasibilities_,sumOfRelaxedDualInfeasibilities_);
    840      if ((moreSpecialOptions_&8388608)!=0) {
    841        if (algorithm_<0) {
    842          bool doneFiddling=false;
    843          // Optimization may make exact test iffy
    844          double testTolerance=minimumPrimalTolerance_+1.0e-15;
    845          while (!doneFiddling) {
    846            doneFiddling=true;
    847            while( !sumOfRelaxedPrimalInfeasibilities_&&
    848                   primalTolerance_>testTolerance) {
    849              // feasible - adjust tolerance
    850              double saveTolerance=primalTolerance_;
    851              primalTolerance_=CoinMax(0.25*primalTolerance_,
    852                                       minimumPrimalTolerance_);
    853              printf("Resetting primal tolerance from %g to %g\n",
    854                     saveTolerance,primalTolerance_);
    855              dblParam_[ClpPrimalTolerance]=primalTolerance_;
    856              moreSpecialOptions_ &= ~8388608;
    857              // redo with switch off
    858              returnCode=gutsOfSolution ( givenDuals,givenPrimals,valuesPass);
    859            }
    860            if(primalTolerance_>testTolerance)
    861              moreSpecialOptions_ |= 8388608; // back on
    862            if ((moreSpecialOptions_&8388608)!=0) {
    863              assert( numberPrimalInfeasibilities_);
    864              // average infeasibility
    865              double average=sumPrimalInfeasibilities_/numberPrimalInfeasibilities_;
    866              double minimum=COIN_DBL_MAX;
    867              double averageTotal=average;
    868              bool firstTime=averageInfeasibility_[0]==COIN_DBL_MAX;
    869              for (int i=0;i<CLP_INFEAS_SAVE-1;i++) {
    870                double value = averageInfeasibility_[i+1];
    871                averageTotal+=value;
    872                averageInfeasibility_[i]=value;
    873                minimum=CoinMin(minimum,value);
    874              }
    875              averageInfeasibility_[CLP_INFEAS_SAVE-1]=average;
    876              averageTotal /= CLP_INFEAS_SAVE;
    877              double oldTolerance=primalTolerance_;
    878              if (averageInfeasibility_[0]!=COIN_DBL_MAX) {
    879                if (firstTime) {
    880                  primalTolerance_=CoinMin(0.1,0.1*averageTotal);
    881                  primalTolerance_ = CoinMin(primalTolerance_,average);
    882                } else if (primalTolerance_>0.1*minimum) {
    883                  primalTolerance_=0.1*minimum;
    884                }
    885                primalTolerance_=
    886                  CoinMax(primalTolerance_,minimumPrimalTolerance_);
    887              }
    888              if (primalTolerance_!=oldTolerance) {
    889                printf("Changing primal tolerance from %g to %g\n",
    890                       oldTolerance,primalTolerance_);
    891                moreSpecialOptions_ &= ~8388608;
    892                // redo with switch off
    893                returnCode=gutsOfSolution ( givenDuals,givenPrimals,valuesPass);
    894                if(primalTolerance_>testTolerance)
    895                  moreSpecialOptions_ |= 8388608|4194304;
    896                if( !sumOfRelaxedPrimalInfeasibilities_)
    897                  doneFiddling=false; // over done it
    898              }
    899            }
    900          }
    901        }
    902      }
    903      return returnCode;
    904 }
    905 void
    906 ClpSimplex::computePrimals ( const double * rowActivities,
    907                              const double * columnActivities)
    908 {
    909 
    910      //work space
    911      CoinIndexedVector  * workSpace = rowArray_[0];
    912 
    913      CoinIndexedVector * arrayVector = rowArray_[1];
    914      arrayVector->clear();
    915      CoinIndexedVector * previousVector = rowArray_[2];
    916      previousVector->clear();
    917      // accumulate non basic stuff
    918 
    919      int iRow;
    920      // order is this way for scaling
    921      if (columnActivities != columnActivityWork_)
    922           ClpDisjointCopyN(columnActivities, numberColumns_, columnActivityWork_);
    923      if (rowActivities != rowActivityWork_)
    924           ClpDisjointCopyN(rowActivities, numberRows_, rowActivityWork_);
    925      double * array = arrayVector->denseVector();
    926      int * index = arrayVector->getIndices();
    927      int number = 0;
    928      const double * rhsOffset = matrix_->rhsOffset(this, false, true);
    929      if (!rhsOffset) {
    930           // Use whole matrix every time to make it easier for ClpMatrixBase
    931           // So zero out basic
    932           for (iRow = 0; iRow < numberRows_; iRow++) {
    933                int iPivot = pivotVariable_[iRow];
    934                assert (iPivot >= 0);
    935                solution_[iPivot] = 0.0;
     811    printf("Changing pivot tolerance from %g to %g - inverting too often\n",
     812      pivotTolerance, factorization_->pivotTolerance());
     813#endif
     814  }
     815  // Switch off false values pass indicator
     816  if (!valuesPass && algorithm_ > 0)
     817    firstFree_ = -1;
     818  if (handler_->logLevel() == 63)
     819    printf("end getsolution algorithm %d status %d npinf %d sum,relaxed %g,%g ndinf %d sum,relaxed %g,%g\n",
     820      algorithm_, problemStatus_,
     821      numberPrimalInfeasibilities_, sumPrimalInfeasibilities_, sumOfRelaxedPrimalInfeasibilities_,
     822      numberDualInfeasibilities_, sumDualInfeasibilities_, sumOfRelaxedDualInfeasibilities_);
     823  if ((moreSpecialOptions_ & 8388608) != 0) {
     824    if (algorithm_ < 0) {
     825      bool doneFiddling = false;
     826      // Optimization may make exact test iffy
     827      double testTolerance = minimumPrimalTolerance_ + 1.0e-15;
     828      while (!doneFiddling) {
     829        doneFiddling = true;
     830        while (!sumOfRelaxedPrimalInfeasibilities_ && primalTolerance_ > testTolerance) {
     831          // feasible - adjust tolerance
     832          double saveTolerance = primalTolerance_;
     833          primalTolerance_ = CoinMax(0.25 * primalTolerance_,
     834            minimumPrimalTolerance_);
     835          printf("Resetting primal tolerance from %g to %g\n",
     836            saveTolerance, primalTolerance_);
     837          dblParam_[ClpPrimalTolerance] = primalTolerance_;
     838          moreSpecialOptions_ &= ~8388608;
     839          // redo with switch off
     840          returnCode = gutsOfSolution(givenDuals, givenPrimals, valuesPass);
     841        }
     842        if (primalTolerance_ > testTolerance)
     843          moreSpecialOptions_ |= 8388608; // back on
     844        if ((moreSpecialOptions_ & 8388608) != 0) {
     845          assert(numberPrimalInfeasibilities_);
     846          // average infeasibility
     847          double average = sumPrimalInfeasibilities_ / numberPrimalInfeasibilities_;
     848          double minimum = COIN_DBL_MAX;
     849          double averageTotal = average;
     850          bool firstTime = averageInfeasibility_[0] == COIN_DBL_MAX;
     851          for (int i = 0; i < CLP_INFEAS_SAVE - 1; i++) {
     852            double value = averageInfeasibility_[i + 1];
     853            averageTotal += value;
     854            averageInfeasibility_[i] = value;
     855            minimum = CoinMin(minimum, value);
     856          }
     857          averageInfeasibility_[CLP_INFEAS_SAVE - 1] = average;
     858          averageTotal /= CLP_INFEAS_SAVE;
     859          double oldTolerance = primalTolerance_;
     860          if (averageInfeasibility_[0] != COIN_DBL_MAX) {
     861            if (firstTime) {
     862              primalTolerance_ = CoinMin(0.1, 0.1 * averageTotal);
     863              primalTolerance_ = CoinMin(primalTolerance_, average);
     864            } else if (primalTolerance_ > 0.1 * minimum) {
     865              primalTolerance_ = 0.1 * minimum;
     866            }
     867            primalTolerance_ = CoinMax(primalTolerance_, minimumPrimalTolerance_);
     868          }
     869          if (primalTolerance_ != oldTolerance) {
     870            printf("Changing primal tolerance from %g to %g\n",
     871              oldTolerance, primalTolerance_);
     872            moreSpecialOptions_ &= ~8388608;
     873            // redo with switch off
     874            returnCode = gutsOfSolution(givenDuals, givenPrimals, valuesPass);
     875            if (primalTolerance_ > testTolerance)
     876              moreSpecialOptions_ |= 8388608 | 4194304;
     877            if (!sumOfRelaxedPrimalInfeasibilities_)
     878              doneFiddling = false; // over done it
     879          }
     880        }
     881      }
     882    }
     883  }
     884  return returnCode;
     885}
     886void ClpSimplex::computePrimals(const double *rowActivities,
     887  const double *columnActivities)
     888{
     889
     890  //work space
     891  CoinIndexedVector *workSpace = rowArray_[0];
     892
     893  CoinIndexedVector *arrayVector = rowArray_[1];
     894  arrayVector->clear();
     895  CoinIndexedVector *previousVector = rowArray_[2];
     896  previousVector->clear();
     897  // accumulate non basic stuff
     898
     899  int iRow;
     900  // order is this way for scaling
     901  if (columnActivities != columnActivityWork_)
     902    ClpDisjointCopyN(columnActivities, numberColumns_, columnActivityWork_);
     903  if (rowActivities != rowActivityWork_)
     904    ClpDisjointCopyN(rowActivities, numberRows_, rowActivityWork_);
     905  double *array = arrayVector->denseVector();
     906  int *index = arrayVector->getIndices();
     907  int number = 0;
     908  const double *rhsOffset = matrix_->rhsOffset(this, false, true);
     909  if (!rhsOffset) {
     910    // Use whole matrix every time to make it easier for ClpMatrixBase
     911    // So zero out basic
     912    for (iRow = 0; iRow < numberRows_; iRow++) {
     913      int iPivot = pivotVariable_[iRow];
     914      assert(iPivot >= 0);
     915      solution_[iPivot] = 0.0;
    936916#ifdef CLP_INVESTIGATE
    937                assert (getStatus(iPivot) == basic);
    938 #endif
    939           }
    940           // Extended solution before "update"
    941           matrix_->primalExpanded(this, 0);
    942           times(-1.0, columnActivityWork_, array);
    943           for (iRow = 0; iRow < numberRows_; iRow++) {
    944                double value = array[iRow] + rowActivityWork_[iRow];
    945                if (value) {
    946                     array[iRow] = value;
    947                     index[number++] = iRow;
    948                } else {
    949                     array[iRow] = 0.0;
    950                }
    951           }
    952      } else {
    953           // we have an effective rhs lying around
    954           // zero out basic (really just for slacks)
    955           for (iRow = 0; iRow < numberRows_; iRow++) {
    956                int iPivot = pivotVariable_[iRow];
    957                solution_[iPivot] = 0.0;
    958           }
    959           for (iRow = 0; iRow < numberRows_; iRow++) {
    960                double value = rhsOffset[iRow] + rowActivityWork_[iRow];
    961                if (value) {
    962                     array[iRow] = value;
    963                     index[number++] = iRow;
    964                } else {
    965                     array[iRow] = 0.0;
    966                }
    967           }
    968      }
    969      arrayVector->setNumElements(number);
     917      assert(getStatus(iPivot) == basic);
     918#endif
     919    }
     920    // Extended solution before "update"
     921    matrix_->primalExpanded(this, 0);
     922    times(-1.0, columnActivityWork_, array);
     923    for (iRow = 0; iRow < numberRows_; iRow++) {
     924      double value = array[iRow] + rowActivityWork_[iRow];
     925      if (value) {
     926        array[iRow] = value;
     927        index[number++] = iRow;
     928      } else {
     929        array[iRow] = 0.0;
     930      }
     931    }
     932  } else {
     933    // we have an effective rhs lying around
     934    // zero out basic (really just for slacks)
     935    for (iRow = 0; iRow < numberRows_; iRow++) {
     936      int iPivot = pivotVariable_[iRow];
     937      solution_[iPivot] = 0.0;
     938    }
     939    for (iRow = 0; iRow < numberRows_; iRow++) {
     940      double value = rhsOffset[iRow] + rowActivityWork_[iRow];
     941      if (value) {
     942        array[iRow] = value;
     943        index[number++] = iRow;
     944      } else {
     945        array[iRow] = 0.0;
     946      }
     947    }
     948  }
     949  arrayVector->setNumElements(number);
    970950#ifdef CLP_DEBUG
    971      if (numberIterations_ == -3840) {
    972           int i;
    973           for (i = 0; i < numberRows_ + numberColumns_; i++)
    974                printf("%d status %d\n", i, status_[i]);
    975           printf("xxxxx1\n");
    976           for (i = 0; i < numberRows_; i++)
    977                if (array[i])
    978                     printf("%d rhs %g\n", i, array[i]);
    979           printf("xxxxx2\n");
    980           for (i = 0; i < numberRows_ + numberColumns_; i++)
    981                if (getStatus(i) != basic)
    982                     printf("%d non basic %g %g %g\n", i, lower_[i], solution_[i], upper_[i]);
    983           printf("xxxxx3\n");
    984      }
    985 #endif
    986      // Ftran adjusted RHS and iterate to improve accuracy
    987      double lastError = COIN_DBL_MAX;
    988      int iRefine;
    989      CoinIndexedVector * thisVector = arrayVector;
    990      CoinIndexedVector * lastVector = previousVector;
     951  if (numberIterations_ == -3840) {
     952    int i;
     953    for (i = 0; i < numberRows_ + numberColumns_; i++)
     954      printf("%d status %d\n", i, status_[i]);
     955    printf("xxxxx1\n");
     956    for (i = 0; i < numberRows_; i++)
     957      if (array[i])
     958        printf("%d rhs %g\n", i, array[i]);
     959    printf("xxxxx2\n");
     960    for (i = 0; i < numberRows_ + numberColumns_; i++)
     961      if (getStatus(i) != basic)
     962        printf("%d non basic %g %g %g\n", i, lower_[i], solution_[i], upper_[i]);
     963    printf("xxxxx3\n");
     964  }
     965#endif
     966  // Ftran adjusted RHS and iterate to improve accuracy
     967  double lastError = COIN_DBL_MAX;
     968  int iRefine;
     969  CoinIndexedVector *thisVector = arrayVector;
     970  CoinIndexedVector *lastVector = previousVector;
    991971#if 0
    992972     static double * xsave=NULL;
     
    10321012     }
    10331013#endif
    1034      //printf ("ZZ0 n before %d",number);
    1035      if (number)
    1036           factorization_->updateColumn(workSpace, thisVector);
    1037      //printf(" - after %d\n",thisVector->getNumElements());
    1038      double * work = workSpace->denseVector();
     1014  //printf ("ZZ0 n before %d",number);
     1015  if (number)
     1016    factorization_->updateColumn(workSpace, thisVector);
     1017  //printf(" - after %d\n",thisVector->getNumElements());
     1018  double *work = workSpace->denseVector();
    10391019#ifdef CLP_DEBUG
    1040      if (numberIterations_ == -3840) {
    1041           int i;
    1042           for (i = 0; i < numberRows_; i++)
    1043                if (array[i])
    1044                     printf("%d after rhs %g\n", i, array[i]);
    1045           printf("xxxxx4\n");
    1046      }
    1047 #endif
    1048      bool goodSolution = true;
    1049      for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
    1050 
    1051           int numberIn = thisVector->getNumElements();
    1052           int * indexIn = thisVector->getIndices();
    1053           double * arrayIn = thisVector->denseVector();
    1054           // put solution in correct place
    1055           if (!rhsOffset) {
    1056                int j;
    1057                for (j = 0; j < numberIn; j++) {
    1058                     iRow = indexIn[j];
    1059                     int iPivot = pivotVariable_[iRow];
    1060                     solution_[iPivot] = arrayIn[iRow];
    1061                     //assert (fabs(solution_[iPivot])<1.0e100);
    1062                }
     1020  if (numberIterations_ == -3840) {
     1021    int i;
     1022    for (i = 0; i < numberRows_; i++)
     1023      if (array[i])
     1024        printf("%d after rhs %g\n", i, array[i]);
     1025    printf("xxxxx4\n");
     1026  }
     1027#endif
     1028  bool goodSolution = true;
     1029  for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
     1030
     1031    int numberIn = thisVector->getNumElements();
     1032    int *indexIn = thisVector->getIndices();
     1033    double *arrayIn = thisVector->denseVector();
     1034    // put solution in correct place
     1035    if (!rhsOffset) {
     1036      int j;
     1037      for (j = 0; j < numberIn; j++) {
     1038        iRow = indexIn[j];
     1039        int iPivot = pivotVariable_[iRow];
     1040        solution_[iPivot] = arrayIn[iRow];
     1041        //assert (fabs(solution_[iPivot])<1.0e100);
     1042      }
     1043    } else {
     1044      for (iRow = 0; iRow < numberRows_; iRow++) {
     1045        int iPivot = pivotVariable_[iRow];
     1046        solution_[iPivot] = arrayIn[iRow];
     1047        //assert (fabs(solution_[iPivot])<1.0e100);
     1048      }
     1049    }
     1050    // Extended solution after "update"
     1051    matrix_->primalExpanded(this, 1);
     1052    // check Ax == b  (for all)
     1053    // signal column generated matrix to just do basic (and gub)
     1054    unsigned int saveOptions = specialOptions();
     1055    setSpecialOptions(16);
     1056    times(-1.0, columnActivityWork_, work);
     1057    setSpecialOptions(saveOptions);
     1058    largestPrimalError_ = 0.0;
     1059    double multiplier = 131072.0;
     1060    for (iRow = 0; iRow < numberRows_; iRow++) {
     1061      double value = work[iRow] + rowActivityWork_[iRow];
     1062      work[iRow] = value * multiplier;
     1063      if (fabs(value) > largestPrimalError_) {
     1064        largestPrimalError_ = fabs(value);
     1065      }
     1066    }
     1067    if (largestPrimalError_ >= lastError) {
     1068      // restore
     1069      CoinIndexedVector *temp = thisVector;
     1070      thisVector = lastVector;
     1071      lastVector = temp;
     1072      goodSolution = false;
     1073      break;
     1074    }
     1075    if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
     1076      // try and make better
     1077      // save this
     1078      CoinIndexedVector *temp = thisVector;
     1079      thisVector = lastVector;
     1080      lastVector = temp;
     1081      int *indexOut = thisVector->getIndices();
     1082      int number = 0;
     1083      array = thisVector->denseVector();
     1084      thisVector->clear();
     1085      for (iRow = 0; iRow < numberRows_; iRow++) {
     1086        double value = work[iRow];
     1087        if (value) {
     1088          array[iRow] = value;
     1089          indexOut[number++] = iRow;
     1090          work[iRow] = 0.0;
     1091        }
     1092      }
     1093      thisVector->setNumElements(number);
     1094      lastError = largestPrimalError_;
     1095      //printf ("ZZ%d n before %d",iRefine+1,number);
     1096      factorization_->updateColumn(workSpace, thisVector);
     1097      //printf(" - after %d\n",thisVector->getNumElements());
     1098      multiplier = 1.0 / multiplier;
     1099      double *previous = lastVector->denseVector();
     1100      number = 0;
     1101      for (iRow = 0; iRow < numberRows_; iRow++) {
     1102        double value = previous[iRow] + multiplier * array[iRow];
     1103        if (value) {
     1104          array[iRow] = value;
     1105          indexOut[number++] = iRow;
     1106        } else {
     1107          array[iRow] = 0.0;
     1108        }
     1109      }
     1110      thisVector->setNumElements(number);
     1111    } else {
     1112      break;
     1113    }
     1114  }
     1115
     1116  // solution as accurate as we are going to get
     1117  ClpFillN(work, numberRows_, 0.0);
     1118  if (!goodSolution) {
     1119    array = thisVector->denseVector();
     1120    // put solution in correct place
     1121    for (iRow = 0; iRow < numberRows_; iRow++) {
     1122      int iPivot = pivotVariable_[iRow];
     1123      solution_[iPivot] = array[iRow];
     1124      //assert (fabs(solution_[iPivot])<1.0e100);
     1125    }
     1126  }
     1127  arrayVector->clear();
     1128  previousVector->clear();
     1129#ifdef CLP_DEBUG
     1130  if (numberIterations_ == -3840) {
     1131    exit(77);
     1132  }
     1133#endif
     1134}
     1135// now dual side
     1136void ClpSimplex::computeDuals(double *givenDjs)
     1137{
     1138#ifndef SLIM_CLP
     1139  if (objective_->type() == 1 || !objective_->activated()) {
     1140#endif
     1141    // Linear
     1142    //work space
     1143    CoinIndexedVector *workSpace = rowArray_[0];
     1144
     1145    CoinIndexedVector *arrayVector = rowArray_[1];
     1146    arrayVector->clear();
     1147    CoinIndexedVector *previousVector = rowArray_[2];
     1148    previousVector->clear();
     1149    int iRow;
     1150#ifdef CLP_DEBUG
     1151    workSpace->checkClear();
     1152#endif
     1153    double *array = arrayVector->denseVector();
     1154    int *index = arrayVector->getIndices();
     1155    int number = 0;
     1156    if (!givenDjs) {
     1157      for (iRow = 0; iRow < numberRows_; iRow++) {
     1158        int iPivot = pivotVariable_[iRow];
     1159        double value = cost_[iPivot];
     1160        if (value) {
     1161          array[iRow] = value;
     1162          index[number++] = iRow;
     1163        }
     1164      }
     1165    } else {
     1166      // dual values pass - djs may not be zero
     1167      for (iRow = 0; iRow < numberRows_; iRow++) {
     1168        int iPivot = pivotVariable_[iRow];
     1169        // make sure zero if done
     1170        if (!pivoted(iPivot))
     1171          givenDjs[iPivot] = 0.0;
     1172        double value = cost_[iPivot] - givenDjs[iPivot];
     1173        if (value) {
     1174          array[iRow] = value;
     1175          index[number++] = iRow;
     1176        }
     1177      }
     1178    }
     1179    arrayVector->setNumElements(number);
     1180    // Extended duals before "updateTranspose"
     1181    matrix_->dualExpanded(this, arrayVector, givenDjs, 0);
     1182
     1183    // Btran basic costs and get as accurate as possible
     1184    double lastError = COIN_DBL_MAX;
     1185    int iRefine;
     1186    double *work = workSpace->denseVector();
     1187    CoinIndexedVector *thisVector = arrayVector;
     1188    CoinIndexedVector *lastVector = previousVector;
     1189    factorization_->updateColumnTranspose(workSpace, thisVector);
     1190
     1191    for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
     1192      // check basic reduced costs zero
     1193      largestDualError_ = 0.0;
     1194      if (!numberExtraRows_) {
     1195        // Just basic
     1196        int *index2 = workSpace->getIndices();
     1197        // use reduced costs for slacks as work array
     1198        double *work2 = reducedCostWork_ + numberColumns_;
     1199        int numberStructurals = 0;
     1200        for (iRow = 0; iRow < numberRows_; iRow++) {
     1201          int iPivot = pivotVariable_[iRow];
     1202          if (iPivot < numberColumns_)
     1203            index2[numberStructurals++] = iPivot;
     1204        }
     1205        matrix_->listTransposeTimes(this, array, index2, numberStructurals, work2);
     1206        numberStructurals = 0;
     1207        if (!givenDjs) {
     1208          for (iRow = 0; iRow < numberRows_; iRow++) {
     1209            int iPivot = pivotVariable_[iRow];
     1210            double value;
     1211            if (iPivot >= numberColumns_) {
     1212              // slack
     1213              value = rowObjectiveWork_[iPivot - numberColumns_]
     1214                + array[iPivot - numberColumns_];
     1215            } else {
     1216              // column
     1217              value = objectiveWork_[iPivot] - work2[numberStructurals++];
     1218            }
     1219            work[iRow] = value;
     1220            if (fabs(value) > largestDualError_) {
     1221              largestDualError_ = fabs(value);
     1222            }
     1223          }
     1224        } else {
     1225          for (iRow = 0; iRow < numberRows_; iRow++) {
     1226            int iPivot = pivotVariable_[iRow];
     1227            if (iPivot >= numberColumns_) {
     1228              // slack
     1229              work[iRow] = rowObjectiveWork_[iPivot - numberColumns_]
     1230                + array[iPivot - numberColumns_] - givenDjs[iPivot];
     1231            } else {
     1232              // column
     1233              work[iRow] = objectiveWork_[iPivot] - work2[numberStructurals++]
     1234                - givenDjs[iPivot];
     1235            }
     1236            if (fabs(work[iRow]) > largestDualError_) {
     1237              largestDualError_ = fabs(work[iRow]);
     1238              //assert (largestDualError_<1.0e-7);
     1239              //if (largestDualError_>1.0e-7)
     1240              //printf("large dual error %g\n",largestDualError_);
     1241            }
     1242          }
     1243        }
     1244      } else {
     1245        // extra rows - be more careful
     1246#if 1
     1247        // would be faster to do just for basic but this reduces code
     1248        ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
     1249        transposeTimes(-1.0, array, reducedCostWork_);
     1250#else
     1251      // Just basic
     1252      int *index2 = workSpace->getIndices();
     1253      int numberStructurals = 0;
     1254      for (iRow = 0; iRow < numberRows_; iRow++) {
     1255        int iPivot = pivotVariable_[iRow];
     1256        if (iPivot < numberColumns_)
     1257          index2[numberStructurals++] = iPivot;
     1258      }
     1259      matrix_->listTransposeTimes(this, array, index2, numberStructurals, work);
     1260      for (iRow = 0; iRow < numberStructurals; iRow++) {
     1261        int iPivot = index2[iRow];
     1262        reducedCostWork_[iPivot] = objectiveWork_[iPivot] - work[iRow];
     1263      }
     1264#endif
     1265        // update by duals on sets
     1266        matrix_->dualExpanded(this, NULL, NULL, 1);
     1267        if (!givenDjs) {
     1268          for (iRow = 0; iRow < numberRows_; iRow++) {
     1269            int iPivot = pivotVariable_[iRow];
     1270            double value;
     1271            if (iPivot >= numberColumns_) {
     1272              // slack
     1273              value = rowObjectiveWork_[iPivot - numberColumns_]
     1274                + array[iPivot - numberColumns_];
     1275            } else {
     1276              // column
     1277              value = reducedCostWork_[iPivot];
     1278            }
     1279            work[iRow] = value;
     1280            if (fabs(value) > largestDualError_) {
     1281              largestDualError_ = fabs(value);
     1282            }
     1283          }
     1284        } else {
     1285          for (iRow = 0; iRow < numberRows_; iRow++) {
     1286            int iPivot = pivotVariable_[iRow];
     1287            if (iPivot >= numberColumns_) {
     1288              // slack
     1289              work[iRow] = rowObjectiveWork_[iPivot - numberColumns_]
     1290                + array[iPivot - numberColumns_] - givenDjs[iPivot];
     1291            } else {
     1292              // column
     1293              work[iRow] = reducedCostWork_[iPivot] - givenDjs[iPivot];
     1294            }
     1295            if (fabs(work[iRow]) > largestDualError_) {
     1296              largestDualError_ = fabs(work[iRow]);
     1297              //assert (largestDualError_<1.0e-7);
     1298              //if (largestDualError_>1.0e-7)
     1299              //printf("large dual error %g\n",largestDualError_);
     1300            }
     1301          }
     1302        }
     1303      }
     1304      if (largestDualError_ >= lastError) {
     1305        // restore
     1306        CoinIndexedVector *temp = thisVector;
     1307        thisVector = lastVector;
     1308        lastVector = temp;
     1309        break;
     1310      }
     1311      if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
     1312        && !givenDjs) {
     1313        // try and make better
     1314        // save this
     1315        CoinIndexedVector *temp = thisVector;
     1316        thisVector = lastVector;
     1317        lastVector = temp;
     1318        int *indexOut = thisVector->getIndices();
     1319        int number = 0;
     1320        array = thisVector->denseVector();
     1321        thisVector->clear();
     1322        double multiplier = 131072.0;
     1323        for (iRow = 0; iRow < numberRows_; iRow++) {
     1324          double value = multiplier * work[iRow];
     1325          if (value) {
     1326            array[iRow] = value;
     1327            indexOut[number++] = iRow;
     1328            work[iRow] = 0.0;
     1329          }
     1330          work[iRow] = 0.0;
     1331        }
     1332        thisVector->setNumElements(number);
     1333        lastError = largestDualError_;
     1334        factorization_->updateColumnTranspose(workSpace, thisVector);
     1335        multiplier = 1.0 / multiplier;
     1336        double *previous = lastVector->denseVector();
     1337        number = 0;
     1338        for (iRow = 0; iRow < numberRows_; iRow++) {
     1339          double value = previous[iRow] + multiplier * array[iRow];
     1340          if (value) {
     1341            array[iRow] = value;
     1342            indexOut[number++] = iRow;
    10631343          } else {
    1064                for (iRow = 0; iRow < numberRows_; iRow++) {
    1065                     int iPivot = pivotVariable_[iRow];
    1066                     solution_[iPivot] = arrayIn[iRow];
    1067                     //assert (fabs(solution_[iPivot])<1.0e100);
    1068                }
    1069           }
    1070           // Extended solution after "update"
    1071           matrix_->primalExpanded(this, 1);
    1072           // check Ax == b  (for all)
    1073           // signal column generated matrix to just do basic (and gub)
    1074           unsigned int saveOptions = specialOptions();
    1075           setSpecialOptions(16);
    1076           times(-1.0, columnActivityWork_, work);
    1077           setSpecialOptions(saveOptions);
    1078           largestPrimalError_ = 0.0;
    1079           double multiplier = 131072.0;
    1080           for (iRow = 0; iRow < numberRows_; iRow++) {
    1081                double value = work[iRow] + rowActivityWork_[iRow];
    1082                work[iRow] = value * multiplier;
    1083                if (fabs(value) > largestPrimalError_) {
    1084                     largestPrimalError_ = fabs(value);
    1085                }
    1086           }
    1087           if (largestPrimalError_ >= lastError) {
    1088                // restore
    1089                CoinIndexedVector * temp = thisVector;
    1090                thisVector = lastVector;
    1091                lastVector = temp;
    1092                goodSolution = false;
    1093                break;
    1094           }
    1095           if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
    1096                // try and make better
    1097                // save this
    1098                CoinIndexedVector * temp = thisVector;
    1099                thisVector = lastVector;
    1100                lastVector = temp;
    1101                int * indexOut = thisVector->getIndices();
    1102                int number = 0;
    1103                array = thisVector->denseVector();
    1104                thisVector->clear();
    1105                for (iRow = 0; iRow < numberRows_; iRow++) {
    1106                     double value = work[iRow];
    1107                     if (value) {
    1108                          array[iRow] = value;
    1109                          indexOut[number++] = iRow;
    1110                          work[iRow] = 0.0;
    1111                     }
    1112                }
    1113                thisVector->setNumElements(number);
    1114                lastError = largestPrimalError_;
    1115                //printf ("ZZ%d n before %d",iRefine+1,number);
    1116                factorization_->updateColumn(workSpace, thisVector);
    1117                //printf(" - after %d\n",thisVector->getNumElements());
    1118                multiplier = 1.0 / multiplier;
    1119                double * previous = lastVector->denseVector();
    1120                number = 0;
    1121                for (iRow = 0; iRow < numberRows_; iRow++) {
    1122                     double value = previous[iRow] + multiplier * array[iRow];
    1123                     if (value) {
    1124                          array[iRow] = value;
    1125                          indexOut[number++] = iRow;
    1126                     } else {
    1127                          array[iRow] = 0.0;
    1128                     }
    1129                }
    1130                thisVector->setNumElements(number);
    1131           } else {
    1132                break;
    1133           }
    1134      }
    1135 
    1136      // solution as accurate as we are going to get
    1137      ClpFillN(work, numberRows_, 0.0);
    1138      if (!goodSolution) {
    1139           array = thisVector->denseVector();
    1140           // put solution in correct place
    1141           for (iRow = 0; iRow < numberRows_; iRow++) {
    1142                int iPivot = pivotVariable_[iRow];
    1143                solution_[iPivot] = array[iRow];
    1144                //assert (fabs(solution_[iPivot])<1.0e100);
    1145           }
    1146      }
    1147      arrayVector->clear();
    1148      previousVector->clear();
    1149 #ifdef CLP_DEBUG
    1150      if (numberIterations_ == -3840) {
    1151           exit(77);
    1152      }
    1153 #endif
    1154 }
    1155 // now dual side
    1156 void
    1157 ClpSimplex::computeDuals(double * givenDjs)
    1158 {
     1344            array[iRow] = 0.0;
     1345          }
     1346        }
     1347        thisVector->setNumElements(number);
     1348      } else {
     1349        break;
     1350      }
     1351    }
     1352    // now look at dual solution
     1353    array = thisVector->denseVector();
     1354    for (iRow = 0; iRow < numberRows_; iRow++) {
     1355      // slack
     1356      double value = array[iRow];
     1357      dual_[iRow] = value;
     1358      value += rowObjectiveWork_[iRow];
     1359      rowReducedCost_[iRow] = value;
     1360    }
     1361    // can use work if problem scaled (for better cache)
     1362    ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix_);
     1363    double *saveRowScale = rowScale_;
     1364    //double * saveColumnScale = columnScale_;
     1365    if (scaledMatrix_) {
     1366      rowScale_ = NULL;
     1367      clpMatrix = scaledMatrix_;
     1368    }
     1369    if (clpMatrix && (clpMatrix->flags() & 2) == 0) {
     1370      CoinIndexedVector *cVector = columnArray_[0];
     1371      int *whichColumn = cVector->getIndices();
     1372      assert(!cVector->getNumElements());
     1373      int n = 0;
     1374      for (int i = 0; i < numberColumns_; i++) {
     1375        if (getColumnStatus(i) != basic) {
     1376          whichColumn[n++] = i;
     1377          reducedCostWork_[i] = objectiveWork_[i];
     1378        } else {
     1379          reducedCostWork_[i] = 0.0;
     1380        }
     1381      }
     1382      if (numberRows_ > 4000)
     1383        clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
     1384          rowScale_, columnScale_, work);
     1385      else
     1386        clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
     1387          rowScale_, columnScale_, NULL);
     1388    } else {
     1389      ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
     1390      if (numberRows_ > 4000)
     1391        matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
     1392          rowScale_, columnScale_, work);
     1393      else
     1394        matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
     1395          rowScale_, columnScale_, NULL);
     1396    }
     1397    rowScale_ = saveRowScale;
     1398    //columnScale_ = saveColumnScale;
     1399    ClpFillN(work, numberRows_, 0.0);
     1400    // Extended duals and check dual infeasibility
     1401    if (!matrix_->skipDualCheck() || algorithm_ < 0 || problemStatus_ != -2)
     1402      matrix_->dualExpanded(this, NULL, NULL, 2);
     1403    // If necessary - override results
     1404    if (givenDjs) {
     1405      // restore accurate duals
     1406      CoinMemcpyN(dj_, (numberRows_ + numberColumns_), givenDjs);
     1407    }
     1408    arrayVector->clear();
     1409    previousVector->clear();
    11591410#ifndef SLIM_CLP
    1160      if (objective_->type() == 1 || !objective_->activated()) {
    1161 #endif
    1162           // Linear
    1163           //work space
    1164           CoinIndexedVector  * workSpace = rowArray_[0];
    1165 
    1166           CoinIndexedVector * arrayVector = rowArray_[1];
    1167           arrayVector->clear();
    1168           CoinIndexedVector * previousVector = rowArray_[2];
    1169           previousVector->clear();
    1170           int iRow;
    1171 #ifdef CLP_DEBUG
    1172           workSpace->checkClear();
    1173 #endif
    1174           double * array = arrayVector->denseVector();
    1175           int * index = arrayVector->getIndices();
    1176           int number = 0;
    1177           if (!givenDjs) {
    1178                for (iRow = 0; iRow < numberRows_; iRow++) {
    1179                     int iPivot = pivotVariable_[iRow];
    1180                     double value = cost_[iPivot];
    1181                     if (value) {
    1182                          array[iRow] = value;
    1183                          index[number++] = iRow;
    1184                     }
    1185                }
    1186           } else {
    1187                // dual values pass - djs may not be zero
    1188                for (iRow = 0; iRow < numberRows_; iRow++) {
    1189                     int iPivot = pivotVariable_[iRow];
    1190                     // make sure zero if done
    1191                     if (!pivoted(iPivot))
    1192                          givenDjs[iPivot] = 0.0;
    1193                     double value = cost_[iPivot] - givenDjs[iPivot];
    1194                     if (value) {
    1195                          array[iRow] = value;
    1196                          index[number++] = iRow;
    1197                     }
    1198                }
    1199           }
    1200           arrayVector->setNumElements(number);
    1201           // Extended duals before "updateTranspose"
    1202           matrix_->dualExpanded(this, arrayVector, givenDjs, 0);
    1203 
    1204           // Btran basic costs and get as accurate as possible
    1205           double lastError = COIN_DBL_MAX;
    1206           int iRefine;
    1207           double * work = workSpace->denseVector();
    1208           CoinIndexedVector * thisVector = arrayVector;
    1209           CoinIndexedVector * lastVector = previousVector;
    1210           factorization_->updateColumnTranspose(workSpace, thisVector);
    1211 
    1212           for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
    1213                // check basic reduced costs zero
    1214                largestDualError_ = 0.0;
    1215                if (!numberExtraRows_) {
    1216                     // Just basic
    1217                     int * index2 = workSpace->getIndices();
    1218                     // use reduced costs for slacks as work array
    1219                     double * work2 = reducedCostWork_ + numberColumns_;
    1220                     int numberStructurals = 0;
    1221                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1222                          int iPivot = pivotVariable_[iRow];
    1223                          if (iPivot < numberColumns_)
    1224                               index2[numberStructurals++] = iPivot;
    1225                     }
    1226                     matrix_->listTransposeTimes(this, array, index2, numberStructurals, work2);
    1227                     numberStructurals = 0;
    1228                     if (!givenDjs) {
    1229                          for (iRow = 0; iRow < numberRows_; iRow++) {
    1230                               int iPivot = pivotVariable_[iRow];
    1231                               double value;
    1232                               if (iPivot >= numberColumns_) {
    1233                                    // slack
    1234                                    value = rowObjectiveWork_[iPivot-numberColumns_]
    1235                                            + array[iPivot-numberColumns_];
    1236                               } else {
    1237                                    // column
    1238                                    value = objectiveWork_[iPivot] - work2[numberStructurals++];
    1239                               }
    1240                               work[iRow] = value;
    1241                               if (fabs(value) > largestDualError_) {
    1242                                    largestDualError_ = fabs(value);
    1243                               }
    1244                          }
    1245                     } else {
    1246                          for (iRow = 0; iRow < numberRows_; iRow++) {
    1247                               int iPivot = pivotVariable_[iRow];
    1248                               if (iPivot >= numberColumns_) {
    1249                                    // slack
    1250                                    work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
    1251                                                 + array[iPivot-numberColumns_] - givenDjs[iPivot];
    1252                               } else {
    1253                                    // column
    1254                                    work[iRow] = objectiveWork_[iPivot] - work2[numberStructurals++]
    1255                                                 - givenDjs[iPivot];
    1256                               }
    1257                               if (fabs(work[iRow]) > largestDualError_) {
    1258                                    largestDualError_ = fabs(work[iRow]);
    1259                                    //assert (largestDualError_<1.0e-7);
    1260                                    //if (largestDualError_>1.0e-7)
    1261                                    //printf("large dual error %g\n",largestDualError_);
    1262                               }
    1263                          }
    1264                     }
    1265                } else {
    1266                     // extra rows - be more careful
    1267 #if 1
    1268                     // would be faster to do just for basic but this reduces code
    1269                     ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
    1270                     transposeTimes(-1.0, array, reducedCostWork_);
    1271 #else
    1272                     // Just basic
    1273                     int * index2 = workSpace->getIndices();
    1274                     int numberStructurals = 0;
    1275                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1276                          int iPivot = pivotVariable_[iRow];
    1277                          if (iPivot < numberColumns_)
    1278                               index2[numberStructurals++] = iPivot;
    1279                     }
    1280                     matrix_->listTransposeTimes(this, array, index2, numberStructurals, work);
    1281                     for (iRow = 0; iRow < numberStructurals; iRow++) {
    1282                          int iPivot = index2[iRow];
    1283                          reducedCostWork_[iPivot] = objectiveWork_[iPivot] - work[iRow];
    1284                     }
    1285 #endif
    1286                     // update by duals on sets
    1287                     matrix_->dualExpanded(this, NULL, NULL, 1);
    1288                     if (!givenDjs) {
    1289                          for (iRow = 0; iRow < numberRows_; iRow++) {
    1290                               int iPivot = pivotVariable_[iRow];
    1291                               double value;
    1292                               if (iPivot >= numberColumns_) {
    1293                                    // slack
    1294                                    value = rowObjectiveWork_[iPivot-numberColumns_]
    1295                                            + array[iPivot-numberColumns_];
    1296                               } else {
    1297                                    // column
    1298                                    value = reducedCostWork_[iPivot];
    1299                               }
    1300                               work[iRow] = value;
    1301                               if (fabs(value) > largestDualError_) {
    1302                                    largestDualError_ = fabs(value);
    1303                               }
    1304                          }
    1305                     } else {
    1306                          for (iRow = 0; iRow < numberRows_; iRow++) {
    1307                               int iPivot = pivotVariable_[iRow];
    1308                               if (iPivot >= numberColumns_) {
    1309                                    // slack
    1310                                    work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
    1311                                                 + array[iPivot-numberColumns_] - givenDjs[iPivot];
    1312                               } else {
    1313                                    // column
    1314                                    work[iRow] = reducedCostWork_[iPivot] - givenDjs[iPivot];
    1315                               }
    1316                               if (fabs(work[iRow]) > largestDualError_) {
    1317                                    largestDualError_ = fabs(work[iRow]);
    1318                                    //assert (largestDualError_<1.0e-7);
    1319                                    //if (largestDualError_>1.0e-7)
    1320                                    //printf("large dual error %g\n",largestDualError_);
    1321                               }
    1322                          }
    1323                     }
    1324                }
    1325                if (largestDualError_ >= lastError) {
    1326                     // restore
    1327                     CoinIndexedVector * temp = thisVector;
    1328                     thisVector = lastVector;
    1329                     lastVector = temp;
    1330                     break;
    1331                }
    1332                if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
    1333                          && !givenDjs) {
    1334                     // try and make better
    1335                     // save this
    1336                     CoinIndexedVector * temp = thisVector;
    1337                     thisVector = lastVector;
    1338                     lastVector = temp;
    1339                     int * indexOut = thisVector->getIndices();
    1340                     int number = 0;
    1341                     array = thisVector->denseVector();
    1342                     thisVector->clear();
    1343                     double multiplier = 131072.0;
    1344                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1345                          double value = multiplier * work[iRow];
    1346                          if (value) {
    1347                               array[iRow] = value;
    1348                               indexOut[number++] = iRow;
    1349                               work[iRow] = 0.0;
    1350                          }
    1351                          work[iRow] = 0.0;
    1352                     }
    1353                     thisVector->setNumElements(number);
    1354                     lastError = largestDualError_;
    1355                     factorization_->updateColumnTranspose(workSpace, thisVector);
    1356                     multiplier = 1.0 / multiplier;
    1357                     double * previous = lastVector->denseVector();
    1358                     number = 0;
    1359                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1360                          double value = previous[iRow] + multiplier * array[iRow];
    1361                          if (value) {
    1362                               array[iRow] = value;
    1363                               indexOut[number++] = iRow;
    1364                          } else {
    1365                               array[iRow] = 0.0;
    1366                          }
    1367                     }
    1368                     thisVector->setNumElements(number);
    1369                } else {
    1370                     break;
    1371                }
    1372           }
    1373           // now look at dual solution
    1374           array = thisVector->denseVector();
    1375           for (iRow = 0; iRow < numberRows_; iRow++) {
    1376                // slack
    1377                double value = array[iRow];
    1378                dual_[iRow] = value;
    1379                value += rowObjectiveWork_[iRow];
    1380                rowReducedCost_[iRow] = value;
    1381           }
    1382           // can use work if problem scaled (for better cache)
    1383           ClpPackedMatrix* clpMatrix =
    1384                dynamic_cast< ClpPackedMatrix*>(matrix_);
    1385           double * saveRowScale = rowScale_;
    1386           //double * saveColumnScale = columnScale_;
    1387           if (scaledMatrix_) {
    1388                rowScale_ = NULL;
    1389                clpMatrix = scaledMatrix_;
    1390           }
    1391           if (clpMatrix && (clpMatrix->flags() & 2) == 0) {
    1392                CoinIndexedVector * cVector = columnArray_[0];
    1393                int * whichColumn = cVector->getIndices();
    1394                assert (!cVector->getNumElements());
    1395                int n = 0;
    1396                for (int i = 0; i < numberColumns_; i++) {
    1397                     if (getColumnStatus(i) != basic) {
    1398                          whichColumn[n++] = i;
    1399                          reducedCostWork_[i] = objectiveWork_[i];
    1400                     } else {
    1401                          reducedCostWork_[i] = 0.0;
    1402                     }
    1403                }
    1404                if (numberRows_ > 4000)
    1405                     clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
    1406                                                     rowScale_, columnScale_, work);
    1407                else
    1408                     clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
    1409                                                     rowScale_, columnScale_, NULL);
    1410           } else {
    1411                ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
    1412                if (numberRows_ > 4000)
    1413                     matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
    1414                                             rowScale_, columnScale_, work);
    1415                else
    1416                     matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
    1417                                             rowScale_, columnScale_, NULL);
    1418           }
    1419           rowScale_ = saveRowScale;
    1420           //columnScale_ = saveColumnScale;
    1421           ClpFillN(work, numberRows_, 0.0);
    1422           // Extended duals and check dual infeasibility
    1423           if (!matrix_->skipDualCheck() || algorithm_ < 0 || problemStatus_ != -2)
    1424                matrix_->dualExpanded(this, NULL, NULL, 2);
    1425           // If necessary - override results
    1426           if (givenDjs) {
    1427                // restore accurate duals
    1428                CoinMemcpyN(dj_, (numberRows_ + numberColumns_), givenDjs);
    1429           }
    1430           arrayVector->clear();
    1431           previousVector->clear();
    1432 #ifndef SLIM_CLP
    1433      } else {
    1434           // Nonlinear
    1435           objective_->reducedGradient(this, dj_, false);
    1436           // get dual_ by moving from reduced costs for slacks
    1437           CoinMemcpyN(dj_ + numberColumns_, numberRows_, dual_);
    1438      }
     1411  } else {
     1412    // Nonlinear
     1413    objective_->reducedGradient(this, dj_, false);
     1414    // get dual_ by moving from reduced costs for slacks
     1415    CoinMemcpyN(dj_ + numberColumns_, numberRows_, dual_);
     1416  }
    14391417#endif
    14401418}
     
    14421420   primal and dual solutions.  Uses input arrays for variables at
    14431421   bounds.  Returns feasibility states */
    1444 int ClpSimplex::getSolution ( const double * /*rowActivities*/,
    1445                               const double * /*columnActivities*/)
    1446 {
    1447      if (!factorization_->status()) {
    1448           // put in standard form
    1449           createRim(7 + 8 + 16 + 32, false, -1);
    1450           if (pivotVariable_[0] < 0)
    1451                internalFactorize(0);
    1452           // do work
    1453           gutsOfSolution ( NULL, NULL);
    1454           // release extra memory
    1455           deleteRim(0);
    1456      }
    1457      return factorization_->status();
     1422int ClpSimplex::getSolution(const double * /*rowActivities*/,
     1423  const double * /*columnActivities*/)
     1424{
     1425  if (!factorization_->status()) {
     1426    // put in standard form
     1427    createRim(7 + 8 + 16 + 32, false, -1);
     1428    if (pivotVariable_[0] < 0)
     1429      internalFactorize(0);
     1430    // do work
     1431    gutsOfSolution(NULL, NULL);
     1432    // release extra memory
     1433    deleteRim(0);
     1434  }
     1435  return factorization_->status();
    14581436}
    14591437/* Given an existing factorization computes and checks
    14601438   primal and dual solutions.  Uses current problem arrays for
    14611439   bounds.  Returns feasibility states */
    1462 int ClpSimplex::getSolution ( )
    1463 {
    1464      double * rowActivities = new double[numberRows_];
    1465      double * columnActivities = new double[numberColumns_];
    1466      ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
    1467      ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
    1468      int status = getSolution( rowActivities, columnActivities);
    1469      delete [] rowActivities;
    1470      delete [] columnActivities;
    1471      return status;
     1440int ClpSimplex::getSolution()
     1441{
     1442  double *rowActivities = new double[numberRows_];
     1443  double *columnActivities = new double[numberColumns_];
     1444  ClpDisjointCopyN(rowActivityWork_, numberRows_, rowActivities);
     1445  ClpDisjointCopyN(columnActivityWork_, numberColumns_, columnActivities);
     1446  int status = getSolution(rowActivities, columnActivities);
     1447  delete[] rowActivities;
     1448  delete[] columnActivities;
     1449  return status;
    14721450}
    14731451// Factorizes using current basis.  This is for external use
    14741452// Return codes are as from ClpFactorization
    1475 int ClpSimplex::factorize ()
    1476 {
    1477      // put in standard form
    1478      createRim(7 + 8 + 16 + 32, false);
    1479      // do work
    1480      int status = internalFactorize(-1);
    1481      // release extra memory
    1482      deleteRim(0);
    1483 
    1484      return status;
     1453int ClpSimplex::factorize()
     1454{
     1455  // put in standard form
     1456  createRim(7 + 8 + 16 + 32, false);
     1457  // do work
     1458  int status = internalFactorize(-1);
     1459  // release extra memory
     1460  deleteRim(0);
     1461
     1462  return status;
    14851463}
    14861464// Clean up status
    1487 void
    1488 ClpSimplex::cleanStatus()
    1489 {
    1490      int iRow, iColumn;
    1491      int numberBasic = 0;
    1492      // make row activities correct
    1493      memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
    1494      times(1.0, columnActivityWork_, rowActivityWork_);
    1495      if (!status_)
    1496           createStatus();
    1497      for (iRow = 0; iRow < numberRows_; iRow++) {
    1498           if (getRowStatus(iRow) == basic)
    1499                numberBasic++;
    1500           else {
    1501                setRowStatus(iRow, superBasic);
    1502                // but put to bound if close
    1503                if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
    1504                          <= primalTolerance_) {
    1505                     rowActivityWork_[iRow] = rowLowerWork_[iRow];
    1506                     setRowStatus(iRow, atLowerBound);
    1507                } else if (fabs(rowActivityWork_[iRow] - rowUpperWork_[iRow])
    1508                           <= primalTolerance_) {
    1509                     rowActivityWork_[iRow] = rowUpperWork_[iRow];
    1510                     setRowStatus(iRow, atUpperBound);
    1511                }
    1512           }
    1513      }
    1514      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1515           if (getColumnStatus(iColumn) == basic) {
    1516                if (numberBasic == numberRows_) {
    1517                     // take out of basis
    1518                     setColumnStatus(iColumn, superBasic);
    1519                     // but put to bound if close
    1520                     if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
    1521                               <= primalTolerance_) {
    1522                          columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1523                          setColumnStatus(iColumn, atLowerBound);
    1524                     } else if (fabs(columnActivityWork_[iColumn]
    1525                                     - columnUpperWork_[iColumn])
    1526                                <= primalTolerance_) {
    1527                          columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1528                          setColumnStatus(iColumn, atUpperBound);
    1529                     }
    1530                } else
    1531                     numberBasic++;
    1532           } else {
    1533                setColumnStatus(iColumn, superBasic);
    1534                // but put to bound if close
    1535                if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
    1536                          <= primalTolerance_) {
    1537                     columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1538                     setColumnStatus(iColumn, atLowerBound);
    1539                } else if (fabs(columnActivityWork_[iColumn]
    1540                                - columnUpperWork_[iColumn])
    1541                           <= primalTolerance_) {
    1542                     columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1543                     setColumnStatus(iColumn, atUpperBound);
    1544                }
    1545           }
    1546      }
     1465void ClpSimplex::cleanStatus()
     1466{
     1467  int iRow, iColumn;
     1468  int numberBasic = 0;
     1469  // make row activities correct
     1470  memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
     1471  times(1.0, columnActivityWork_, rowActivityWork_);
     1472  if (!status_)
     1473    createStatus();
     1474  for (iRow = 0; iRow < numberRows_; iRow++) {
     1475    if (getRowStatus(iRow) == basic)
     1476      numberBasic++;
     1477    else {
     1478      setRowStatus(iRow, superBasic);
     1479      // but put to bound if close
     1480      if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
     1481        <= primalTolerance_) {
     1482        rowActivityWork_[iRow] = rowLowerWork_[iRow];
     1483        setRowStatus(iRow, atLowerBound);
     1484      } else if (fabs(rowActivityWork_[iRow] - rowUpperWork_[iRow])
     1485        <= primalTolerance_) {
     1486        rowActivityWork_[iRow] = rowUpperWork_[iRow];
     1487        setRowStatus(iRow, atUpperBound);
     1488      }
     1489    }
     1490  }
     1491  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1492    if (getColumnStatus(iColumn) == basic) {
     1493      if (numberBasic == numberRows_) {
     1494        // take out of basis
     1495        setColumnStatus(iColumn, superBasic);
     1496        // but put to bound if close
     1497        if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
     1498          <= primalTolerance_) {
     1499          columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1500          setColumnStatus(iColumn, atLowerBound);
     1501        } else if (fabs(columnActivityWork_[iColumn]
     1502                     - columnUpperWork_[iColumn])
     1503          <= primalTolerance_) {
     1504          columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1505          setColumnStatus(iColumn, atUpperBound);
     1506        }
     1507      } else
     1508        numberBasic++;
     1509    } else {
     1510      setColumnStatus(iColumn, superBasic);
     1511      // but put to bound if close
     1512      if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
     1513        <= primalTolerance_) {
     1514        columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1515        setColumnStatus(iColumn, atLowerBound);
     1516      } else if (fabs(columnActivityWork_[iColumn]
     1517                   - columnUpperWork_[iColumn])
     1518        <= primalTolerance_) {
     1519        columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1520        setColumnStatus(iColumn, atUpperBound);
     1521      }
     1522    }
     1523  }
    15471524}
    15481525
     
    15551532   Special case is numberRows_+1 -> all slack basis.
    15561533*/
    1557 int ClpSimplex::internalFactorize ( int solveType)
    1558 {
    1559      int iRow, iColumn;
    1560      int totalSlacks = numberRows_;
    1561      if (!status_)
    1562           createStatus();
    1563 
    1564      bool valuesPass = false;
    1565      if (solveType >= 10) {
    1566           valuesPass = true;
    1567           solveType -= 10;
    1568      }
     1534int ClpSimplex::internalFactorize(int solveType)
     1535{
     1536  int iRow, iColumn;
     1537  int totalSlacks = numberRows_;
     1538  if (!status_)
     1539    createStatus();
     1540
     1541  bool valuesPass = false;
     1542  if (solveType >= 10) {
     1543    valuesPass = true;
     1544    solveType -= 10;
     1545  }
    15691546#ifdef CLP_DEBUG
    1570      if (solveType > 0) {
    1571           int numberFreeIn = 0, numberFreeOut = 0;
    1572           double biggestDj = 0.0;
    1573           for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1574                switch(getColumnStatus(iColumn)) {
    1575 
    1576                case basic:
    1577                     if (columnLower_[iColumn] < -largeValue_
    1578                               && columnUpper_[iColumn] > largeValue_)
    1579                          numberFreeIn++;
    1580                     break;
    1581                default:
    1582                     if (columnLower_[iColumn] < -largeValue_
    1583                               && columnUpper_[iColumn] > largeValue_) {
    1584                          numberFreeOut++;
    1585                          biggestDj = CoinMax(fabs(dj_[iColumn]), biggestDj);
    1586                     }
    1587                     break;
    1588                }
    1589           }
    1590           if (numberFreeIn + numberFreeOut)
    1591                printf("%d in basis, %d out - largest dj %g\n",
    1592                       numberFreeIn, numberFreeOut, biggestDj);
    1593      }
    1594 #endif
    1595      if (solveType <= 0) {
    1596           // Make sure everything is clean
    1597           for (iRow = 0; iRow < numberRows_; iRow++) {
    1598                if(getRowStatus(iRow) == isFixed) {
    1599                     // double check fixed
    1600                     if (rowUpperWork_[iRow] > rowLowerWork_[iRow])
    1601                          setRowStatus(iRow, atLowerBound);
    1602                } else if (getRowStatus(iRow) == isFree) {
    1603                     // may not be free after all
    1604                     if (rowLowerWork_[iRow] > -largeValue_ || rowUpperWork_[iRow] < largeValue_)
    1605                          setRowStatus(iRow, superBasic);
    1606                }
    1607           }
    1608           for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1609                if(getColumnStatus(iColumn) == isFixed) {
    1610                     // double check fixed
    1611                     if (columnUpperWork_[iColumn] > columnLowerWork_[iColumn])
    1612                          setColumnStatus(iColumn, atLowerBound);
    1613                } else if (getColumnStatus(iColumn) == isFree) {
    1614                     // may not be free after all
    1615                     if (columnLowerWork_[iColumn] > -largeValue_ || columnUpperWork_[iColumn] < largeValue_)
    1616                          setColumnStatus(iColumn, superBasic);
    1617                }
    1618           }
    1619           if (!valuesPass) {
    1620                // not values pass so set to bounds
    1621                bool allSlack = true;
    1622                if (status_) {
    1623                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1624                          if (getRowStatus(iRow) != basic) {
    1625                               allSlack = false;
    1626                               break;
    1627                          }
    1628                     }
    1629                }
    1630                if (!allSlack) {
    1631                     //#define CLP_INVESTIGATE2
     1547  if (solveType > 0) {
     1548    int numberFreeIn = 0, numberFreeOut = 0;
     1549    double biggestDj = 0.0;
     1550    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1551      switch (getColumnStatus(iColumn)) {
     1552
     1553      case basic:
     1554        if (columnLower_[iColumn] < -largeValue_
     1555          && columnUpper_[iColumn] > largeValue_)
     1556          numberFreeIn++;
     1557        break;
     1558      default:
     1559        if (columnLower_[iColumn] < -largeValue_
     1560          && columnUpper_[iColumn] > largeValue_) {
     1561          numberFreeOut++;
     1562          biggestDj = CoinMax(fabs(dj_[iColumn]), biggestDj);
     1563        }
     1564        break;
     1565      }
     1566    }
     1567    if (numberFreeIn + numberFreeOut)
     1568      printf("%d in basis, %d out - largest dj %g\n",
     1569        numberFreeIn, numberFreeOut, biggestDj);
     1570  }
     1571#endif
     1572  if (solveType <= 0) {
     1573    // Make sure everything is clean
     1574    for (iRow = 0; iRow < numberRows_; iRow++) {
     1575      if (getRowStatus(iRow) == isFixed) {
     1576        // double check fixed
     1577        if (rowUpperWork_[iRow] > rowLowerWork_[iRow])
     1578          setRowStatus(iRow, atLowerBound);
     1579      } else if (getRowStatus(iRow) == isFree) {
     1580        // may not be free after all
     1581        if (rowLowerWork_[iRow] > -largeValue_ || rowUpperWork_[iRow] < largeValue_)
     1582          setRowStatus(iRow, superBasic);
     1583      }
     1584    }
     1585    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1586      if (getColumnStatus(iColumn) == isFixed) {
     1587        // double check fixed
     1588        if (columnUpperWork_[iColumn] > columnLowerWork_[iColumn])
     1589          setColumnStatus(iColumn, atLowerBound);
     1590      } else if (getColumnStatus(iColumn) == isFree) {
     1591        // may not be free after all
     1592        if (columnLowerWork_[iColumn] > -largeValue_ || columnUpperWork_[iColumn] < largeValue_)
     1593          setColumnStatus(iColumn, superBasic);
     1594      }
     1595    }
     1596    if (!valuesPass) {
     1597      // not values pass so set to bounds
     1598      bool allSlack = true;
     1599      if (status_) {
     1600        for (iRow = 0; iRow < numberRows_; iRow++) {
     1601          if (getRowStatus(iRow) != basic) {
     1602            allSlack = false;
     1603            break;
     1604          }
     1605        }
     1606      }
     1607      if (!allSlack) {
     1608        //#define CLP_INVESTIGATE2
    16321609#ifdef CLP_INVESTIGATE3
    1633                     int numberTotal = numberRows_ + numberColumns_;
    1634                     double * saveSol = valuesPass ?
    1635                                        CoinCopyOfArray(solution_, numberTotal) : NULL;
    1636 #endif
    1637                     // set values from warm start (if sensible)
    1638                     int numberBasic = 0;
    1639                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1640                          switch(getRowStatus(iRow)) {
    1641 
    1642                          case basic:
    1643                               numberBasic++;
    1644                               break;
    1645                          case atUpperBound:
    1646                               rowActivityWork_[iRow] = rowUpperWork_[iRow];
    1647                               if (rowActivityWork_[iRow] > largeValue_) {
    1648                                    if (rowLowerWork_[iRow] > -largeValue_) {
    1649                                         rowActivityWork_[iRow] = rowLowerWork_[iRow];
    1650                                         setRowStatus(iRow, atLowerBound);
    1651                                    } else {
    1652                                         // say free
    1653                                         setRowStatus(iRow, isFree);
    1654                                         rowActivityWork_[iRow] = 0.0;
    1655                                    }
    1656                               }
    1657                               break;
    1658                          case ClpSimplex::isFixed:
    1659                          case atLowerBound:
    1660                               rowActivityWork_[iRow] = rowLowerWork_[iRow];
    1661                               if (rowActivityWork_[iRow] < -largeValue_) {
    1662                                    if (rowUpperWork_[iRow] < largeValue_) {
    1663                                         rowActivityWork_[iRow] = rowUpperWork_[iRow];
    1664                                         setRowStatus(iRow, atUpperBound);
    1665                                    } else {
    1666                                         // say free
    1667                                         setRowStatus(iRow, isFree);
    1668                                         rowActivityWork_[iRow] = 0.0;
    1669                                    }
    1670                               }
    1671                               break;
    1672                          case isFree:
    1673                               break;
    1674                               // not really free - fall through to superbasic
    1675                          case superBasic:
    1676                               if (rowUpperWork_[iRow] > largeValue_) {
    1677                                    if (rowLowerWork_[iRow] > -largeValue_) {
    1678                                         rowActivityWork_[iRow] = rowLowerWork_[iRow];
    1679                                         setRowStatus(iRow, atLowerBound);
    1680                                    } else {
    1681                                         // say free
    1682                                         setRowStatus(iRow, isFree);
    1683                                         rowActivityWork_[iRow] = 0.0;
    1684                                    }
    1685                               } else {
    1686                                    if (rowLowerWork_[iRow] > -largeValue_) {
    1687                                         // set to nearest
    1688                                         if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
    1689                                                   < fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])) {
    1690                                              rowActivityWork_[iRow] = rowLowerWork_[iRow];
    1691                                              setRowStatus(iRow, atLowerBound);
    1692                                         } else {
    1693                                              rowActivityWork_[iRow] = rowUpperWork_[iRow];
    1694                                              setRowStatus(iRow, atUpperBound);
    1695                                         }
    1696                                    } else {
    1697                                         rowActivityWork_[iRow] = rowUpperWork_[iRow];
    1698                                         setRowStatus(iRow, atUpperBound);
    1699                                    }
    1700                               }
    1701                               break;
    1702                          }
    1703                     }
    1704                     totalSlacks = numberBasic;
    1705 
    1706                     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1707                          switch(getColumnStatus(iColumn)) {
    1708 
    1709                          case basic:
    1710                               if (numberBasic == maximumBasic_) {
    1711                                    // take out of basis
    1712                                    if (columnLowerWork_[iColumn] > -largeValue_) {
    1713                                         if (columnActivityWork_[iColumn] - columnLowerWork_[iColumn] <
    1714                                                   columnUpperWork_[iColumn] - columnActivityWork_[iColumn]) {
    1715                                              columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1716                                              setColumnStatus(iColumn, atLowerBound);
    1717                                         } else {
    1718                                              columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1719                                              setColumnStatus(iColumn, atUpperBound);
    1720                                         }
    1721                                    } else if (columnUpperWork_[iColumn] < largeValue_) {
    1722                                         columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1723                                         setColumnStatus(iColumn, atUpperBound);
    1724                                    } else {
    1725                                         columnActivityWork_[iColumn] = 0.0;
    1726                                         setColumnStatus(iColumn, isFree);
    1727                                    }
    1728                               } else {
    1729                                    numberBasic++;
    1730                               }
    1731                               break;
    1732                          case atUpperBound:
    1733                               columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1734                               if (columnActivityWork_[iColumn] > largeValue_) {
    1735                                    if (columnLowerWork_[iColumn] < -largeValue_) {
    1736                                         columnActivityWork_[iColumn] = 0.0;
    1737                                         setColumnStatus(iColumn, isFree);
    1738                                    } else {
    1739                                         columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1740                                         setColumnStatus(iColumn, atLowerBound);
    1741                                    }
    1742                               }
    1743                               break;
    1744                          case isFixed:
    1745                          case atLowerBound:
    1746                               columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1747                               if (columnActivityWork_[iColumn] < -largeValue_) {
    1748                                    if (columnUpperWork_[iColumn] > largeValue_) {
    1749                                         columnActivityWork_[iColumn] = 0.0;
    1750                                         setColumnStatus(iColumn, isFree);
    1751                                    } else {
    1752                                         columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1753                                         setColumnStatus(iColumn, atUpperBound);
    1754                                    }
    1755                               }
    1756                               break;
    1757                          case isFree:
    1758                               break;
    1759                               // not really free - fall through to superbasic
    1760                          case superBasic:
    1761                               if (columnUpperWork_[iColumn] > largeValue_) {
    1762                                    if (columnLowerWork_[iColumn] > -largeValue_) {
    1763                                         columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1764                                         setColumnStatus(iColumn, atLowerBound);
    1765                                    } else {
    1766                                         // say free
    1767                                         setColumnStatus(iColumn, isFree);
    1768                                         columnActivityWork_[iColumn] = 0.0;
    1769                                    }
    1770                               } else {
    1771                                    if (columnLowerWork_[iColumn] > -largeValue_) {
    1772                                         // set to nearest
    1773                                         if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
    1774                                                   < fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])) {
    1775                                              columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1776                                              setColumnStatus(iColumn, atLowerBound);
    1777                                         } else {
    1778                                              columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1779                                              setColumnStatus(iColumn, atUpperBound);
    1780                                         }
    1781                                    } else {
    1782                                         columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1783                                         setColumnStatus(iColumn, atUpperBound);
    1784                                    }
    1785                               }
    1786                               break;
    1787                          }
    1788                     }
     1610        int numberTotal = numberRows_ + numberColumns_;
     1611        double *saveSol = valuesPass ? CoinCopyOfArray(solution_, numberTotal) : NULL;
     1612#endif
     1613        // set values from warm start (if sensible)
     1614        int numberBasic = 0;
     1615        for (iRow = 0; iRow < numberRows_; iRow++) {
     1616          switch (getRowStatus(iRow)) {
     1617
     1618          case basic:
     1619            numberBasic++;
     1620            break;
     1621          case atUpperBound:
     1622            rowActivityWork_[iRow] = rowUpperWork_[iRow];
     1623            if (rowActivityWork_[iRow] > largeValue_) {
     1624              if (rowLowerWork_[iRow] > -largeValue_) {
     1625                rowActivityWork_[iRow] = rowLowerWork_[iRow];
     1626                setRowStatus(iRow, atLowerBound);
     1627              } else {
     1628                // say free
     1629                setRowStatus(iRow, isFree);
     1630                rowActivityWork_[iRow] = 0.0;
     1631              }
     1632            }
     1633            break;
     1634          case ClpSimplex::isFixed:
     1635          case atLowerBound:
     1636            rowActivityWork_[iRow] = rowLowerWork_[iRow];
     1637            if (rowActivityWork_[iRow] < -largeValue_) {
     1638              if (rowUpperWork_[iRow] < largeValue_) {
     1639                rowActivityWork_[iRow] = rowUpperWork_[iRow];
     1640                setRowStatus(iRow, atUpperBound);
     1641              } else {
     1642                // say free
     1643                setRowStatus(iRow, isFree);
     1644                rowActivityWork_[iRow] = 0.0;
     1645              }
     1646            }
     1647            break;
     1648          case isFree:
     1649            break;
     1650            // not really free - fall through to superbasic
     1651          case superBasic:
     1652            if (rowUpperWork_[iRow] > largeValue_) {
     1653              if (rowLowerWork_[iRow] > -largeValue_) {
     1654                rowActivityWork_[iRow] = rowLowerWork_[iRow];
     1655                setRowStatus(iRow, atLowerBound);
     1656              } else {
     1657                // say free
     1658                setRowStatus(iRow, isFree);
     1659                rowActivityWork_[iRow] = 0.0;
     1660              }
     1661            } else {
     1662              if (rowLowerWork_[iRow] > -largeValue_) {
     1663                // set to nearest
     1664                if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
     1665                  < fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])) {
     1666                  rowActivityWork_[iRow] = rowLowerWork_[iRow];
     1667                  setRowStatus(iRow, atLowerBound);
     1668                } else {
     1669                  rowActivityWork_[iRow] = rowUpperWork_[iRow];
     1670                  setRowStatus(iRow, atUpperBound);
     1671                }
     1672              } else {
     1673                rowActivityWork_[iRow] = rowUpperWork_[iRow];
     1674                setRowStatus(iRow, atUpperBound);
     1675              }
     1676            }
     1677            break;
     1678          }
     1679        }
     1680        totalSlacks = numberBasic;
     1681
     1682        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1683          switch (getColumnStatus(iColumn)) {
     1684
     1685          case basic:
     1686            if (numberBasic == maximumBasic_) {
     1687              // take out of basis
     1688              if (columnLowerWork_[iColumn] > -largeValue_) {
     1689                if (columnActivityWork_[iColumn] - columnLowerWork_[iColumn] < columnUpperWork_[iColumn] - columnActivityWork_[iColumn]) {
     1690                  columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1691                  setColumnStatus(iColumn, atLowerBound);
     1692                } else {
     1693                  columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1694                  setColumnStatus(iColumn, atUpperBound);
     1695                }
     1696              } else if (columnUpperWork_[iColumn] < largeValue_) {
     1697                columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1698                setColumnStatus(iColumn, atUpperBound);
     1699              } else {
     1700                columnActivityWork_[iColumn] = 0.0;
     1701                setColumnStatus(iColumn, isFree);
     1702              }
     1703            } else {
     1704              numberBasic++;
     1705            }
     1706            break;
     1707          case atUpperBound:
     1708            columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1709            if (columnActivityWork_[iColumn] > largeValue_) {
     1710              if (columnLowerWork_[iColumn] < -largeValue_) {
     1711                columnActivityWork_[iColumn] = 0.0;
     1712                setColumnStatus(iColumn, isFree);
     1713              } else {
     1714                columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1715                setColumnStatus(iColumn, atLowerBound);
     1716              }
     1717            }
     1718            break;
     1719          case isFixed:
     1720          case atLowerBound:
     1721            columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1722            if (columnActivityWork_[iColumn] < -largeValue_) {
     1723              if (columnUpperWork_[iColumn] > largeValue_) {
     1724                columnActivityWork_[iColumn] = 0.0;
     1725                setColumnStatus(iColumn, isFree);
     1726              } else {
     1727                columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1728                setColumnStatus(iColumn, atUpperBound);
     1729              }
     1730            }
     1731            break;
     1732          case isFree:
     1733            break;
     1734            // not really free - fall through to superbasic
     1735          case superBasic:
     1736            if (columnUpperWork_[iColumn] > largeValue_) {
     1737              if (columnLowerWork_[iColumn] > -largeValue_) {
     1738                columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1739                setColumnStatus(iColumn, atLowerBound);
     1740              } else {
     1741                // say free
     1742                setColumnStatus(iColumn, isFree);
     1743                columnActivityWork_[iColumn] = 0.0;
     1744              }
     1745            } else {
     1746              if (columnLowerWork_[iColumn] > -largeValue_) {
     1747                // set to nearest
     1748                if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
     1749                  < fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])) {
     1750                  columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1751                  setColumnStatus(iColumn, atLowerBound);
     1752                } else {
     1753                  columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1754                  setColumnStatus(iColumn, atUpperBound);
     1755                }
     1756              } else {
     1757                columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1758                setColumnStatus(iColumn, atUpperBound);
     1759              }
     1760            }
     1761            break;
     1762          }
     1763        }
    17891764#ifdef CLP_INVESTIGATE3
    1790                     if (saveSol) {
    1791                          int numberChanged = 0;
    1792                          double largestChanged = 0.0;
    1793                          for (int i = 0; i < numberTotal; i++) {
    1794                               double difference = fabs(solution_[i] - saveSol[i]);
    1795                               if (difference > 1.0e-7) {
    1796                                    numberChanged++;
    1797                                    if (difference > largestChanged)
    1798                                         largestChanged = difference;
    1799                               }
    1800                          }
    1801                          if (numberChanged)
    1802                               printf("%d changed, largest %g\n", numberChanged, largestChanged);
    1803                          delete [] saveSol;
    1804                     }
     1765        if (saveSol) {
     1766          int numberChanged = 0;
     1767          double largestChanged = 0.0;
     1768          for (int i = 0; i < numberTotal; i++) {
     1769            double difference = fabs(solution_[i] - saveSol[i]);
     1770            if (difference > 1.0e-7) {
     1771              numberChanged++;
     1772              if (difference > largestChanged)
     1773                largestChanged = difference;
     1774            }
     1775          }
     1776          if (numberChanged)
     1777            printf("%d changed, largest %g\n", numberChanged, largestChanged);
     1778          delete[] saveSol;
     1779        }
    18051780#endif
    18061781#if 0
     
    18241799                    }
    18251800#endif
    1826                } else {
    1827                     // all slack basis
    1828                     int numberBasic = 0;
    1829                     if (!status_) {
    1830                          createStatus();
    1831                     }
    1832                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1833                          double lower = rowLowerWork_[iRow];
    1834                          double upper = rowUpperWork_[iRow];
    1835                          if (lower > -largeValue_ || upper < largeValue_) {
    1836                               if (fabs(lower) <= fabs(upper)) {
    1837                                    rowActivityWork_[iRow] = lower;
    1838                               } else {
    1839                                    rowActivityWork_[iRow] = upper;
    1840                               }
    1841                          } else {
    1842                               rowActivityWork_[iRow] = 0.0;
    1843                          }
    1844                          setRowStatus(iRow, basic);
    1845                          numberBasic++;
    1846                     }
    1847                     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1848                          double lower = columnLowerWork_[iColumn];
    1849                          double upper = columnUpperWork_[iColumn];
    1850                          double big_bound = largeValue_;
    1851                          double value = columnActivityWork_[iColumn];
    1852                          if (lower > -big_bound || upper < big_bound) {
    1853                               if ((getColumnStatus(iColumn) == atLowerBound &&
    1854                                         columnActivityWork_[iColumn] == lower) ||
    1855                                         (getColumnStatus(iColumn) == atUpperBound &&
    1856                                          columnActivityWork_[iColumn] == upper)) {
    1857                                    // status looks plausible
    1858                               } else {
    1859                                    // set to sensible
    1860                                 if (getColumnStatus(iColumn) == atUpperBound
    1861                                     && upper < 1.0e20) {
    1862                                   columnActivityWork_[iColumn] = upper;
    1863                                 } else if (getColumnStatus(iColumn) == atLowerBound
    1864                                            && lower > -1.0e20) {
    1865                                   columnActivityWork_[iColumn] = lower;
    1866                                 } else {
    1867                                   if (fabs(lower) <= fabs(upper)) {
    1868                                     setColumnStatus(iColumn, atLowerBound);
    1869                                     columnActivityWork_[iColumn] = lower;
    1870                                   } else {
    1871                                     setColumnStatus(iColumn, atUpperBound);
    1872                                     columnActivityWork_[iColumn] = upper;
    1873                                   }
    1874                                 }
    1875                               }
    1876                          } else {
    1877                               setColumnStatus(iColumn, isFree);
    1878                               columnActivityWork_[iColumn] = 0.0;
    1879                          }
    1880                     }
    1881                }
     1801      } else {
     1802        // all slack basis
     1803        int numberBasic = 0;
     1804        if (!status_) {
     1805          createStatus();
     1806        }
     1807        for (iRow = 0; iRow < numberRows_; iRow++) {
     1808          double lower = rowLowerWork_[iRow];
     1809          double upper = rowUpperWork_[iRow];
     1810          if (lower > -largeValue_ || upper < largeValue_) {
     1811            if (fabs(lower) <= fabs(upper)) {
     1812              rowActivityWork_[iRow] = lower;
     1813            } else {
     1814              rowActivityWork_[iRow] = upper;
     1815            }
    18821816          } else {
    1883                // values pass has less coding
    1884                // make row activities correct and clean basis a bit
    1885                cleanStatus();
    1886                if (status_) {
    1887                     int numberBasic = 0;
    1888                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1889                          if (getRowStatus(iRow) == basic)
    1890                               numberBasic++;
    1891                     }
    1892                     totalSlacks = numberBasic;
     1817            rowActivityWork_[iRow] = 0.0;
     1818          }
     1819          setRowStatus(iRow, basic);
     1820          numberBasic++;
     1821        }
     1822        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1823          double lower = columnLowerWork_[iColumn];
     1824          double upper = columnUpperWork_[iColumn];
     1825          double big_bound = largeValue_;
     1826          double value = columnActivityWork_[iColumn];
     1827          if (lower > -big_bound || upper < big_bound) {
     1828            if ((getColumnStatus(iColumn) == atLowerBound && columnActivityWork_[iColumn] == lower) || (getColumnStatus(iColumn) == atUpperBound && columnActivityWork_[iColumn] == upper)) {
     1829              // status looks plausible
     1830            } else {
     1831              // set to sensible
     1832              if (getColumnStatus(iColumn) == atUpperBound
     1833                && upper < 1.0e20) {
     1834                columnActivityWork_[iColumn] = upper;
     1835              } else if (getColumnStatus(iColumn) == atLowerBound
     1836                && lower > -1.0e20) {
     1837                columnActivityWork_[iColumn] = lower;
     1838              } else {
     1839                if (fabs(lower) <= fabs(upper)) {
     1840                  setColumnStatus(iColumn, atLowerBound);
     1841                  columnActivityWork_[iColumn] = lower;
     1842                } else {
     1843                  setColumnStatus(iColumn, atUpperBound);
     1844                  columnActivityWork_[iColumn] = upper;
     1845                }
     1846              }
     1847            }
     1848          } else {
     1849            setColumnStatus(iColumn, isFree);
     1850            columnActivityWork_[iColumn] = 0.0;
     1851          }
     1852        }
     1853      }
     1854    } else {
     1855      // values pass has less coding
     1856      // make row activities correct and clean basis a bit
     1857      cleanStatus();
     1858      if (status_) {
     1859        int numberBasic = 0;
     1860        for (iRow = 0; iRow < numberRows_; iRow++) {
     1861          if (getRowStatus(iRow) == basic)
     1862            numberBasic++;
     1863        }
     1864        totalSlacks = numberBasic;
    18931865#if 0
    18941866                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     
    18971869                    }
    18981870#endif
    1899                } else {
    1900                     // all slack basis
    1901                     int numberBasic = 0;
    1902                     if (!status_) {
    1903                          createStatus();
    1904                     }
    1905                     for (iRow = 0; iRow < numberRows_; iRow++) {
    1906                          setRowStatus(iRow, basic);
    1907                          numberBasic++;
    1908                     }
    1909                     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1910                          setColumnStatus(iColumn, superBasic);
    1911                          // but put to bound if close
    1912                          if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
    1913                                    <= primalTolerance_) {
    1914                               columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1915                               setColumnStatus(iColumn, atLowerBound);
    1916                          } else if (fabs(columnActivityWork_[iColumn]
    1917                                          - columnUpperWork_[iColumn])
    1918                                     <= primalTolerance_) {
    1919                               columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
    1920                               setColumnStatus(iColumn, atUpperBound);
    1921                          }
    1922                     }
    1923                }
    1924           }
    1925           numberRefinements_ = 1;
    1926           // set fixed if they are
    1927           for (iRow = 0; iRow < numberRows_; iRow++) {
    1928                if (getRowStatus(iRow) != basic ) {
    1929                     if (rowLowerWork_[iRow] == rowUpperWork_[iRow]) {
    1930                          rowActivityWork_[iRow] = rowLowerWork_[iRow];
    1931                          setRowStatus(iRow, isFixed);
    1932                     }
    1933                }
    1934           }
    1935           for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    1936                if (getColumnStatus(iColumn) != basic ) {
    1937                     if (columnLowerWork_[iColumn] == columnUpperWork_[iColumn]) {
    1938                          columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
    1939                          setColumnStatus(iColumn, isFixed);
    1940                     }
    1941                }
    1942           }
    1943      }
    1944      //for (iRow=0;iRow<numberRows_+numberColumns_;iRow++) {
    1945      //if (fabs(solution_[iRow])>1.0e10) {
    1946      //  printf("large %g at %d - status %d\n",
    1947      //         solution_[iRow],iRow,status_[iRow]);
    1948      //}
    1949      //}
    1950 #    if 0 //ndef _MSC_VER
    1951          // The local static var k is a problem when trying to build a DLL. Since this is
    1952          // just for debugging (likely done on *nix), just hide it from Windows
    1953         // -- lh, 101016 --
     1871      } else {
     1872        // all slack basis
     1873        int numberBasic = 0;
     1874        if (!status_) {
     1875          createStatus();
     1876        }
     1877        for (iRow = 0; iRow < numberRows_; iRow++) {
     1878          setRowStatus(iRow, basic);
     1879          numberBasic++;
     1880        }
     1881        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1882          setColumnStatus(iColumn, superBasic);
     1883          // but put to bound if close
     1884          if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
     1885            <= primalTolerance_) {
     1886            columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1887            setColumnStatus(iColumn, atLowerBound);
     1888          } else if (fabs(columnActivityWork_[iColumn]
     1889                       - columnUpperWork_[iColumn])
     1890            <= primalTolerance_) {
     1891            columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
     1892            setColumnStatus(iColumn, atUpperBound);
     1893          }
     1894        }
     1895      }
     1896    }
     1897    numberRefinements_ = 1;
     1898    // set fixed if they are
     1899    for (iRow = 0; iRow < numberRows_; iRow++) {
     1900      if (getRowStatus(iRow) != basic) {
     1901        if (rowLowerWork_[iRow] == rowUpperWork_[iRow]) {
     1902          rowActivityWork_[iRow] = rowLowerWork_[iRow];
     1903          setRowStatus(iRow, isFixed);
     1904        }
     1905      }
     1906    }
     1907    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1908      if (getColumnStatus(iColumn) != basic) {
     1909        if (columnLowerWork_[iColumn] == columnUpperWork_[iColumn]) {
     1910          columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
     1911          setColumnStatus(iColumn, isFixed);
     1912        }
     1913      }
     1914    }
     1915  }
     1916  //for (iRow=0;iRow<numberRows_+numberColumns_;iRow++) {
     1917  //if (fabs(solution_[iRow])>1.0e10) {
     1918  //  printf("large %g at %d - status %d\n",
     1919  //         solution_[iRow],iRow,status_[iRow]);
     1920  //}
     1921  //}
     1922#if 0 //ndef _MSC_VER                                                              \
     1923  // The local static var k is a problem when trying to build a DLL. Since this is \
     1924  // just for debugging (likely done on *nix), just hide it from Windows           \
     1925  // -- lh, 101016 --
    19541926     if (0)  {
    19551927          static int k = 0;
     
    19651937          k++;
    19661938     }
    1967 #    endif
     1939#endif
    19681940#if 0 //ndef NDEBUG
    19691941  // Make sure everything is clean
     
    20191991           numberInside,sumInside,numberInsideLarge);
    20201992#endif
    2021      int status = factorization_->factorize(this, solveType, valuesPass);
    2022      if (status) {
    2023           handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
    2024                     << status
    2025                     << CoinMessageEol;
     1993  int status = factorization_->factorize(this, solveType, valuesPass);
     1994  if (status) {
     1995    handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
     1996      << status
     1997      << CoinMessageEol;
    20261998#ifdef CLP_USEFUL_PRINTOUT
    2027           printf("Basis singular - pivot tolerance %g\n",
    2028                 factorization_->pivotTolerance());
    2029 #endif
    2030           return -1;
    2031      } else if (!solveType) {
    2032           // Initial basis - return number of singularities
    2033           int numberSlacks = 0;
    2034           for (iRow = 0; iRow < numberRows_; iRow++) {
    2035                if (getRowStatus(iRow) == basic)
    2036                     numberSlacks++;
    2037           }
    2038           status = CoinMax(numberSlacks - totalSlacks, 0);
    2039           // special case if all slack
    2040           if (numberSlacks == numberRows_) {
    2041                status = numberRows_ + 1;
    2042           }
    2043      }
    2044 
    2045      // sparse methods
    2046      //if (factorization_->sparseThreshold()) {
    2047      // get default value
    2048      factorization_->sparseThreshold(0);
    2049      if (!(moreSpecialOptions_&1024))
    2050        factorization_->goSparse();
    2051      //}
    2052 
    2053      return status;
     1999    printf("Basis singular - pivot tolerance %g\n",
     2000      factorization_->pivotTolerance());
     2001#endif
     2002    return -1;
     2003  } else if (!solveType) {
     2004    // Initial basis - return number of singularities
     2005    int numberSlacks = 0;
     2006    for (iRow = 0; iRow < numberRows_; iRow++) {
     2007      if (getRowStatus(iRow) == basic)
     2008        numberSlacks++;
     2009    }
     2010    status = CoinMax(numberSlacks - totalSlacks, 0);
     2011    // special case if all slack
     2012    if (numberSlacks == numberRows_) {
     2013      status = numberRows_ + 1;
     2014    }
     2015  }
     2016
     2017  // sparse methods
     2018  //if (factorization_->sparseThreshold()) {
     2019  // get default value
     2020  factorization_->sparseThreshold(0);
     2021  if (!(moreSpecialOptions_ & 1024))
     2022    factorization_->goSparse();
     2023  //}
     2024
     2025  return status;
    20542026}
    20552027/*
     
    20572029   Can also decide to re-factorize
    20582030*/
    2059 int
    2060 ClpSimplex::housekeeping(double objectiveChange)
    2061 {
    2062      // save value of incoming and outgoing
    2063      double oldIn = solution_[sequenceIn_];
    2064      double oldOut = solution_[sequenceOut_];
    2065      numberIterations_++;
    2066      changeMade_++; // something has happened
    2067      // incoming variable
    2068      if (handler_->logLevel() > 7) {
    2069           //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
    2070           handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
    2071                     << directionOut_
    2072                     << directionIn_ << theta_
    2073                     << dualOut_ << dualIn_ << alpha_
    2074                     << CoinMessageEol;
    2075           if (getStatus(sequenceIn_) == isFree) {
    2076                handler_->message(CLP_SIMPLEX_FREEIN, messages_)
    2077                          << sequenceIn_
    2078                          << CoinMessageEol;
    2079           }
    2080      }
     2031int ClpSimplex::housekeeping(double objectiveChange)
     2032{
     2033  // save value of incoming and outgoing
     2034  double oldIn = solution_[sequenceIn_];
     2035  double oldOut = solution_[sequenceOut_];
     2036  numberIterations_++;
     2037  changeMade_++; // something has happened
     2038  // incoming variable
     2039  if (handler_->logLevel() > 7) {
     2040    //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
     2041    handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
     2042      << directionOut_
     2043      << directionIn_ << theta_
     2044      << dualOut_ << dualIn_ << alpha_
     2045      << CoinMessageEol;
     2046    if (getStatus(sequenceIn_) == isFree) {
     2047      handler_->message(CLP_SIMPLEX_FREEIN, messages_)
     2048        << sequenceIn_
     2049        << CoinMessageEol;
     2050    }
     2051  }
    20812052#if 0
    20822053     printf("h1 %d %d %g %g %g %g",
     
    20852056            , dualOut_, dualIn_, alpha_);
    20862057#endif
    2087      // change of incoming
    2088      char rowcol[] = {'R', 'C'};
    2089      if (pivotRow_ >= 0)
    2090           pivotVariable_[pivotRow_] = sequenceIn();
    2091      if (upper_[sequenceIn_] > 1.0e20 && lower_[sequenceIn_] < -1.0e20)
    2092           progressFlag_ |= 2; // making real progress
    2093      solution_[sequenceIn_] = valueIn_;
    2094      if (upper_[sequenceOut_] - lower_[sequenceOut_] < 1.0e-12)
    2095           progressFlag_ |= 1; // making real progress
    2096      if (sequenceIn_ != sequenceOut_) {
    2097           if (alphaAccuracy_ > 0.0) {
    2098                double value = fabs(alpha_);
    2099                if (value > 1.0)
    2100                     alphaAccuracy_ *= value;
    2101                else
    2102                     alphaAccuracy_ /= value;
    2103           }
    2104           //assert( getStatus(sequenceOut_)== basic);
    2105           setStatus(sequenceIn_, basic);
    2106           if (upper_[sequenceOut_] - lower_[sequenceOut_] > 0) {
    2107                // As Nonlinear costs may have moved bounds (to more feasible)
    2108                // Redo using value
    2109                if (fabs(valueOut_ - lower_[sequenceOut_]) < fabs(valueOut_ - upper_[sequenceOut_])) {
    2110                     // going to lower
    2111                     setStatus(sequenceOut_, atLowerBound);
    2112                     oldOut = lower_[sequenceOut_];
    2113                } else {
    2114                     // going to upper
    2115                     setStatus(sequenceOut_, atUpperBound);
    2116                     oldOut = upper_[sequenceOut_];
    2117                }
    2118           } else {
    2119                // fixed
    2120                setStatus(sequenceOut_, isFixed);
    2121           }
    2122           solution_[sequenceOut_] = valueOut_;
    2123      } else {
    2124           //if (objective_->type()<2)
    2125           //assert (fabs(theta_)>1.0e-13);
    2126           // flip from bound to bound
    2127           // As Nonlinear costs may have moved bounds (to more feasible)
    2128           // Redo using value
    2129           if (fabs(valueIn_ - lower_[sequenceIn_]) < fabs(valueIn_ - upper_[sequenceIn_])) {
    2130                // as if from upper bound
    2131                setStatus(sequenceIn_, atLowerBound);
    2132           } else {
    2133                // as if from lower bound
    2134                setStatus(sequenceIn_, atUpperBound);
    2135           }
    2136      }
    2137 
    2138      // Update hidden stuff e.g. effective RHS and gub
    2139      int invertNow=matrix_->updatePivot(this, oldIn, oldOut);
    2140      objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
    2141      if (handler_->logLevel() > 7) {
    2142           //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
    2143           handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
    2144                     << numberIterations_ << objectiveValue()
    2145                     << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
    2146                     << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
    2147           handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
    2148           handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
    2149           handler_->message() << CoinMessageEol;
    2150      }
     2058  // change of incoming
     2059  char rowcol[] = { 'R', 'C' };
     2060  if (pivotRow_ >= 0)
     2061    pivotVariable_[pivotRow_] = sequenceIn();
     2062  if (upper_[sequenceIn_] > 1.0e20 && lower_[sequenceIn_] < -1.0e20)
     2063    progressFlag_ |= 2; // making real progress
     2064  solution_[sequenceIn_] = valueIn_;
     2065  if (upper_[sequenceOut_] - lower_[sequenceOut_] < 1.0e-12)
     2066    progressFlag_ |= 1; // making real progress
     2067  if (sequenceIn_ != sequenceOut_) {
     2068    if (alphaAccuracy_ > 0.0) {
     2069      double value = fabs(alpha_);
     2070      if (value > 1.0)
     2071        alphaAccuracy_ *= value;
     2072      else
     2073        alphaAccuracy_ /= value;
     2074    }
     2075    //assert( getStatus(sequenceOut_)== basic);
     2076    setStatus(sequenceIn_, basic);
     2077    if (upper_[sequenceOut_] - lower_[sequenceOut_] > 0) {
     2078      // As Nonlinear costs may have moved bounds (to more feasible)
     2079      // Redo using value
     2080      if (fabs(valueOut_ - lower_[sequenceOut_]) < fabs(valueOut_ - upper_[sequenceOut_])) {
     2081        // going to lower
     2082        setStatus(sequenceOut_, atLowerBound);
     2083        oldOut = lower_[sequenceOut_];
     2084      } else {
     2085        // going to upper
     2086        setStatus(sequenceOut_, atUpperBound);
     2087        oldOut = upper_[sequenceOut_];
     2088      }
     2089    } else {
     2090      // fixed
     2091      setStatus(sequenceOut_, isFixed);
     2092    }
     2093    solution_[sequenceOut_] = valueOut_;
     2094  } else {
     2095    //if (objective_->type()<2)
     2096    //assert (fabs(theta_)>1.0e-13);
     2097    // flip from bound to bound
     2098    // As Nonlinear costs may have moved bounds (to more feasible)
     2099    // Redo using value
     2100    if (fabs(valueIn_ - lower_[sequenceIn_]) < fabs(valueIn_ - upper_[sequenceIn_])) {
     2101      // as if from upper bound
     2102      setStatus(sequenceIn_, atLowerBound);
     2103    } else {
     2104      // as if from lower bound
     2105      setStatus(sequenceIn_, atUpperBound);
     2106    }
     2107  }
     2108
     2109  // Update hidden stuff e.g. effective RHS and gub
     2110  int invertNow = matrix_->updatePivot(this, oldIn, oldOut);
     2111  objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
     2112  if (handler_->logLevel() > 7) {
     2113    //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
     2114    handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
     2115      << numberIterations_ << objectiveValue()
     2116      << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
     2117      << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
     2118    handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
     2119    handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
     2120    handler_->message() << CoinMessageEol;
     2121  }
    21512122#if 0
    21522123     if (numberIterations_ > 10000)
     
    21562127                 , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
    21572128#endif
    2158      if (trustedUserPointer_ && trustedUserPointer_->typeStruct == 1) {
    2159           if (algorithm_ > 0 && integerType_ && !nonLinearCost_->numberInfeasibilities()) {
    2160                if (fabs(theta_) > 1.0e-6 || !numberIterations_) {
    2161                     // For saving solutions
    2162                     typedef struct {
    2163                          int numberSolutions;
    2164                          int maximumSolutions;
    2165                          int numberColumns;
    2166                          double ** solution;
    2167                          int * numberUnsatisfied;
    2168                     } clpSolution;
    2169                     clpSolution * solution = reinterpret_cast<clpSolution *> (trustedUserPointer_->data);
    2170                     if (solution->numberSolutions == solution->maximumSolutions) {
    2171                          int n =  solution->maximumSolutions;
    2172                          int n2 = (n * 3) / 2 + 10;
    2173                          solution->maximumSolutions = n2;
    2174                          double ** temp = new double * [n2];
    2175                          for (int i = 0; i < n; i++)
    2176                               temp[i] = solution->solution[i];
    2177                          delete [] solution->solution;
    2178                          solution->solution = temp;
    2179                          int * tempN = new int [n2];
    2180                          for (int i = 0; i < n; i++)
    2181                               tempN[i] = solution->numberUnsatisfied[i];
    2182                          delete [] solution->numberUnsatisfied;
    2183                          solution->numberUnsatisfied = tempN;
    2184                     }
    2185                     assert (numberColumns_ == solution->numberColumns);
    2186                     double * sol = new double [numberColumns_];
    2187                     solution->solution[solution->numberSolutions] = sol;
    2188                     int numberFixed = 0;
    2189                     int numberUnsat = 0;
    2190                     int numberSat = 0;
    2191                     double sumUnsat = 0.0;
    2192                     double tolerance = 10.0 * primalTolerance_;
    2193                     double mostAway = 0.0;
    2194                     for (int i = 0; i < numberColumns_; i++) {
    2195                          // Save anyway
    2196                          sol[i] = columnScale_ ? solution_[i] * columnScale_[i] : solution_[i];
    2197                          // rest is optional
    2198                          if (upper_[i] > lower_[i]) {
    2199                               double value = solution_[i];
    2200                               if (value > lower_[i] + tolerance &&
    2201                                         value < upper_[i] - tolerance && integerType_[i]) {
    2202                                    // may have to modify value if scaled
    2203                                    if (columnScale_)
    2204                                         value *= columnScale_[i];
    2205                                    double closest = floor(value + 0.5);
    2206                                    // problem may be perturbed so relax test
    2207                                    if (fabs(value - closest) > 1.0e-4) {
    2208                                         numberUnsat++;
    2209                                         sumUnsat += fabs(value - closest);
    2210                                         if (mostAway < fabs(value - closest)) {
    2211                                              mostAway = fabs(value - closest);
    2212                                         }
    2213                                    } else {
    2214                                         numberSat++;
    2215                                    }
    2216                               } else {
    2217                                    numberSat++;
    2218                               }
    2219                          } else {
    2220                               numberFixed++;
    2221                          }
    2222                     }
    2223                     solution->numberUnsatisfied[solution->numberSolutions++] = numberUnsat;
    2224                     COIN_DETAIL_PRINT(printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n",
    2225                                              numberIterations_, numberUnsat, sumUnsat, mostAway, numberFixed, numberSat));
    2226                }
    2227           }
    2228      }
    2229      if (hitMaximumIterations())
    2230           return 2;
     2129  if (trustedUserPointer_ && trustedUserPointer_->typeStruct == 1) {
     2130    if (algorithm_ > 0 && integerType_ && !nonLinearCost_->numberInfeasibilities()) {
     2131      if (fabs(theta_) > 1.0e-6 || !numberIterations_) {
     2132        // For saving solutions
     2133        typedef struct {
     2134          int numberSolutions;
     2135          int maximumSolutions;
     2136          int numberColumns;
     2137          double **solution;
     2138          int *numberUnsatisfied;
     2139        } clpSolution;
     2140        clpSolution *solution = reinterpret_cast< clpSolution * >(trustedUserPointer_->data);
     2141        if (solution->numberSolutions == solution->maximumSolutions) {
     2142          int n = solution->maximumSolutions;
     2143          int n2 = (n * 3) / 2 + 10;
     2144          solution->maximumSolutions = n2;
     2145          double **temp = new double *[n2];
     2146          for (int i = 0; i < n; i++)
     2147            temp[i] = solution->solution[i];
     2148          delete[] solution->solution;
     2149          solution->solution = temp;
     2150          int *tempN = new int[n2];
     2151          for (int i = 0; i < n; i++)
     2152            tempN[i] = solution->numberUnsatisfied[i];
     2153          delete[] solution->numberUnsatisfied;
     2154          solution->numberUnsatisfied = tempN;
     2155        }
     2156        assert(numberColumns_ == solution->numberColumns);
     2157        double *sol = new double[numberColumns_];
     2158        solution->solution[solution->numberSolutions] = sol;
     2159        int numberFixed = 0;
     2160        int numberUnsat = 0;
     2161        int numberSat = 0;
     2162        double sumUnsat = 0.0;
     2163        double tolerance = 10.0 * primalTolerance_;
     2164        double mostAway = 0.0;
     2165        for (int i = 0; i < numberColumns_; i++) {
     2166          // Save anyway
     2167          sol[i] = columnScale_ ? solution_[i] * columnScale_[i] : solution_[i];
     2168          // rest is optional
     2169          if (upper_[i] > lower_[i]) {
     2170            double value = solution_[i];
     2171            if (value > lower_[i] + tolerance && value < upper_[i] - tolerance && integerType_[i]) {
     2172              // may have to modify value if scaled
     2173              if (columnScale_)
     2174                value *= columnScale_[i];
     2175              double closest = floor(value + 0.5);
     2176              // problem may be perturbed so relax test
     2177              if (fabs(value - closest) > 1.0e-4) {
     2178                numberUnsat++;
     2179                sumUnsat += fabs(value - closest);
     2180                if (mostAway < fabs(value - closest)) {
     2181                  mostAway = fabs(value - closest);
     2182                }
     2183              } else {
     2184                numberSat++;
     2185              }
     2186            } else {
     2187              numberSat++;
     2188            }
     2189          } else {
     2190            numberFixed++;
     2191          }
     2192        }
     2193        solution->numberUnsatisfied[solution->numberSolutions++] = numberUnsat;
     2194        COIN_DETAIL_PRINT(printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n",
     2195          numberIterations_, numberUnsat, sumUnsat, mostAway, numberFixed, numberSat));
     2196      }
     2197    }
     2198  }
     2199  if (hitMaximumIterations())
     2200    return 2;
    22312201#if 1
    2232      //if (numberIterations_>14000)
    2233      //handler_->setLogLevel(63);
    2234      //if (numberIterations_>24000)
    2235      //exit(77);
    2236      // check for small cycles
    2237      int in = sequenceIn_;
    2238      int out = sequenceOut_;
    2239      matrix_->correctSequence(this, in, out);
    2240      int cycle = progress_.cycle(in, out,
    2241                                  directionIn_, directionOut_);
    2242      if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
    2243           //if (cycle>0) {
    2244           if (handler_->logLevel() >= 63)
    2245                printf("Cycle of %d\n", cycle);
    2246           // reset
    2247           progress_.startCheck();
    2248           double random = randomNumberGenerator_.randomDouble();
    2249           int extra = static_cast<int> (9.999 * random);
    2250           int off[] = {1, 1, 1, 1, 2, 2, 2, 3, 3, 4};
    2251           if (factorization_->pivots() > cycle) {
    2252                forceFactorization_ = CoinMax(1, cycle - off[extra]);
    2253           } else {
    2254             /* need to reject something
     2202  //if (numberIterations_>14000)
     2203  //handler_->setLogLevel(63);
     2204  //if (numberIterations_>24000)
     2205  //exit(77);
     2206  // check for small cycles
     2207  int in = sequenceIn_;
     2208  int out = sequenceOut_;
     2209  matrix_->correctSequence(this, in, out);
     2210  int cycle = progress_.cycle(in, out,
     2211    directionIn_, directionOut_);
     2212  if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
     2213    //if (cycle>0) {
     2214    if (handler_->logLevel() >= 63)
     2215      printf("Cycle of %d\n", cycle);
     2216    // reset
     2217    progress_.startCheck();
     2218    double random = randomNumberGenerator_.randomDouble();
     2219    int extra = static_cast< int >(9.999 * random);
     2220    int off[] = { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 };
     2221    if (factorization_->pivots() > cycle) {
     2222      forceFactorization_ = CoinMax(1, cycle - off[extra]);
     2223    } else {
     2224      /* need to reject something
    22552225               should be better if don't reject incoming
    22562226               as it is in basis */
    2257                int iSequence;
    2258                //if (algorithm_ > 0)
    2259                //   iSequence = sequenceIn_;
    2260                //else
    2261                     iSequence = sequenceOut_;
    2262                char x = isColumn(iSequence) ? 'C' : 'R';
    2263                if (handler_->logLevel() >= 63)
    2264                     handler_->message(CLP_SIMPLEX_FLAG, messages_)
    2265                               << x << sequenceWithin(iSequence)
    2266                               << CoinMessageEol;
    2267                setFlagged(iSequence);
    2268                //printf("flagging %d\n",iSequence);
    2269           }
    2270           return 1;
    2271      }
    2272 #endif
    2273      // only time to re-factorize if one before real time
    2274      // this is so user won't be surprised that maximumPivots has exact meaning
    2275      int numberPivots = factorization_->pivots();
    2276      int maximumPivots = factorization_->maximumPivots();
    2277      int numberDense = factorization_->numberDense();
    2278      bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 >
    2279                         2 * maximumIterations());
    2280      if (numberPivots == maximumPivots ||
    2281                maximumPivots < 2) {
    2282           // If dense then increase
    2283           if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots
    2284               && false) {
    2285                factorization_->maximumPivots(numberDense);
    2286                dualRowPivot_->maximumPivotsChanged();
    2287                primalColumnPivot_->maximumPivotsChanged();
    2288                // and redo arrays
    2289                for (int iRow = 0; iRow < 4; iRow++) {
    2290                     int length = rowArray_[iRow]->capacity() + numberDense - maximumPivots;
    2291                     rowArray_[iRow]->reserve(length);
    2292                }
    2293           }
    2294 #if CLP_FACTORIZATION_NEW_TIMING>1
    2295           factorization_->statsRefactor('M');
    2296 #endif
    2297           return 1;
    2298      } else if ((factorization_->timeToRefactorize() && !dontInvert)
    2299                 ||invertNow) {
    2300           //printf("ret after %d pivots\n",factorization_->pivots());
    2301 #if CLP_FACTORIZATION_NEW_TIMING>1
    2302           factorization_->statsRefactor('T');
    2303 #endif
    2304           return 1;
    2305      } else if (forceFactorization_ > 0 &&
    2306                 factorization_->pivots() == forceFactorization_) {
    2307           // relax
    2308           forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
    2309           if (forceFactorization_ > factorization_->maximumPivots())
    2310                forceFactorization_ = -1; //off
    2311 #if CLP_FACTORIZATION_NEW_TIMING>1
    2312           factorization_->statsRefactor('F');
    2313 #endif
    2314           return 1;
    2315      } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
    2316           double random = randomNumberGenerator_.randomDouble();
    2317           while (random<0.45)
    2318             random *= 2.0;
    2319           int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
    2320           if (factorization_->pivots() >= random * maxNumber) {
    2321                return 1;
    2322           } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) &&
    2323                      numberIterations_ < 1001000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
    2324                return 1;
    2325           } else {
    2326                // carry on iterating
    2327                return 0;
    2328           }
    2329      } else {
    2330           // carry on iterating
    2331           return 0;
    2332      }
     2227      int iSequence;
     2228      //if (algorithm_ > 0)
     2229      //   iSequence = sequenceIn_;
     2230      //else
     2231      iSequence = sequenceOut_;
     2232      char x = isColumn(iSequence) ? 'C' : 'R';
     2233      if (handler_->logLevel() >= 63)
     2234        handler_->message(CLP_SIMPLEX_FLAG, messages_)
     2235          << x << sequenceWithin(iSequence)
     2236          << CoinMessageEol;
     2237      setFlagged(iSequence);
     2238      //printf("flagging %d\n",iSequence);
     2239    }
     2240    return 1;
     2241  }
     2242#endif
     2243  // only time to re-factorize if one before real time
     2244  // this is so user won't be surprised that maximumPivots has exact meaning
     2245  int numberPivots = factorization_->pivots();
     2246  int maximumPivots = factorization_->maximumPivots();
     2247  int numberDense = factorization_->numberDense();
     2248  bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 > 2 * maximumIterations());
     2249  if (numberPivots == maximumPivots || maximumPivots < 2) {
     2250    // If dense then increase
     2251    if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots
     2252      && false) {
     2253      factorization_->maximumPivots(numberDense);
     2254      dualRowPivot_->maximumPivotsChanged();
     2255      primalColumnPivot_->maximumPivotsChanged();
     2256      // and redo arrays
     2257      for (int iRow = 0; iRow < 4; iRow++) {
     2258        int length = rowArray_[iRow]->capacity() + numberDense - maximumPivots;
     2259        rowArray_[iRow]->reserve(length);
     2260      }
     2261    }
     2262#if CLP_FACTORIZATION_NEW_TIMING > 1
     2263    factorization_->statsRefactor('M');
     2264#endif
     2265    return 1;
     2266  } else if ((factorization_->timeToRefactorize() && !dontInvert)
     2267    || invertNow) {
     2268    //printf("ret after %d pivots\n",factorization_->pivots());
     2269#if CLP_FACTORIZATION_NEW_TIMING > 1
     2270    factorization_->statsRefactor('T');
     2271#endif
     2272    return 1;
     2273  } else if (forceFactorization_ > 0 && factorization_->pivots() == forceFactorization_) {
     2274    // relax
     2275    forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
     2276    if (forceFactorization_ > factorization_->maximumPivots())
     2277      forceFactorization_ = -1; //off
     2278#if CLP_FACTORIZATION_NEW_TIMING > 1
     2279    factorization_->statsRefactor('F');
     2280#endif
     2281    return 1;
     2282  } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type() < 15) {
     2283    double random = randomNumberGenerator_.randomDouble();
     2284    while (random < 0.45)
     2285      random *= 2.0;
     2286    int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
     2287    if (factorization_->pivots() >= random * maxNumber) {
     2288      return 1;
     2289    } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && numberIterations_ < 1001000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
     2290      return 1;
     2291    } else {
     2292      // carry on iterating
     2293      return 0;
     2294    }
     2295  } else {
     2296    // carry on iterating
     2297    return 0;
     2298  }
    23332299}
    23342300// Copy constructor.
    2335 ClpSimplex::ClpSimplex(const ClpSimplex &rhs, int scalingMode) :
    2336      ClpModel(rhs, scalingMode),
    2337      bestPossibleImprovement_(0.0),
    2338      zeroTolerance_(1.0e-13),
    2339      columnPrimalSequence_(-2),
    2340      rowPrimalSequence_(-2),
    2341      bestObjectiveValue_(rhs.bestObjectiveValue_),
    2342      moreSpecialOptions_(2),
    2343      baseIteration_(0),
    2344      vectorMode_(0),
    2345      primalToleranceToGetOptimal_(-1.0),
    2346      largeValue_(1.0e15),
    2347      largestPrimalError_(0.0),
    2348      largestDualError_(0.0),
    2349      alphaAccuracy_(-1.0),
    2350      dualBound_(1.0e10),
    2351      alpha_(0.0),
    2352      theta_(0.0),
    2353      lowerIn_(0.0),
    2354      valueIn_(0.0),
    2355      upperIn_(-COIN_DBL_MAX),
    2356      dualIn_(0.0),
    2357      lowerOut_(-1),
    2358      valueOut_(-1),
    2359      upperOut_(-1),
    2360      dualOut_(-1),
    2361      dualTolerance_(1.0e-7),
    2362      primalTolerance_(1.0e-7),
    2363      sumDualInfeasibilities_(0.0),
    2364      sumPrimalInfeasibilities_(0.0),
    2365      infeasibilityCost_(1.0e10),
    2366      sumOfRelaxedDualInfeasibilities_(0.0),
    2367      sumOfRelaxedPrimalInfeasibilities_(0.0),
    2368      acceptablePivot_(1.0e-8),
    2369      lower_(NULL),
    2370      rowLowerWork_(NULL),
    2371      columnLowerWork_(NULL),
    2372      upper_(NULL),
    2373      rowUpperWork_(NULL),
    2374      columnUpperWork_(NULL),
    2375      cost_(NULL),
    2376      rowObjectiveWork_(NULL),
    2377      objectiveWork_(NULL),
    2378      sequenceIn_(-1),
    2379      directionIn_(-1),
    2380      sequenceOut_(-1),
    2381      directionOut_(-1),
    2382      pivotRow_(-1),
    2383      lastGoodIteration_(-100),
    2384      dj_(NULL),
    2385      rowReducedCost_(NULL),
    2386      reducedCostWork_(NULL),
    2387      solution_(NULL),
    2388      rowActivityWork_(NULL),
    2389      columnActivityWork_(NULL),
    2390      numberDualInfeasibilities_(0),
    2391      numberDualInfeasibilitiesWithoutFree_(0),
    2392      numberPrimalInfeasibilities_(100),
    2393      numberRefinements_(0),
    2394      pivotVariable_(NULL),
    2395      factorization_(NULL),
    2396      savedSolution_(NULL),
    2397      numberTimesOptimal_(0),
    2398      disasterArea_(NULL),
    2399      changeMade_(1),
    2400      algorithm_(0),
    2401      forceFactorization_(-1),
    2402      perturbation_(100),
    2403      nonLinearCost_(NULL),
    2404      lastBadIteration_(-999999),
    2405      lastFlaggedIteration_(-999999),
    2406      numberFake_(0),
    2407      numberChanged_(0),
    2408      progressFlag_(0),
    2409      firstFree_(-1),
    2410      numberExtraRows_(0),
    2411      maximumBasic_(0),
    2412      dontFactorizePivots_(0),
    2413      incomingInfeasibility_(1.0),
    2414      allowedInfeasibility_(10.0),
    2415      automaticScale_(0),
    2416      maximumPerturbationSize_(0),
    2417      perturbationArray_(NULL),
    2418     baseModel_(NULL)
     2301ClpSimplex::ClpSimplex(const ClpSimplex &rhs, int scalingMode)
     2302  : ClpModel(rhs, scalingMode)
     2303  , bestPossibleImprovement_(0.0)
     2304  , zeroTolerance_(1.0e-13)
     2305  , columnPrimalSequence_(-2)
     2306  , rowPrimalSequence_(-2)
     2307  , bestObjectiveValue_(rhs.bestObjectiveValue_)
     2308  , moreSpecialOptions_(2)
     2309  , baseIteration_(0)
     2310  , vectorMode_(0)
     2311  , primalToleranceToGetOptimal_(-1.0)
     2312  , largeValue_(1.0e15)
     2313  , largestPrimalError_(0.0)
     2314  , largestDualError_(0.0)
     2315  , alphaAccuracy_(-1.0)
     2316  , dualBound_(1.0e10)
     2317  , alpha_(0.0)
     2318  , theta_(0.0)
     2319  , lowerIn_(0.0)
     2320  , valueIn_(0.0)
     2321  , upperIn_(-COIN_DBL_MAX)
     2322  , dualIn_(0.0)
     2323  , lowerOut_(-1)
     2324  , valueOut_(-1)
     2325  , upperOut_(-1)
     2326  , dualOut_(-1)
     2327  , dualTolerance_(1.0e-7)
     2328  , primalTolerance_(1.0e-7)
     2329  , sumDualInfeasibilities_(0.0)
     2330  , sumPrimalInfeasibilities_(0.0)
     2331  , infeasibilityCost_(1.0e10)
     2332  , sumOfRelaxedDualInfeasibilities_(0.0)
     2333  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     2334  , acceptablePivot_(1.0e-8)
     2335  , lower_(NULL)
     2336  , rowLowerWork_(NULL)
     2337  , columnLowerWork_(NULL)
     2338  , upper_(NULL)
     2339  , rowUpperWork_(NULL)
     2340  , columnUpperWork_(NULL)
     2341  , cost_(NULL)
     2342  , rowObjectiveWork_(NULL)
     2343  , objectiveWork_(NULL)
     2344  , sequenceIn_(-1)
     2345  , directionIn_(-1)
     2346  , sequenceOut_(-1)
     2347  , directionOut_(-1)
     2348  , pivotRow_(-1)
     2349  , lastGoodIteration_(-100)
     2350  , dj_(NULL)
     2351  , rowReducedCost_(NULL)
     2352  , reducedCostWork_(NULL)
     2353  , solution_(NULL)
     2354  , rowActivityWork_(NULL)
     2355  , columnActivityWork_(NULL)
     2356  , numberDualInfeasibilities_(0)
     2357  , numberDualInfeasibilitiesWithoutFree_(0)
     2358  , numberPrimalInfeasibilities_(100)
     2359  , numberRefinements_(0)
     2360  , pivotVariable_(NULL)
     2361  , factorization_(NULL)
     2362  , savedSolution_(NULL)
     2363  , numberTimesOptimal_(0)
     2364  , disasterArea_(NULL)
     2365  , changeMade_(1)
     2366  , algorithm_(0)
     2367  , forceFactorization_(-1)
     2368  , perturbation_(100)
     2369  , nonLinearCost_(NULL)
     2370  , lastBadIteration_(-999999)
     2371  , lastFlaggedIteration_(-999999)
     2372  , numberFake_(0)
     2373  , numberChanged_(0)
     2374  , progressFlag_(0)
     2375  , firstFree_(-1)
     2376  , numberExtraRows_(0)
     2377  , maximumBasic_(0)
     2378  , dontFactorizePivots_(0)
     2379  , incomingInfeasibility_(1.0)
     2380  , allowedInfeasibility_(10.0)
     2381  , automaticScale_(0)
     2382  , maximumPerturbationSize_(0)
     2383  , perturbationArray_(NULL)
     2384  , baseModel_(NULL)
    24192385#ifdef ABC_INHERIT
    2420      ,abcSimplex_(NULL),
    2421     abcState_(0)
    2422 #endif
    2423 {
    2424      int i;
    2425      for (i = 0; i < 6; i++) {
    2426           rowArray_[i] = NULL;
    2427           columnArray_[i] = NULL;
    2428      }
    2429      for (i = 0; i < 4; i++) {
    2430           spareIntArray_[i] = 0;
    2431           spareDoubleArray_[i] = 0.0;
    2432      }
    2433      saveStatus_ = NULL;
    2434      factorization_ = NULL;
    2435      dualRowPivot_ = NULL;
    2436      primalColumnPivot_ = NULL;
    2437      gutsOfDelete(0);
    2438      delete nonLinearCost_;
    2439      nonLinearCost_ = NULL;
    2440      gutsOfCopy(rhs);
    2441      solveType_ = 1; // say simplex based life form
     2386  , abcSimplex_(NULL)
     2387  , abcState_(0)
     2388#endif
     2389{
     2390  int i;
     2391  for (i = 0; i < 6; i++) {
     2392    rowArray_[i] = NULL;
     2393    columnArray_[i] = NULL;
     2394  }
     2395  for (i = 0; i < 4; i++) {
     2396    spareIntArray_[i] = 0;
     2397    spareDoubleArray_[i] = 0.0;
     2398  }
     2399  saveStatus_ = NULL;
     2400  factorization_ = NULL;
     2401  dualRowPivot_ = NULL;
     2402  primalColumnPivot_ = NULL;
     2403  gutsOfDelete(0);
     2404  delete nonLinearCost_;
     2405  nonLinearCost_ = NULL;
     2406  gutsOfCopy(rhs);
     2407  solveType_ = 1; // say simplex based life form
    24422408}
    24432409// Copy constructor from model
    2444 ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
    2445      ClpModel(rhs, scalingMode),
    2446      bestPossibleImprovement_(0.0),
    2447      zeroTolerance_(1.0e-13),
    2448      columnPrimalSequence_(-2),
    2449      rowPrimalSequence_(-2),
    2450      bestObjectiveValue_(-COIN_DBL_MAX),
    2451      moreSpecialOptions_(2),
    2452      baseIteration_(0),
    2453      vectorMode_(0),
    2454      primalToleranceToGetOptimal_(-1.0),
    2455      largeValue_(1.0e15),
    2456      largestPrimalError_(0.0),
    2457      largestDualError_(0.0),
    2458      alphaAccuracy_(-1.0),
    2459      dualBound_(1.0e10),
    2460      alpha_(0.0),
    2461      theta_(0.0),
    2462      lowerIn_(0.0),
    2463      valueIn_(0.0),
    2464      upperIn_(-COIN_DBL_MAX),
    2465      dualIn_(0.0),
    2466      lowerOut_(-1),
    2467      valueOut_(-1),
    2468      upperOut_(-1),
    2469      dualOut_(-1),
    2470      dualTolerance_(1.0e-7),
    2471      primalTolerance_(1.0e-7),
    2472      sumDualInfeasibilities_(0.0),
    2473      sumPrimalInfeasibilities_(0.0),
    2474      infeasibilityCost_(1.0e10),
    2475      sumOfRelaxedDualInfeasibilities_(0.0),
    2476      sumOfRelaxedPrimalInfeasibilities_(0.0),
    2477      acceptablePivot_(1.0e-8),
    2478      lower_(NULL),
    2479      rowLowerWork_(NULL),
    2480      columnLowerWork_(NULL),
    2481      upper_(NULL),
    2482      rowUpperWork_(NULL),
    2483      columnUpperWork_(NULL),
    2484      cost_(NULL),
    2485      rowObjectiveWork_(NULL),
    2486      objectiveWork_(NULL),
    2487      sequenceIn_(-1),
    2488      directionIn_(-1),
    2489      sequenceOut_(-1),
    2490      directionOut_(-1),
    2491      pivotRow_(-1),
    2492      lastGoodIteration_(-100),
    2493      dj_(NULL),
    2494      rowReducedCost_(NULL),
    2495      reducedCostWork_(NULL),
    2496      solution_(NULL),
    2497      rowActivityWork_(NULL),
    2498      columnActivityWork_(NULL),
    2499      numberDualInfeasibilities_(0),
    2500      numberDualInfeasibilitiesWithoutFree_(0),
    2501      numberPrimalInfeasibilities_(100),
    2502      numberRefinements_(0),
    2503      pivotVariable_(NULL),
    2504      factorization_(NULL),
    2505      savedSolution_(NULL),
    2506      numberTimesOptimal_(0),
    2507      disasterArea_(NULL),
    2508      changeMade_(1),
    2509      algorithm_(0),
    2510      forceFactorization_(-1),
    2511      perturbation_(100),
    2512      nonLinearCost_(NULL),
    2513      lastBadIteration_(-999999),
    2514      lastFlaggedIteration_(-999999),
    2515      numberFake_(0),
    2516      numberChanged_(0),
    2517      progressFlag_(0),
    2518      firstFree_(-1),
    2519      numberExtraRows_(0),
    2520      maximumBasic_(0),
    2521      dontFactorizePivots_(0),
    2522      incomingInfeasibility_(1.0),
    2523      allowedInfeasibility_(10.0),
    2524      automaticScale_(0),
    2525      maximumPerturbationSize_(0),
    2526      perturbationArray_(NULL),
    2527     baseModel_(NULL)
     2410ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode)
     2411  : ClpModel(rhs, scalingMode)
     2412  , bestPossibleImprovement_(0.0)
     2413  , zeroTolerance_(1.0e-13)
     2414  , columnPrimalSequence_(-2)
     2415  , rowPrimalSequence_(-2)
     2416  , bestObjectiveValue_(-COIN_DBL_MAX)
     2417  , moreSpecialOptions_(2)
     2418  , baseIteration_(0)
     2419  , vectorMode_(0)
     2420  , primalToleranceToGetOptimal_(-1.0)
     2421  , largeValue_(1.0e15)
     2422  , largestPrimalError_(0.0)
     2423  , largestDualError_(0.0)
     2424  , alphaAccuracy_(-1.0)
     2425  , dualBound_(1.0e10)
     2426  , alpha_(0.0)
     2427  , theta_(0.0)
     2428  , lowerIn_(0.0)
     2429  , valueIn_(0.0)
     2430  , upperIn_(-COIN_DBL_MAX)
     2431  , dualIn_(0.0)
     2432  , lowerOut_(-1)
     2433  , valueOut_(-1)
     2434  , upperOut_(-1)
     2435  , dualOut_(-1)
     2436  , dualTolerance_(1.0e-7)
     2437  , primalTolerance_(1.0e-7)
     2438  , sumDualInfeasibilities_(0.0)
     2439  , sumPrimalInfeasibilities_(0.0)
     2440  , infeasibilityCost_(1.0e10)
     2441  , sumOfRelaxedDualInfeasibilities_(0.0)
     2442  , sumOfRelaxedPrimalInfeasibilities_(0.0)
     2443  , acceptablePivot_(1.0e-8)
     2444  , lower_(NULL)
     2445  , rowLowerWork_(NULL)
     2446  , columnLowerWork_(NULL)
     2447  , upper_(NULL)
     2448  , rowUpperWork_(NULL)
     2449  , columnUpperWork_(NULL)
     2450  , cost_(NULL)
     2451  , rowObjectiveWork_(NULL)
     2452  , objectiveWork_(NULL)
     2453  , sequenceIn_(-1)
     2454  , directionIn_(-1)
     2455  , sequenceOut_(-1)
     2456  , directionOut_(-1)
     2457  , pivotRow_(-1)
     2458  , lastGoodIteration_(-100)
     2459  , dj_(NULL)
     2460  , rowReducedCost_(NULL)
     2461  , reducedCostWork_(NULL)
     2462  , solution_(NULL)
     2463  , rowActivityWork_(NULL)
     2464  , columnActivityWork_(NULL)
     2465  , numberDualInfeasibilities_(0)
     2466  , numberDualInfeasibilitiesWithoutFree_(0)
     2467  , numberPrimalInfeasibilities_(100)
     2468  , numberRefinements_(0)
     2469  , pivotVariable_(NULL)
     2470  , factorization_(NULL)
     2471  , savedSolution_(NULL)
     2472  , numberTimesOptimal_(0)
     2473  , disasterArea_(NULL)
     2474  , changeMade_(1)
     2475  , algorithm_(0)
     2476  , forceFactorization_(-1)
     2477  , perturbation_(100)
     2478  , nonLinearCost_(NULL)
     2479  , lastBadIteration_(-999999)
     2480  , lastFlaggedIteration_(-999999)
     2481  , numberFake_(0)
     2482  , numberChanged_(0)
     2483  , progressFlag_(0)
     2484  , firstFree_(-1)
     2485  , numberExtraRows_(0)
     2486  , maximumBasic_(0)
     2487  , dontFactorizePivots_(0)
     2488  , incomingInfeasibility_(1.0)
     2489  , allowedInfeasibility_(10.0)
     2490  , automaticScale_(0)
     2491  , maximumPerturbationSize_(0)
     2492  , perturbationArray_(NULL)
     2493  , baseModel_(NULL)
    25282494#ifdef ABC_INHERIT
    2529      ,abcSimplex_(NULL),
    2530      abcState_(0)
    2531 #endif
    2532 {
    2533      int i;
    2534      for (i = 0; i < 6; i++) {
    2535           rowArray_[i] = NULL;
    2536           columnArray_[i] = NULL;
    2537      }
    2538      for (i = 0; i < 4; i++) {
    2539           spareIntArray_[i] = 0;
    2540           spareDoubleArray_[i] = 0.0;
    2541      }
    2542      saveStatus_ = NULL;
    2543      // get an empty factorization so we can set tolerances etc
    2544      getEmptyFactorization();
    2545      // say Steepest pricing
    2546      dualRowPivot_ = new ClpDualRowSteepest();
    2547      // say Steepest pricing
    2548      primalColumnPivot_ = new ClpPrimalColumnSteepest();
    2549      solveType_ = 1; // say simplex based life form
    2550 
     2495  , abcSimplex_(NULL)
     2496  , abcState_(0)
     2497#endif
     2498{
     2499  int i;
     2500  for (i = 0; i < 6; i++) {
     2501    rowArray_[i] = NULL;
     2502    columnArray_[i] = NULL;
     2503  }
     2504  for (i = 0; i < 4; i++) {
     2505    spareIntArray_[i] = 0;
     2506    spareDoubleArray_[i] = 0.0;
     2507  }
     2508  saveStatus_ = NULL;
     2509  // get an empty factorization so we can set tolerances etc
     2510  getEmptyFactorization();
     2511  // say Steepest pricing
     2512  dualRowPivot_ = new ClpDualRowSteepest();
     2513  // say Steepest pricing
     2514  primalColumnPivot_ = new ClpPrimalColumnSteepest();
     2515  solveType_ = 1; // say simplex based life form
    25512516}
    25522517// Assignment operator. This copies the data
    25532518ClpSimplex &
    2554 ClpSimplex::operator=(const ClpSimplex & rhs)
    2555 {
    2556      if (this != &rhs) {
    2557           gutsOfDelete(0);
    2558           delete nonLinearCost_;
    2559           nonLinearCost_ = NULL;
    2560           ClpModel::operator=(rhs);
    2561           gutsOfCopy(rhs);
    2562      }
    2563      return *this;
    2564 }
    2565 void
    2566 ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
    2567 {
    2568      assert (numberRows_ == rhs.numberRows_);
    2569      assert (numberColumns_ == rhs.numberColumns_);
    2570      numberExtraRows_ = rhs.numberExtraRows_;
    2571      maximumBasic_ = rhs.maximumBasic_;
    2572      dontFactorizePivots_ = rhs.dontFactorizePivots_;
    2573      int numberRows2 = numberRows_ + numberExtraRows_;
    2574      moreSpecialOptions_ = rhs.moreSpecialOptions_;
    2575      if ((whatsChanged_ & 1) != 0) {
    2576           int numberTotal = numberColumns_ + numberRows2;
    2577           if ((specialOptions_ & 65536) != 0 && maximumRows_ >= 0) {
    2578                assert (maximumInternalRows_ >= numberRows2);
    2579                assert (maximumInternalColumns_ >= numberColumns_);
    2580                numberTotal = 2 * (maximumInternalColumns_ + maximumInternalRows_);
    2581           }
    2582           lower_ = ClpCopyOfArray(rhs.lower_, numberTotal);
    2583           rowLowerWork_ = lower_ + numberColumns_;
    2584           columnLowerWork_ = lower_;
    2585           upper_ = ClpCopyOfArray(rhs.upper_, numberTotal);
    2586           rowUpperWork_ = upper_ + numberColumns_;
    2587           columnUpperWork_ = upper_;
    2588           cost_ = ClpCopyOfArray(rhs.cost_, numberTotal);
    2589           objectiveWork_ = cost_;
    2590           rowObjectiveWork_ = cost_ + numberColumns_;
    2591           dj_ = ClpCopyOfArray(rhs.dj_, numberTotal);
    2592           if (dj_) {
    2593                reducedCostWork_ = dj_;
    2594                rowReducedCost_ = dj_ + numberColumns_;
    2595           }
    2596           solution_ = ClpCopyOfArray(rhs.solution_, numberTotal);
    2597           if (solution_) {
    2598                columnActivityWork_ = solution_;
    2599                rowActivityWork_ = solution_ + numberColumns_;
    2600           }
    2601           if (rhs.pivotVariable_) {
    2602                pivotVariable_ = new int[numberRows2];
    2603                CoinMemcpyN ( rhs.pivotVariable_, numberRows2 , pivotVariable_);
    2604           } else {
    2605                pivotVariable_ = NULL;
    2606           }
    2607           savedSolution_ = ClpCopyOfArray(rhs.savedSolution_, numberTotal);
    2608           int i;
    2609           for (i = 0; i < 6; i++) {
    2610                rowArray_[i] = NULL;
    2611                if (rhs.rowArray_[i])
    2612                     rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
    2613                columnArray_[i] = NULL;
    2614                if (rhs.columnArray_[i])
    2615                     columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
    2616           }
    2617           if (rhs.saveStatus_) {
    2618                saveStatus_ = ClpCopyOfArray( rhs.saveStatus_, numberTotal);
    2619           }
    2620      } else {
    2621           lower_ = NULL;
    2622           rowLowerWork_ = NULL;
    2623           columnLowerWork_ = NULL;
    2624           upper_ = NULL;
    2625           rowUpperWork_ = NULL;
    2626           columnUpperWork_ = NULL;
    2627           cost_ = NULL;
    2628           objectiveWork_ = NULL;
    2629           rowObjectiveWork_ = NULL;
    2630           dj_ = NULL;
    2631           reducedCostWork_ = NULL;
    2632           rowReducedCost_ = NULL;
    2633           solution_ = NULL;
    2634           columnActivityWork_ = NULL;
    2635           rowActivityWork_ = NULL;
    2636           pivotVariable_ = NULL;
    2637           savedSolution_ = NULL;
    2638           int i;
    2639           for (i = 0; i < 6; i++) {
    2640                rowArray_[i] = NULL;
    2641                columnArray_[i] = NULL;
    2642           }
    2643           saveStatus_ = NULL;
    2644      }
    2645      if (rhs.factorization_) {
    2646           setFactorization(*rhs.factorization_);
    2647      } else {
    2648           delete factorization_;
    2649           factorization_ = NULL;
    2650      }
    2651      bestPossibleImprovement_ = rhs.bestPossibleImprovement_;
    2652      columnPrimalSequence_ = rhs.columnPrimalSequence_;
    2653      zeroTolerance_ = rhs.zeroTolerance_;
    2654      rowPrimalSequence_ = rhs.rowPrimalSequence_;
    2655      bestObjectiveValue_ = rhs.bestObjectiveValue_;
    2656      baseIteration_ = rhs.baseIteration_;
    2657      vectorMode_ = rhs.vectorMode_;
    2658      primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
    2659      largeValue_ = rhs.largeValue_;
    2660      largestPrimalError_ = rhs.largestPrimalError_;
    2661      largestDualError_ = rhs.largestDualError_;
    2662      alphaAccuracy_ = rhs.alphaAccuracy_;
    2663      dualBound_ = rhs.dualBound_;
    2664      alpha_ = rhs.alpha_;
    2665      theta_ = rhs.theta_;
    2666      lowerIn_ = rhs.lowerIn_;
    2667      valueIn_ = rhs.valueIn_;
    2668      upperIn_ = rhs.upperIn_;
    2669      dualIn_ = rhs.dualIn_;
    2670      sequenceIn_ = rhs.sequenceIn_;
    2671      directionIn_ = rhs.directionIn_;
    2672      lowerOut_ = rhs.lowerOut_;
    2673      valueOut_ = rhs.valueOut_;
    2674      upperOut_ = rhs.upperOut_;
    2675      dualOut_ = rhs.dualOut_;
    2676      sequenceOut_ = rhs.sequenceOut_;
    2677      directionOut_ = rhs.directionOut_;
    2678      pivotRow_ = rhs.pivotRow_;
    2679      lastGoodIteration_ = rhs.lastGoodIteration_;
    2680      numberRefinements_ = rhs.numberRefinements_;
    2681      dualTolerance_ = rhs.dualTolerance_;
    2682      primalTolerance_ = rhs.primalTolerance_;
    2683      sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
    2684      numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    2685      numberDualInfeasibilitiesWithoutFree_ =
    2686           rhs.numberDualInfeasibilitiesWithoutFree_;
    2687      sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    2688      numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    2689      dualRowPivot_ = rhs.dualRowPivot_->clone(true);
    2690      dualRowPivot_->setModel(this);
    2691      primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
    2692      primalColumnPivot_->setModel(this);
    2693      numberTimesOptimal_ = rhs.numberTimesOptimal_;
    2694      disasterArea_ = NULL;
    2695      changeMade_ = rhs.changeMade_;
    2696      algorithm_ = rhs.algorithm_;
    2697      forceFactorization_ = rhs.forceFactorization_;
    2698      perturbation_ = rhs.perturbation_;
    2699      infeasibilityCost_ = rhs.infeasibilityCost_;
    2700      lastBadIteration_ = rhs.lastBadIteration_;
    2701      lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
    2702      numberFake_ = rhs.numberFake_;
    2703      numberChanged_ = rhs.numberChanged_;
    2704      progressFlag_ = rhs.progressFlag_;
    2705      firstFree_ = rhs.firstFree_;
    2706      incomingInfeasibility_ = rhs.incomingInfeasibility_;
    2707      allowedInfeasibility_ = rhs.allowedInfeasibility_;
    2708      automaticScale_ = rhs.automaticScale_;
     2519ClpSimplex::operator=(const ClpSimplex &rhs)
     2520{
     2521  if (this != &rhs) {
     2522    gutsOfDelete(0);
     2523    delete nonLinearCost_;
     2524    nonLinearCost_ = NULL;
     2525    ClpModel::operator=(rhs);
     2526    gutsOfCopy(rhs);
     2527  }
     2528  return *this;
     2529}
     2530void ClpSimplex::gutsOfCopy(const ClpSimplex &rhs)
     2531{
     2532  assert(numberRows_ == rhs.numberRows_);
     2533  assert(numberColumns_ == rhs.numberColumns_);
     2534  numberExtraRows_ = rhs.numberExtraRows_;
     2535  maximumBasic_ = rhs.maximumBasic_;
     2536  dontFactorizePivots_ = rhs.dontFactorizePivots_;
     2537  int numberRows2 = numberRows_ + numberExtraRows_;
     2538  moreSpecialOptions_ = rhs.moreSpecialOptions_;
     2539  if ((whatsChanged_ & 1) != 0) {
     2540    int numberTotal = numberColumns_ + numberRows2;
     2541    if ((specialOptions_ & 65536) != 0 && maximumRows_ >= 0) {
     2542      assert(maximumInternalRows_ >= numberRows2);
     2543      assert(maximumInternalColumns_ >= numberColumns_);
     2544      numberTotal = 2 * (maximumInternalColumns_ + maximumInternalRows_);
     2545    }
     2546    lower_ = ClpCopyOfArray(rhs.lower_, numberTotal);
     2547    rowLowerWork_ = lower_ + numberColumns_;
     2548    columnLowerWork_ = lower_;
     2549    upper_ = ClpCopyOfArray(rhs.upper_, numberTotal);
     2550    rowUpperWork_ = upper_ + numberColumns_;
     2551    columnUpperWork_ = upper_;
     2552    cost_ = ClpCopyOfArray(rhs.cost_, numberTotal);
     2553    objectiveWork_ = cost_;
     2554    rowObjectiveWork_ = cost_ + numberColumns_;
     2555    dj_ = ClpCopyOfArray(rhs.dj_, numberTotal);
     2556    if (dj_) {
     2557      reducedCostWork_ = dj_;
     2558      rowReducedCost_ = dj_ + numberColumns_;
     2559    }
     2560    solution_ = ClpCopyOfArray(rhs.solution_, numberTotal);
     2561    if (solution_) {
     2562      columnActivityWork_ = solution_;
     2563      rowActivityWork_ = solution_ + numberColumns_;
     2564    }
     2565    if (rhs.pivotVariable_) {
     2566      pivotVariable_ = new int[numberRows2];
     2567      CoinMemcpyN(rhs.pivotVariable_, numberRows2, pivotVariable_);
     2568    } else {
     2569      pivotVariable_ = NULL;
     2570    }
     2571    savedSolution_ = ClpCopyOfArray(rhs.savedSolution_, numberTotal);
     2572    int i;
     2573    for (i = 0; i < 6; i++) {
     2574      rowArray_[i] = NULL;
     2575      if (rhs.rowArray_[i])
     2576        rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
     2577      columnArray_[i] = NULL;
     2578      if (rhs.columnArray_[i])
     2579        columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
     2580    }
     2581    if (rhs.saveStatus_) {
     2582      saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberTotal);
     2583    }
     2584  } else {
     2585    lower_ = NULL;
     2586    rowLowerWork_ = NULL;
     2587    columnLowerWork_ = NULL;
     2588    upper_ = NULL;
     2589    rowUpperWork_ = NULL;
     2590    columnUpperWork_ = NULL;
     2591    cost_ = NULL;
     2592    objectiveWork_ = NULL;
     2593    rowObjectiveWork_ = NULL;
     2594    dj_ = NULL;
     2595    reducedCostWork_ = NULL;
     2596    rowReducedCost_ = NULL;
     2597    solution_ = NULL;
     2598    columnActivityWork_ = NULL;
     2599    rowActivityWork_ = NULL;
     2600    pivotVariable_ = NULL;
     2601    savedSolution_ = NULL;
     2602    int i;
     2603    for (i = 0; i < 6; i++) {
     2604      rowArray_[i] = NULL;
     2605      columnArray_[i] = NULL;
     2606    }
     2607    saveStatus_ = NULL;
     2608  }
     2609  if (rhs.factorization_) {
     2610    setFactorization(*rhs.factorization_);
     2611  } else {
     2612    delete factorization_;
     2613    factorization_ = NULL;
     2614  }
     2615  bestPossibleImprovement_ = rhs.bestPossibleImprovement_;
     2616  columnPrimalSequence_ = rhs.columnPrimalSequence_;
     2617  zeroTolerance_ = rhs.zeroTolerance_;
     2618  rowPrimalSequence_ = rhs.rowPrimalSequence_;
     2619  bestObjectiveValue_ = rhs.bestObjectiveValue_;
     2620  baseIteration_ = rhs.baseIteration_;
     2621  vectorMode_ = rhs.vectorMode_;
     2622  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
     2623  largeValue_ = rhs.largeValue_;
     2624  largestPrimalError_ = rhs.largestPrimalError_;
     2625  largestDualError_ = rhs.largestDualError_;
     2626  alphaAccuracy_ = rhs.alphaAccuracy_;
     2627  dualBound_ = rhs.dualBound_;
     2628  alpha_ = rhs.alpha_;
     2629  theta_ = rhs.theta_;
     2630  lowerIn_ = rhs.lowerIn_;
     2631  valueIn_ = rhs.valueIn_;
     2632  upperIn_ = rhs.upperIn_;
     2633  dualIn_ = rhs.dualIn_;
     2634  sequenceIn_ = rhs.sequenceIn_;
     2635  directionIn_ = rhs.directionIn_;
     2636  lowerOut_ = rhs.lowerOut_;
     2637  valueOut_ = rhs.valueOut_;
     2638  upperOut_ = rhs.upperOut_;
     2639  dualOut_ = rhs.dualOut_;
     2640  sequenceOut_ = rhs.sequenceOut_;
     2641  directionOut_ = rhs.directionOut_;
     2642  pivotRow_ = rhs.pivotRow_;
     2643  lastGoodIteration_ = rhs.lastGoodIteration_;
     2644  numberRefinements_ = rhs.numberRefinements_;
     2645  dualTolerance_ = rhs.dualTolerance_;
     2646  primalTolerance_ = rhs.primalTolerance_;
     2647  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
     2648  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
     2649  numberDualInfeasibilitiesWithoutFree_ = rhs.numberDualInfeasibilitiesWithoutFree_;
     2650  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
     2651  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
     2652  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
     2653  dualRowPivot_->setModel(this);
     2654  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
     2655  primalColumnPivot_->setModel(this);
     2656  numberTimesOptimal_ = rhs.numberTimesOptimal_;
     2657  disasterArea_ = NULL;
     2658  changeMade_ = rhs.changeMade_;
     2659  algorithm_ = rhs.algorithm_;
     2660  forceFactorization_ = rhs.forceFactorization_;
     2661  perturbation_ = rhs.perturbation_;
     2662  infeasibilityCost_ = rhs.infeasibilityCost_;
     2663  lastBadIteration_ = rhs.lastBadIteration_;
     2664  lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
     2665  numberFake_ = rhs.numberFake_;
     2666  numberChanged_ = rhs.numberChanged_;
     2667  progressFlag_ = rhs.progressFlag_;
     2668  firstFree_ = rhs.firstFree_;
     2669  incomingInfeasibility_ = rhs.incomingInfeasibility_;
     2670  allowedInfeasibility_ = rhs.allowedInfeasibility_;
     2671  automaticScale_ = rhs.automaticScale_;
    27092672#ifdef ABC_INHERIT
    2710      abcSimplex_ = NULL;
    2711      abcState_ = rhs.abcState_;
    2712 #endif
    2713      maximumPerturbationSize_ = rhs.maximumPerturbationSize_;
    2714      if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) {
    2715           perturbationArray_ = CoinCopyOfArray(rhs.perturbationArray_,
    2716                                                maximumPerturbationSize_);
    2717      } else {
    2718           maximumPerturbationSize_ = 0;
    2719           perturbationArray_ = NULL;
    2720      }
    2721      if (rhs.baseModel_) {
    2722           baseModel_ = new ClpSimplex(*rhs.baseModel_);
    2723      } else {
    2724           baseModel_ = NULL;
    2725      }
    2726      progress_ = rhs.progress_;
    2727      for (int i = 0; i < 4; i++) {
    2728           spareIntArray_[i] = rhs.spareIntArray_[i];
    2729           spareDoubleArray_[i] = rhs.spareDoubleArray_[i];
    2730      }
    2731      sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    2732      sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    2733      acceptablePivot_ = rhs.acceptablePivot_;
    2734      if (rhs.nonLinearCost_ != NULL)
    2735           nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
    2736      else
    2737           nonLinearCost_ = NULL;
    2738      solveType_ = rhs.solveType_;
    2739      eventHandler_->setSimplex(this);
     2673  abcSimplex_ = NULL;
     2674  abcState_ = rhs.abcState_;
     2675#endif
     2676  maximumPerturbationSize_ = rhs.maximumPerturbationSize_;
     2677  if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) {
     2678    perturbationArray_ = CoinCopyOfArray(rhs.perturbationArray_,
     2679      maximumPerturbationSize_);
     2680  } else {
     2681    maximumPerturbationSize_ = 0;
     2682    perturbationArray_ = NULL;
     2683  }
     2684  if (rhs.baseModel_) {
     2685    baseModel_ = new ClpSimplex(*rhs.baseModel_);
     2686  } else {
     2687    baseModel_ = NULL;
     2688  }
     2689  progress_ = rhs.progress_;
     2690  for (int i = 0; i < 4; i++) {
     2691    spareIntArray_[i] = rhs.spareIntArray_[i];
     2692    spareDoubleArray_[i] = rhs.spareDoubleArray_[i];
     2693  }
     2694  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
     2695  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
     2696  acceptablePivot_ = rhs.acceptablePivot_;
     2697  if (rhs.nonLinearCost_ != NULL)
     2698    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
     2699  else
     2700    nonLinearCost_ = NULL;
     2701  solveType_ = rhs.solveType_;
     2702  eventHandler_->setSimplex(this);
    27402703}
    27412704// type == 0 do everything, most + pivot data, 2 factorization data as well
    2742 void
    2743 ClpSimplex::gutsOfDelete(int type)
    2744 {
    2745      if (!type || (specialOptions_ & 65536) == 0) {
    2746           maximumInternalColumns_ = -1;
    2747           maximumInternalRows_ = -1;
    2748           delete [] lower_;
    2749           lower_ = NULL;
    2750           rowLowerWork_ = NULL;
    2751           columnLowerWork_ = NULL;
    2752           delete [] upper_;
    2753           upper_ = NULL;
    2754           rowUpperWork_ = NULL;
    2755           columnUpperWork_ = NULL;
    2756           delete [] cost_;
    2757           cost_ = NULL;
    2758           objectiveWork_ = NULL;
    2759           rowObjectiveWork_ = NULL;
    2760           delete [] dj_;
    2761           dj_ = NULL;
    2762           reducedCostWork_ = NULL;
    2763           rowReducedCost_ = NULL;
    2764           delete [] solution_;
    2765           solution_ = NULL;
    2766           rowActivityWork_ = NULL;
    2767           columnActivityWork_ = NULL;
    2768           delete [] savedSolution_;
    2769           savedSolution_ = NULL;
    2770      }
    2771      if ((specialOptions_ & 2) == 0) {
    2772           delete nonLinearCost_;
    2773           nonLinearCost_ = NULL;
    2774      }
    2775      int i;
    2776      if ((specialOptions_ & 65536) == 0) {
    2777           for (i = 0; i < 6; i++) {
    2778                delete rowArray_[i];
    2779                rowArray_[i] = NULL;
    2780                delete columnArray_[i];
    2781                columnArray_[i] = NULL;
    2782           }
    2783      }
    2784      delete [] saveStatus_;
    2785      saveStatus_ = NULL;
    2786      if (type != 1) {
    2787           delete rowCopy_;
    2788           rowCopy_ = NULL;
    2789      }
    2790      if (!type) {
    2791           // delete everything
    2792           setEmptyFactorization();
    2793           delete [] pivotVariable_;
    2794           pivotVariable_ = NULL;
    2795           delete dualRowPivot_;
    2796           dualRowPivot_ = NULL;
    2797           delete primalColumnPivot_;
    2798           primalColumnPivot_ = NULL;
    2799           delete baseModel_;
    2800           baseModel_ = NULL;
    2801           delete [] perturbationArray_;
    2802           perturbationArray_ = NULL;
    2803           maximumPerturbationSize_ = 0;
    2804      } else {
    2805           // delete any size information in methods
    2806           if (type > 1) {
    2807                //assert (factorization_);
    2808                if (factorization_)
    2809                     factorization_->clearArrays();
    2810                delete [] pivotVariable_;
    2811                pivotVariable_ = NULL;
    2812           }
    2813           dualRowPivot_->clearArrays();
    2814           primalColumnPivot_->clearArrays();
    2815      }
     2705void ClpSimplex::gutsOfDelete(int type)
     2706{
     2707  if (!type || (specialOptions_ & 65536) == 0) {
     2708    maximumInternalColumns_ = -1;
     2709    maximumInternalRows_ = -1;
     2710    delete[] lower_;
     2711    lower_ = NULL;
     2712    rowLowerWork_ = NULL;
     2713    columnLowerWork_ = NULL;
     2714    delete[] upper_;
     2715    upper_ = NULL;
     2716    rowUpperWork_ = NULL;
     2717    columnUpperWork_ = NULL;
     2718    delete[] cost_;
     2719    cost_ = NULL;
     2720    objectiveWork_ = NULL;
     2721    rowObjectiveWork_ = NULL;
     2722    delete[] dj_;
     2723    dj_ = NULL;
     2724    reducedCostWork_ = NULL;
     2725    rowReducedCost_ = NULL;
     2726    delete[] solution_;
     2727    solution_ = NULL;
     2728    rowActivityWork_ = NULL;
     2729    columnActivityWork_ = NULL;
     2730    delete[] savedSolution_;
     2731    savedSolution_ = NULL;
     2732  }
     2733  if ((specialOptions_ & 2) == 0) {
     2734    delete nonLinearCost_;
     2735    nonLinearCost_ = NULL;
     2736  }
     2737  int i;
     2738  if ((specialOptions_ & 65536) == 0) {
     2739    for (i = 0; i < 6; i++) {
     2740      delete rowArray_[i];
     2741      rowArray_[i] = NULL;
     2742      delete columnArray_[i];
     2743      columnArray_[i] = NULL;
     2744    }
     2745  }
     2746  delete[] saveStatus_;
     2747  saveStatus_ = NULL;
     2748  if (type != 1) {
     2749    delete rowCopy_;
     2750    rowCopy_ = NULL;
     2751  }
     2752  if (!type) {
     2753    // delete everything
     2754    setEmptyFactorization();
     2755    delete[] pivotVariable_;
     2756    pivotVariable_ = NULL;
     2757    delete dualRowPivot_;
     2758    dualRowPivot_ = NULL;
     2759    delete primalColumnPivot_;
     2760    primalColumnPivot_ = NULL;
     2761    delete baseModel_;
     2762    baseModel_ = NULL;
     2763    delete[] perturbationArray_;
     2764    perturbationArray_ = NULL;
     2765    maximumPerturbationSize_ = 0;
     2766  } else {
     2767    // delete any size information in methods
     2768    if (type > 1) {
     2769      //assert (factorization_);
     2770      if (factorization_)
     2771        factorization_->clearArrays();
     2772      delete[] pivotVariable_;
     2773      pivotVariable_ = NULL;
     2774    }
     2775    dualRowPivot_->clearArrays();
     2776    primalColumnPivot_->clearArrays();
     2777  }
    28162778}
    28172779// This sets largest infeasibility and most infeasible
    2818 void
    2819 ClpSimplex::checkPrimalSolution(const double * rowActivities,
    2820                                 const double * columnActivities)
    2821 {
    2822      double * solution;
    2823      int iRow, iColumn;
    2824 
    2825      objectiveValue_ = 0.0;
    2826      // now look at primal solution
    2827      solution = rowActivityWork_;
    2828      sumPrimalInfeasibilities_ = 0.0;
    2829      numberPrimalInfeasibilities_ = 0;
    2830      double primalTolerance = primalTolerance_;
    2831      double relaxedTolerance = primalTolerance_;
    2832      // we can't really trust infeasibilities if there is primal error
    2833      double error = CoinMin(1.0e-2, largestPrimalError_);
    2834      // allow tolerance at least slightly bigger than standard
    2835      relaxedTolerance = relaxedTolerance +  error;
    2836      sumOfRelaxedPrimalInfeasibilities_ = 0.0;
    2837      for (iRow = 0; iRow < numberRows_; iRow++) {
    2838           //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
    2839           double infeasibility = 0.0;
    2840           objectiveValue_ += solution[iRow] * rowObjectiveWork_[iRow];
    2841           if (solution[iRow] > rowUpperWork_[iRow]) {
    2842                infeasibility = solution[iRow] - rowUpperWork_[iRow];
    2843           } else if (solution[iRow] < rowLowerWork_[iRow]) {
    2844                infeasibility = rowLowerWork_[iRow] - solution[iRow];
    2845           }
    2846           if (infeasibility > primalTolerance) {
    2847                sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
    2848                if (infeasibility > relaxedTolerance)
    2849                     sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
    2850                numberPrimalInfeasibilities_ ++;
    2851           }
    2852           infeasibility = fabs(rowActivities[iRow] - solution[iRow]);
    2853      }
    2854      // Check any infeasibilities from dynamic rows
    2855      matrix_->primalExpanded(this, 2);
    2856      solution = columnActivityWork_;
    2857      if (!matrix_->rhsOffset(this)) {
    2858           for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2859                //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
    2860                double infeasibility = 0.0;
    2861                objectiveValue_ += objectiveWork_[iColumn] * solution[iColumn];
    2862                if (solution[iColumn] > columnUpperWork_[iColumn]) {
    2863                     infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
    2864                } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
    2865                     infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
    2866                }
    2867                if (infeasibility > primalTolerance) {
    2868                     sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
    2869                     if (infeasibility > relaxedTolerance)
    2870                          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
    2871                     numberPrimalInfeasibilities_ ++;
    2872                }
    2873                infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
    2874           }
    2875      } else {
    2876           // as we are using effective rhs we only check basics
    2877           // But we do need to get objective
    2878           objectiveValue_ += innerProduct(objectiveWork_, numberColumns_, solution);
    2879           for (int j = 0; j < numberRows_; j++) {
    2880                int iColumn = pivotVariable_[j];
    2881                //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
    2882                double infeasibility = 0.0;
    2883                if (solution[iColumn] > columnUpperWork_[iColumn]) {
    2884                     infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
    2885                } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
    2886                     infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
    2887                }
    2888                if (infeasibility > primalTolerance) {
    2889                     sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
    2890                     if (infeasibility > relaxedTolerance)
    2891                          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
    2892                     numberPrimalInfeasibilities_ ++;
    2893                }
    2894                infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
    2895           }
    2896      }
    2897      objectiveValue_ += objective_->nonlinearOffset();
    2898      objectiveValue_ /= (objectiveScale_ * rhsScale_);
    2899 }
    2900 void
    2901 ClpSimplex::checkDualSolution()
    2902 {
    2903 
    2904      int iRow, iColumn;
    2905      sumDualInfeasibilities_ = 0.0;
    2906      numberDualInfeasibilities_ = 0;
    2907      numberDualInfeasibilitiesWithoutFree_ = 0;
    2908      if (matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) {
    2909           // pretend we found dual infeasibilities
    2910           sumOfRelaxedDualInfeasibilities_ = 1.0;
    2911           sumDualInfeasibilities_ = 1.0;
    2912           numberDualInfeasibilities_ = 1;
    2913           return;
    2914      }
    2915      int firstFreePrimal = -1;
    2916      int firstFreeDual = -1;
    2917      int numberSuperBasicWithDj = 0;
    2918      bestPossibleImprovement_ = 0.0;
    2919      // we can't really trust infeasibilities if there is dual error
    2920      double error = CoinMin(1.0e-2, largestDualError_);
    2921      // allow tolerance at least slightly bigger than standard
    2922      double relaxedTolerance = dualTolerance_ +  error;
    2923      // allow bigger tolerance for possible improvement
    2924      double possTolerance = 5.0 * relaxedTolerance;
    2925      sumOfRelaxedDualInfeasibilities_ = 0.0;
    2926 
    2927      // Check any djs from dynamic rows
    2928      matrix_->dualExpanded(this, NULL, NULL, 3);
    2929      numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_;
    2930      objectiveValue_ = 0.0;
    2931      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    2932           objectiveValue_ += objectiveWork_[iColumn] * columnActivityWork_[iColumn];
    2933           if (getColumnStatus(iColumn) != basic && !flagged(iColumn)) {
    2934                // not basic
    2935                double distanceUp = columnUpperWork_[iColumn] -
    2936                                    columnActivityWork_[iColumn];
    2937                double distanceDown = columnActivityWork_[iColumn] -
    2938                                      columnLowerWork_[iColumn];
    2939                if (distanceUp > primalTolerance_) {
    2940                     double value = reducedCostWork_[iColumn];
    2941                     // Check if "free"
    2942                     if (distanceDown > primalTolerance_) {
    2943                          if (fabs(value) > 1.0e2 * relaxedTolerance) {
    2944                               numberSuperBasicWithDj++;
    2945                               if (firstFreeDual < 0)
    2946                                    firstFreeDual = iColumn;
    2947                          }
    2948                          if (firstFreePrimal < 0)
    2949                               firstFreePrimal = iColumn;
    2950                     }
    2951                     // should not be negative
    2952                     if (value < 0.0) {
    2953                          value = - value;
    2954                          if (value > dualTolerance_) {
    2955                               if (getColumnStatus(iColumn) != isFree) {
    2956                                    numberDualInfeasibilitiesWithoutFree_ ++;
    2957                                    sumDualInfeasibilities_ += value - dualTolerance_;
    2958                                    if (value > possTolerance)
    2959                                         bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
    2960                                    if (value > relaxedTolerance)
    2961                                         sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
    2962                                    numberDualInfeasibilities_ ++;
    2963                               } else {
    2964                                    // free so relax a lot
    2965                                    value *= 0.01;
    2966                                    if (value > dualTolerance_) {
    2967                                         sumDualInfeasibilities_ += value - dualTolerance_;
    2968                                         if (value > possTolerance)
    2969                                              bestPossibleImprovement_ = 1.0e100;
    2970                                         if (value > relaxedTolerance)
    2971                                              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
    2972                                         numberDualInfeasibilities_ ++;
    2973                                    }
    2974                               }
    2975                          }
    2976                     }
    2977                }
    2978                if (distanceDown > primalTolerance_) {
    2979                     double value = reducedCostWork_[iColumn];
    2980                     // should not be positive
    2981                     if (value > 0.0) {
    2982                          if (value > dualTolerance_) {
    2983                               sumDualInfeasibilities_ += value - dualTolerance_;
    2984                               if (value > possTolerance)
    2985                                    bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
    2986                               if (value > relaxedTolerance)
    2987                                    sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
    2988                               numberDualInfeasibilities_ ++;
    2989                               if (getColumnStatus(iColumn) != isFree)
    2990                                    numberDualInfeasibilitiesWithoutFree_ ++;
    2991                               // maybe we can make feasible by increasing tolerance
    2992                          }
    2993                     }
    2994                }
    2995           }
    2996      }
    2997      for (iRow = 0; iRow < numberRows_; iRow++) {
    2998           objectiveValue_ += rowActivityWork_[iRow] * rowObjectiveWork_[iRow];
    2999           if (getRowStatus(iRow) != basic && !flagged(iRow + numberColumns_)) {
    3000                // not basic
    3001                double distanceUp = rowUpperWork_[iRow] - rowActivityWork_[iRow];
    3002                double distanceDown = rowActivityWork_[iRow] - rowLowerWork_[iRow];
    3003                if (distanceUp > primalTolerance_) {
    3004                     double value = rowReducedCost_[iRow];
    3005                     // Check if "free"
    3006                     if (distanceDown > primalTolerance_) {
    3007                          if (fabs(value) > 1.0e2 * relaxedTolerance) {
    3008                               numberSuperBasicWithDj++;
    3009                               if (firstFreeDual < 0)
    3010                                    firstFreeDual = iRow + numberColumns_;
    3011                          }
    3012                          if (firstFreePrimal < 0)
    3013                               firstFreePrimal = iRow + numberColumns_;
    3014                     }
    3015                     // should not be negative
    3016                     if (value < 0.0) {
    3017                          value = - value;
    3018                          if (value > dualTolerance_) {
    3019                               sumDualInfeasibilities_ += value - dualTolerance_;
    3020                               if (value > possTolerance)
    3021                                    bestPossibleImprovement_ += value * CoinMin(distanceUp, 1.0e10);
    3022                               if (value > relaxedTolerance)
    3023                                    sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
    3024                               numberDualInfeasibilities_ ++;
    3025                               if (getRowStatus(iRow) != isFree)
    3026                                    numberDualInfeasibilitiesWithoutFree_ ++;
    3027                          }
    3028                     }
    3029                }
    3030                if (distanceDown > primalTolerance_) {
    3031                     double value = rowReducedCost_[iRow];
    3032                     // should not be positive
    3033                     if (value > 0.0) {
    3034                          if (value > dualTolerance_) {
    3035                               sumDualInfeasibilities_ += value - dualTolerance_;
    3036                               if (value > possTolerance)
    3037                                    bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
    3038                               if (value > relaxedTolerance)
    3039                                    sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
    3040                               numberDualInfeasibilities_ ++;
    3041                               if (getRowStatus(iRow) != isFree)
    3042                                    numberDualInfeasibilitiesWithoutFree_ ++;
    3043                               // maybe we can make feasible by increasing tolerance
    3044                          }
    3045                     }
    3046                }
    3047           }
    3048      }
    3049      if (algorithm_ < 0 && firstFreeDual >= 0) {
    3050           // dual
    3051           firstFree_ = firstFreeDual;
    3052      } else if (numberSuperBasicWithDj ||
    3053                 (progress_.lastIterationNumber(0) <= 0)) {
    3054           firstFree_ = firstFreePrimal;
    3055      }
    3056      objectiveValue_ += objective_->nonlinearOffset();
    3057      objectiveValue_ /= (objectiveScale_ * rhsScale_);
     2780void ClpSimplex::checkPrimalSolution(const double *rowActivities,
     2781  const double *columnActivities)
     2782{
     2783  double *solution;
     2784  int iRow, iColumn;
     2785
     2786  objectiveValue_ = 0.0;
     2787  // now look at primal solution
     2788  solution = rowActivityWork_;
     2789  sumPrimalInfeasibilities_ = 0.0;
     2790  numberPrimalInfeasibilities_ = 0;
     2791  double primalTolerance = primalTolerance_;
     2792  double relaxedTolerance = primalTolerance_;
     2793  // we can't really trust infeasibilities if there is primal error
     2794  double error = CoinMin(1.0e-2, largestPrimalError_);
     2795  // allow tolerance at least slightly bigger than standard
     2796  relaxedTolerance = relaxedTolerance + error;
     2797  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
     2798  for (iRow = 0; iRow < numberRows_; iRow++) {
     2799    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
     2800    double infeasibility = 0.0;
     2801    objectiveValue_ += solution[iRow] * rowObjectiveWork_[iRow];
     2802    if (solution[iRow] > rowUpperWork_[iRow]) {
     2803      infeasibility = solution[iRow] - rowUpperWork_[iRow];
     2804    } else if (solution[iRow] < rowLowerWork_[iRow]) {
     2805      infeasibility = rowLowerWork_[iRow] - solution[iRow];
     2806    }
     2807    if (infeasibility > primalTolerance) {
     2808      sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
     2809      if (infeasibility > relaxedTolerance)
     2810        sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
     2811      numberPrimalInfeasibilities_++;
     2812    }
     2813    infeasibility = fabs(rowActivities[iRow] - solution[iRow]);
     2814  }
     2815  // Check any infeasibilities from dynamic rows
     2816  matrix_->primalExpanded(this, 2);
     2817  solution = columnActivityWork_;
     2818  if (!matrix_->rhsOffset(this)) {
     2819    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2820      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
     2821      double infeasibility = 0.0;
     2822      objectiveValue_ += objectiveWork_[iColumn] * solution[iColumn];
     2823      if (solution[iColumn] > columnUpperWork_[iColumn]) {
     2824        infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
     2825      } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
     2826        infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
     2827      }
     2828      if (infeasibility > primalTolerance) {
     2829        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
     2830        if (infeasibility > relaxedTolerance)
     2831          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
     2832        numberPrimalInfeasibilities_++;
     2833      }
     2834      infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
     2835    }
     2836  } else {
     2837    // as we are using effective rhs we only check basics
     2838    // But we do need to get objective
     2839    objectiveValue_ += innerProduct(objectiveWork_, numberColumns_, solution);
     2840    for (int j = 0; j < numberRows_; j++) {
     2841      int iColumn = pivotVariable_[j];
     2842      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
     2843      double infeasibility = 0.0;
     2844      if (solution[iColumn] > columnUpperWork_[iColumn]) {
     2845        infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
     2846      } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
     2847        infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
     2848      }
     2849      if (infeasibility > primalTolerance) {
     2850        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
     2851        if (infeasibility > relaxedTolerance)
     2852          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
     2853        numberPrimalInfeasibilities_++;
     2854      }
     2855      infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
     2856    }
     2857  }
     2858  objectiveValue_ += objective_->nonlinearOffset();
     2859  objectiveValue_ /= (objectiveScale_ * rhsScale_);
     2860}
     2861void ClpSimplex::checkDualSolution()
     2862{
     2863
     2864  int iRow, iColumn;
     2865  sumDualInfeasibilities_ = 0.0;
     2866  numberDualInfeasibilities_ = 0;
     2867  numberDualInfeasibilitiesWithoutFree_ = 0;
     2868  if (matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) {
     2869    // pretend we found dual infeasibilities
     2870    sumOfRelaxedDualInfeasibilities_ = 1.0;
     2871    sumDualInfeasibilities_ = 1.0;
     2872    numberDualInfeasibilities_ = 1;
     2873    return;
     2874  }
     2875  int firstFreePrimal = -1;
     2876  int firstFreeDual = -1;
     2877  int numberSuperBasicWithDj = 0;
     2878  bestPossibleImprovement_ = 0.0;
     2879  // we can't really trust infeasibilities if there is dual error
     2880  double error = CoinMin(1.0e-2, largestDualError_);
     2881  // allow tolerance at least slightly bigger than standard
     2882  double relaxedTolerance = dualTolerance_ + error;
     2883  // allow bigger tolerance for possible improvement
     2884  double possTolerance = 5.0 * relaxedTolerance;
     2885  sumOfRelaxedDualInfeasibilities_ = 0.0;
     2886
     2887  // Check any djs from dynamic rows
     2888  matrix_->dualExpanded(this, NULL, NULL, 3);
     2889  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_;
     2890  objectiveValue_ = 0.0;
     2891  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2892    objectiveValue_ += objectiveWork_[iColumn] * columnActivityWork_[iColumn];
     2893    if (getColumnStatus(iColumn) != basic && !flagged(iColumn)) {
     2894      // not basic
     2895      double distanceUp = columnUpperWork_[iColumn] - columnActivityWork_[iColumn];
     2896      double distanceDown = columnActivityWork_[iColumn] - columnLowerWork_[iColumn];
     2897      if (distanceUp > primalTolerance_) {
     2898        double value = reducedCostWork_[iColumn];
     2899        // Check if "free"
     2900        if (distanceDown > primalTolerance_) {
     2901          if (fabs(value) > 1.0e2 * relaxedTolerance) {
     2902            numberSuperBasicWithDj++;
     2903            if (firstFreeDual < 0)
     2904              firstFreeDual = iColumn;
     2905          }
     2906          if (firstFreePrimal < 0)
     2907            firstFreePrimal = iColumn;
     2908        }
     2909        // should not be negative
     2910        if (value < 0.0) {
     2911          value = -value;
     2912          if (value > dualTolerance_) {
     2913            if (getColumnStatus(iColumn) != isFree) {
     2914              numberDualInfeasibilitiesWithoutFree_++;
     2915              sumDualInfeasibilities_ += value - dualTolerance_;
     2916              if (value > possTolerance)
     2917                bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
     2918              if (value > relaxedTolerance)
     2919                sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
     2920              numberDualInfeasibilities_++;
     2921            } else {
     2922              // free so relax a lot
     2923              value *= 0.01;
     2924              if (value > dualTolerance_) {
     2925                sumDualInfeasibilities_ += value - dualTolerance_;
     2926                if (value > possTolerance)
     2927                  bestPossibleImprovement_ = 1.0e100;
     2928                if (value > relaxedTolerance)
     2929                  sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
     2930                numberDualInfeasibilities_++;
     2931              }
     2932            }
     2933          }
     2934        }
     2935      }
     2936      if (distanceDown > primalTolerance_) {
     2937        double value = reducedCostWork_[iColumn];
     2938        // should not be positive
     2939        if (value > 0.0) {
     2940          if (value > dualTolerance_) {
     2941            sumDualInfeasibilities_ += value - dualTolerance_;
     2942            if (value > possTolerance)
     2943              bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
     2944            if (value > relaxedTolerance)
     2945              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
     2946            numberDualInfeasibilities_++;
     2947            if (getColumnStatus(iColumn) != isFree)
     2948              numberDualInfeasibilitiesWithoutFree_++;
     2949            // maybe we can make feasible by increasing tolerance
     2950          }
     2951        }
     2952      }
     2953    }
     2954  }
     2955  for (iRow = 0; iRow < numberRows_; iRow++) {
     2956    objectiveValue_ += rowActivityWork_[iRow] * rowObjectiveWork_[iRow];
     2957    if (getRowStatus(iRow) != basic && !flagged(iRow + numberColumns_)) {
     2958      // not basic
     2959      double distanceUp = rowUpperWork_[iRow] - rowActivityWork_[iRow];
     2960      double distanceDown = rowActivityWork_[iRow] - rowLowerWork_[iRow];
     2961      if (distanceUp > primalTolerance_) {
     2962        double value = rowReducedCost_[iRow];
     2963        // Check if "free"
     2964        if (distanceDown > primalTolerance_) {
     2965          if (fabs(value) > 1.0e2 * relaxedTolerance) {
     2966            numberSuperBasicWithDj++;
     2967            if (firstFreeDual < 0)
     2968              firstFreeDual = iRow + numberColumns_;
     2969          }
     2970          if (firstFreePrimal < 0)
     2971            firstFreePrimal = iRow + numberColumns_;
     2972        }
     2973        // should not be negative
     2974        if (value < 0.0) {
     2975          value = -value;
     2976          if (value > dualTolerance_) {
     2977            sumDualInfeasibilities_ += value - dualTolerance_;
     2978            if (value > possTolerance)
     2979              bestPossibleImprovement_ += value * CoinMin(distanceUp, 1.0e10);
     2980            if (value > relaxedTolerance)
     2981              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
     2982            numberDualInfeasibilities_++;
     2983            if (getRowStatus(iRow) != isFree)
     2984              numberDualInfeasibilitiesWithoutFree_++;
     2985          }
     2986        }
     2987      }
     2988      if (distanceDown > primalTolerance_) {
     2989        double value = rowReducedCost_[iRow];
     2990        // should not be positive
     2991        if (value > 0.0) {
     2992          if (value > dualTolerance_) {
     2993            sumDualInfeasibilities_ += value - dualTolerance_;
     2994            if (value > possTolerance)
     2995              bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
     2996            if (value > relaxedTolerance)
     2997              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
     2998            numberDualInfeasibilities_++;
     2999            if (getRowStatus(iRow) != isFree)
     3000              numberDualInfeasibilitiesWithoutFree_++;
     3001            // maybe we can make feasible by increasing tolerance
     3002          }
     3003        }
     3004      }
     3005    }
     3006  }
     3007  if (algorithm_ < 0 && firstFreeDual >= 0) {
     3008    // dual
     3009    firstFree_ = firstFreeDual;
     3010  } else if (numberSuperBasicWithDj || (progress_.lastIterationNumber(0) <= 0)) {
     3011    firstFree_ = firstFreePrimal;
     3012  }
     3013  objectiveValue_ += objective_->nonlinearOffset();
     3014  objectiveValue_ /= (objectiveScale_ * rhsScale_);
    30583015}
    30593016/* This sets sum and number of infeasibilities (Dual and Primal) */
    3060 void
    3061 ClpSimplex::checkBothSolutions()
    3062 {
    3063      if ((matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) ||
    3064                matrix_->rhsOffset(this)) {
    3065           // Say may be free or superbasic
     3017void ClpSimplex::checkBothSolutions()
     3018{
     3019  if ((matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) || matrix_->rhsOffset(this)) {
     3020    // Say may be free or superbasic
     3021    moreSpecialOptions_ &= ~8;
     3022    // old way
     3023    checkPrimalSolution(rowActivityWork_, columnActivityWork_);
     3024    checkDualSolution();
     3025    return;
     3026  }
     3027  int iSequence;
     3028  assert(dualTolerance_ > 0.0 && dualTolerance_ < 1.0e10);
     3029  assert(primalTolerance_ > 0.0 && primalTolerance_ < 1.0e10);
     3030  objectiveValue_ = 0.0;
     3031  sumPrimalInfeasibilities_ = 0.0;
     3032  numberPrimalInfeasibilities_ = 0;
     3033  double primalTolerance = primalTolerance_;
     3034  double relaxedToleranceP = primalTolerance_;
     3035  // we can't really trust infeasibilities if there is primal error
     3036  double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_, 0.0 * primalTolerance_));
     3037  // allow tolerance at least slightly bigger than standard
     3038  relaxedToleranceP = relaxedToleranceP + error;
     3039  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
     3040  sumDualInfeasibilities_ = 0.0;
     3041  numberDualInfeasibilities_ = 0;
     3042  double dualTolerance = dualTolerance_;
     3043  double relaxedToleranceD = dualTolerance;
     3044  // we can't really trust infeasibilities if there is dual error
     3045  error = CoinMin(1.0e-2, CoinMax(largestDualError_, 5.0 * dualTolerance_));
     3046  // allow tolerance at least slightly bigger than standard
     3047  relaxedToleranceD = relaxedToleranceD + error;
     3048  // allow bigger tolerance for possible improvement
     3049  double possTolerance = 5.0 * relaxedToleranceD;
     3050  sumOfRelaxedDualInfeasibilities_ = 0.0;
     3051  bestPossibleImprovement_ = 0.0;
     3052
     3053  // Check any infeasibilities from dynamic rows
     3054  matrix_->primalExpanded(this, 2);
     3055  // Check any djs from dynamic rows
     3056  matrix_->dualExpanded(this, NULL, NULL, 3);
     3057  int numberDualInfeasibilitiesFree = 0;
     3058  int firstFreePrimal = -1;
     3059  int firstFreeDual = -1;
     3060  int numberSuperBasicWithDj = 0;
     3061
     3062  int numberTotal = numberRows_ + numberColumns_;
     3063  // Say no free or superbasic
     3064  moreSpecialOptions_ |= 8;
     3065  //#define PRINT_INFEAS
     3066#ifdef PRINT_INFEAS
     3067  int seqInf[10];
     3068#endif
     3069  for (iSequence = 0; iSequence < numberTotal; iSequence++) {
     3070    //#define LIKELY_SUPERBASIC
     3071#ifdef LIKELY_SUPERBASIC
     3072    if (getStatus(iSequence) == isFree || getStatus(iSequence) == superBasic)
     3073      moreSpecialOptions_ &= ~8; // Say superbasic variables exist
     3074#endif
     3075    double value = solution_[iSequence];
     3076#ifdef COIN_DEBUG
     3077    if (fabs(value) > 1.0e20)
     3078      printf("%d values %g %g %g - status %d\n", iSequence, lower_[iSequence],
     3079        solution_[iSequence], upper_[iSequence], status_[iSequence]);
     3080#endif
     3081    objectiveValue_ += value * cost_[iSequence];
     3082    double distanceUp = upper_[iSequence] - value;
     3083    double distanceDown = value - lower_[iSequence];
     3084    if (distanceUp < -primalTolerance) {
     3085      double infeasibility = -distanceUp;
     3086#ifndef LIKELY_SUPERBASIC
     3087      if (getStatus(iSequence) != basic)
     3088        moreSpecialOptions_ &= ~8; // Say superbasic variables exist
     3089#endif
     3090      sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
     3091      if (infeasibility > relaxedToleranceP)
     3092        sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
     3093#ifdef PRINT_INFEAS
     3094      if (numberPrimalInfeasibilities_ < 10) {
     3095        seqInf[numberPrimalInfeasibilities_] = iSequence;
     3096      }
     3097#endif
     3098      numberPrimalInfeasibilities_++;
     3099    } else if (distanceDown < -primalTolerance) {
     3100      double infeasibility = -distanceDown;
     3101#ifndef LIKELY_SUPERBASIC
     3102      if (getStatus(iSequence) != basic)
     3103        moreSpecialOptions_ &= ~8; // Say superbasic variables exist
     3104#endif
     3105      sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
     3106      if (infeasibility > relaxedToleranceP)
     3107        sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
     3108#ifdef PRINT_INFEAS
     3109      if (numberPrimalInfeasibilities_ < 10) {
     3110        seqInf[numberPrimalInfeasibilities_] = iSequence;
     3111      }
     3112#endif
     3113      numberPrimalInfeasibilities_++;
     3114    } else {
     3115      // feasible (so could be free)
     3116      if (getStatus(iSequence) != basic && !flagged(iSequence)) {
     3117        // not basic
     3118        double djValue = dj_[iSequence];
     3119        if (distanceDown < primalTolerance) {
     3120          if (distanceUp > primalTolerance && djValue < -dualTolerance) {
     3121            sumDualInfeasibilities_ -= djValue + dualTolerance;
     3122            if (djValue < -possTolerance)
     3123              bestPossibleImprovement_ -= distanceUp * djValue;
     3124            if (djValue < -relaxedToleranceD)
     3125              sumOfRelaxedDualInfeasibilities_ -= djValue + relaxedToleranceD;
     3126            numberDualInfeasibilities_++;
     3127          }
     3128        } else if (distanceUp < primalTolerance) {
     3129          if (djValue > dualTolerance) {
     3130            sumDualInfeasibilities_ += djValue - dualTolerance;
     3131            if (djValue > possTolerance)
     3132              bestPossibleImprovement_ += distanceDown * djValue;
     3133            if (djValue > relaxedToleranceD)
     3134              sumOfRelaxedDualInfeasibilities_ += djValue - relaxedToleranceD;
     3135            numberDualInfeasibilities_++;
     3136          }
     3137        } else {
     3138          // may be free
     3139          // Say free or superbasic