source: trunk/Clp/src/ClpSimplex.cpp @ 2393

Last change on this file since 2393 was 2393, checked in by forrest, 4 months ago

change problem status of 4 to 1 for bound conflict

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 436.0 KB
Line 
1/* $Id: ClpSimplex.cpp 2393 2019-02-21 12:01:58Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6//#undef NDEBUG
7
8#include "ClpConfig.h"
9
10#include "CoinPragma.hpp"
11#include <math.h>
12
13#if SLIM_CLP == 2
14#define SLIM_NOIO
15#endif
16#include "CoinHelperFunctions.hpp"
17#include "CoinFloatEqual.hpp"
18#include "ClpSimplex.hpp"
19#include "ClpFactorization.hpp"
20#include "ClpPackedMatrix.hpp"
21#include "CoinIndexedVector.hpp"
22#include "ClpDualRowDantzig.hpp"
23#include "ClpDualRowSteepest.hpp"
24#include "ClpPrimalColumnDantzig.hpp"
25#include "ClpPrimalColumnSteepest.hpp"
26#include "ClpPEDualRowSteepest.hpp"
27#include "ClpPEPrimalColumnSteepest.hpp"
28#include "ClpNonLinearCost.hpp"
29#include "ClpMessage.hpp"
30#include "ClpEventHandler.hpp"
31#include "ClpLinearObjective.hpp"
32#include "ClpHelperFunctions.hpp"
33#include "CoinModel.hpp"
34#include "CoinLpIO.hpp"
35#include <cfloat>
36#if CLP_HAS_ABC
37#include "CoinAbcCommon.hpp"
38#endif
39
40#include <string>
41#include <stdio.h>
42#include <iostream>
43//#############################################################################
44
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)
131#ifdef ABC_INHERIT
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
155}
156
157// Subproblem constructor
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)
246#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  }
309}
310// Subproblem constructor
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)
399#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 = 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  }
486}
487// Puts solution back small model
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_);
516#if 0
517     if (!problemStatus_) {
518          ClpDisjointCopyN(smallModel.objective(), smallModel.numberColumns_, smallModel.reducedCost_);
519          smallModel.matrix_->transposeTimes(-1.0, smallModel.dual_, smallModel.reducedCost_);
520          for (int i = 0; i < smallModel.numberColumns_; i++) {
521               if (smallModel.getColumnStatus(i) == basic)
522                    assert (fabs(smallModel.reducedCost_[i]) < 1.0e-5);
523          }
524          ClpDisjointCopyN(objective(), numberColumns_, reducedCost_);
525          matrix_->transposeTimes(-1.0, dual_, reducedCost_);
526          for (int i = 0; i < numberColumns_; i++) {
527               if (getColumnStatus(i) == basic)
528                    assert (fabs(reducedCost_[i]) < 1.0e-5);
529          }
530     }
531#endif
532  matrix()->times(columnActivity_, rowActivity_);
533}
534
535//-----------------------------------------------------------------------------
536
537ClpSimplex::~ClpSimplex()
538{
539  setPersistenceFlag(0);
540  gutsOfDelete(0);
541  delete nonLinearCost_;
542}
543//#############################################################################
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) {
583#ifdef CLP_USER_DRIVEN
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();
604#ifdef CLP_DEBUG
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);
655#if 0
656                    allSlackBasis(true);
657                    CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
658#else
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_);
739#ifndef COIN_REUSE_RANDOM
740    memset(dj_, 0, (numberRows_ + numberColumns_) * sizeof(double));
741#else
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();
790#ifdef CLP_USEFUL_PRINTOUT
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);
810#ifdef CLP_USEFUL_PRINTOUT
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;
916#ifdef CLP_INVESTIGATE
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);
950#ifdef CLP_DEBUG
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;
971#if 0
972     static double * xsave=NULL;
973     {
974       double * xx = thisVector->denseVector();
975       double largest=0.0;
976       int iLargest=-1;
977       for (int i=0;i<numberRows_;i++) {
978         if (fabs(xx[i])>largest) {
979           largest=fabs(xx[i]);
980           iLargest=i;
981         }
982       }
983       printf("largest incoming rhs %g on row %d\n",largest,iLargest);
984     }
985     if (numberIterations_<-40722) {
986       double * xx = thisVector->denseVector();
987       if (xsave) {
988         double * sol = xsave+numberRows_;
989         double largest=0.0;
990         int iLargest=-1;
991         for (int i=0;i<numberRows_;i++) {
992           if (fabs(xx[i]-xsave[i])>largest) {
993             largest=fabs(xx[i]-xsave[i]);
994             iLargest=i;
995           }
996         }
997         printf("error %g on row %d\n",largest,iLargest);
998         largest=0.0;
999         iLargest=-1;
1000         for (int i=0;i<numberColumns_;i++) {
1001           if (fabs(solution_[i]-sol[i])>largest) {
1002             largest=fabs(solution_[i]-sol[i]);
1003             iLargest=i;
1004           }
1005         }
1006         printf("error %g on col %d\n",largest,iLargest);
1007       } else {
1008         xsave=new double[2*numberRows_+numberColumns_];
1009       }
1010       memcpy(xsave,xx,numberRows_*sizeof(double));
1011       memcpy(xsave+numberRows_,solution_,(numberRows_+numberColumns_)*sizeof(double));
1012     }
1013#endif
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();
1019#ifdef CLP_DEBUG
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;
1343          } else {
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();
1410#ifndef SLIM_CLP
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  }
1417#endif
1418}
1419/* Given an existing factorization computes and checks
1420   primal and dual solutions.  Uses input arrays for variables at
1421   bounds.  Returns feasibility states */
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();
1436}
1437/* Given an existing factorization computes and checks
1438   primal and dual solutions.  Uses current problem arrays for
1439   bounds.  Returns feasibility states */
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;
1450}
1451// Factorizes using current basis.  This is for external use
1452// Return codes are as from ClpFactorization
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;
1463}
1464// Clean up status
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  }
1524}
1525
1526/* Factorizes using current basis.
1527   solveType - 1 iterating, 0 initial, -1 external
1528   - 2 then iterating but can throw out of basis
1529   If 10 added then in primal values pass
1530   Return codes are as from ClpFactorization unless initial factorization
1531   when total number of singularities is returned.
1532   Special case is numberRows_+1 -> all slack basis.
1533*/
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  }
1546#ifdef CLP_DEBUG
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
1609#ifdef CLP_INVESTIGATE3
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        }
1764#ifdef CLP_INVESTIGATE3
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        }
1780#endif
1781#if 0
1782                    if (numberBasic < numberRows_) {
1783                         // add some slacks in case odd warmstart
1784#ifdef CLP_INVESTIGATE
1785                         printf("BAD %d basic, %d rows %d slacks\n",
1786                                numberBasic, numberRows_, totalSlacks);
1787#endif
1788                         int iRow = numberRows_ - 1;
1789                         while (numberBasic < numberRows_) {
1790                              if (getRowStatus(iRow) != basic) {
1791                                   setRowStatus(iRow, basic);
1792                                   numberBasic++;
1793                                   totalSlacks++;
1794                                   iRow--;
1795                              } else {
1796                                   break;
1797                              }
1798                         }
1799                    }
1800#endif
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            }
1816          } else {
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;
1865#if 0
1866                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1867                         if (getColumnStatus(iColumn) == basic)
1868                              numberBasic++;
1869                    }
1870#endif
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 --
1926     if (0)  {
1927          static int k = 0;
1928          printf("start basis\n");
1929          int i;
1930          for (i = 0; i < numberRows_; i++)
1931               printf ("xx %d %d\n", i, pivotVariable_[i]);
1932          for (i = 0; i < numberRows_ + numberColumns_; i++)
1933               if (getColumnStatus(i) == basic)
1934                    printf ("yy %d basic\n", i);
1935          if (k > 20)
1936               exit(0);
1937          k++;
1938     }
1939#endif
1940#if 0 //ndef NDEBUG
1941  // Make sure everything is clean
1942  double sumOutside=0.0;
1943  int numberOutside=0;
1944  //double sumOutsideLarge=0.0;
1945  int numberOutsideLarge=0;
1946  double sumInside=0.0;
1947  int numberInside=0;
1948  //double sumInsideLarge=0.0;
1949  int numberInsideLarge=0;
1950  int numberTotal=numberRows_+numberColumns_;
1951  for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
1952    if(getStatus(iSequence) == isFixed) {
1953      // double check fixed
1954      assert (upper_[iSequence] == lower_[iSequence]);
1955      assert (fabs(solution_[iSequence]-lower_[iSequence])<primalTolerance_);
1956    } else if (getStatus(iSequence) == isFree) {
1957      assert (upper_[iSequence] == COIN_DBL_MAX && lower_[iSequence]==-COIN_DBL_MAX);
1958    } else if (getStatus(iSequence) == atLowerBound) {
1959      assert (fabs(solution_[iSequence]-lower_[iSequence])<1000.0*primalTolerance_);
1960      if (solution_[iSequence]<lower_[iSequence]) {
1961        numberOutside++;
1962        sumOutside-=solution_[iSequence]-lower_[iSequence];
1963        if (solution_[iSequence]<lower_[iSequence]-primalTolerance_)
1964          numberOutsideLarge++;
1965      } else if (solution_[iSequence]>lower_[iSequence]) {
1966        numberInside++;
1967        sumInside+=solution_[iSequence]-lower_[iSequence];
1968        if (solution_[iSequence]>lower_[iSequence]+primalTolerance_)
1969          numberInsideLarge++;
1970      }
1971    } else if (getStatus(iSequence) == atUpperBound) {
1972      assert (fabs(solution_[iSequence]-upper_[iSequence])<1000.0*primalTolerance_);
1973      if (solution_[iSequence]>upper_[iSequence]) {
1974        numberOutside++;
1975        sumOutside+=solution_[iSequence]-upper_[iSequence];
1976        if (solution_[iSequence]>upper_[iSequence]+primalTolerance_)
1977          numberOutsideLarge++;
1978      } else if (solution_[iSequence]<upper_[iSequence]) {
1979        numberInside++;
1980        sumInside-=solution_[iSequence]-upper_[iSequence];
1981        if (solution_[iSequence]<upper_[iSequence]-primalTolerance_)
1982          numberInsideLarge++;
1983      }
1984    } else if (getStatus(iSequence) == superBasic) {
1985      //assert (!valuesPass);
1986    }
1987  }
1988  if (numberInside+numberOutside)
1989    printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
1990           numberOutside,sumOutside,numberOutsideLarge,
1991           numberInside,sumInside,numberInsideLarge);
1992#endif
1993  int status = factorization_->factorize(this, solveType, valuesPass);
1994  if (status) {
1995    handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
1996      << status
1997      << CoinMessageEol;
1998#ifdef CLP_USEFUL_PRINTOUT
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;
2026}
2027/*
2028   This does basis housekeeping and does values for in/out variables.
2029   Can also decide to re-factorize
2030*/
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  }
2052#if 0
2053     printf("h1 %d %d %g %g %g %g",
2054            directionOut_
2055            , directionIn_, theta_
2056            , dualOut_, dualIn_, alpha_);
2057#endif
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  }
2122#if 0
2123     if (numberIterations_ > 10000)
2124          printf(" it %d %g %c%d %c%d\n"
2125                 , numberIterations_, objectiveValue()
2126                 , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
2127                 , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
2128#endif
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;
2201#if 1
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
2225               should be better if don't reject incoming
2226               as it is in basis */
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  }
2299}
2300// Copy constructor.
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)
2385#ifdef ABC_INHERIT
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
2408}
2409// Copy constructor from model
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)
2494#ifdef ABC_INHERIT
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
2516}
2517// Assignment operator. This copies the data
2518ClpSimplex &
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_;
2672#ifdef ABC_INHERIT
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);
2703}
2704// type == 0 do everything, most + pivot data, 2 factorization data as well
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  }
2778}
2779// This sets largest infeasibility and most infeasible
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_);
3015}
3016/* This sets sum and number of infeasibilities (Dual and Primal) */
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
3140          moreSpecialOptions_ &= ~8;
3141          djValue *= 0.01;
3142          if (fabs(djValue) > dualTolerance) {
3143            if (getStatus(iSequence) == isFree)
3144              numberDualInfeasibilitiesFree++;
3145            sumDualInfeasibilities_ += fabs(djValue) - dualTolerance;
3146            bestPossibleImprovement_ = 1.0e100;
3147            numberDualInfeasibilities_++;
3148            if (fabs(djValue) > relaxedToleranceD) {
3149              sumOfRelaxedDualInfeasibilities_ += value - relaxedToleranceD;
3150              numberSuperBasicWithDj++;
3151              if (firstFreeDual < 0)
3152                firstFreeDual = iSequence;
3153            }
3154          }
3155          if (firstFreePrimal < 0)
3156            firstFreePrimal = iSequence;
3157        }
3158      }
3159    }
3160  }
3161  objectiveValue_ += objective_->nonlinearOffset();
3162  objectiveValue_ /= (objectiveScale_ * rhsScale_);
3163  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberDualInfeasibilitiesFree;
3164#ifdef PRINT_INFEAS
3165  if (numberPrimalInfeasibilities_ <= 10) {
3166    printf("---------------start-----------\n");
3167    if (!rowScale_) {
3168      for (int i = 0; i < numberPrimalInfeasibilities_; i++) {
3169        int iSeq = seqInf[i];
3170        double infeas;
3171        if (solution_[iSeq] < lower_[iSeq])
3172          infeas = lower_[iSeq] - solution_[iSeq];
3173        else
3174          infeas = solution_[iSeq] - upper_[iSeq];
3175        if (iSeq < numberColumns_) {
3176          printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g\n",
3177            iSeq, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas);
3178        } else {
3179          printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g\n",
3180            iSeq - numberColumns_, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas);
3181        }
3182      }
3183    } else {
3184      for (int i = 0; i < numberPrimalInfeasibilities_; i++) {
3185        int iSeq = seqInf[i];
3186        double infeas;
3187        if (solution_[iSeq] < lower_[iSeq])
3188          infeas = lower_[iSeq] - solution_[iSeq];
3189        else
3190          infeas = solution_[iSeq] - upper_[iSeq];
3191        double unscaled = infeas;
3192        if (iSeq < numberColumns_) {
3193          unscaled *= columnScale_[iSeq];
3194          printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3195            iSeq, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas, unscaled);
3196        } else {
3197          unscaled /= rowScale_[iSeq - numberColumns_];
3198          printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3199            iSeq - numberColumns_, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas, unscaled);
3200        }
3201      }
3202    }
3203  }
3204#endif
3205  if (algorithm_ < 0 && firstFreeDual >= 0) {
3206    // dual
3207    firstFree_ = firstFreeDual;
3208  } else if (numberSuperBasicWithDj || (progress_.lastIterationNumber(0) <= 0)) {
3209    firstFree_ = firstFreePrimal;
3210  }
3211}
3212/* Adds multiple of a column into an array */
3213void ClpSimplex::add(double *array,
3214  int sequence, double multiplier) const
3215{
3216  if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3217    //slack
3218    array[sequence - numberColumns_] -= multiplier;
3219  } else {
3220    // column
3221    matrix_->add(this, array, sequence, multiplier);
3222  }
3223}
3224/*
3225  Unpacks one column of the matrix into indexed array
3226*/
3227void ClpSimplex::unpack(CoinIndexedVector *rowArray) const
3228{
3229  rowArray->clear();
3230  if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3231    //slack
3232    rowArray->insert(sequenceIn_ - numberColumns_, -1.0);
3233  } else {
3234    // column
3235    matrix_->unpack(this, rowArray, sequenceIn_);
3236  }
3237}
3238void ClpSimplex::unpack(CoinIndexedVector *rowArray, int sequence) const
3239{
3240  rowArray->clear();
3241  if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3242    //slack
3243    rowArray->insert(sequence - numberColumns_, -1.0);
3244  } else {
3245    // column
3246    matrix_->unpack(this, rowArray, sequence);
3247  }
3248}
3249/*
3250  Unpacks one column of the matrix into indexed array
3251*/
3252void ClpSimplex::unpackPacked(CoinIndexedVector *rowArray)
3253{
3254  rowArray->clear();
3255  if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3256    //slack
3257    int *index = rowArray->getIndices();
3258    double *array = rowArray->denseVector();
3259    array[0] = -1.0;
3260    index[0] = sequenceIn_ - numberColumns_;
3261    rowArray->setNumElements(1);
3262    rowArray->setPackedMode(true);
3263  } else {
3264    // column
3265    matrix_->unpackPacked(this, rowArray, sequenceIn_);
3266  }
3267}
3268void ClpSimplex::unpackPacked(CoinIndexedVector *rowArray, int sequence)
3269{
3270  rowArray->clear();
3271  if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3272    //slack
3273    int *index = rowArray->getIndices();
3274    double *array = rowArray->denseVector();
3275    array[0] = -1.0;
3276    index[0] = sequence - numberColumns_;
3277    rowArray->setNumElements(1);
3278    rowArray->setPackedMode(true);
3279  } else {
3280    // column
3281    matrix_->unpackPacked(this, rowArray, sequence);
3282  }
3283}
3284//static int x_gaps[4]={0,0,0,0};
3285//static int scale_times[]={0,0,0,0};
3286bool ClpSimplex::createRim(int what, bool makeRowCopy, int startFinishOptions)
3287{
3288  bool goodMatrix = true;
3289  int saveLevel = handler_->logLevel();
3290  spareIntArray_[0] = 0;
3291  if (!matrix_->canGetRowCopy())
3292    makeRowCopy = false; // switch off row copy if can't produce
3293  // Arrays will be there and correct size unless what is 63
3294  bool newArrays = (what == 63);
3295  // We may be restarting with same size
3296  bool keepPivots = false;
3297  if (startFinishOptions == -1) {
3298    startFinishOptions = 0;
3299    keepPivots = true;
3300  }
3301  bool oldMatrix = ((startFinishOptions & 4) != 0 && (whatsChanged_ & 1) != 0);
3302  if (what == 63) {
3303    pivotRow_ = -1;
3304    if (!status_)
3305      createStatus();
3306    if (oldMatrix)
3307      newArrays = false;
3308    if (problemStatus_ == 10) {
3309      handler_->setLogLevel(0); // switch off messages
3310      if (rowArray_[0]) {
3311        // stuff is still there
3312        oldMatrix = true;
3313        newArrays = false;
3314        keepPivots = true;
3315        for (int iRow = 0; iRow < 4; iRow++) {
3316          rowArray_[iRow]->clear();
3317        }
3318        for (int iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
3319          columnArray_[iColumn]->clear();
3320        }
3321      }
3322    } else if (factorization_) {
3323      // match up factorization messages
3324      if (handler_->logLevel() < 3)
3325        factorization_->messageLevel(0);
3326      else
3327        factorization_->messageLevel(CoinMax(3, factorization_->messageLevel()));
3328      /* Faster to keep pivots rather than re-scan matrix.  Matrix may have changed
3329                  i.e. oldMatrix false but okay as long as same number rows and status array exists
3330               */
3331      if ((startFinishOptions & 2) != 0 && factorization_->numberRows() == numberRows_ && status_)
3332        keepPivots = true;
3333    }
3334    numberExtraRows_ = matrix_->generalExpanded(this, 2, maximumBasic_);
3335    if (numberExtraRows_ && newArrays) {
3336      // make sure status array large enough
3337      assert(status_);
3338      int numberOld = numberRows_ + numberColumns_;
3339      int numberNew = numberRows_ + numberColumns_ + numberExtraRows_;
3340      unsigned char *newStatus = new unsigned char[numberNew];
3341      memset(newStatus + numberOld, 0, numberExtraRows_);
3342      CoinMemcpyN(status_, numberOld, newStatus);
3343      delete[] status_;
3344      status_ = newStatus;
3345    }
3346  }
3347  int numberRows2 = numberRows_ + numberExtraRows_;
3348  int numberTotal = numberRows2 + numberColumns_;
3349  if ((specialOptions_ & 65536) != 0) {
3350    assert(!numberExtraRows_);
3351    if (!cost_ || numberRows2 > maximumInternalRows_ || numberColumns_ > maximumInternalColumns_) {
3352      newArrays = true;
3353      keepPivots = false;
3354      COIN_DETAIL_PRINT(printf("createrim a %d rows, %d maximum rows %d maxinternal\n",
3355        numberRows_, maximumRows_, maximumInternalRows_));
3356      int oldMaximumRows = maximumInternalRows_;
3357      int oldMaximumColumns = maximumInternalColumns_;
3358      if (cost_) {
3359        if (numberRows2 > maximumInternalRows_)
3360          maximumInternalRows_ = numberRows2;
3361        if (numberColumns_ > maximumInternalColumns_)
3362          maximumInternalColumns_ = numberColumns_;
3363      } else {
3364        maximumInternalRows_ = numberRows2;
3365        maximumInternalColumns_ = numberColumns_;
3366      }
3367      //maximumRows_=CoinMax(maximumInternalRows_,maximumRows_);
3368      //maximumColumns_=CoinMax(maximumInternalColumns_,maximumColumns_);
3369      assert(maximumInternalRows_ == maximumRows_);
3370      assert(maximumInternalColumns_ == maximumColumns_);
3371      COIN_DETAIL_PRINT(printf("createrim b %d rows, %d maximum rows, %d maxinternal\n",
3372        numberRows_, maximumRows_, maximumInternalRows_));
3373      int numberTotal2 = (maximumInternalRows_ + maximumInternalColumns_) * 2;
3374      delete[] cost_;
3375      cost_ = new double[numberTotal2];
3376      delete[] lower_;
3377      delete[] upper_;
3378      lower_ = new double[numberTotal2];
3379      upper_ = new double[numberTotal2];
3380      delete[] dj_;
3381      dj_ = new double[numberTotal2];
3382      delete[] solution_;
3383      solution_ = new double[numberTotal2];
3384      // ***** should be non NULL but seems to be too much
3385      //printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_);
3386      if (savedRowScale_) {
3387        assert(oldMaximumRows > 0);
3388        double *temp;
3389        temp = new double[4 * maximumRows_];
3390        CoinFillN(temp, 4 * maximumRows_, 1.0);
3391        CoinMemcpyN(savedRowScale_, numberRows_, temp);
3392        CoinMemcpyN(savedRowScale_ + oldMaximumRows, numberRows_, temp + maximumRows_);
3393        CoinMemcpyN(savedRowScale_ + 2 * oldMaximumRows, numberRows_, temp + 2 * maximumRows_);
3394        CoinMemcpyN(savedRowScale_ + 3 * oldMaximumRows, numberRows_, temp + 3 * maximumRows_);
3395        delete[] savedRowScale_;
3396        savedRowScale_ = temp;
3397        temp = new double[4 * maximumColumns_];
3398        CoinFillN(temp, 4 * maximumColumns_, 1.0);
3399        CoinMemcpyN(savedColumnScale_, numberColumns_, temp);
3400        CoinMemcpyN(savedColumnScale_ + oldMaximumColumns, numberColumns_, temp + maximumColumns_);
3401        CoinMemcpyN(savedColumnScale_ + 2 * oldMaximumColumns, numberColumns_, temp + 2 * maximumColumns_);
3402        CoinMemcpyN(savedColumnScale_ + 3 * oldMaximumColumns, numberColumns_, temp + 3 * maximumColumns_);
3403        delete[] savedColumnScale_;
3404        savedColumnScale_ = temp;
3405      }
3406    }
3407  }
3408  int i;
3409  bool doSanityCheck = true;
3410  if (what == 63) {
3411    // We may want to switch stuff off for speed
3412    if ((specialOptions_ & 256) != 0)
3413      makeRowCopy = false; // no row copy
3414    if ((specialOptions_ & 128) != 0)
3415      doSanityCheck = false; // no sanity check
3416    //check matrix
3417    if (!matrix_)
3418      matrix_ = new ClpPackedMatrix();
3419    int checkType = (doSanityCheck) ? 15 : 14;
3420    if (oldMatrix)
3421      checkType = 14;
3422    bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
3423    if (inCbcOrOther)
3424      checkType -= 4; // don't check for duplicates
3425    if (!matrix_->allElementsInRange(this, smallElement_, 1.0e20, checkType)) {
3426      problemStatus_ = 4;
3427      secondaryStatus_ = 8;
3428      //goodMatrix= false;
3429      return false;
3430    }
3431    bool rowCopyIsScaled;
3432    if (makeRowCopy) {
3433      if (!oldMatrix || !rowCopy_) {
3434        delete rowCopy_;
3435        // may return NULL if can't give row copy
3436        rowCopy_ = matrix_->reverseOrderedCopy();
3437        rowCopyIsScaled = false;
3438      } else {
3439        rowCopyIsScaled = true;
3440      }
3441    }
3442#if 0
3443          if (what == 63) {
3444               int k = rowScale_ ? 1 : 0;
3445               if (oldMatrix)
3446                    k += 2;
3447               scale_times[k]++;
3448               if ((scale_times[0] + scale_times[1] + scale_times[2] + scale_times[3]) % 1000 == 0)
3449                    printf("scale counts %d %d %d %d\n",
3450                           scale_times[0], scale_times[1], scale_times[2], scale_times[3]);
3451          }
3452#endif
3453    // do scaling if needed
3454    if (!oldMatrix && scalingFlag_ < 0) {
3455      if (scalingFlag_ < 0 && rowScale_) {
3456        //if (handler_->logLevel()>0)
3457        printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
3458          scalingFlag_);
3459        scalingFlag_ = 0;
3460      }
3461      delete[] rowScale_;
3462      delete[] columnScale_;
3463      rowScale_ = NULL;
3464      columnScale_ = NULL;
3465    }
3466    inverseRowScale_ = NULL;
3467    inverseColumnScale_ = NULL;
3468    if (scalingFlag_ > 0 && (specialOptions_ & 65536) != 0 && rowScale_ && rowScale_ == savedRowScale_)
3469      rowScale_ = NULL;
3470    if (scalingFlag_ > 0 && !rowScale_) {
3471      if ((specialOptions_ & 65536) != 0) {
3472        assert(!rowScale_);
3473        rowScale_ = savedRowScale_;
3474        columnScale_ = savedColumnScale_;
3475        // put back original
3476        if (savedRowScale_) {
3477          inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
3478          inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
3479          CoinMemcpyN(savedRowScale_ + 2 * maximumInternalRows_,
3480            numberRows2, savedRowScale_);
3481          CoinMemcpyN(savedRowScale_ + 3 * maximumInternalRows_,
3482            numberRows2, inverseRowScale_);
3483          CoinMemcpyN(savedColumnScale_ + 2 * maximumColumns_,
3484            numberColumns_, savedColumnScale_);
3485          CoinMemcpyN(savedColumnScale_ + 3 * maximumColumns_,
3486            numberColumns_, inverseColumnScale_);
3487        }
3488      }
3489      if (matrix_->scale(this, this))
3490        scalingFlag_ = -scalingFlag_; // not scaled after all
3491      if (rowScale_ && automaticScale_) {
3492        if (!savedRowScale_) {
3493          inverseRowScale_ = rowScale_ + numberRows2;
3494          inverseColumnScale_ = columnScale_ + numberColumns_;
3495        }
3496        // try automatic scaling
3497        double smallestObj = 1.0e100;
3498        double largestObj = 0.0;
3499        double largestRhs = 0.0;
3500        const double *obj = objective();
3501        for (i = 0; i < numberColumns_; i++) {
3502          double value = fabs(obj[i]);
3503          value *= columnScale_[i];
3504          if (value && columnLower_[i] != columnUpper_[i]) {
3505            smallestObj = CoinMin(smallestObj, value);
3506            largestObj = CoinMax(largestObj, value);
3507          }
3508          if (columnLower_[i] > 0.0 || columnUpper_[i] < 0.0) {
3509            double scale = 1.0 * inverseColumnScale_[i];
3510            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3511            if (columnLower_[i] > 0)
3512              largestRhs = CoinMax(largestRhs, columnLower_[i] * scale);
3513            if (columnUpper_[i] < 0.0)
3514              largestRhs = CoinMax(largestRhs, -columnUpper_[i] * scale);
3515          }
3516        }
3517        for (i = 0; i < numberRows_; i++) {
3518          if (rowLower_[i] > 0.0 || rowUpper_[i] < 0.0) {
3519            double scale = rowScale_[i];
3520            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3521            if (rowLower_[i] > 0)
3522              largestRhs = CoinMax(largestRhs, rowLower_[i] * scale);
3523            if (rowUpper_[i] < 0.0)
3524              largestRhs = CoinMax(largestRhs, -rowUpper_[i] * scale);
3525          }
3526        }
3527        COIN_DETAIL_PRINT(printf("small obj %g, large %g - rhs %g\n", smallestObj, largestObj, largestRhs));
3528        bool scalingDone = false;
3529        // look at element range
3530        double smallestNegative;
3531        double largestNegative;
3532        double smallestPositive;
3533        double largestPositive;
3534        matrix_->rangeOfElements(smallestNegative, largestNegative,
3535          smallestPositive, largestPositive);
3536        smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
3537        largestPositive = CoinMax(fabs(largestNegative), largestPositive);
3538        if (largestObj) {
3539          double ratio = largestObj / smallestObj;
3540          double scale = 1.0;
3541          if (ratio < 1.0e8) {
3542            // reasonable
3543            if (smallestObj < 1.0e-4) {
3544              // may as well scale up
3545              scalingDone = true;
3546              scale = 1.0e-3 / smallestObj;
3547            } else if (largestObj < 1.0e6 || (algorithm_ > 0 && largestObj < 1.0e-4 * infeasibilityCost_)) {
3548              //done=true;
3549            } else {
3550              scalingDone = true;
3551              if (algorithm_ < 0) {
3552                scale = 1.0e6 / largestObj;
3553              } else {
3554                scale = CoinMax(1.0e6, 1.0e-4 * infeasibilityCost_) / largestObj;
3555              }
3556            }
3557          } else if (ratio < 1.0e12) {
3558            // not so good
3559            if (smallestObj < 1.0e-7) {
3560              // may as well scale up
3561              scalingDone = true;
3562              scale = 1.0e-6 / smallestObj;
3563            } else if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3564              //done=true;
3565            } else {
3566              scalingDone = true;
3567              if (algorithm_ < 0) {
3568                scale = 1.0e7 / largestObj;
3569              } else {
3570                scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3571              }
3572            }
3573          } else {
3574            // Really nasty problem
3575            if (smallestObj < 1.0e-8) {
3576              // may as well scale up
3577              scalingDone = true;
3578              scale = 1.0e-7 / smallestObj;
3579              largestObj *= scale;
3580            }
3581            if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3582              //done=true;
3583            } else {
3584              scalingDone = true;
3585              if (algorithm_ < 0) {
3586                scale = 1.0e7 / largestObj;
3587              } else {
3588                scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3589              }
3590            }
3591          }
3592          objectiveScale_ = scale;
3593        }
3594        if (largestRhs > 1.0e12) {
3595          scalingDone = true;
3596          rhsScale_ = 1.0e9 / largestRhs;
3597        } else if (largestPositive > 1.0e-14 * smallestPositive && largestRhs > 1.0e6) {
3598          scalingDone = true;
3599          rhsScale_ = 1.0e6 / largestRhs;
3600        } else {
3601          rhsScale_ = 1.0;
3602        }
3603        if (scalingDone) {
3604          handler_->message(CLP_RIM_SCALE, messages_)
3605            << objectiveScale_ << rhsScale_
3606            << CoinMessageEol;
3607        }
3608      }
3609    } else if (makeRowCopy && scalingFlag_ > 0 && !rowCopyIsScaled) {
3610      matrix_->scaleRowCopy(this);
3611    }
3612    if (rowScale_ && !savedRowScale_) {
3613      inverseRowScale_ = rowScale_ + numberRows2;
3614      inverseColumnScale_ = columnScale_ + numberColumns_;
3615    }
3616    // See if we can try for faster row copy
3617    if (makeRowCopy && !oldMatrix) {
3618      ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix_);
3619      if (clpMatrix && numberThreads_)
3620        clpMatrix->specialRowCopy(this, rowCopy_);
3621      if (clpMatrix)
3622        clpMatrix->specialColumnCopy(this);
3623    }
3624  }
3625  if (what == 63) {
3626#if 0
3627          {
3628               x_gaps[0]++;
3629               ClpPackedMatrix* clpMatrix =
3630               dynamic_cast< ClpPackedMatrix*>(matrix_);
3631               if (clpMatrix) {
3632                    if (!clpMatrix->getPackedMatrix()->hasGaps())
3633                         x_gaps[1]++;
3634                    if ((clpMatrix->flags() & 2) == 0)
3635                         x_gaps[3]++;
3636               } else {
3637                    x_gaps[2]++;
3638               }
3639               if ((x_gaps[0] % 1000) == 0)
3640                    printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n",
3641                           x_gaps[0], x_gaps[1], x_gaps[2], x_gaps[3]);
3642          }
3643#endif
3644    if (newArrays && (specialOptions_ & 65536) == 0) {
3645      delete[] cost_;
3646      cost_ = new double[2 * numberTotal];
3647      delete[] lower_;
3648      delete[] upper_;
3649      lower_ = new double[numberTotal];
3650      upper_ = new double[numberTotal];
3651      delete[] dj_;
3652      dj_ = new double[numberTotal];
3653      delete[] solution_;
3654      solution_ = new double[numberTotal];
3655    }
3656    reducedCostWork_ = dj_;
3657    rowReducedCost_ = dj_ + numberColumns_;
3658    columnActivityWork_ = solution_;
3659    rowActivityWork_ = solution_ + numberColumns_;
3660    objectiveWork_ = cost_;
3661    rowObjectiveWork_ = cost_ + numberColumns_;
3662    rowLowerWork_ = lower_ + numberColumns_;
3663    columnLowerWork_ = lower_;
3664    rowUpperWork_ = upper_ + numberColumns_;
3665    columnUpperWork_ = upper_;
3666  }
3667  if ((what & 4) != 0) {
3668    double direction = optimizationDirection_ * objectiveScale_;
3669    const double *obj = objective();
3670    const double *rowScale = rowScale_;
3671    const double *columnScale = columnScale_;
3672    // and also scale by scale factors
3673    if (rowScale) {
3674      if (rowObjective_) {
3675        for (i = 0; i < numberRows_; i++)
3676          rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
3677      } else {
3678        memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3679      }
3680      // If scaled then do all columns later in one loop
3681      if (what != 63) {
3682        for (i = 0; i < numberColumns_; i++) {
3683          CoinAssert(fabs(obj[i]) < 1.0e25);
3684          objectiveWork_[i] = obj[i] * direction * columnScale[i];
3685        }
3686      }
3687    } else {
3688      if (rowObjective_) {
3689        for (i = 0; i < numberRows_; i++)
3690          rowObjectiveWork_[i] = rowObjective_[i] * direction;
3691      } else {
3692        memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3693      }
3694      for (i = 0; i < numberColumns_; i++) {
3695        CoinAssert(fabs(obj[i]) < 1.0e25);
3696        objectiveWork_[i] = obj[i] * direction;
3697      }
3698    }
3699  }
3700  if ((what & 1) != 0) {
3701    const double *rowScale = rowScale_;
3702    // clean up any mismatches on infinity
3703    // and fix any variables with tiny gaps
3704    double primalTolerance = dblParam_[ClpPrimalTolerance];
3705    if (rowScale) {
3706      // If scaled then do all columns later in one loop
3707      if (what != 63) {
3708        const double *inverseScale = inverseColumnScale_;
3709        for (i = 0; i < numberColumns_; i++) {
3710          double multiplier = rhsScale_ * inverseScale[i];
3711          double lowerValue = columnLower_[i];
3712          double upperValue = columnUpper_[i];
3713          if (lowerValue > -1.0e20) {
3714            columnLowerWork_[i] = lowerValue * multiplier;
3715            if (upperValue >= 1.0e20) {
3716              columnUpperWork_[i] = COIN_DBL_MAX;
3717            } else {
3718              columnUpperWork_[i] = upperValue * multiplier;
3719              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3720                if (columnLowerWork_[i] >= 0.0) {
3721                  columnUpperWork_[i] = columnLowerWork_[i];
3722                } else if (columnUpperWork_[i] <= 0.0) {
3723                  columnLowerWork_[i] = columnUpperWork_[i];
3724                } else {
3725                  columnUpperWork_[i] = 0.0;
3726                  columnLowerWork_[i] = 0.0;
3727                }
3728              }
3729            }
3730          } else if (upperValue < 1.0e20) {
3731            columnLowerWork_[i] = -COIN_DBL_MAX;
3732            columnUpperWork_[i] = upperValue * multiplier;
3733          } else {
3734            // free
3735            columnLowerWork_[i] = -COIN_DBL_MAX;
3736            columnUpperWork_[i] = COIN_DBL_MAX;
3737          }
3738        }
3739      }
3740      for (i = 0; i < numberRows_; i++) {
3741        double multiplier = rhsScale_ * rowScale[i];
3742        double lowerValue = rowLower_[i];
3743        double upperValue = rowUpper_[i];
3744        if (lowerValue > -1.0e20) {
3745          rowLowerWork_[i] = lowerValue * multiplier;
3746          if (upperValue >= 1.0e20) {
3747            rowUpperWork_[i] = COIN_DBL_MAX;
3748          } else {
3749            rowUpperWork_[i] = upperValue * multiplier;
3750            if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3751              if (rowLowerWork_[i] >= 0.0) {
3752                rowUpperWork_[i] = rowLowerWork_[i];
3753              } else if (rowUpperWork_[i] <= 0.0) {
3754                rowLowerWork_[i] = rowUpperWork_[i];
3755              } else {
3756                rowUpperWork_[i] = 0.0;
3757                rowLowerWork_[i] = 0.0;
3758              }
3759            }
3760          }
3761        } else if (upperValue < 1.0e20) {
3762          rowLowerWork_[i] = -COIN_DBL_MAX;
3763          rowUpperWork_[i] = upperValue * multiplier;
3764        } else {
3765          // free
3766          rowLowerWork_[i] = -COIN_DBL_MAX;
3767          rowUpperWork_[i] = COIN_DBL_MAX;
3768        }
3769      }
3770    } else if (rhsScale_ != 1.0) {
3771      for (i = 0; i < numberColumns_; i++) {
3772        double lowerValue = columnLower_[i];
3773        double upperValue = columnUpper_[i];
3774        if (lowerValue > -1.0e20) {
3775          columnLowerWork_[i] = lowerValue * rhsScale_;
3776          if (upperValue >= 1.0e20) {
3777            columnUpperWork_[i] = COIN_DBL_MAX;
3778          } else {
3779            columnUpperWork_[i] = upperValue * rhsScale_;
3780            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3781              if (columnLowerWork_[i] >= 0.0) {
3782                columnUpperWork_[i] = columnLowerWork_[i];
3783              } else if (columnUpperWork_[i] <= 0.0) {
3784                columnLowerWork_[i] = columnUpperWork_[i];
3785              } else {
3786                columnUpperWork_[i] = 0.0;
3787                columnLowerWork_[i] = 0.0;
3788              }
3789            }
3790          }
3791        } else if (upperValue < 1.0e20) {
3792          columnLowerWork_[i] = -COIN_DBL_MAX;
3793          columnUpperWork_[i] = upperValue * rhsScale_;
3794        } else {
3795          // free
3796          columnLowerWork_[i] = -COIN_DBL_MAX;
3797          columnUpperWork_[i] = COIN_DBL_MAX;
3798        }
3799      }
3800      for (i = 0; i < numberRows_; i++) {
3801        double lowerValue = rowLower_[i];
3802        double upperValue = rowUpper_[i];
3803        if (lowerValue > -1.0e20) {
3804          rowLowerWork_[i] = lowerValue * rhsScale_;
3805          if (upperValue >= 1.0e20) {
3806            rowUpperWork_[i] = COIN_DBL_MAX;
3807          } else {
3808            rowUpperWork_[i] = upperValue * rhsScale_;
3809            if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3810              if (rowLowerWork_[i] >= 0.0) {
3811                rowUpperWork_[i] = rowLowerWork_[i];
3812              } else if (rowUpperWork_[i] <= 0.0) {
3813                rowLowerWork_[i] = rowUpperWork_[i];
3814              } else {
3815                rowUpperWork_[i] = 0.0;
3816                rowLowerWork_[i] = 0.0;
3817              }
3818            }
3819          }
3820        } else if (upperValue < 1.0e20) {
3821          rowLowerWork_[i] = -COIN_DBL_MAX;
3822          rowUpperWork_[i] = upperValue * rhsScale_;
3823        } else {
3824          // free
3825          rowLowerWork_[i] = -COIN_DBL_MAX;
3826          rowUpperWork_[i] = COIN_DBL_MAX;
3827        }
3828      }
3829    } else {
3830      for (i = 0; i < numberColumns_; i++) {
3831        double lowerValue = columnLower_[i];
3832        double upperValue = columnUpper_[i];
3833        if (lowerValue > -1.0e20) {
3834          columnLowerWork_[i] = lowerValue;
3835          if (upperValue >= 1.0e20) {
3836            columnUpperWork_[i] = COIN_DBL_MAX;
3837          } else {
3838            columnUpperWork_[i] = upperValue;
3839            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3840              if (columnLowerWork_[i] >= 0.0) {
3841                columnUpperWork_[i] = columnLowerWork_[i];
3842              } else if (columnUpperWork_[i] <= 0.0) {
3843                columnLowerWork_[i] = columnUpperWork_[i];
3844              } else {
3845                columnUpperWork_[i] = 0.0;
3846                columnLowerWork_[i] = 0.0;
3847              }
3848            }
3849          }
3850        } else if (upperValue < 1.0e20) {
3851          columnLowerWork_[i] = -COIN_DBL_MAX;
3852          columnUpperWork_[i] = upperValue;
3853        } else {
3854          // free
3855          columnLowerWork_[i] = -COIN_DBL_MAX;
3856          columnUpperWork_[i] = COIN_DBL_MAX;
3857        }
3858      }
3859      for (i = 0; i < numberRows_; i++) {
3860        double lowerValue = rowLower_[i];
3861        double upperValue = rowUpper_[i];
3862        if (lowerValue > -1.0e20) {
3863          rowLowerWork_[i] = lowerValue;
3864          if (upperValue >= 1.0e20) {
3865            rowUpperWork_[i] = COIN_DBL_MAX;
3866          } else {
3867            rowUpperWork_[i] = upperValue;
3868            if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3869              if (rowLowerWork_[i] >= 0.0) {
3870                rowUpperWork_[i] = rowLowerWork_[i];
3871              } else if (rowUpperWork_[i] <= 0.0) {
3872                rowLowerWork_[i] = rowUpperWork_[i];
3873              } else {
3874                rowUpperWork_[i] = 0.0;
3875                rowLowerWork_[i] = 0.0;
3876              }
3877            }
3878          }
3879        } else if (upperValue < 1.0e20) {
3880          rowLowerWork_[i] = -COIN_DBL_MAX;
3881          rowUpperWork_[i] = upperValue;
3882        } else {
3883          // free
3884          rowLowerWork_[i] = -COIN_DBL_MAX;
3885          rowUpperWork_[i] = COIN_DBL_MAX;
3886        }
3887      }
3888    }
3889  }
3890  if (what == 63) {
3891    // move information to work arrays
3892    double direction = optimizationDirection_;
3893    // direction is actually scale out not scale in
3894    if (direction)
3895      direction = 1.0 / direction;
3896    if (direction != 1.0) {
3897      // reverse all dual signs
3898      for (i = 0; i < numberColumns_; i++)
3899        reducedCost_[i] *= direction;
3900      for (i = 0; i < numberRows_; i++)
3901        dual_[i] *= direction;
3902    }
3903    for (i = 0; i < numberRows_ + numberColumns_; i++) {
3904      setFakeBound(i, noFake);
3905    }
3906    if (rowScale_) {
3907      const double *obj = objective();
3908      double direction = optimizationDirection_ * objectiveScale_;
3909      // clean up any mismatches on infinity
3910      // and fix any variables with tiny gaps
3911      double primalTolerance = dblParam_[ClpPrimalTolerance];
3912      // on entry
3913      const double *inverseScale = inverseColumnScale_;
3914      for (i = 0; i < numberColumns_; i++) {
3915        CoinAssert(fabs(obj[i]) < 1.0e25);
3916        double scaleFactor = columnScale_[i];
3917        double multiplier = rhsScale_ * inverseScale[i];
3918        scaleFactor *= direction;
3919        objectiveWork_[i] = obj[i] * scaleFactor;
3920        reducedCostWork_[i] = reducedCost_[i] * scaleFactor;
3921        double lowerValue = columnLower_[i];
3922        double upperValue = columnUpper_[i];
3923        if (lowerValue > -1.0e20) {
3924          columnLowerWork_[i] = lowerValue * multiplier;
3925          if (upperValue >= 1.0e20) {
3926            columnUpperWork_[i] = COIN_DBL_MAX;
3927          } else {
3928            columnUpperWork_[i] = upperValue * multiplier;
3929            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3930              if (columnLowerWork_[i] >= 0.0) {
3931                columnUpperWork_[i] = columnLowerWork_[i];
3932              } else if (columnUpperWork_[i] <= 0.0) {
3933                columnLowerWork_[i] = columnUpperWork_[i];
3934              } else {
3935                columnUpperWork_[i] = 0.0;
3936                columnLowerWork_[i] = 0.0;
3937              }
3938            }
3939          }
3940        } else if (upperValue < 1.0e20) {
3941          columnLowerWork_[i] = -COIN_DBL_MAX;
3942          columnUpperWork_[i] = upperValue * multiplier;
3943        } else {
3944          // free
3945          columnLowerWork_[i] = -COIN_DBL_MAX;
3946          columnUpperWork_[i] = COIN_DBL_MAX;
3947        }
3948        double value = columnActivity_[i] * multiplier;
3949        if (fabs(value) > 1.0e20) {
3950          //printf("bad value of %g for column %d\n",value,i);
3951          setColumnStatus(i, superBasic);
3952          if (columnUpperWork_[i] < 0.0) {
3953            value = columnUpperWork_[i];
3954          } else if (columnLowerWork_[i] > 0.0) {
3955            value = columnLowerWork_[i];
3956          } else {
3957            value = 0.0;
3958          }
3959        }
3960        columnActivityWork_[i] = value;
3961      }
3962      inverseScale = inverseRowScale_;
3963      for (i = 0; i < numberRows_; i++) {
3964        dual_[i] *= inverseScale[i];
3965        dual_[i] *= objectiveScale_;
3966        rowReducedCost_[i] = dual_[i];
3967        double multiplier = rhsScale_ * rowScale_[i];
3968        double value = rowActivity_[i] * multiplier;
3969        if (fabs(value) > 1.0e20) {
3970          //printf("bad value of %g for row %d\n",value,i);
3971          setRowStatus(i, superBasic);
3972          if (rowUpperWork_[i] < 0.0) {
3973            value = rowUpperWork_[i];
3974          } else if (rowLowerWork_[i] > 0.0) {
3975            value = rowLowerWork_[i];
3976          } else {
3977            value = 0.0;
3978          }
3979        }
3980        rowActivityWork_[i] = value;
3981      }
3982    } else if (objectiveScale_ != 1.0 || rhsScale_ != 1.0) {
3983      // on entry
3984      for (i = 0; i < numberColumns_; i++) {
3985        double value = columnActivity_[i];
3986        value *= rhsScale_;
3987        if (fabs(value) > 1.0e20) {
3988          //printf("bad value of %g for column %d\n",value,i);
3989          setColumnStatus(i, superBasic);
3990          if (columnUpperWork_[i] < 0.0) {
3991            value = columnUpperWork_[i];
3992          } else if (columnLowerWork_[i] > 0.0) {
3993            value = columnLowerWork_[i];
3994          } else {
3995            value = 0.0;
3996          }
3997        }
3998        columnActivityWork_[i] = value;
3999        reducedCostWork_[i] = reducedCost_[i] * objectiveScale_;
4000      }
4001      for (i = 0; i < numberRows_; i++) {
4002        double value = rowActivity_[i];
4003        value *= rhsScale_;
4004        if (fabs(value) > 1.0e20) {
4005          //printf("bad value of %g for row %d\n",value,i);
4006          setRowStatus(i, superBasic);
4007          if (rowUpperWork_[i] < 0.0) {
4008            value = rowUpperWork_[i];
4009          } else if (rowLowerWork_[i] > 0.0) {
4010            value = rowLowerWork_[i];
4011          } else {
4012            value = 0.0;
4013          }
4014        }
4015        rowActivityWork_[i] = value;
4016        dual_[i] *= objectiveScale_;
4017        rowReducedCost_[i] = dual_[i];
4018      }
4019    } else {
4020      // on entry
4021      for (i = 0; i < numberColumns_; i++) {
4022        double value = columnActivity_[i];
4023        if (fabs(value) > 1.0e20) {
4024          //printf("bad value of %g for column %d\n",value,i);
4025          setColumnStatus(i, superBasic);
4026          if (columnUpperWork_[i] < 0.0) {
4027            value = columnUpperWork_[i];
4028          } else if (columnLowerWork_[i] > 0.0) {
4029            value = columnLowerWork_[i];
4030          } else {
4031            value = 0.0;
4032          }
4033        }
4034        columnActivityWork_[i] = value;
4035        reducedCostWork_[i] = reducedCost_[i];
4036      }
4037      for (i = 0; i < numberRows_; i++) {
4038        double value = rowActivity_[i];
4039        if (fabs(value) > 1.0e20) {
4040          //printf("bad value of %g for row %d\n",value,i);
4041          setRowStatus(i, superBasic);
4042          if (rowUpperWork_[i] < 0.0) {
4043            value = rowUpperWork_[i];
4044          } else if (rowLowerWork_[i] > 0.0) {
4045            value = rowLowerWork_[i];
4046          } else {
4047            value = 0.0;
4048          }
4049        }
4050        rowActivityWork_[i] = value;
4051        rowReducedCost_[i] = dual_[i];
4052      }
4053    }
4054  }
4055
4056  if (what == 63 && doSanityCheck) {
4057    // check rim of problem okay
4058    if (!sanityCheck())
4059      goodMatrix = false;
4060  }
4061  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
4062  // maybe we need to move scales to SimplexModel for factorization?
4063  if ((what == 63 && !pivotVariable_) || (newArrays && !keepPivots)) {
4064    delete[] pivotVariable_;
4065    pivotVariable_ = new int[numberRows2 + 1];
4066    for (int i = 0; i < numberRows2 + 1; i++)
4067      pivotVariable_[i] = -1;
4068  } else if (what == 63 && !keepPivots) {
4069    // just reset
4070    for (int i = 0; i < numberRows2; i++)
4071      pivotVariable_[i] = -1;
4072  } else if (what == 63) {
4073    // check pivots
4074    for (int i = 0; i < numberRows2; i++) {
4075      int iSequence = pivotVariable_[i];
4076      if (iSequence < 0 || (iSequence < numberRows_ + numberColumns_ && getStatus(iSequence) != basic)) {
4077        keepPivots = false;
4078        break;
4079      }
4080    }
4081    if (!keepPivots) {
4082      // reset
4083      for (int i = 0; i < numberRows2; i++)
4084        pivotVariable_[i] = -1;
4085    } else {
4086      // clean
4087      for (int i = 0; i < numberColumns_ + numberRows_; i++) {
4088        Status status = getStatus(i);
4089        if (status != basic) {
4090          if (upper_[i] == lower_[i]) {
4091            setStatus(i, isFixed);
4092            solution_[i] = lower_[i];
4093          } else if (status == atLowerBound) {
4094            if (lower_[i] > -1.0e20) {
4095              solution_[i] = lower_[i];
4096            } else {
4097              //printf("seq %d at lower of %g\n",i,lower_[i]);
4098              if (upper_[i] < 1.0e20) {
4099                solution_[i] = upper_[i];
4100                setStatus(i, atUpperBound);
4101              } else {
4102                setStatus(i, isFree);
4103              }
4104            }
4105          } else if (status == atUpperBound) {
4106            if (upper_[i] < 1.0e20) {
4107              solution_[i] = upper_[i];
4108            } else {
4109              //printf("seq %d at upper of %g\n",i,upper_[i]);
4110              if (lower_[i] > -1.0e20) {
4111                solution_[i] = lower_[i];
4112                setStatus(i, atLowerBound);
4113              } else {
4114                setStatus(i, isFree);
4115              }
4116            }
4117          } else if (status == isFixed && upper_[i] > lower_[i]) {
4118            // was fixed - not now
4119            if (solution_[i] <= lower_[i] + primalTolerance_) {
4120              setStatus(i, atLowerBound);
4121            } else if (solution_[i] >= upper_[i] - primalTolerance_) {
4122              setStatus(i, atUpperBound);
4123            } else {
4124              setStatus(i, superBasic);
4125            }
4126          }
4127        }
4128      }
4129    }
4130  }
4131
4132  if (what == 63) {
4133    // Safer to get new arrays anyway so test on newArrays removed
4134    // get some arrays
4135    int iRow, iColumn;
4136    // these are "indexed" arrays so we always know where nonzeros are
4137    /**********************************************************
4138               rowArray_[3] is long enough for rows+columns (2 also maybe)
4139               rowArray_[1] is long enough for max(rows,columns)
4140          *********************************************************/
4141    for (iRow = 0; iRow < 4; iRow++) {
4142      int length = numberRows2 + factorization_->maximumPivots();
4143      if (iRow > SHORT_REGION || objective_->type() > 1)
4144        length += numberColumns_;
4145      else if (iRow == 1)
4146        length = CoinMax(length, numberColumns_);
4147      if ((specialOptions_ & 65536) == 0 || !rowArray_[iRow]) {
4148        delete rowArray_[iRow];
4149        rowArray_[iRow] = new CoinIndexedVector();
4150      }
4151      rowArray_[iRow]->reserve(length);
4152    }
4153
4154    for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
4155      if ((specialOptions_ & 65536) == 0 || !columnArray_[iColumn]) {
4156        delete columnArray_[iColumn];
4157        columnArray_[iColumn] = new CoinIndexedVector();
4158      }
4159      columnArray_[iColumn]->reserve(numberColumns_ + numberRows2);
4160    }
4161  }
4162  if (problemStatus_ == 10) {
4163    problemStatus_ = -1;
4164    handler_->setLogLevel(saveLevel); // switch back messages
4165  }
4166  if ((what & 5) != 0)
4167    matrix_->generalExpanded(this, 9, what); // update costs and bounds if necessary
4168  if (goodMatrix && (specialOptions_ & 65536) != 0) {
4169    int save = maximumColumns_ + maximumRows_;
4170    CoinMemcpyN(cost_, numberTotal, cost_ + save);
4171    CoinMemcpyN(lower_, numberTotal, lower_ + save);
4172    CoinMemcpyN(upper_, numberTotal, upper_ + save);
4173    CoinMemcpyN(dj_, numberTotal, dj_ + save);
4174    CoinMemcpyN(solution_, numberTotal, solution_ + save);
4175    if (rowScale_ && !savedRowScale_) {
4176      double *temp;
4177      temp = new double[4 * maximumRows_];
4178      CoinFillN(temp, 4 * maximumRows_, 1.0);
4179      CoinMemcpyN(rowScale_, numberRows2, temp);
4180      CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + maximumRows_);
4181      CoinMemcpyN(rowScale_, numberRows2, temp + 2 * maximumRows_);
4182      CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + 3 * maximumRows_);
4183      delete[] rowScale_;
4184      savedRowScale_ = temp;
4185      rowScale_ = savedRowScale_;
4186      inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
4187      temp = new double[4 * maximumColumns_];
4188      CoinFillN(temp, 4 * maximumColumns_, 1.0);
4189      CoinMemcpyN(columnScale_, numberColumns_, temp);
4190      CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + maximumColumns_);
4191      CoinMemcpyN(columnScale_, numberColumns_, temp + 2 * maximumColumns_);
4192      CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + 3 * maximumColumns_);
4193      delete[] columnScale_;
4194      savedColumnScale_ = temp;
4195      columnScale_ = savedColumnScale_;
4196      inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
4197    }
4198  }
4199#ifdef CLP_USER_DRIVEN
4200  eventHandler_->event(ClpEventHandler::endOfCreateRim);
4201#endif
4202  return goodMatrix;
4203}
4204// Does rows and columns
4205void ClpSimplex::createRim1(bool initial)
4206{
4207  int i;
4208  int numberRows2 = numberRows_ + numberExtraRows_;
4209  int numberTotal = numberRows2 + numberColumns_;
4210  if ((specialOptions_ & 65536) != 0 && true) {
4211    assert(!initial);
4212    int save = maximumColumns_ + maximumRows_;
4213    CoinMemcpyN(lower_ + save, numberTotal, lower_);
4214    CoinMemcpyN(upper_ + save, numberTotal, upper_);
4215    return;
4216  }
4217  const double *rowScale = rowScale_;
4218  // clean up any mismatches on infinity
4219  // and fix any variables with tiny gaps
4220  double primalTolerance = dblParam_[ClpPrimalTolerance];
4221  if (rowScale) {
4222    // If scaled then do all columns later in one loop
4223    if (!initial) {
4224      const double *inverseScale = inverseColumnScale_;
4225      for (i = 0; i < numberColumns_; i++) {
4226        double multiplier = rhsScale_ * inverseScale[i];
4227        double lowerValue = columnLower_[i];
4228        double upperValue = columnUpper_[i];
4229        if (lowerValue > -1.0e20) {
4230          columnLowerWork_[i] = lowerValue * multiplier;
4231          if (upperValue >= 1.0e20) {
4232            columnUpperWork_[i] = COIN_DBL_MAX;
4233          } else {
4234            columnUpperWork_[i] = upperValue * multiplier;
4235            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4236              if (columnLowerWork_[i] >= 0.0) {
4237                columnUpperWork_[i] = columnLowerWork_[i];
4238              } else if (columnUpperWork_[i] <= 0.0) {
4239                columnLowerWork_[i] = columnUpperWork_[i];
4240              } else {
4241                columnUpperWork_[i] = 0.0;
4242                columnLowerWork_[i] = 0.0;
4243              }
4244            }
4245          }
4246        } else if (upperValue < 1.0e20) {
4247          columnLowerWork_[i] = -COIN_DBL_MAX;
4248          columnUpperWork_[i] = upperValue * multiplier;
4249        } else {
4250          // free
4251          columnLowerWork_[i] = -COIN_DBL_MAX;
4252          columnUpperWork_[i] = COIN_DBL_MAX;
4253        }
4254      }
4255    }
4256    for (i = 0; i < numberRows_; i++) {
4257      double multiplier = rhsScale_ * rowScale[i];
4258      double lowerValue = rowLower_[i];
4259      double upperValue = rowUpper_[i];
4260      if (lowerValue > -1.0e20) {
4261        rowLowerWork_[i] = lowerValue * multiplier;
4262        if (upperValue >= 1.0e20) {
4263          rowUpperWork_[i] = COIN_DBL_MAX;
4264        } else {
4265          rowUpperWork_[i] = upperValue * multiplier;
4266          if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4267            if (rowLowerWork_[i] >= 0.0) {
4268              rowUpperWork_[i] = rowLowerWork_[i];
4269            } else if (rowUpperWork_[i] <= 0.0) {
4270              rowLowerWork_[i] = rowUpperWork_[i];
4271            } else {
4272              rowUpperWork_[i] = 0.0;
4273              rowLowerWork_[i] = 0.0;
4274            }
4275          }
4276        }
4277      } else if (upperValue < 1.0e20) {
4278        rowLowerWork_[i] = -COIN_DBL_MAX;
4279        rowUpperWork_[i] = upperValue * multiplier;
4280      } else {
4281        // free
4282        rowLowerWork_[i] = -COIN_DBL_MAX;
4283        rowUpperWork_[i] = COIN_DBL_MAX;
4284      }
4285    }
4286  } else if (rhsScale_ != 1.0) {
4287    for (i = 0; i < numberColumns_; i++) {
4288      double lowerValue = columnLower_[i];
4289      double upperValue = columnUpper_[i];
4290      if (lowerValue > -1.0e20) {
4291        columnLowerWork_[i] = lowerValue * rhsScale_;
4292        if (upperValue >= 1.0e20) {
4293          columnUpperWork_[i] = COIN_DBL_MAX;
4294        } else {
4295          columnUpperWork_[i] = upperValue * rhsScale_;
4296          if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4297            if (columnLowerWork_[i] >= 0.0) {
4298              columnUpperWork_[i] = columnLowerWork_[i];
4299            } else if (columnUpperWork_[i] <= 0.0) {
4300              columnLowerWork_[i] = columnUpperWork_[i];
4301            } else {
4302              columnUpperWork_[i] = 0.0;
4303              columnLowerWork_[i] = 0.0;
4304            }
4305          }
4306        }
4307      } else if (upperValue < 1.0e20) {
4308        columnLowerWork_[i] = -COIN_DBL_MAX;
4309        columnUpperWork_[i] = upperValue * rhsScale_;
4310      } else {
4311        // free
4312        columnLowerWork_[i] = -COIN_DBL_MAX;
4313        columnUpperWork_[i] = COIN_DBL_MAX;
4314      }
4315    }
4316    for (i = 0; i < numberRows_; i++) {
4317      double lowerValue = rowLower_[i];
4318      double upperValue = rowUpper_[i];
4319      if (lowerValue > -1.0e20) {
4320        rowLowerWork_[i] = lowerValue * rhsScale_;
4321        if (upperValue >= 1.0e20) {
4322          rowUpperWork_[i] = COIN_DBL_MAX;
4323        } else {
4324          rowUpperWork_[i] = upperValue * rhsScale_;
4325          if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4326            if (rowLowerWork_[i] >= 0.0) {
4327              rowUpperWork_[i] = rowLowerWork_[i];
4328            } else if (rowUpperWork_[i] <= 0.0) {
4329              rowLowerWork_[i] = rowUpperWork_[i];
4330            } else {
4331              rowUpperWork_[i] = 0.0;
4332              rowLowerWork_[i] = 0.0;
4333            }
4334          }
4335        }
4336      } else if (upperValue < 1.0e20) {
4337        rowLowerWork_[i] = -COIN_DBL_MAX;
4338        rowUpperWork_[i] = upperValue * rhsScale_;
4339      } else {
4340        // free
4341        rowLowerWork_[i] = -COIN_DBL_MAX;
4342        rowUpperWork_[i] = COIN_DBL_MAX;
4343      }
4344    }
4345  } else {
4346    for (i = 0; i < numberColumns_; i++) {
4347      double lowerValue = columnLower_[i];
4348      double upperValue = columnUpper_[i];
4349      if (lowerValue > -1.0e20) {
4350        columnLowerWork_[i] = lowerValue;
4351        if (upperValue >= 1.0e20) {
4352          columnUpperWork_[i] = COIN_DBL_MAX;
4353        } else {
4354          columnUpperWork_[i] = upperValue;
4355          if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4356            if (columnLowerWork_[i] >= 0.0) {
4357              columnUpperWork_[i] = columnLowerWork_[i];
4358            } else if (columnUpperWork_[i] <= 0.0) {
4359              columnLowerWork_[i] = columnUpperWork_[i];
4360            } else {
4361              columnUpperWork_[i] = 0.0;
4362              columnLowerWork_[i] = 0.0;
4363            }
4364          }
4365        }
4366      } else if (upperValue < 1.0e20) {
4367        columnLowerWork_[i] = -COIN_DBL_MAX;
4368        columnUpperWork_[i] = upperValue;
4369      } else {
4370        // free
4371        columnLowerWork_[i] = -COIN_DBL_MAX;
4372        columnUpperWork_[i] = COIN_DBL_MAX;
4373      }
4374    }
4375    for (i = 0; i < numberRows_; i++) {
4376      double lowerValue = rowLower_[i];
4377      double upperValue = rowUpper_[i];
4378      if (lowerValue > -1.0e20) {
4379        rowLowerWork_[i] = lowerValue;
4380        if (upperValue >= 1.0e20) {
4381          rowUpperWork_[i] = COIN_DBL_MAX;
4382        } else {
4383          rowUpperWork_[i] = upperValue;
4384          if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4385            if (rowLowerWork_[i] >= 0.0) {
4386              rowUpperWork_[i] = rowLowerWork_[i];
4387            } else if (rowUpperWork_[i] <= 0.0) {
4388              rowLowerWork_[i] = rowUpperWork_[i];
4389            } else {
4390              rowUpperWork_[i] = 0.0;
4391              rowLowerWork_[i] = 0.0;
4392            }
4393          }
4394        }
4395      } else if (upperValue < 1.0e20) {
4396        rowLowerWork_[i] = -COIN_DBL_MAX;
4397        rowUpperWork_[i] = upperValue;
4398      } else {
4399        // free
4400        rowLowerWork_[i] = -COIN_DBL_MAX;
4401        rowUpperWork_[i] = COIN_DBL_MAX;
4402      }
4403    }
4404  }
4405#ifndef NDEBUG
4406  if ((specialOptions_ & 65536) != 0 && false) {
4407    assert(!initial);
4408    int save = maximumColumns_ + maximumRows_;
4409    for (int i = 0; i < numberTotal; i++) {
4410      assert(fabs(lower_[i] - lower_[i + save]) < 1.0e-5);
4411      assert(fabs(upper_[i] - upper_[i + save]) < 1.0e-5);
4412    }
4413  }
4414#endif
4415}
4416// Does objective
4417void ClpSimplex::createRim4(bool initial)
4418{
4419  int i;
4420  int numberRows2 = numberRows_ + numberExtraRows_;
4421  int numberTotal = numberRows2 + numberColumns_;
4422  if ((specialOptions_ & 65536) != 0 && true) {
4423    assert(!initial);
4424    int save = maximumColumns_ + maximumRows_;
4425    CoinMemcpyN(cost_ + save, numberTotal, cost_);
4426    return;
4427  }
4428  double direction = optimizationDirection_ * objectiveScale_;
4429  const double *obj = objective();
4430  const double *rowScale = rowScale_;
4431  const double *columnScale = columnScale_;
4432  // and also scale by scale factors
4433  if (rowScale) {
4434    if (rowObjective_) {
4435      for (i = 0; i < numberRows_; i++)
4436        rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
4437    } else {
4438      memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4439    }
4440    // If scaled then do all columns later in one loop
4441    if (!initial) {
4442      for (i = 0; i < numberColumns_; i++) {
4443        CoinAssert(fabs(obj[i]) < 1.0e25);
4444        objectiveWork_[i] = obj[i] * direction * columnScale[i];
4445      }
4446    }
4447  } else {
4448    if (rowObjective_) {
4449      for (i = 0; i < numberRows_; i++)
4450        rowObjectiveWork_[i] = rowObjective_[i] * direction;
4451    } else {
4452      memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4453    }
4454    for (i = 0; i < numberColumns_; i++) {
4455      CoinAssert(fabs(obj[i]) < 1.0e25);
4456      objectiveWork_[i] = obj[i] * direction;
4457    }
4458  }
4459}
4460// Does rows and columns and objective
4461void ClpSimplex::createRim5(bool initial)
4462{
4463  createRim4(initial);
4464  createRim1(initial);
4465}
4466void ClpSimplex::deleteRim(int getRidOfFactorizationData)
4467{
4468#ifdef CLP_USER_DRIVEN
4469  eventHandler_->event(ClpEventHandler::beforeDeleteRim);
4470#endif
4471  // Just possible empty problem
4472  int numberRows = numberRows_;
4473  int numberColumns = numberColumns_;
4474  if (!numberRows || !numberColumns) {
4475    numberRows = 0;
4476    if (objective_->type() < 2)
4477      numberColumns = 0;
4478  }
4479  int i;
4480  if (problemStatus_ != 1 && problemStatus_ != 2) {
4481    delete[] ray_;
4482    ray_ = NULL;
4483  }
4484  // set upperOut_ to furthest away from bound so can use in dual for dualBound_
4485  upperOut_ = 1.0;
4486#if 0
4487     {
4488          int nBad = 0;
4489          for (i = 0; i < numberColumns; i++) {
4490               if (lower_[i] == upper_[i] && getColumnStatus(i) == basic)
4491                    nBad++;
4492          }
4493          if (nBad)
4494               printf("yy %d basic fixed\n", nBad);
4495     }
4496#endif
4497  if ((moreSpecialOptions_ & 4194304) != 0) {
4498    // preset tolerances were changed
4499    moreSpecialOptions_ &= ~4194304;
4500    primalTolerance_ = 1.0e-7;
4501    dblParam_[ClpPrimalTolerance] = primalTolerance_;
4502    dualTolerance_ = 1.0e-7;
4503    dblParam_[ClpDualTolerance] = dualTolerance_;
4504  }
4505  // ray may be null if in branch and bound
4506  if (rowScale_ && solution_) {
4507    // Collect infeasibilities
4508    int numberPrimalScaled = 0;
4509    int numberPrimalUnscaled = 0;
4510    int numberDualScaled = 0;
4511    int numberDualUnscaled = 0;
4512    double scaleC = 1.0 / objectiveScale_;
4513    double scaleR = 1.0 / rhsScale_;
4514    const double *inverseScale = inverseColumnScale_;
4515    for (i = 0; i < numberColumns; i++) {
4516      double scaleFactor = columnScale_[i];
4517      double valueScaled = columnActivityWork_[i];
4518      double lowerScaled = columnLowerWork_[i];
4519      double upperScaled = columnUpperWork_[i];
4520      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4521        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4522          numberPrimalScaled++;
4523        else
4524          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4525      }
4526      columnActivity_[i] = valueScaled * scaleFactor * scaleR;
4527      double value = columnActivity_[i];
4528      if (value < columnLower_[i] - primalTolerance_)
4529        numberPrimalUnscaled++;
4530      else if (value > columnUpper_[i] + primalTolerance_)
4531        numberPrimalUnscaled++;
4532      double valueScaledDual = reducedCostWork_[i];
4533      if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4534        numberDualScaled++;
4535      if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4536        numberDualScaled++;
4537      reducedCost_[i] = (valueScaledDual * scaleC) * inverseScale[i];
4538      double valueDual = reducedCost_[i];
4539      if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4540        numberDualUnscaled++;
4541      if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4542        numberDualUnscaled++;
4543    }
4544    inverseScale = inverseRowScale_;
4545    for (i = 0; i < numberRows; i++) {
4546      double scaleFactor = rowScale_[i];
4547      double valueScaled = rowActivityWork_[i];
4548      double lowerScaled = rowLowerWork_[i];
4549      double upperScaled = rowUpperWork_[i];
4550      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4551        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4552          numberPrimalScaled++;
4553        else
4554          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4555      }
4556      rowActivity_[i] = (valueScaled * scaleR) * inverseScale[i];
4557      double value = rowActivity_[i];
4558      if (value < rowLower_[i] - primalTolerance_)
4559        numberPrimalUnscaled++;
4560      else if (value > rowUpper_[i] + primalTolerance_)
4561        numberPrimalUnscaled++;
4562      double valueScaledDual = dual_[i] + rowObjectiveWork_[i];
4563      ;
4564      if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4565        numberDualScaled++;
4566      if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4567        numberDualScaled++;
4568      dual_[i] *= scaleFactor * scaleC;
4569      double valueDual = dual_[i];
4570      if (rowObjective_)
4571        valueDual += rowObjective_[i];
4572      if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4573        numberDualUnscaled++;
4574      if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4575        numberDualUnscaled++;
4576    }
4577    if (!problemStatus_ && !secondaryStatus_) {
4578      // See if we need to set secondary status
4579      if (numberPrimalUnscaled) {
4580        if (numberDualUnscaled)
4581          secondaryStatus_ = 4;
4582        else
4583          secondaryStatus_ = 2;
4584      } else {
4585        if (numberDualUnscaled)
4586          secondaryStatus_ = 3;
4587      }
4588    }
4589    if (problemStatus_ == 2 && ray_) {
4590      for (i = 0; i < numberColumns; i++) {
4591        ray_[i] *= columnScale_[i];
4592      }
4593    } else if (problemStatus_ == 1 && ray_) {
4594      for (i = 0; i < numberRows; i++) {
4595        ray_[i] *= rowScale_[i];
4596      }
4597    }
4598  } else if (rhsScale_ != 1.0 || objectiveScale_ != 1.0) {
4599    // Collect infeasibilities
4600    int numberPrimalScaled = 0;
4601    int numberPrimalUnscaled = 0;
4602    int numberDualScaled = 0;
4603    int numberDualUnscaled = 0;
4604    double scaleC = 1.0 / objectiveScale_;
4605    double scaleR = 1.0 / rhsScale_;
4606    for (i = 0; i < numberColumns; i++) {
4607      double valueScaled = columnActivityWork_[i];
4608      double lowerScaled = columnLowerWork_[i];
4609      double upperScaled = columnUpperWork_[i];
4610      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4611        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4612          numberPrimalScaled++;
4613        else
4614          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4615      }
4616      columnActivity_[i] = valueScaled * scaleR;
4617      double value = columnActivity_[i];
4618      if (value < columnLower_[i] - primalTolerance_)
4619        numberPrimalUnscaled++;
4620      else if (value > columnUpper_[i] + primalTolerance_)
4621        numberPrimalUnscaled++;
4622      double valueScaledDual = reducedCostWork_[i];
4623      if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4624        numberDualScaled++;
4625      if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4626        numberDualScaled++;
4627      reducedCost_[i] = valueScaledDual * scaleC;
4628      double valueDual = reducedCost_[i];
4629      if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4630        numberDualUnscaled++;
4631      if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4632        numberDualUnscaled++;
4633    }
4634    for (i = 0; i < numberRows; i++) {
4635      double valueScaled = rowActivityWork_[i];
4636      double lowerScaled = rowLowerWork_[i];
4637      double upperScaled = rowUpperWork_[i];
4638      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4639        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4640          numberPrimalScaled++;
4641        else
4642          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4643      }
4644      rowActivity_[i] = valueScaled * scaleR;
4645      double value = rowActivity_[i];
4646      if (value < rowLower_[i] - primalTolerance_)
4647        numberPrimalUnscaled++;
4648      else if (value > rowUpper_[i] + primalTolerance_)
4649        numberPrimalUnscaled++;
4650      double valueScaledDual = dual_[i] + rowObjectiveWork_[i];
4651      ;
4652      if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4653        numberDualScaled++;
4654      if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4655        numberDualScaled++;
4656      dual_[i] *= scaleC;
4657      double valueDual = dual_[i];
4658      if (rowObjective_)
4659        valueDual += rowObjective_[i];
4660      if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4661        numberDualUnscaled++;
4662      if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4663        numberDualUnscaled++;
4664    }
4665    if (!problemStatus_ && !secondaryStatus_) {
4666      // See if we need to set secondary status
4667      if (numberPrimalUnscaled) {
4668        if (numberDualUnscaled)
4669          secondaryStatus_ = 4;
4670        else
4671          secondaryStatus_ = 2;
4672      } else {
4673        if (numberDualUnscaled)
4674          secondaryStatus_ = 3;
4675      }
4676    }
4677  } else {
4678    if (columnActivityWork_) {
4679      for (i = 0; i < numberColumns; i++) {
4680        double value = columnActivityWork_[i];
4681        double lower = columnLowerWork_[i];
4682        double upper = columnUpperWork_[i];
4683        if (lower > -1.0e20 || upper < 1.0e20) {
4684          if (value > lower && value < upper)
4685            upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4686        }
4687        columnActivity_[i] = columnActivityWork_[i];
4688        reducedCost_[i] = reducedCostWork_[i];
4689      }
4690      for (i = 0; i < numberRows; i++) {
4691        double value = rowActivityWork_[i];
4692        double lower = rowLowerWork_[i];
4693        double upper = rowUpperWork_[i];
4694        if (lower > -1.0e20 || upper < 1.0e20) {
4695          if (value > lower && value < upper)
4696            upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4697        }
4698        rowActivity_[i] = rowActivityWork_[i];
4699      }
4700    }
4701  }
4702  // switch off scalefactor if auto
4703  if (automaticScale_) {
4704    rhsScale_ = 1.0;
4705    objectiveScale_ = 1.0;
4706  }
4707  if (optimizationDirection_ != 1.0) {
4708    // and modify all dual signs
4709    for (i = 0; i < numberColumns; i++)
4710      reducedCost_[i] *= optimizationDirection_;
4711    for (i = 0; i < numberRows; i++)
4712      dual_[i] *= optimizationDirection_;
4713  }
4714  // scaling may have been turned off
4715  scalingFlag_ = abs(scalingFlag_);
4716  if (getRidOfFactorizationData > 0) {
4717    gutsOfDelete(getRidOfFactorizationData + 1);
4718  } else {
4719    // at least get rid of nonLinearCost_
4720    delete nonLinearCost_;
4721    nonLinearCost_ = NULL;
4722  }
4723  if (!rowObjective_ && problemStatus_ == 0 && objective_->type() == 1 && numberRows && numberColumns) {
4724    // Redo objective value
4725    double objectiveValue = 0.0;
4726    const double *cost = objective();
4727    for (int i = 0; i < numberColumns; i++) {
4728      double value = columnActivity_[i];
4729      objectiveValue += value * cost[i];
4730    }
4731    //if (fabs(objectiveValue_ -objectiveValue*optimizationDirection())>1.0e-5)
4732    //printf("old obj %g new %g\n",objectiveValue_, objectiveValue*optimizationDirection());
4733    objectiveValue_ = objectiveValue * optimizationDirection();
4734  }
4735  // get rid of data
4736  matrix_->generalExpanded(this, 13, scalingFlag_);
4737}
4738void ClpSimplex::setDualBound(double value)
4739{
4740  if (value > 0.0)
4741    dualBound_ = value;
4742}
4743void ClpSimplex::setInfeasibilityCost(double value)
4744{
4745  if (value > 0.0)
4746    infeasibilityCost_ = value;
4747}
4748void ClpSimplex::setNumberRefinements(int value)
4749{
4750  if (value >= 0 && value < 10)
4751    numberRefinements_ = value;
4752}
4753// Sets row pivot choice algorithm in dual
4754void ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot &choice)
4755{
4756  delete dualRowPivot_;
4757  dualRowPivot_ = choice.clone(true);
4758  dualRowPivot_->setModel(this);
4759}
4760// Sets column pivot choice algorithm in primal
4761void ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot &choice)
4762{
4763  delete primalColumnPivot_;
4764  primalColumnPivot_ = choice.clone(true);
4765  primalColumnPivot_->setModel(this);
4766}
4767void ClpSimplex::setFactorization(ClpFactorization &factorization)
4768{
4769  if (factorization_)
4770    factorization_->setFactorization(factorization);
4771  else
4772    factorization_ = new ClpFactorization(factorization,
4773      numberRows_);
4774}
4775
4776// Swaps factorization
4777ClpFactorization *
4778ClpSimplex::swapFactorization(ClpFactorization *factorization)
4779{
4780  ClpFactorization *swap = factorization_;
4781  factorization_ = factorization;
4782  return swap;
4783}
4784// Copies in factorization to existing one
4785void ClpSimplex::copyFactorization(ClpFactorization &factorization)
4786{
4787  *factorization_ = factorization;
4788}
4789/* Perturbation:
4790   -50 to +50 - perturb by this power of ten (-6 sounds good)
4791   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
4792   101 - we are perturbed
4793   102 - don't try perturbing again
4794   default is 100
4795*/
4796void ClpSimplex::setPerturbation(int value)
4797{
4798  if (value <= 102 && value >= -1000) {
4799    perturbation_ = value;
4800  }
4801}
4802// Sparsity on or off
4803bool ClpSimplex::sparseFactorization() const
4804{
4805  return factorization_->sparseThreshold() != 0;
4806}
4807void ClpSimplex::setSparseFactorization(bool value)
4808{
4809  if (value) {
4810    if (!factorization_->sparseThreshold())
4811      factorization_->goSparse();
4812  } else {
4813    factorization_->sparseThreshold(0);
4814  }
4815}
4816void checkCorrect(ClpSimplex * /*model*/, int iRow,
4817  const double *element, const int *rowStart, const int *rowLength,
4818  const int *column,
4819  const double *columnLower_, const double *columnUpper_,
4820  int /*infiniteUpperC*/,
4821  int /*infiniteLowerC*/,
4822  double &maximumUpC,
4823  double &maximumDownC)
4824{
4825  int infiniteUpper = 0;
4826  int infiniteLower = 0;
4827  double maximumUp = 0.0;
4828  double maximumDown = 0.0;
4829  CoinBigIndex rStart = rowStart[iRow];
4830  CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4831  CoinBigIndex j;
4832  double large = 1.0e15;
4833  int iColumn;
4834  // Compute possible lower and upper ranges
4835
4836  for (j = rStart; j < rEnd; ++j) {
4837    double value = element[j];
4838    iColumn = column[j];
4839    if (value > 0.0) {
4840      if (columnUpper_[iColumn] >= large) {
4841        ++infiniteUpper;
4842      } else {
4843        maximumUp += columnUpper_[iColumn] * value;
4844      }
4845      if (columnLower_[iColumn] <= -large) {
4846        ++infiniteLower;
4847      } else {
4848        maximumDown += columnLower_[iColumn] * value;
4849      }
4850    } else if (value < 0.0) {
4851      if (columnUpper_[iColumn] >= large) {
4852        ++infiniteLower;
4853      } else {
4854        maximumDown += columnUpper_[iColumn] * value;
4855      }
4856      if (columnLower_[iColumn] <= -large) {
4857        ++infiniteUpper;
4858      } else {
4859        maximumUp += columnLower_[iColumn] * value;
4860      }
4861    }
4862  }
4863  //assert (infiniteLowerC==infiniteLower);
4864  //assert (infiniteUpperC==infiniteUpper);
4865  if (fabs(maximumUp - maximumUpC) > 1.0e-12 * CoinMax(fabs(maximumUp), fabs(maximumUpC)))
4866    COIN_DETAIL_PRINT(printf("row %d comp up %g, true up %g\n", iRow,
4867      maximumUpC, maximumUp));
4868  if (fabs(maximumDown - maximumDownC) > 1.0e-12 * CoinMax(fabs(maximumDown), fabs(maximumDownC)))
4869    COIN_DETAIL_PRINT(printf("row %d comp down %g, true down %g\n", iRow,
4870      maximumDownC, maximumDown));
4871  maximumUpC = maximumUp;
4872  maximumDownC = maximumDown;
4873}
4874
4875/* Tightens primal bounds to make dual faster.  Unless
4876   fixed, bounds are slightly looser than they could be.
4877   This is to make dual go faster and is probably not needed
4878   with a presolve.  Returns non-zero if problem infeasible
4879
4880   Fudge for branch and bound - put bounds on columns of factor *
4881   largest value (at continuous) - should improve stability
4882   in branch and bound on infeasible branches (0.0 is off)
4883*/
4884int ClpSimplex::tightenPrimalBounds(double factor, int doTight, bool tightIntegers)
4885{
4886
4887  // Get a row copy in standard format
4888  CoinPackedMatrix copy;
4889  copy.setExtraGap(0.0);
4890  copy.setExtraMajor(0.0);
4891  copy.reverseOrderedCopyOf(*matrix());
4892  // Matrix may have been created so get rid of it
4893  matrix_->releasePackedMatrix();
4894  // get matrix data pointers
4895  const int *column = copy.getIndices();
4896  const CoinBigIndex *rowStart = copy.getVectorStarts();
4897  const int *rowLength = copy.getVectorLengths();
4898  const double *element = copy.getElements();
4899  int numberChanged = 1, iPass = 0;
4900  double large = largeValue(); // treat bounds > this as infinite
4901#ifndef NDEBUG
4902  double large2 = 1.0e10 * large;
4903#endif
4904  int numberInfeasible = 0;
4905  int totalTightened = 0;
4906
4907  double tolerance = primalTolerance();
4908
4909  // Save column bounds
4910  double *saveLower = new double[numberColumns_];
4911  CoinMemcpyN(columnLower_, numberColumns_, saveLower);
4912  double *saveUpper = new double[numberColumns_];
4913  CoinMemcpyN(columnUpper_, numberColumns_, saveUpper);
4914
4915  int iRow, iColumn;
4916  // If wanted compute a reasonable dualBound_
4917  if (factor == COIN_DBL_MAX) {
4918    factor = 0.0;
4919    if (dualBound_ == 1.0e10) {
4920      // get largest scaled away from bound
4921      double largest = 1.0e-12;
4922      double largestScaled = 1.0e-12;
4923      int iRow;
4924      for (iRow = 0; iRow < numberRows_; iRow++) {
4925        double value = rowActivity_[iRow];
4926        double above = value - rowLower_[iRow];
4927        double below = rowUpper_[iRow] - value;
4928        if (above < 1.0e12) {
4929          largest = CoinMax(largest, above);
4930        }
4931        if (below < 1.0e12) {
4932          largest = CoinMax(largest, below);
4933        }
4934        if (rowScale_) {
4935          double multiplier = rowScale_[iRow];
4936          above *= multiplier;
4937          below *= multiplier;
4938        }
4939        if (above < 1.0e12) {
4940          largestScaled = CoinMax(largestScaled, above);
4941        }
4942        if (below < 1.0e12) {
4943          largestScaled = CoinMax(largestScaled, below);
4944        }
4945      }
4946
4947      int iColumn;
4948      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4949        double value = columnActivity_[iColumn];
4950        double above = value - columnLower_[iColumn];
4951        double below = columnUpper_[iColumn] - value;
4952        if (above < 1.0e12) {
4953          largest = CoinMax(largest, above);
4954        }
4955        if (below < 1.0e12) {
4956          largest = CoinMax(largest, below);
4957        }
4958        if (columnScale_) {
4959          double multiplier = 1.0 / columnScale_[iColumn];
4960          above *= multiplier;
4961          below *= multiplier;
4962        }
4963        if (above < 1.0e12) {
4964          largestScaled = CoinMax(largestScaled, above);
4965        }
4966        if (below < 1.0e12) {
4967          largestScaled = CoinMax(largestScaled, below);
4968        }
4969      }
4970      std::cout << "Largest (scaled) away from bound " << largestScaled
4971                << " unscaled " << largest << std::endl;
4972      dualBound_ = CoinMax(1.0001e7, CoinMin(100.0 * largest, 1.00001e10));
4973    }
4974  }
4975
4976  // If wanted - tighten column bounds using solution
4977  if (factor) {
4978    double largest = 0.0;
4979    if (factor > 0.0) {
4980      assert(factor > 1.0);
4981      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4982        if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
4983          largest = CoinMax(largest, fabs(columnActivity_[iColumn]));
4984        }
4985      }
4986      largest *= factor;
4987    } else {
4988      // absolute
4989      largest = -factor;
4990    }
4991    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4992      if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
4993        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn], largest);
4994        columnLower_[iColumn] = CoinMax(columnLower_[iColumn], -largest);
4995      }
4996    }
4997  }
4998#define MAXPASS 10
4999
5000  // Loop round seeing if we can tighten bounds
5001  // Would be faster to have a stack of possible rows
5002  // and we put altered rows back on stack
5003  int numberCheck = -1;
5004  while (numberChanged > numberCheck) {
5005
5006    numberChanged = 0; // Bounds tightened this pass
5007
5008    if (iPass == MAXPASS)
5009      break;
5010    iPass++;
5011
5012    for (iRow = 0; iRow < numberRows_; iRow++) {
5013
5014      if (rowLower_[iRow] > -large || rowUpper_[iRow] < large) {
5015
5016        // possible row
5017        int infiniteUpper = 0;
5018        int infiniteLower = 0;
5019        double maximumUp = 0.0;
5020        double maximumDown = 0.0;
5021        double newBound;
5022        CoinBigIndex rStart = rowStart[iRow];
5023        CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
5024        CoinBigIndex j;
5025        // Compute possible lower and upper ranges
5026
5027        for (j = rStart; j < rEnd; ++j) {
5028          double value = element[j];
5029          iColumn = column[j];
5030          if (value > 0.0) {
5031            if (columnUpper_[iColumn] >= large) {
5032              ++infiniteUpper;
5033            } else {
5034              maximumUp += columnUpper_[iColumn] * value;
5035            }
5036            if (columnLower_[iColumn] <= -large) {
5037              ++infiniteLower;
5038            } else {
5039              maximumDown += columnLower_[iColumn] * value;
5040            }
5041          } else if (value < 0.0) {
5042            if (columnUpper_[iColumn] >= large) {
5043              ++infiniteLower;
5044            } else {
5045              maximumDown += columnUpper_[iColumn] * value;
5046            }
5047            if (columnLower_[iColumn] <= -large) {
5048              ++infiniteUpper;
5049            } else {
5050              maximumUp += columnLower_[iColumn] * value;
5051            }
5052          }
5053        }
5054        // Build in a margin of error
5055        maximumUp += 1.0e-8 * fabs(maximumUp);
5056        maximumDown -= 1.0e-8 * fabs(maximumDown);
5057        double maxUp = maximumUp + infiniteUpper * 1.0e31;
5058        double maxDown = maximumDown - infiniteLower * 1.0e31;
5059        if (maxUp <= rowUpper_[iRow] + tolerance && maxDown >= rowLower_[iRow] - tolerance) {
5060
5061          // Row is redundant - make totally free
5062          // NO - can't do this for postsolve
5063          // rowLower_[iRow]=-COIN_DBL_MAX;
5064          // rowUpper_[iRow]=COIN_DBL_MAX;
5065          //printf("Redundant row in presolveX %d\n",iRow);
5066
5067        } else {
5068          if (maxUp < rowLower_[iRow] - 100.0 * tolerance || maxDown > rowUpper_[iRow] + 100.0 * tolerance) {
5069            // problem is infeasible - exit at once
5070            numberInfeasible++;
5071            break;
5072          }
5073          double lower = rowLower_[iRow];
5074          double upper = rowUpper_[iRow];
5075          for (j = rStart; j < rEnd; ++j) {
5076            double value = element[j];
5077            iColumn = column[j];
5078            double nowLower = columnLower_[iColumn];
5079            double nowUpper = columnUpper_[iColumn];
5080            if (value > 0.0) {
5081              // positive value
5082              if (lower > -large) {
5083                if (!infiniteUpper) {
5084                  assert(nowUpper < large2);
5085                  newBound = nowUpper + (lower - maximumUp) / value;
5086                  // relax if original was large
5087                  if (fabs(maximumUp) > 1.0e8)
5088                    newBound -= 1.0e-12 * fabs(maximumUp);
5089                } else if (infiniteUpper == 1 && nowUpper > large) {
5090                  newBound = (lower - maximumUp) / value;
5091                  // relax if original was large
5092                  if (fabs(maximumUp) > 1.0e8)
5093                    newBound -= 1.0e-12 * fabs(maximumUp);
5094                } else {
5095                  newBound = -COIN_DBL_MAX;
5096                }
5097                if (newBound > nowLower + 1.0e-12 && newBound > -large) {
5098                  // Tighten the lower bound
5099                  numberChanged++;
5100                  // check infeasible (relaxed)
5101                  if (nowUpper < newBound) {
5102                    if (nowUpper - newBound < -100.0 * tolerance)
5103                      numberInfeasible++;
5104                    else
5105                      newBound = nowUpper;
5106                  }
5107                  columnLower_[iColumn] = newBound;
5108                  // adjust
5109                  double now;
5110                  if (nowLower < -large) {
5111                    now = 0.0;
5112                    infiniteLower--;
5113                  } else {
5114                    now = nowLower;
5115                  }
5116                  maximumDown += (newBound - now) * value;
5117                  nowLower = newBound;
5118#ifdef DEBUG
5119                  checkCorrect(this, iRow,
5120                    element, rowStart, rowLength,
5121                    column,
5122                    columnLower_, columnUpper_,
5123                    infiniteUpper,
5124                    infiniteLower,
5125                    maximumUp,
5126                    maximumDown);
5127#endif
5128                }
5129              }
5130              if (upper < large) {
5131                if (!infiniteLower) {
5132                  assert(nowLower > -large2);
5133                  newBound = nowLower + (upper - maximumDown) / value;
5134                  // relax if original was large
5135                  if (fabs(maximumDown) > 1.0e8)
5136                    newBound += 1.0e-12 * fabs(maximumDown);
5137                } else if (infiniteLower == 1 && nowLower < -large) {
5138                  newBound = (upper - maximumDown) / value;
5139                  // relax if original was large
5140                  if (fabs(maximumDown) > 1.0e8)
5141                    newBound += 1.0e-12 * fabs(maximumDown);
5142                } else {
5143                  newBound = COIN_DBL_MAX;
5144                }
5145                if (newBound < nowUpper - 1.0e-12 && newBound < large) {
5146                  // Tighten the upper bound
5147                  numberChanged++;
5148                  // check infeasible (relaxed)
5149                  if (nowLower > newBound) {
5150                    if (newBound - nowLower < -100.0 * tolerance)
5151                      numberInfeasible++;
5152                    else
5153                      newBound = nowLower;
5154                  }
5155                  columnUpper_[iColumn] = newBound;
5156                  // adjust
5157                  double now;
5158                  if (nowUpper > large) {
5159                    now = 0.0;
5160                    infiniteUpper--;
5161                  } else {
5162                    now = nowUpper;
5163                  }
5164                  maximumUp += (newBound - now) * value;
5165                  nowUpper = newBound;
5166#ifdef DEBUG
5167                  checkCorrect(this, iRow,
5168                    element, rowStart, rowLength,
5169                    column,
5170                    columnLower_, columnUpper_,
5171                    infiniteUpper,
5172                    infiniteLower,
5173                    maximumUp,
5174                    maximumDown);
5175#endif
5176                }
5177              }
5178            } else {
5179              // negative value
5180              if (lower > -large) {
5181                if (!infiniteUpper) {
5182                  assert(nowLower < large2);
5183                  newBound = nowLower + (lower - maximumUp) / value;
5184                  // relax if original was large
5185                  if (fabs(maximumUp) > 1.0e8)
5186                    newBound += 1.0e-12 * fabs(maximumUp);
5187                } else if (infiniteUpper == 1 && nowLower < -large) {
5188                  newBound = (lower - maximumUp) / value;
5189                  // relax if original was large
5190                  if (fabs(maximumUp) > 1.0e8)
5191                    newBound += 1.0e-12 * fabs(maximumUp);
5192                } else {
5193                  newBound = COIN_DBL_MAX;
5194                }
5195                if (newBound < nowUpper - 1.0e-12 && newBound < large) {
5196                  // Tighten the upper bound
5197                  numberChanged++;
5198                  // check infeasible (relaxed)
5199                  if (nowLower > newBound) {
5200                    if (newBound - nowLower < -100.0 * tolerance)
5201                      numberInfeasible++;
5202                    else
5203                      newBound = nowLower;
5204                  }
5205                  columnUpper_[iColumn] = newBound;
5206                  // adjust
5207                  double now;
5208                  if (nowUpper > large) {
5209                    now = 0.0;
5210                    infiniteLower--;
5211                  } else {
5212                    now = nowUpper;
5213                  }
5214                  maximumDown += (newBound - now) * value;
5215                  nowUpper = newBound;
5216#ifdef DEBUG
5217                  checkCorrect(this, iRow,
5218                    element, rowStart, rowLength,
5219                    column,
5220                    columnLower_, columnUpper_,
5221                    infiniteUpper,
5222                    infiniteLower,
5223                    maximumUp,
5224                    maximumDown);
5225#endif
5226                }
5227              }
5228              if (upper < large) {
5229                if (!infiniteLower) {
5230                  assert(nowUpper < large2);
5231                  newBound = nowUpper + (upper - maximumDown) / value;
5232                  // relax if original was large
5233                  if (fabs(maximumDown) > 1.0e8)
5234                    newBound -= 1.0e-12 * fabs(maximumDown);
5235                } else if (infiniteLower == 1 && nowUpper > large) {
5236                  newBound = (upper - maximumDown) / value;
5237                  // relax if original was large
5238                  if (fabs(maximumDown) > 1.0e8)
5239                    newBound -= 1.0e-12 * fabs(maximumDown);
5240                } else {
5241                  newBound = -COIN_DBL_MAX;
5242                }
5243                if (newBound > nowLower + 1.0e-12 && newBound > -large) {
5244                  // Tighten the lower bound
5245                  numberChanged++;
5246                  // check infeasible (relaxed)
5247                  if (nowUpper < newBound) {
5248                    if (nowUpper - newBound < -100.0 * tolerance)
5249                      numberInfeasible++;
5250                    else
5251                      newBound = nowUpper;
5252                  }
5253                  columnLower_[iColumn] = newBound;
5254                  // adjust
5255                  double now;
5256                  if (nowLower < -large) {
5257                    now = 0.0;
5258                    infiniteUpper--;
5259                  } else {
5260                    now = nowLower;
5261                  }
5262                  maximumUp += (newBound - now) * value;
5263                  nowLower = newBound;
5264#ifdef DEBUG
5265                  checkCorrect(this, iRow,
5266                    element, rowStart, rowLength,
5267                    column,
5268                    columnLower_, columnUpper_,
5269                    infiniteUpper,
5270                    infiniteLower,
5271                    maximumUp,
5272                    maximumDown);
5273#endif
5274                }
5275              }
5276            }
5277          }
5278        }
5279      }
5280    }
5281    totalTightened += numberChanged;
5282    if (iPass == 1)
5283      numberCheck = numberChanged >> 4;
5284    if (numberInfeasible)
5285      break;
5286  }
5287  if (!numberInfeasible) {
5288    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN, messages_)
5289      << totalTightened
5290      << CoinMessageEol;
5291    // Set bounds slightly loose
5292    double useTolerance = 1.0e-3;
5293    if (doTight > 0) {
5294      if (doTight > 10) {
5295        useTolerance = 0.0;
5296      } else {
5297        while (doTight) {
5298          useTolerance *= 0.1;
5299          doTight--;
5300        }
5301      }
5302    }
5303    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5304      if (saveUpper[iColumn] > saveLower[iColumn] + tolerance) {
5305        // Make large bounds stay infinite
5306        if (saveUpper[iColumn] > 1.0e30 && columnUpper_[iColumn] > 1.0e10) {
5307          columnUpper_[iColumn] = COIN_DBL_MAX;
5308        }
5309        if (saveLower[iColumn] < -1.0e30 && columnLower_[iColumn] < -1.0e10) {
5310          columnLower_[iColumn] = -COIN_DBL_MAX;
5311        }
5312#ifdef KEEP_GOING_IF_FIXED
5313        double multiplier = 5.0e-3 * floor(100.0 * randomNumberGenerator_.randomDouble()) + 1.0;
5314        multiplier *= 100.0;
5315#else
5316        double multiplier = 100.0;
5317#endif
5318        if (columnUpper_[iColumn] - columnLower_[iColumn] < useTolerance + 1.0e-8) {
5319          // relax enough so will have correct dj
5320#if 1
5321          columnLower_[iColumn] = CoinMax(saveLower[iColumn],
5322            columnLower_[iColumn] - multiplier * useTolerance);
5323          columnUpper_[iColumn] = CoinMin(saveUpper[iColumn],
5324            columnUpper_[iColumn] + multiplier * useTolerance);
5325#else
5326          if (fabs(columnUpper_[iColumn]) < fabs(columnLower_[iColumn])) {
5327            if (columnUpper_[iColumn] - multiplier * useTolerance > saveLower[iColumn]) {
5328              columnLower_[iColumn] = columnUpper_[iColumn] - multiplier * useTolerance;
5329            } else {
5330              columnLower_[iColumn] = saveLower[iColumn];
5331              columnUpper_[iColumn] = CoinMin(saveUpper[iColumn],
5332                saveLower[iColumn] + multiplier * useTolerance);
5333            }
5334          } else {
5335            if (columnLower_[iColumn] + multiplier * useTolerance < saveUpper[iColumn]) {
5336              columnUpper_[iColumn] = columnLower_[iColumn] + multiplier * useTolerance;
5337            } else {
5338              columnUpper_[iColumn] = saveUpper[iColumn];
5339              columnLower_[iColumn] = CoinMax(saveLower[iColumn],
5340                saveUpper[iColumn] - multiplier * useTolerance);
5341            }
5342          }
5343#endif
5344        } else {
5345          if (columnUpper_[iColumn] < saveUpper[iColumn]) {
5346            // relax a bit
5347            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn] + multiplier * useTolerance,
5348              saveUpper[iColumn]);
5349          }
5350          if (columnLower_[iColumn] > saveLower[iColumn]) {
5351            // relax a bit
5352            columnLower_[iColumn] = CoinMax(columnLower_[iColumn] - multiplier * useTolerance,
5353              saveLower[iColumn]);
5354          }
5355        }
5356      }
5357    }
5358    if (tightIntegers && integerType_) {
5359      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5360        if (integerType_[iColumn]) {
5361          double value;
5362          value = floor(columnLower_[iColumn] + 0.5);
5363          if (fabs(value - columnLower_[iColumn]) > primalTolerance_)
5364            value = ceil(columnLower_[iColumn]);
5365          columnLower_[iColumn] = value;
5366          value = floor(columnUpper_[iColumn] + 0.5);
5367          if (fabs(value - columnUpper_[iColumn]) > primalTolerance_)
5368            value = floor(columnUpper_[iColumn]);
5369          columnUpper_[iColumn] = value;
5370          if (columnLower_[iColumn] > columnUpper_[iColumn])
5371            numberInfeasible++;
5372        }
5373      }
5374      if (numberInfeasible) {
5375        handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
5376          << numberInfeasible
5377          << CoinMessageEol;
5378        // restore column bounds
5379        CoinMemcpyN(saveLower, numberColumns_, columnLower_);
5380        CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
5381      }
5382    }
5383  } else {
5384    handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
5385      << numberInfeasible
5386      << CoinMessageEol;
5387    // restore column bounds
5388    CoinMemcpyN(saveLower, numberColumns_, columnLower_);
5389    CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
5390  }
5391  delete[] saveLower;
5392  delete[] saveUpper;
5393  return (numberInfeasible);
5394}
5395//#define SAVE_AND_RESTORE
5396// dual
5397#include "ClpSimplexDual.hpp"
5398#include "ClpSimplexPrimal.hpp"
5399#ifndef SAVE_AND_RESTORE
5400int ClpSimplex::dual(int ifValuesPass, int startFinishOptions)
5401#else
5402int ClpSimplex::dual(int ifValuesPass, int startFinishOptions)
5403{
5404  // May be empty problem
5405  if (numberRows_ && numberColumns_) {
5406    // Save on file for debug
5407    int returnCode;
5408    returnCode = saveModel("debug.sav");
5409    if (returnCode) {
5410