source: stable/1.17/Clp/src/ClpSimplex.cpp @ 2437

Last change on this file since 2437 was 2437, checked in by forrest, 3 months ago

try and fix eventHandler pointing to correct model

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 436.2 KB
Line 
1/* $Id: ClpSimplex.cpp 2437 2019-03-22 18:40:08Z 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  eventHandler_->setSimplex(this);
156}
157
158// Subproblem constructor
159ClpSimplex::ClpSimplex(const ClpModel *rhs,
160  int numberRows, const int *whichRow,
161  int numberColumns, const int *whichColumn,
162  bool dropNames, bool dropIntegers, bool fixOthers)
163  : ClpModel(rhs, numberRows, whichRow,
164      numberColumns, whichColumn, dropNames, dropIntegers)
165  , bestPossibleImprovement_(0.0)
166  , zeroTolerance_(1.0e-13)
167  , columnPrimalSequence_(-2)
168  , rowPrimalSequence_(-2)
169  , bestObjectiveValue_(-COIN_DBL_MAX)
170  , moreSpecialOptions_(2)
171  , baseIteration_(0)
172  , vectorMode_(0)
173  , primalToleranceToGetOptimal_(-1.0)
174  , largeValue_(1.0e15)
175  , largestPrimalError_(0.0)
176  , largestDualError_(0.0)
177  , alphaAccuracy_(-1.0)
178  , dualBound_(1.0e10)
179  , alpha_(0.0)
180  , theta_(0.0)
181  , lowerIn_(0.0)
182  , valueIn_(0.0)
183  , upperIn_(-COIN_DBL_MAX)
184  , dualIn_(0.0)
185  , lowerOut_(-1)
186  , valueOut_(-1)
187  , upperOut_(-1)
188  , dualOut_(-1)
189  , dualTolerance_(1.0e-7)
190  , primalTolerance_(1.0e-7)
191  , sumDualInfeasibilities_(0.0)
192  , sumPrimalInfeasibilities_(0.0)
193  , infeasibilityCost_(1.0e10)
194  , sumOfRelaxedDualInfeasibilities_(0.0)
195  , sumOfRelaxedPrimalInfeasibilities_(0.0)
196  , acceptablePivot_(1.0e-8)
197  , lower_(NULL)
198  , rowLowerWork_(NULL)
199  , columnLowerWork_(NULL)
200  , upper_(NULL)
201  , rowUpperWork_(NULL)
202  , columnUpperWork_(NULL)
203  , cost_(NULL)
204  , rowObjectiveWork_(NULL)
205  , objectiveWork_(NULL)
206  , sequenceIn_(-1)
207  , directionIn_(-1)
208  , sequenceOut_(-1)
209  , directionOut_(-1)
210  , pivotRow_(-1)
211  , lastGoodIteration_(-100)
212  , dj_(NULL)
213  , rowReducedCost_(NULL)
214  , reducedCostWork_(NULL)
215  , solution_(NULL)
216  , rowActivityWork_(NULL)
217  , columnActivityWork_(NULL)
218  , numberDualInfeasibilities_(0)
219  , numberDualInfeasibilitiesWithoutFree_(0)
220  , numberPrimalInfeasibilities_(100)
221  , numberRefinements_(0)
222  , pivotVariable_(NULL)
223  , factorization_(NULL)
224  , savedSolution_(NULL)
225  , numberTimesOptimal_(0)
226  , disasterArea_(NULL)
227  , changeMade_(1)
228  , algorithm_(0)
229  , forceFactorization_(-1)
230  , perturbation_(100)
231  , nonLinearCost_(NULL)
232  , lastBadIteration_(-999999)
233  , lastFlaggedIteration_(-999999)
234  , numberFake_(0)
235  , numberChanged_(0)
236  , progressFlag_(0)
237  , firstFree_(-1)
238  , numberExtraRows_(0)
239  , maximumBasic_(0)
240  , dontFactorizePivots_(0)
241  , incomingInfeasibility_(1.0)
242  , allowedInfeasibility_(10.0)
243  , automaticScale_(0)
244  , maximumPerturbationSize_(0)
245  , perturbationArray_(NULL)
246  , baseModel_(NULL)
247#ifdef ABC_INHERIT
248  , abcSimplex_(NULL)
249  , abcState_(0)
250#endif
251{
252  int i;
253  for (i = 0; i < 6; i++) {
254    rowArray_[i] = NULL;
255    columnArray_[i] = NULL;
256  }
257  for (i = 0; i < 4; i++) {
258    spareIntArray_[i] = 0;
259    spareDoubleArray_[i] = 0.0;
260  }
261  saveStatus_ = NULL;
262  // get an empty factorization so we can set tolerances etc
263  getEmptyFactorization();
264  // say Steepest pricing
265  dualRowPivot_ = new ClpDualRowSteepest();
266  // say Steepest pricing
267  primalColumnPivot_ = new ClpPrimalColumnSteepest();
268  solveType_ = 1; // say simplex based life form
269  eventHandler_->setSimplex(this);
270  if (fixOthers) {
271    int numberOtherColumns = rhs->numberColumns();
272    int numberOtherRows = rhs->numberRows();
273    double *solution = new double[numberOtherColumns];
274    CoinZeroN(solution, numberOtherColumns);
275    int i;
276    for (i = 0; i < numberColumns; i++) {
277      int iColumn = whichColumn[i];
278      if (solution[iColumn])
279        fixOthers = false; // duplicates
280      solution[iColumn] = 1.0;
281    }
282    if (fixOthers) {
283      const double *otherSolution = rhs->primalColumnSolution();
284      const double *objective = rhs->objective();
285      double offset = 0.0;
286      for (i = 0; i < numberOtherColumns; i++) {
287        if (solution[i]) {
288          solution[i] = 0.0; // in
289        } else {
290          solution[i] = otherSolution[i];
291          offset += objective[i] * otherSolution[i];
292        }
293      }
294      double *rhsModification = new double[numberOtherRows];
295      CoinZeroN(rhsModification, numberOtherRows);
296      rhs->matrix()->times(solution, rhsModification);
297      for (i = 0; i < numberRows; i++) {
298        int iRow = whichRow[i];
299        if (rowLower_[i] > -1.0e20)
300          rowLower_[i] -= rhsModification[iRow];
301        if (rowUpper_[i] < 1.0e20)
302          rowUpper_[i] -= rhsModification[iRow];
303      }
304      delete[] rhsModification;
305      setObjectiveOffset(rhs->objectiveOffset() - offset);
306      // And set objective value to match
307      setObjectiveValue(rhs->objectiveValue());
308    }
309    delete[] solution;
310  }
311}
312// Subproblem constructor
313ClpSimplex::ClpSimplex(const ClpSimplex *rhs,
314  int numberRows, const int *whichRow,
315  int numberColumns, const int *whichColumn,
316  bool dropNames, bool dropIntegers, bool fixOthers)
317  : ClpModel(rhs, numberRows, whichRow,
318      numberColumns, whichColumn, dropNames, dropIntegers)
319  , bestPossibleImprovement_(0.0)
320  , zeroTolerance_(1.0e-13)
321  , columnPrimalSequence_(-2)
322  , rowPrimalSequence_(-2)
323  , bestObjectiveValue_(-COIN_DBL_MAX)
324  , moreSpecialOptions_(2)
325  , baseIteration_(0)
326  , vectorMode_(0)
327  , primalToleranceToGetOptimal_(-1.0)
328  , largeValue_(1.0e15)
329  , largestPrimalError_(0.0)
330  , largestDualError_(0.0)
331  , alphaAccuracy_(-1.0)
332  , dualBound_(1.0e10)
333  , alpha_(0.0)
334  , theta_(0.0)
335  , lowerIn_(0.0)
336  , valueIn_(0.0)
337  , upperIn_(-COIN_DBL_MAX)
338  , dualIn_(0.0)
339  , lowerOut_(-1)
340  , valueOut_(-1)
341  , upperOut_(-1)
342  , dualOut_(-1)
343  , dualTolerance_(rhs->dualTolerance_)
344  , primalTolerance_(rhs->primalTolerance_)
345  , sumDualInfeasibilities_(0.0)
346  , sumPrimalInfeasibilities_(0.0)
347  , infeasibilityCost_(1.0e10)
348  , sumOfRelaxedDualInfeasibilities_(0.0)
349  , sumOfRelaxedPrimalInfeasibilities_(0.0)
350  , acceptablePivot_(1.0e-8)
351  , lower_(NULL)
352  , rowLowerWork_(NULL)
353  , columnLowerWork_(NULL)
354  , upper_(NULL)
355  , rowUpperWork_(NULL)
356  , columnUpperWork_(NULL)
357  , cost_(NULL)
358  , rowObjectiveWork_(NULL)
359  , objectiveWork_(NULL)
360  , sequenceIn_(-1)
361  , directionIn_(-1)
362  , sequenceOut_(-1)
363  , directionOut_(-1)
364  , pivotRow_(-1)
365  , lastGoodIteration_(-100)
366  , dj_(NULL)
367  , rowReducedCost_(NULL)
368  , reducedCostWork_(NULL)
369  , solution_(NULL)
370  , rowActivityWork_(NULL)
371  , columnActivityWork_(NULL)
372  , numberDualInfeasibilities_(0)
373  , numberDualInfeasibilitiesWithoutFree_(0)
374  , numberPrimalInfeasibilities_(100)
375  , numberRefinements_(0)
376  , pivotVariable_(NULL)
377  , factorization_(NULL)
378  , savedSolution_(NULL)
379  , numberTimesOptimal_(0)
380  , disasterArea_(NULL)
381  , changeMade_(1)
382  , algorithm_(0)
383  , forceFactorization_(-1)
384  , perturbation_(100)
385  , nonLinearCost_(NULL)
386  , lastBadIteration_(-999999)
387  , lastFlaggedIteration_(-999999)
388  , numberFake_(0)
389  , numberChanged_(0)
390  , progressFlag_(0)
391  , firstFree_(-1)
392  , numberExtraRows_(0)
393  , maximumBasic_(0)
394  , dontFactorizePivots_(0)
395  , incomingInfeasibility_(1.0)
396  , allowedInfeasibility_(10.0)
397  , automaticScale_(0)
398  , maximumPerturbationSize_(0)
399  , perturbationArray_(NULL)
400  , baseModel_(NULL)
401#ifdef ABC_INHERIT
402  , abcSimplex_(NULL)
403  , abcState_(rhs->abcState_)
404#endif
405{
406  int i;
407  for (i = 0; i < 6; i++) {
408    rowArray_[i] = NULL;
409    columnArray_[i] = NULL;
410  }
411  for (i = 0; i < 4; i++) {
412    spareIntArray_[i] = 0;
413    spareDoubleArray_[i] = 0.0;
414  }
415  saveStatus_ = NULL;
416  eventHandler_->setSimplex(this);
417  factorization_ = new ClpFactorization(*rhs->factorization_, -numberRows_);
418  //factorization_ = new ClpFactorization(*rhs->factorization_,
419  //                            rhs->factorization_->goDenseThreshold());
420  ClpPEDualRowSteepest *pivotDualPE = dynamic_cast< ClpPEDualRowSteepest * >(rhs->dualRowPivot_);
421  if (pivotDualPE) {
422    dualRowPivot_ = new ClpPEDualRowSteepest(pivotDualPE->psi());
423  } else {
424    ClpDualRowDantzig *pivot = dynamic_cast< ClpDualRowDantzig * >(rhs->dualRowPivot_);
425    // say Steepest pricing
426    if (!pivot)
427      dualRowPivot_ = new ClpDualRowSteepest();
428    else
429      dualRowPivot_ = new ClpDualRowDantzig();
430  }
431  ClpPEPrimalColumnSteepest *pivotPrimalPE = dynamic_cast< ClpPEPrimalColumnSteepest * >(rhs->primalColumnPivot_);
432  if (pivotPrimalPE) {
433    primalColumnPivot_ = new ClpPEPrimalColumnSteepest(pivotPrimalPE->psi());
434  } else {
435    // say Steepest pricing
436    primalColumnPivot_ = new ClpPrimalColumnSteepest();
437  }
438  solveType_ = 1; // say simplex based life form
439  if (fixOthers) {
440    int numberOtherColumns = rhs->numberColumns();
441    int numberOtherRows = rhs->numberRows();
442    double *solution = new double[numberOtherColumns];
443    CoinZeroN(solution, numberOtherColumns);
444    int i;
445    for (i = 0; i < numberColumns; i++) {
446      int iColumn = whichColumn[i];
447      if (solution[iColumn])
448        fixOthers = false; // duplicates
449      solution[iColumn] = 1.0;
450    }
451    if (fixOthers) {
452      const double *otherSolution = rhs->primalColumnSolution();
453      const double *objective = rhs->objective();
454      double offset = 0.0;
455      for (i = 0; i < numberOtherColumns; i++) {
456        if (solution[i]) {
457          solution[i] = 0.0; // in
458        } else {
459          solution[i] = otherSolution[i];
460          offset += objective[i] * otherSolution[i];
461        }
462      }
463      double *rhsModification = new double[numberOtherRows];
464      CoinZeroN(rhsModification, numberOtherRows);
465      rhs->matrix()->times(solution, rhsModification);
466      for (i = 0; i < numberRows; i++) {
467        int iRow = whichRow[i];
468        if (rowLower_[i] > -1.0e20)
469          rowLower_[i] -= rhsModification[iRow];
470        if (rowUpper_[i] < 1.0e20)
471          rowUpper_[i] -= rhsModification[iRow];
472      }
473      delete[] rhsModification;
474      setObjectiveOffset(rhs->objectiveOffset() - offset);
475      // And set objective value to match
476      setObjectiveValue(rhs->objectiveValue());
477    }
478    delete[] solution;
479  }
480  if (rhs->maximumPerturbationSize_) {
481    maximumPerturbationSize_ = 2 * numberColumns;
482    perturbationArray_ = new double[maximumPerturbationSize_];
483    for (i = 0; i < numberColumns; i++) {
484      int iColumn = whichColumn[i];
485      perturbationArray_[2 * i] = rhs->perturbationArray_[2 * iColumn];
486      perturbationArray_[2 * i + 1] = rhs->perturbationArray_[2 * iColumn + 1];
487    }
488  }
489}
490// Puts solution back small model
491void ClpSimplex::getbackSolution(const ClpSimplex &smallModel, const int *whichRow, const int *whichColumn)
492{
493  setSumDualInfeasibilities(smallModel.sumDualInfeasibilities());
494  setNumberDualInfeasibilities(smallModel.numberDualInfeasibilities());
495  setSumPrimalInfeasibilities(smallModel.sumPrimalInfeasibilities());
496  setNumberPrimalInfeasibilities(smallModel.numberPrimalInfeasibilities());
497  setNumberIterations(smallModel.numberIterations());
498  setProblemStatus(smallModel.status());
499  setObjectiveValue(smallModel.objectiveValue());
500  const double *solution2 = smallModel.primalColumnSolution();
501  int i;
502  int numberRows2 = smallModel.numberRows();
503  int numberColumns2 = smallModel.numberColumns();
504  const double *dj2 = smallModel.dualColumnSolution();
505  for (i = 0; i < numberColumns2; i++) {
506    int iColumn = whichColumn[i];
507    columnActivity_[iColumn] = solution2[i];
508    reducedCost_[iColumn] = dj2[i];
509    setStatus(iColumn, smallModel.getStatus(i));
510  }
511  const double *dual2 = smallModel.dualRowSolution();
512  memset(dual_, 0, numberRows_ * sizeof(double));
513  for (i = 0; i < numberRows2; i++) {
514    int iRow = whichRow[i];
515    setRowStatus(iRow, smallModel.getRowStatus(i));
516    dual_[iRow] = dual2[i];
517  }
518  CoinZeroN(rowActivity_, numberRows_);
519#if 0
520     if (!problemStatus_) {
521          ClpDisjointCopyN(smallModel.objective(), smallModel.numberColumns_, smallModel.reducedCost_);
522          smallModel.matrix_->transposeTimes(-1.0, smallModel.dual_, smallModel.reducedCost_);
523          for (int i = 0; i < smallModel.numberColumns_; i++) {
524               if (smallModel.getColumnStatus(i) == basic)
525                    assert (fabs(smallModel.reducedCost_[i]) < 1.0e-5);
526          }
527          ClpDisjointCopyN(objective(), numberColumns_, reducedCost_);
528          matrix_->transposeTimes(-1.0, dual_, reducedCost_);
529          for (int i = 0; i < numberColumns_; i++) {
530               if (getColumnStatus(i) == basic)
531                    assert (fabs(reducedCost_[i]) < 1.0e-5);
532          }
533     }
534#endif
535  matrix()->times(columnActivity_, rowActivity_);
536}
537
538//-----------------------------------------------------------------------------
539
540ClpSimplex::~ClpSimplex()
541{
542  setPersistenceFlag(0);
543  gutsOfDelete(0);
544  delete nonLinearCost_;
545}
546//#############################################################################
547void ClpSimplex::setLargeValue(double value)
548{
549  if (value > 0.0 && value < COIN_DBL_MAX)
550    largeValue_ = value;
551}
552int ClpSimplex::gutsOfSolution(double *givenDuals,
553  const double *givenPrimals,
554  bool valuesPass)
555{
556
557  // if values pass, save values of basic variables
558  double *save = NULL;
559  double oldValue = 0.0;
560  double oldLargestPrimalError = largestPrimalError_;
561  double oldLargestDualError = largestDualError_;
562  if (valuesPass) {
563    assert(algorithm_ > 0); // only primal at present
564    assert(nonLinearCost_);
565    int iRow;
566    checkPrimalSolution(rowActivityWork_, columnActivityWork_);
567    // get correct bounds on all variables
568    nonLinearCost_->checkInfeasibilities(primalTolerance_);
569    oldValue = nonLinearCost_->largestInfeasibility();
570    save = new double[numberRows_];
571    for (iRow = 0; iRow < numberRows_; iRow++) {
572      int iPivot = pivotVariable_[iRow];
573      save[iRow] = solution_[iPivot];
574    }
575  }
576  // do work
577  computePrimals(rowActivityWork_, columnActivityWork_);
578  // If necessary - override results
579  if (givenPrimals) {
580    CoinMemcpyN(givenPrimals, numberColumns_, columnActivityWork_);
581    memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
582    times(-1.0, columnActivityWork_, rowActivityWork_);
583  }
584  double objectiveModification = 0.0;
585  if (algorithm_ > 0 && nonLinearCost_ != NULL) {
586#ifdef CLP_USER_DRIVEN
587    eventHandler_->eventWithInfo(ClpEventHandler::goodFactorization, NULL);
588#endif
589    // primal algorithm
590    // get correct bounds on all variables
591    // If  4 bit set - Force outgoing variables to exact bound (primal)
592    if ((specialOptions_ & 4) == 0)
593      nonLinearCost_->checkInfeasibilities(primalTolerance_);
594    else
595      nonLinearCost_->checkInfeasibilities(0.0);
596    objectiveModification += nonLinearCost_->changeInCost();
597    if (nonLinearCost_->numberInfeasibilities())
598      if (handler_->detail(CLP_SIMPLEX_NONLINEAR, messages_) < 100) {
599        handler_->message(CLP_SIMPLEX_NONLINEAR, messages_)
600          << nonLinearCost_->changeInCost()
601          << nonLinearCost_->numberInfeasibilities()
602          << CoinMessageEol;
603      }
604  }
605  if (valuesPass) {
606    double badInfeasibility = nonLinearCost_->largestInfeasibility();
607#ifdef CLP_DEBUG
608    std::cout << "Largest given infeasibility " << oldValue
609              << " now " << nonLinearCost_->largestInfeasibility() << std::endl;
610#endif
611    int numberOut = 0;
612    // But may be very large rhs etc
613    double useError = CoinMin(largestPrimalError_,
614      1.0e5 / maximumAbsElement(solution_, numberRows_ + numberColumns_));
615    if ((oldValue < incomingInfeasibility_ || badInfeasibility > (CoinMax(10.0 * allowedInfeasibility_, 100.0 * oldValue)))
616      && (badInfeasibility > CoinMax(incomingInfeasibility_, allowedInfeasibility_) || useError > 1.0e-3)) {
617      if (algorithm_ > 1) {
618        // nonlinear
619        //printf("Original largest infeas %g, now %g, primalError %g\n",
620        //      oldValue,nonLinearCost_->largestInfeasibility(),
621        //      largestPrimalError_);
622        //printf("going to all slack\n");
623        allSlackBasis(true);
624        CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
625        return 1;
626      }
627      //printf("Original largest infeas %g, now %g, primalError %g\n",
628      //     oldValue,nonLinearCost_->largestInfeasibility(),
629      //     largestPrimalError_);
630      // throw out up to 1000 structurals
631      int maxOut = (allowedInfeasibility_ == 10.0) ? 1000 : 100;
632      int iRow;
633      int *sort = new int[numberRows_];
634      // first put back solution and store difference
635      for (iRow = 0; iRow < numberRows_; iRow++) {
636        int iPivot = pivotVariable_[iRow];
637        double difference = fabs(solution_[iPivot] - save[iRow]);
638        solution_[iPivot] = save[iRow];
639        save[iRow] = difference;
640      }
641      int numberBasic = 0;
642      for (iRow = 0; iRow < numberRows_; iRow++) {
643        int iPivot = pivotVariable_[iRow];
644
645        if (iPivot < numberColumns_) {
646          // column
647          double difference = save[iRow];
648          if (difference > 1.0e-4) {
649            sort[numberOut] = iRow;
650            save[numberOut++] = -difference;
651            if (getStatus(iPivot) == basic)
652              numberBasic++;
653          }
654        }
655      }
656      if (!numberBasic) {
657        //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
658#if 0
659                    allSlackBasis(true);
660                    CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
661#else
662        // allow
663        numberOut = 0;
664#endif
665      }
666      CoinSort_2(save, save + numberOut, sort);
667      numberOut = CoinMin(maxOut, numberOut);
668      for (iRow = 0; iRow < numberOut; iRow++) {
669        int jRow = sort[iRow];
670        int iColumn = pivotVariable_[jRow];
671        setColumnStatus(iColumn, superBasic);
672        setRowStatus(jRow, basic);
673        pivotVariable_[jRow] = jRow + numberColumns_;
674        if (fabs(solution_[iColumn]) > 1.0e10) {
675          if (upper_[iColumn] < 0.0) {
676            solution_[iColumn] = upper_[iColumn];
677          } else if (lower_[iColumn] > 0.0) {
678            solution_[iColumn] = lower_[iColumn];
679          } else {
680            solution_[iColumn] = 0.0;
681          }
682        }
683      }
684      delete[] sort;
685    }
686    delete[] save;
687    if (numberOut)
688      return numberOut;
689  }
690  if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
691    //printf("trying feas pump\n");
692    const char *integerType = integerInformation();
693    assert(integerType);
694    assert(perturbationArray_);
695    CoinZeroN(cost_, numberRows_ + numberColumns_);
696    for (int i = 0; i < numberRows_ - numberRows_; i++) {
697      int iSequence = pivotVariable_[i];
698      if (iSequence < numberColumns_ && integerType[iSequence]) {
699        double lower = lower_[iSequence];
700        double upper = upper_[iSequence];
701        double value = solution_[iSequence];
702        if (value >= lower - primalTolerance_ && value <= upper + primalTolerance_) {
703          double sign;
704          if (value - lower < upper - value)
705            sign = 1.0;
706          else
707            sign = -1.0;
708          cost_[iSequence] = sign * perturbationArray_[iSequence];
709        }
710      }
711    }
712  }
713#if CAN_HAVE_ZERO_OBJ > 1
714  if ((specialOptions_ & 16777216) == 0) {
715#endif
716    computeDuals(givenDuals);
717    if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
718      const char *integerType = integerInformation();
719      // Need to do columns and rows to stay dual feasible
720      for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
721        if (integerType[iSequence] && getStatus(iSequence) != basic) {
722          double djValue = dj_[iSequence];
723          double change = 0.0;
724          if (getStatus(iSequence) == atLowerBound)
725            change = CoinMax(-djValue, 10.0 * perturbationArray_[iSequence]);
726          else if (getStatus(iSequence) == atUpperBound)
727            change = CoinMin(-djValue, -10.0 * perturbationArray_[iSequence]);
728          cost_[iSequence] = change;
729          dj_[iSequence] += change;
730        }
731      }
732    }
733
734    // now check solutions
735    //checkPrimalSolution( rowActivityWork_, columnActivityWork_);
736    //checkDualSolution();
737    checkBothSolutions();
738    objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_);
739#if CAN_HAVE_ZERO_OBJ > 1
740  } else {
741    checkPrimalSolution(rowActivityWork_, columnActivityWork_);
742#ifndef COIN_REUSE_RANDOM
743    memset(dj_, 0, (numberRows_ + numberColumns_) * sizeof(double));
744#else
745    for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
746      double value;
747      switch (getStatus(iSequence)) {
748      case atLowerBound:
749        value = 1.0e-9 * (1.0 + randomNumberGenerator_.randomDouble());
750        break;
751      case atUpperBound:
752        value = -1.0e-9 * (1.0 + randomNumberGenerator_.randomDouble());
753        break;
754      default:
755        value = 0.0;
756        break;
757      }
758    }
759#endif
760    objectiveValue_ = 0.0;
761  }
762#endif
763  if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 || largestDualError_ > 1.0e-2))
764    handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
765      << largestPrimalError_
766      << largestDualError_
767      << CoinMessageEol;
768  if (largestPrimalError_ > 1.0e-1 && numberRows_ > 100 && numberIterations_) {
769    // Change factorization tolerance
770    if (factorization_->zeroTolerance() > 1.0e-18)
771      factorization_->zeroTolerance(1.0e-18);
772  }
773  int returnCode = 0;
774  bool notChanged = true;
775  if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_ < 0 || factorization_->pivotTolerance() < 0.9899999999) && (oldLargestDualError || oldLargestPrimalError)) {
776    double useOldDualError = oldLargestDualError;
777    double useDualError = largestDualError_;
778    if (algorithm_ > 0 && nonLinearCost_ && nonLinearCost_->sumInfeasibilities()) {
779      double factor = CoinMax(1.0, CoinMin(1.0e3, infeasibilityCost_ * 1.0e-6));
780      useOldDualError /= factor;
781      useDualError /= factor;
782    }
783    if ((largestPrimalError_ > 1.0e3 && oldLargestPrimalError * 1.0e2 < largestPrimalError_) || (useDualError > 1.0e3 && useOldDualError * 1.0e2 < useDualError)) {
784      double pivotTolerance = factorization_->pivotTolerance();
785      double factor = (largestPrimalError_ > 1.0e10 || largestDualError_ > 1.0e10)
786        ? 2.0
787        : 1.2;
788      if (pivotTolerance < 0.1)
789        factorization_->pivotTolerance(0.1);
790      else if (pivotTolerance < 0.98999999)
791        factorization_->pivotTolerance(CoinMin(0.99, pivotTolerance * factor));
792      notChanged = pivotTolerance == factorization_->pivotTolerance();
793#ifdef CLP_USEFUL_PRINTOUT
794      if (pivotTolerance < 0.9899999) {
795        printf("Changing pivot tolerance from %g to %g and backtracking\n",
796          pivotTolerance, factorization_->pivotTolerance());
797      }
798      printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n",
799        oldLargestPrimalError, largestPrimalError_,
800        oldLargestDualError, largestDualError_,
801        pivotTolerance);
802#endif
803      if (pivotTolerance < 0.9899999) {
804        largestPrimalError_ = 0.0;
805        largestDualError_ = 0.0;
806        returnCode = -123456789;
807      }
808    }
809  }
810  if (progress_.iterationNumber_[0] > 0 && progress_.iterationNumber_[CLP_PROGRESS - 1] - progress_.iterationNumber_[0] < CLP_PROGRESS * 3 && factorization_->pivotTolerance() < 0.25 && notChanged) {
811    double pivotTolerance = factorization_->pivotTolerance();
812    factorization_->pivotTolerance(pivotTolerance * 1.5);
813#ifdef CLP_USEFUL_PRINTOUT
814    printf("Changing pivot tolerance from %g to %g - inverting too often\n",
815      pivotTolerance, factorization_->pivotTolerance());
816#endif
817  }
818  // Switch off false values pass indicator
819  if (!valuesPass && algorithm_ > 0)
820    firstFree_ = -1;
821  if (handler_->logLevel() == 63)
822    printf("end getsolution algorithm %d status %d npinf %d sum,relaxed %g,%g ndinf %d sum,relaxed %g,%g\n",
823      algorithm_, problemStatus_,
824      numberPrimalInfeasibilities_, sumPrimalInfeasibilities_, sumOfRelaxedPrimalInfeasibilities_,
825      numberDualInfeasibilities_, sumDualInfeasibilities_, sumOfRelaxedDualInfeasibilities_);
826  if ((moreSpecialOptions_ & 8388608) != 0) {
827    if (algorithm_ < 0) {
828      bool doneFiddling = false;
829      // Optimization may make exact test iffy
830      double testTolerance = minimumPrimalTolerance_ + 1.0e-15;
831      while (!doneFiddling) {
832        doneFiddling = true;
833        while (!sumOfRelaxedPrimalInfeasibilities_ && primalTolerance_ > testTolerance) {
834          // feasible - adjust tolerance
835          double saveTolerance = primalTolerance_;
836          primalTolerance_ = CoinMax(0.25 * primalTolerance_,
837            minimumPrimalTolerance_);
838          printf("Resetting primal tolerance from %g to %g\n",
839            saveTolerance, primalTolerance_);
840          dblParam_[ClpPrimalTolerance] = primalTolerance_;
841          moreSpecialOptions_ &= ~8388608;
842          // redo with switch off
843          returnCode = gutsOfSolution(givenDuals, givenPrimals, valuesPass);
844        }
845        if (primalTolerance_ > testTolerance)
846          moreSpecialOptions_ |= 8388608; // back on
847        if ((moreSpecialOptions_ & 8388608) != 0) {
848          assert(numberPrimalInfeasibilities_);
849          // average infeasibility
850          double average = sumPrimalInfeasibilities_ / numberPrimalInfeasibilities_;
851          double minimum = COIN_DBL_MAX;
852          double averageTotal = average;
853          bool firstTime = averageInfeasibility_[0] == COIN_DBL_MAX;
854          for (int i = 0; i < CLP_INFEAS_SAVE - 1; i++) {
855            double value = averageInfeasibility_[i + 1];
856            averageTotal += value;
857            averageInfeasibility_[i] = value;
858            minimum = CoinMin(minimum, value);
859          }
860          averageInfeasibility_[CLP_INFEAS_SAVE - 1] = average;
861          averageTotal /= CLP_INFEAS_SAVE;
862          double oldTolerance = primalTolerance_;
863          if (averageInfeasibility_[0] != COIN_DBL_MAX) {
864            if (firstTime) {
865              primalTolerance_ = CoinMin(0.1, 0.1 * averageTotal);
866              primalTolerance_ = CoinMin(primalTolerance_, average);
867            } else if (primalTolerance_ > 0.1 * minimum) {
868              primalTolerance_ = 0.1 * minimum;
869            }
870            primalTolerance_ = CoinMax(primalTolerance_, minimumPrimalTolerance_);
871          }
872          if (primalTolerance_ != oldTolerance) {
873            printf("Changing primal tolerance from %g to %g\n",
874              oldTolerance, primalTolerance_);
875            moreSpecialOptions_ &= ~8388608;
876            // redo with switch off
877            returnCode = gutsOfSolution(givenDuals, givenPrimals, valuesPass);
878            if (primalTolerance_ > testTolerance)
879              moreSpecialOptions_ |= 8388608 | 4194304;
880            if (!sumOfRelaxedPrimalInfeasibilities_)
881              doneFiddling = false; // over done it
882          }
883        }
884      }
885    }
886  }
887  return returnCode;
888}
889void ClpSimplex::computePrimals(const double *rowActivities,
890  const double *columnActivities)
891{
892
893  //work space
894  CoinIndexedVector *workSpace = rowArray_[0];
895
896  CoinIndexedVector *arrayVector = rowArray_[1];
897  arrayVector->clear();
898  CoinIndexedVector *previousVector = rowArray_[2];
899  previousVector->clear();
900  // accumulate non basic stuff
901
902  int iRow;
903  // order is this way for scaling
904  if (columnActivities != columnActivityWork_)
905    ClpDisjointCopyN(columnActivities, numberColumns_, columnActivityWork_);
906  if (rowActivities != rowActivityWork_)
907    ClpDisjointCopyN(rowActivities, numberRows_, rowActivityWork_);
908  double *array = arrayVector->denseVector();
909  int *index = arrayVector->getIndices();
910  int number = 0;
911  const double *rhsOffset = matrix_->rhsOffset(this, false, true);
912  if (!rhsOffset) {
913    // Use whole matrix every time to make it easier for ClpMatrixBase
914    // So zero out basic
915    for (iRow = 0; iRow < numberRows_; iRow++) {
916      int iPivot = pivotVariable_[iRow];
917      assert(iPivot >= 0);
918      solution_[iPivot] = 0.0;
919#ifdef CLP_INVESTIGATE
920      assert(getStatus(iPivot) == basic);
921#endif
922    }
923    // Extended solution before "update"
924    matrix_->primalExpanded(this, 0);
925    times(-1.0, columnActivityWork_, array);
926    for (iRow = 0; iRow < numberRows_; iRow++) {
927      double value = array[iRow] + rowActivityWork_[iRow];
928      if (value) {
929        array[iRow] = value;
930        index[number++] = iRow;
931      } else {
932        array[iRow] = 0.0;
933      }
934    }
935  } else {
936    // we have an effective rhs lying around
937    // zero out basic (really just for slacks)
938    for (iRow = 0; iRow < numberRows_; iRow++) {
939      int iPivot = pivotVariable_[iRow];
940      solution_[iPivot] = 0.0;
941    }
942    for (iRow = 0; iRow < numberRows_; iRow++) {
943      double value = rhsOffset[iRow] + rowActivityWork_[iRow];
944      if (value) {
945        array[iRow] = value;
946        index[number++] = iRow;
947      } else {
948        array[iRow] = 0.0;
949      }
950    }
951  }
952  arrayVector->setNumElements(number);
953#ifdef CLP_DEBUG
954  if (numberIterations_ == -3840) {
955    int i;
956    for (i = 0; i < numberRows_ + numberColumns_; i++)
957      printf("%d status %d\n", i, status_[i]);
958    printf("xxxxx1\n");
959    for (i = 0; i < numberRows_; i++)
960      if (array[i])
961        printf("%d rhs %g\n", i, array[i]);
962    printf("xxxxx2\n");
963    for (i = 0; i < numberRows_ + numberColumns_; i++)
964      if (getStatus(i) != basic)
965        printf("%d non basic %g %g %g\n", i, lower_[i], solution_[i], upper_[i]);
966    printf("xxxxx3\n");
967  }
968#endif
969  // Ftran adjusted RHS and iterate to improve accuracy
970  double lastError = COIN_DBL_MAX;
971  int iRefine;
972  CoinIndexedVector *thisVector = arrayVector;
973  CoinIndexedVector *lastVector = previousVector;
974#if 0
975     static double * xsave=NULL;
976     {
977       double * xx = thisVector->denseVector();
978       double largest=0.0;
979       int iLargest=-1;
980       for (int i=0;i<numberRows_;i++) {
981         if (fabs(xx[i])>largest) {
982           largest=fabs(xx[i]);
983           iLargest=i;
984         }
985       }
986       printf("largest incoming rhs %g on row %d\n",largest,iLargest);
987     }
988     if (numberIterations_<-40722) {
989       double * xx = thisVector->denseVector();
990       if (xsave) {
991         double * sol = xsave+numberRows_;
992         double largest=0.0;
993         int iLargest=-1;
994         for (int i=0;i<numberRows_;i++) {
995           if (fabs(xx[i]-xsave[i])>largest) {
996             largest=fabs(xx[i]-xsave[i]);
997             iLargest=i;
998           }
999         }
1000         printf("error %g on row %d\n",largest,iLargest);
1001         largest=0.0;
1002         iLargest=-1;
1003         for (int i=0;i<numberColumns_;i++) {
1004           if (fabs(solution_[i]-sol[i])>largest) {
1005             largest=fabs(solution_[i]-sol[i]);
1006             iLargest=i;
1007           }
1008         }
1009         printf("error %g on col %d\n",largest,iLargest);
1010       } else {
1011         xsave=new double[2*numberRows_+numberColumns_];
1012       }
1013       memcpy(xsave,xx,numberRows_*sizeof(double));
1014       memcpy(xsave+numberRows_,solution_,(numberRows_+numberColumns_)*sizeof(double));
1015     }
1016#endif
1017  //printf ("ZZ0 n before %d",number);
1018  if (number)
1019    factorization_->updateColumn(workSpace, thisVector);
1020  //printf(" - after %d\n",thisVector->getNumElements());
1021  double *work = workSpace->denseVector();
1022#ifdef CLP_DEBUG
1023  if (numberIterations_ == -3840) {
1024    int i;
1025    for (i = 0; i < numberRows_; i++)
1026      if (array[i])
1027        printf("%d after rhs %g\n", i, array[i]);
1028    printf("xxxxx4\n");
1029  }
1030#endif
1031  bool goodSolution = true;
1032  for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
1033
1034    int numberIn = thisVector->getNumElements();
1035    int *indexIn = thisVector->getIndices();
1036    double *arrayIn = thisVector->denseVector();
1037    // put solution in correct place
1038    if (!rhsOffset) {
1039      int j;
1040      for (j = 0; j < numberIn; j++) {
1041        iRow = indexIn[j];
1042        int iPivot = pivotVariable_[iRow];
1043        solution_[iPivot] = arrayIn[iRow];
1044        //assert (fabs(solution_[iPivot])<1.0e100);
1045      }
1046    } else {
1047      for (iRow = 0; iRow < numberRows_; iRow++) {
1048        int iPivot = pivotVariable_[iRow];
1049        solution_[iPivot] = arrayIn[iRow];
1050        //assert (fabs(solution_[iPivot])<1.0e100);
1051      }
1052    }
1053    // Extended solution after "update"
1054    matrix_->primalExpanded(this, 1);
1055    // check Ax == b  (for all)
1056    // signal column generated matrix to just do basic (and gub)
1057    unsigned int saveOptions = specialOptions();
1058    setSpecialOptions(16);
1059    times(-1.0, columnActivityWork_, work);
1060    setSpecialOptions(saveOptions);
1061    largestPrimalError_ = 0.0;
1062    double multiplier = 131072.0;
1063    for (iRow = 0; iRow < numberRows_; iRow++) {
1064      double value = work[iRow] + rowActivityWork_[iRow];
1065      work[iRow] = value * multiplier;
1066      if (fabs(value) > largestPrimalError_) {
1067        largestPrimalError_ = fabs(value);
1068      }
1069    }
1070    if (largestPrimalError_ >= lastError) {
1071      // restore
1072      CoinIndexedVector *temp = thisVector;
1073      thisVector = lastVector;
1074      lastVector = temp;
1075      goodSolution = false;
1076      break;
1077    }
1078    if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
1079      // try and make better
1080      // save this
1081      CoinIndexedVector *temp = thisVector;
1082      thisVector = lastVector;
1083      lastVector = temp;
1084      int *indexOut = thisVector->getIndices();
1085      int number = 0;
1086      array = thisVector->denseVector();
1087      thisVector->clear();
1088      for (iRow = 0; iRow < numberRows_; iRow++) {
1089        double value = work[iRow];
1090        if (value) {
1091          array[iRow] = value;
1092          indexOut[number++] = iRow;
1093          work[iRow] = 0.0;
1094        }
1095      }
1096      thisVector->setNumElements(number);
1097      lastError = largestPrimalError_;
1098      //printf ("ZZ%d n before %d",iRefine+1,number);
1099      factorization_->updateColumn(workSpace, thisVector);
1100      //printf(" - after %d\n",thisVector->getNumElements());
1101      multiplier = 1.0 / multiplier;
1102      double *previous = lastVector->denseVector();
1103      number = 0;
1104      for (iRow = 0; iRow < numberRows_; iRow++) {
1105        double value = previous[iRow] + multiplier * array[iRow];
1106        if (value) {
1107          array[iRow] = value;
1108          indexOut[number++] = iRow;
1109        } else {
1110          array[iRow] = 0.0;
1111        }
1112      }
1113      thisVector->setNumElements(number);
1114    } else {
1115      break;
1116    }
1117  }
1118
1119  // solution as accurate as we are going to get
1120  ClpFillN(work, numberRows_, 0.0);
1121  if (!goodSolution) {
1122    array = thisVector->denseVector();
1123    // put solution in correct place
1124    for (iRow = 0; iRow < numberRows_; iRow++) {
1125      int iPivot = pivotVariable_[iRow];
1126      solution_[iPivot] = array[iRow];
1127      //assert (fabs(solution_[iPivot])<1.0e100);
1128    }
1129  }
1130  arrayVector->clear();
1131  previousVector->clear();
1132#ifdef CLP_DEBUG
1133  if (numberIterations_ == -3840) {
1134    exit(77);
1135  }
1136#endif
1137}
1138// now dual side
1139void ClpSimplex::computeDuals(double *givenDjs)
1140{
1141#ifndef SLIM_CLP
1142  if (objective_->type() == 1 || !objective_->activated()) {
1143#endif
1144    // Linear
1145    //work space
1146    CoinIndexedVector *workSpace = rowArray_[0];
1147
1148    CoinIndexedVector *arrayVector = rowArray_[1];
1149    arrayVector->clear();
1150    CoinIndexedVector *previousVector = rowArray_[2];
1151    previousVector->clear();
1152    int iRow;
1153#ifdef CLP_DEBUG
1154    workSpace->checkClear();
1155#endif
1156    double *array = arrayVector->denseVector();
1157    int *index = arrayVector->getIndices();
1158    int number = 0;
1159    if (!givenDjs) {
1160      for (iRow = 0; iRow < numberRows_; iRow++) {
1161        int iPivot = pivotVariable_[iRow];
1162        double value = cost_[iPivot];
1163        if (value) {
1164          array[iRow] = value;
1165          index[number++] = iRow;
1166        }
1167      }
1168    } else {
1169      // dual values pass - djs may not be zero
1170      for (iRow = 0; iRow < numberRows_; iRow++) {
1171        int iPivot = pivotVariable_[iRow];
1172        // make sure zero if done
1173        if (!pivoted(iPivot))
1174          givenDjs[iPivot] = 0.0;
1175        double value = cost_[iPivot] - givenDjs[iPivot];
1176        if (value) {
1177          array[iRow] = value;
1178          index[number++] = iRow;
1179        }
1180      }
1181    }
1182    arrayVector->setNumElements(number);
1183    // Extended duals before "updateTranspose"
1184    matrix_->dualExpanded(this, arrayVector, givenDjs, 0);
1185
1186    // Btran basic costs and get as accurate as possible
1187    double lastError = COIN_DBL_MAX;
1188    int iRefine;
1189    double *work = workSpace->denseVector();
1190    CoinIndexedVector *thisVector = arrayVector;
1191    CoinIndexedVector *lastVector = previousVector;
1192    factorization_->updateColumnTranspose(workSpace, thisVector);
1193
1194    for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
1195      // check basic reduced costs zero
1196      largestDualError_ = 0.0;
1197      if (!numberExtraRows_) {
1198        // Just basic
1199        int *index2 = workSpace->getIndices();
1200        // use reduced costs for slacks as work array
1201        double *work2 = reducedCostWork_ + numberColumns_;
1202        int numberStructurals = 0;
1203        for (iRow = 0; iRow < numberRows_; iRow++) {
1204          int iPivot = pivotVariable_[iRow];
1205          if (iPivot < numberColumns_)
1206            index2[numberStructurals++] = iPivot;
1207        }
1208        matrix_->listTransposeTimes(this, array, index2, numberStructurals, work2);
1209        numberStructurals = 0;
1210        if (!givenDjs) {
1211          for (iRow = 0; iRow < numberRows_; iRow++) {
1212            int iPivot = pivotVariable_[iRow];
1213            double value;
1214            if (iPivot >= numberColumns_) {
1215              // slack
1216              value = rowObjectiveWork_[iPivot - numberColumns_]
1217                + array[iPivot - numberColumns_];
1218            } else {
1219              // column
1220              value = objectiveWork_[iPivot] - work2[numberStructurals++];
1221            }
1222            work[iRow] = value;
1223            if (fabs(value) > largestDualError_) {
1224              largestDualError_ = fabs(value);
1225            }
1226          }
1227        } else {
1228          for (iRow = 0; iRow < numberRows_; iRow++) {
1229            int iPivot = pivotVariable_[iRow];
1230            if (iPivot >= numberColumns_) {
1231              // slack
1232              work[iRow] = rowObjectiveWork_[iPivot - numberColumns_]
1233                + array[iPivot - numberColumns_] - givenDjs[iPivot];
1234            } else {
1235              // column
1236              work[iRow] = objectiveWork_[iPivot] - work2[numberStructurals++]
1237                - givenDjs[iPivot];
1238            }
1239            if (fabs(work[iRow]) > largestDualError_) {
1240              largestDualError_ = fabs(work[iRow]);
1241              //assert (largestDualError_<1.0e-7);
1242              //if (largestDualError_>1.0e-7)
1243              //printf("large dual error %g\n",largestDualError_);
1244            }
1245          }
1246        }
1247      } else {
1248        // extra rows - be more careful
1249#if 1
1250        // would be faster to do just for basic but this reduces code
1251        ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
1252        transposeTimes(-1.0, array, reducedCostWork_);
1253#else
1254      // Just basic
1255      int *index2 = workSpace->getIndices();
1256      int numberStructurals = 0;
1257      for (iRow = 0; iRow < numberRows_; iRow++) {
1258        int iPivot = pivotVariable_[iRow];
1259        if (iPivot < numberColumns_)
1260          index2[numberStructurals++] = iPivot;
1261      }
1262      matrix_->listTransposeTimes(this, array, index2, numberStructurals, work);
1263      for (iRow = 0; iRow < numberStructurals; iRow++) {
1264        int iPivot = index2[iRow];
1265        reducedCostWork_[iPivot] = objectiveWork_[iPivot] - work[iRow];
1266      }
1267#endif
1268        // update by duals on sets
1269        matrix_->dualExpanded(this, NULL, NULL, 1);
1270        if (!givenDjs) {
1271          for (iRow = 0; iRow < numberRows_; iRow++) {
1272            int iPivot = pivotVariable_[iRow];
1273            double value;
1274            if (iPivot >= numberColumns_) {
1275              // slack
1276              value = rowObjectiveWork_[iPivot - numberColumns_]
1277                + array[iPivot - numberColumns_];
1278            } else {
1279              // column
1280              value = reducedCostWork_[iPivot];
1281            }
1282            work[iRow] = value;
1283            if (fabs(value) > largestDualError_) {
1284              largestDualError_ = fabs(value);
1285            }
1286          }
1287        } else {
1288          for (iRow = 0; iRow < numberRows_; iRow++) {
1289            int iPivot = pivotVariable_[iRow];
1290            if (iPivot >= numberColumns_) {
1291              // slack
1292              work[iRow] = rowObjectiveWork_[iPivot - numberColumns_]
1293                + array[iPivot - numberColumns_] - givenDjs[iPivot];
1294            } else {
1295              // column
1296              work[iRow] = reducedCostWork_[iPivot] - givenDjs[iPivot];
1297            }
1298            if (fabs(work[iRow]) > largestDualError_) {
1299              largestDualError_ = fabs(work[iRow]);
1300              //assert (largestDualError_<1.0e-7);
1301              //if (largestDualError_>1.0e-7)
1302              //printf("large dual error %g\n",largestDualError_);
1303            }
1304          }
1305        }
1306      }
1307      if (largestDualError_ >= lastError) {
1308        // restore
1309        CoinIndexedVector *temp = thisVector;
1310        thisVector = lastVector;
1311        lastVector = temp;
1312        break;
1313      }
1314      if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
1315        && !givenDjs) {
1316        // try and make better
1317        // save this
1318        CoinIndexedVector *temp = thisVector;
1319        thisVector = lastVector;
1320        lastVector = temp;
1321        int *indexOut = thisVector->getIndices();
1322        int number = 0;
1323        array = thisVector->denseVector();
1324        thisVector->clear();
1325        double multiplier = 131072.0;
1326        for (iRow = 0; iRow < numberRows_; iRow++) {
1327          double value = multiplier * work[iRow];
1328          if (value) {
1329            array[iRow] = value;
1330            indexOut[number++] = iRow;
1331            work[iRow] = 0.0;
1332          }
1333          work[iRow] = 0.0;
1334        }
1335        thisVector->setNumElements(number);
1336        lastError = largestDualError_;
1337        factorization_->updateColumnTranspose(workSpace, thisVector);
1338        multiplier = 1.0 / multiplier;
1339        double *previous = lastVector->denseVector();
1340        number = 0;
1341        for (iRow = 0; iRow < numberRows_; iRow++) {
1342          double value = previous[iRow] + multiplier * array[iRow];
1343          if (value) {
1344            array[iRow] = value;
1345            indexOut[number++] = iRow;
1346          } else {
1347            array[iRow] = 0.0;
1348          }
1349        }
1350        thisVector->setNumElements(number);
1351      } else {
1352        break;
1353      }
1354    }
1355    // now look at dual solution
1356    array = thisVector->denseVector();
1357    for (iRow = 0; iRow < numberRows_; iRow++) {
1358      // slack
1359      double value = array[iRow];
1360      dual_[iRow] = value;
1361      value += rowObjectiveWork_[iRow];
1362      rowReducedCost_[iRow] = value;
1363    }
1364    // can use work if problem scaled (for better cache)
1365    ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix_);
1366    double *saveRowScale = rowScale_;
1367    //double * saveColumnScale = columnScale_;
1368    if (scaledMatrix_) {
1369      rowScale_ = NULL;
1370      clpMatrix = scaledMatrix_;
1371    }
1372    if (clpMatrix && (clpMatrix->flags() & 2) == 0) {
1373      CoinIndexedVector *cVector = columnArray_[0];
1374      int *whichColumn = cVector->getIndices();
1375      assert(!cVector->getNumElements());
1376      int n = 0;
1377      for (int i = 0; i < numberColumns_; i++) {
1378        if (getColumnStatus(i) != basic) {
1379          whichColumn[n++] = i;
1380          reducedCostWork_[i] = objectiveWork_[i];
1381        } else {
1382          reducedCostWork_[i] = 0.0;
1383        }
1384      }
1385      if (numberRows_ > 4000)
1386        clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
1387          rowScale_, columnScale_, work);
1388      else
1389        clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
1390          rowScale_, columnScale_, NULL);
1391    } else {
1392      ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
1393      if (numberRows_ > 4000)
1394        matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
1395          rowScale_, columnScale_, work);
1396      else
1397        matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
1398          rowScale_, columnScale_, NULL);
1399    }
1400    rowScale_ = saveRowScale;
1401    //columnScale_ = saveColumnScale;
1402    ClpFillN(work, numberRows_, 0.0);
1403    // Extended duals and check dual infeasibility
1404    if (!matrix_->skipDualCheck() || algorithm_ < 0 || problemStatus_ != -2)
1405      matrix_->dualExpanded(this, NULL, NULL, 2);
1406    // If necessary - override results
1407    if (givenDjs) {
1408      // restore accurate duals
1409      CoinMemcpyN(dj_, (numberRows_ + numberColumns_), givenDjs);
1410    }
1411    arrayVector->clear();
1412    previousVector->clear();
1413#ifndef SLIM_CLP
1414  } else {
1415    // Nonlinear
1416    objective_->reducedGradient(this, dj_, false);
1417    // get dual_ by moving from reduced costs for slacks
1418    CoinMemcpyN(dj_ + numberColumns_, numberRows_, dual_);
1419  }
1420#endif
1421}
1422/* Given an existing factorization computes and checks
1423   primal and dual solutions.  Uses input arrays for variables at
1424   bounds.  Returns feasibility states */
1425int ClpSimplex::getSolution(const double * /*rowActivities*/,
1426  const double * /*columnActivities*/)
1427{
1428  if (!factorization_->status()) {
1429    // put in standard form
1430    createRim(7 + 8 + 16 + 32, false, -1);
1431    if (pivotVariable_[0] < 0)
1432      internalFactorize(0);
1433    // do work
1434    gutsOfSolution(NULL, NULL);
1435    // release extra memory
1436    deleteRim(0);
1437  }
1438  return factorization_->status();
1439}
1440/* Given an existing factorization computes and checks
1441   primal and dual solutions.  Uses current problem arrays for
1442   bounds.  Returns feasibility states */
1443int ClpSimplex::getSolution()
1444{
1445  double *rowActivities = new double[numberRows_];
1446  double *columnActivities = new double[numberColumns_];
1447  ClpDisjointCopyN(rowActivityWork_, numberRows_, rowActivities);
1448  ClpDisjointCopyN(columnActivityWork_, numberColumns_, columnActivities);
1449  int status = getSolution(rowActivities, columnActivities);
1450  delete[] rowActivities;
1451  delete[] columnActivities;
1452  return status;
1453}
1454// Factorizes using current basis.  This is for external use
1455// Return codes are as from ClpFactorization
1456int ClpSimplex::factorize()
1457{
1458  // put in standard form
1459  createRim(7 + 8 + 16 + 32, false);
1460  // do work
1461  int status = internalFactorize(-1);
1462  // release extra memory
1463  deleteRim(0);
1464
1465  return status;
1466}
1467// Clean up status
1468void ClpSimplex::cleanStatus()
1469{
1470  int iRow, iColumn;
1471  int numberBasic = 0;
1472  // make row activities correct
1473  memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
1474  times(1.0, columnActivityWork_, rowActivityWork_);
1475  if (!status_)
1476    createStatus();
1477  for (iRow = 0; iRow < numberRows_; iRow++) {
1478    if (getRowStatus(iRow) == basic)
1479      numberBasic++;
1480    else {
1481      setRowStatus(iRow, superBasic);
1482      // but put to bound if close
1483      if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
1484        <= primalTolerance_) {
1485        rowActivityWork_[iRow] = rowLowerWork_[iRow];
1486        setRowStatus(iRow, atLowerBound);
1487      } else if (fabs(rowActivityWork_[iRow] - rowUpperWork_[iRow])
1488        <= primalTolerance_) {
1489        rowActivityWork_[iRow] = rowUpperWork_[iRow];
1490        setRowStatus(iRow, atUpperBound);
1491      }
1492    }
1493  }
1494  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1495    if (getColumnStatus(iColumn) == basic) {
1496      if (numberBasic == numberRows_) {
1497        // take out of basis
1498        setColumnStatus(iColumn, superBasic);
1499        // but put to bound if close
1500        if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1501          <= primalTolerance_) {
1502          columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1503          setColumnStatus(iColumn, atLowerBound);
1504        } else if (fabs(columnActivityWork_[iColumn]
1505                     - columnUpperWork_[iColumn])
1506          <= primalTolerance_) {
1507          columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1508          setColumnStatus(iColumn, atUpperBound);
1509        }
1510      } else
1511        numberBasic++;
1512    } else {
1513      setColumnStatus(iColumn, superBasic);
1514      // but put to bound if close
1515      if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1516        <= primalTolerance_) {
1517        columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1518        setColumnStatus(iColumn, atLowerBound);
1519      } else if (fabs(columnActivityWork_[iColumn]
1520                   - columnUpperWork_[iColumn])
1521        <= primalTolerance_) {
1522        columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1523        setColumnStatus(iColumn, atUpperBound);
1524      }
1525    }
1526  }
1527}
1528
1529/* Factorizes using current basis.
1530   solveType - 1 iterating, 0 initial, -1 external
1531   - 2 then iterating but can throw out of basis
1532   If 10 added then in primal values pass
1533   Return codes are as from ClpFactorization unless initial factorization
1534   when total number of singularities is returned.
1535   Special case is numberRows_+1 -> all slack basis.
1536*/
1537int ClpSimplex::internalFactorize(int solveType)
1538{
1539  int iRow, iColumn;
1540  int totalSlacks = numberRows_;
1541  if (!status_)
1542    createStatus();
1543
1544  bool valuesPass = false;
1545  if (solveType >= 10) {
1546    valuesPass = true;
1547    solveType -= 10;
1548  }
1549#ifdef CLP_DEBUG
1550  if (solveType > 0) {
1551    int numberFreeIn = 0, numberFreeOut = 0;
1552    double biggestDj = 0.0;
1553    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1554      switch (getColumnStatus(iColumn)) {
1555
1556      case basic:
1557        if (columnLower_[iColumn] < -largeValue_
1558          && columnUpper_[iColumn] > largeValue_)
1559          numberFreeIn++;
1560        break;
1561      default:
1562        if (columnLower_[iColumn] < -largeValue_
1563          && columnUpper_[iColumn] > largeValue_) {
1564          numberFreeOut++;
1565          biggestDj = CoinMax(fabs(dj_[iColumn]), biggestDj);
1566        }
1567        break;
1568      }
1569    }
1570    if (numberFreeIn + numberFreeOut)
1571      printf("%d in basis, %d out - largest dj %g\n",
1572        numberFreeIn, numberFreeOut, biggestDj);
1573  }
1574#endif
1575  if (solveType <= 0) {
1576    // Make sure everything is clean
1577    for (iRow = 0; iRow < numberRows_; iRow++) {
1578      if (getRowStatus(iRow) == isFixed) {
1579        // double check fixed
1580        if (rowUpperWork_[iRow] > rowLowerWork_[iRow])
1581          setRowStatus(iRow, atLowerBound);
1582      } else if (getRowStatus(iRow) == isFree) {
1583        // may not be free after all
1584        if (rowLowerWork_[iRow] > -largeValue_ || rowUpperWork_[iRow] < largeValue_)
1585          setRowStatus(iRow, superBasic);
1586      }
1587    }
1588    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1589      if (getColumnStatus(iColumn) == isFixed) {
1590        // double check fixed
1591        if (columnUpperWork_[iColumn] > columnLowerWork_[iColumn])
1592          setColumnStatus(iColumn, atLowerBound);
1593      } else if (getColumnStatus(iColumn) == isFree) {
1594        // may not be free after all
1595        if (columnLowerWork_[iColumn] > -largeValue_ || columnUpperWork_[iColumn] < largeValue_)
1596          setColumnStatus(iColumn, superBasic);
1597      }
1598    }
1599    if (!valuesPass) {
1600      // not values pass so set to bounds
1601      bool allSlack = true;
1602      if (status_) {
1603        for (iRow = 0; iRow < numberRows_; iRow++) {
1604          if (getRowStatus(iRow) != basic) {
1605            allSlack = false;
1606            break;
1607          }
1608        }
1609      }
1610      if (!allSlack) {
1611        //#define CLP_INVESTIGATE2
1612#ifdef CLP_INVESTIGATE3
1613        int numberTotal = numberRows_ + numberColumns_;
1614        double *saveSol = valuesPass ? CoinCopyOfArray(solution_, numberTotal) : NULL;
1615#endif
1616        // set values from warm start (if sensible)
1617        int numberBasic = 0;
1618        for (iRow = 0; iRow < numberRows_; iRow++) {
1619          switch (getRowStatus(iRow)) {
1620
1621          case basic:
1622            numberBasic++;
1623            break;
1624          case atUpperBound:
1625            rowActivityWork_[iRow] = rowUpperWork_[iRow];
1626            if (rowActivityWork_[iRow] > largeValue_) {
1627              if (rowLowerWork_[iRow] > -largeValue_) {
1628                rowActivityWork_[iRow] = rowLowerWork_[iRow];
1629                setRowStatus(iRow, atLowerBound);
1630              } else {
1631                // say free
1632                setRowStatus(iRow, isFree);
1633                rowActivityWork_[iRow] = 0.0;
1634              }
1635            }
1636            break;
1637          case ClpSimplex::isFixed:
1638          case atLowerBound:
1639            rowActivityWork_[iRow] = rowLowerWork_[iRow];
1640            if (rowActivityWork_[iRow] < -largeValue_) {
1641              if (rowUpperWork_[iRow] < largeValue_) {
1642                rowActivityWork_[iRow] = rowUpperWork_[iRow];
1643                setRowStatus(iRow, atUpperBound);
1644              } else {
1645                // say free
1646                setRowStatus(iRow, isFree);
1647                rowActivityWork_[iRow] = 0.0;
1648              }
1649            }
1650            break;
1651          case isFree:
1652            break;
1653            // not really free - fall through to superbasic
1654          case superBasic:
1655            if (rowUpperWork_[iRow] > largeValue_) {
1656              if (rowLowerWork_[iRow] > -largeValue_) {
1657                rowActivityWork_[iRow] = rowLowerWork_[iRow];
1658                setRowStatus(iRow, atLowerBound);
1659              } else {
1660                // say free
1661                setRowStatus(iRow, isFree);
1662                rowActivityWork_[iRow] = 0.0;
1663              }
1664            } else {
1665              if (rowLowerWork_[iRow] > -largeValue_) {
1666                // set to nearest
1667                if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
1668                  < fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])) {
1669                  rowActivityWork_[iRow] = rowLowerWork_[iRow];
1670                  setRowStatus(iRow, atLowerBound);
1671                } else {
1672                  rowActivityWork_[iRow] = rowUpperWork_[iRow];
1673                  setRowStatus(iRow, atUpperBound);
1674                }
1675              } else {
1676                rowActivityWork_[iRow] = rowUpperWork_[iRow];
1677                setRowStatus(iRow, atUpperBound);
1678              }
1679            }
1680            break;
1681          }
1682        }
1683        totalSlacks = numberBasic;
1684
1685        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1686          switch (getColumnStatus(iColumn)) {
1687
1688          case basic:
1689            if (numberBasic == maximumBasic_) {
1690              // take out of basis
1691              if (columnLowerWork_[iColumn] > -largeValue_) {
1692                if (columnActivityWork_[iColumn] - columnLowerWork_[iColumn] < columnUpperWork_[iColumn] - columnActivityWork_[iColumn]) {
1693                  columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1694                  setColumnStatus(iColumn, atLowerBound);
1695                } else {
1696                  columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1697                  setColumnStatus(iColumn, atUpperBound);
1698                }
1699              } else if (columnUpperWork_[iColumn] < largeValue_) {
1700                columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1701                setColumnStatus(iColumn, atUpperBound);
1702              } else {
1703                columnActivityWork_[iColumn] = 0.0;
1704                setColumnStatus(iColumn, isFree);
1705              }
1706            } else {
1707              numberBasic++;
1708            }
1709            break;
1710          case atUpperBound:
1711            columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1712            if (columnActivityWork_[iColumn] > largeValue_) {
1713              if (columnLowerWork_[iColumn] < -largeValue_) {
1714                columnActivityWork_[iColumn] = 0.0;
1715                setColumnStatus(iColumn, isFree);
1716              } else {
1717                columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1718                setColumnStatus(iColumn, atLowerBound);
1719              }
1720            }
1721            break;
1722          case isFixed:
1723          case atLowerBound:
1724            columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1725            if (columnActivityWork_[iColumn] < -largeValue_) {
1726              if (columnUpperWork_[iColumn] > largeValue_) {
1727                columnActivityWork_[iColumn] = 0.0;
1728                setColumnStatus(iColumn, isFree);
1729              } else {
1730                columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1731                setColumnStatus(iColumn, atUpperBound);
1732              }
1733            }
1734            break;
1735          case isFree:
1736            break;
1737            // not really free - fall through to superbasic
1738          case superBasic:
1739            if (columnUpperWork_[iColumn] > largeValue_) {
1740              if (columnLowerWork_[iColumn] > -largeValue_) {
1741                columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1742                setColumnStatus(iColumn, atLowerBound);
1743              } else {
1744                // say free
1745                setColumnStatus(iColumn, isFree);
1746                columnActivityWork_[iColumn] = 0.0;
1747              }
1748            } else {
1749              if (columnLowerWork_[iColumn] > -largeValue_) {
1750                // set to nearest
1751                if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1752                  < fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])) {
1753                  columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1754                  setColumnStatus(iColumn, atLowerBound);
1755                } else {
1756                  columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1757                  setColumnStatus(iColumn, atUpperBound);
1758                }
1759              } else {
1760                columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1761                setColumnStatus(iColumn, atUpperBound);
1762              }
1763            }
1764            break;
1765          }
1766        }
1767#ifdef CLP_INVESTIGATE3
1768        if (saveSol) {
1769          int numberChanged = 0;
1770          double largestChanged = 0.0;
1771          for (int i = 0; i < numberTotal; i++) {
1772            double difference = fabs(solution_[i] - saveSol[i]);
1773            if (difference > 1.0e-7) {
1774              numberChanged++;
1775              if (difference > largestChanged)
1776                largestChanged = difference;
1777            }
1778          }
1779          if (numberChanged)
1780            printf("%d changed, largest %g\n", numberChanged, largestChanged);
1781          delete[] saveSol;
1782        }
1783#endif
1784#if 0
1785                    if (numberBasic < numberRows_) {
1786                         // add some slacks in case odd warmstart
1787#ifdef CLP_INVESTIGATE
1788                         printf("BAD %d basic, %d rows %d slacks\n",
1789                                numberBasic, numberRows_, totalSlacks);
1790#endif
1791                         int iRow = numberRows_ - 1;
1792                         while (numberBasic < numberRows_) {
1793                              if (getRowStatus(iRow) != basic) {
1794                                   setRowStatus(iRow, basic);
1795                                   numberBasic++;
1796                                   totalSlacks++;
1797                                   iRow--;
1798                              } else {
1799                                   break;
1800                              }
1801                         }
1802                    }
1803#endif
1804      } else {
1805        // all slack basis
1806        int numberBasic = 0;
1807        if (!status_) {
1808          createStatus();
1809        }
1810        for (iRow = 0; iRow < numberRows_; iRow++) {
1811          double lower = rowLowerWork_[iRow];
1812          double upper = rowUpperWork_[iRow];
1813          if (lower > -largeValue_ || upper < largeValue_) {
1814            if (fabs(lower) <= fabs(upper)) {
1815              rowActivityWork_[iRow] = lower;
1816            } else {
1817              rowActivityWork_[iRow] = upper;
1818            }
1819          } else {
1820            rowActivityWork_[iRow] = 0.0;
1821          }
1822          setRowStatus(iRow, basic);
1823          numberBasic++;
1824        }
1825        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1826          double lower = columnLowerWork_[iColumn];
1827          double upper = columnUpperWork_[iColumn];
1828          double big_bound = largeValue_;
1829          double value = columnActivityWork_[iColumn];
1830          if (lower > -big_bound || upper < big_bound) {
1831            if ((getColumnStatus(iColumn) == atLowerBound && columnActivityWork_[iColumn] == lower) || (getColumnStatus(iColumn) == atUpperBound && columnActivityWork_[iColumn] == upper)) {
1832              // status looks plausible
1833            } else {
1834              // set to sensible
1835              if (getColumnStatus(iColumn) == atUpperBound
1836                && upper < 1.0e20) {
1837                columnActivityWork_[iColumn] = upper;
1838              } else if (getColumnStatus(iColumn) == atLowerBound
1839                && lower > -1.0e20) {
1840                columnActivityWork_[iColumn] = lower;
1841              } else {
1842                if (fabs(lower) <= fabs(upper)) {
1843                  setColumnStatus(iColumn, atLowerBound);
1844                  columnActivityWork_[iColumn] = lower;
1845                } else {
1846                  setColumnStatus(iColumn, atUpperBound);
1847                  columnActivityWork_[iColumn] = upper;
1848                }
1849              }
1850            }
1851          } else {
1852            setColumnStatus(iColumn, isFree);
1853            columnActivityWork_[iColumn] = 0.0;
1854          }
1855        }
1856      }
1857    } else {
1858      // values pass has less coding
1859      // make row activities correct and clean basis a bit
1860      cleanStatus();
1861      if (status_) {
1862        int numberBasic = 0;
1863        for (iRow = 0; iRow < numberRows_; iRow++) {
1864          if (getRowStatus(iRow) == basic)
1865            numberBasic++;
1866        }
1867        totalSlacks = numberBasic;
1868#if 0
1869                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1870                         if (getColumnStatus(iColumn) == basic)
1871                              numberBasic++;
1872                    }
1873#endif
1874      } else {
1875        // all slack basis
1876        int numberBasic = 0;
1877        if (!status_) {
1878          createStatus();
1879        }
1880        for (iRow = 0; iRow < numberRows_; iRow++) {
1881          setRowStatus(iRow, basic);
1882          numberBasic++;
1883        }
1884        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1885          setColumnStatus(iColumn, superBasic);
1886          // but put to bound if close
1887          if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1888            <= primalTolerance_) {
1889            columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1890            setColumnStatus(iColumn, atLowerBound);
1891          } else if (fabs(columnActivityWork_[iColumn]
1892                       - columnUpperWork_[iColumn])
1893            <= primalTolerance_) {
1894            columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1895            setColumnStatus(iColumn, atUpperBound);
1896          }
1897        }
1898      }
1899    }
1900    numberRefinements_ = 1;
1901    // set fixed if they are
1902    for (iRow = 0; iRow < numberRows_; iRow++) {
1903      if (getRowStatus(iRow) != basic) {
1904        if (rowLowerWork_[iRow] == rowUpperWork_[iRow]) {
1905          rowActivityWork_[iRow] = rowLowerWork_[iRow];
1906          setRowStatus(iRow, isFixed);
1907        }
1908      }
1909    }
1910    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1911      if (getColumnStatus(iColumn) != basic) {
1912        if (columnLowerWork_[iColumn] == columnUpperWork_[iColumn]) {
1913          columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1914          setColumnStatus(iColumn, isFixed);
1915        }
1916      }
1917    }
1918  }
1919  //for (iRow=0;iRow<numberRows_+numberColumns_;iRow++) {
1920  //if (fabs(solution_[iRow])>1.0e10) {
1921  //  printf("large %g at %d - status %d\n",
1922  //         solution_[iRow],iRow,status_[iRow]);
1923  //}
1924  //}
1925#if 0 //ndef _MSC_VER                                                              \
1926  // The local static var k is a problem when trying to build a DLL. Since this is \
1927  // just for debugging (likely done on *nix), just hide it from Windows           \
1928  // -- lh, 101016 --
1929     if (0)  {
1930          static int k = 0;
1931          printf("start basis\n");
1932          int i;
1933          for (i = 0; i < numberRows_; i++)
1934               printf ("xx %d %d\n", i, pivotVariable_[i]);
1935          for (i = 0; i < numberRows_ + numberColumns_; i++)
1936               if (getColumnStatus(i) == basic)
1937                    printf ("yy %d basic\n", i);
1938          if (k > 20)
1939               exit(0);
1940          k++;
1941     }
1942#endif
1943#if 0 //ndef NDEBUG
1944  // Make sure everything is clean
1945  double sumOutside=0.0;
1946  int numberOutside=0;
1947  //double sumOutsideLarge=0.0;
1948  int numberOutsideLarge=0;
1949  double sumInside=0.0;
1950  int numberInside=0;
1951  //double sumInsideLarge=0.0;
1952  int numberInsideLarge=0;
1953  int numberTotal=numberRows_+numberColumns_;
1954  for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
1955    if(getStatus(iSequence) == isFixed) {
1956      // double check fixed
1957      assert (upper_[iSequence] == lower_[iSequence]);
1958      assert (fabs(solution_[iSequence]-lower_[iSequence])<primalTolerance_);
1959    } else if (getStatus(iSequence) == isFree) {
1960      assert (upper_[iSequence] == COIN_DBL_MAX && lower_[iSequence]==-COIN_DBL_MAX);
1961    } else if (getStatus(iSequence) == atLowerBound) {
1962      assert (fabs(solution_[iSequence]-lower_[iSequence])<1000.0*primalTolerance_);
1963      if (solution_[iSequence]<lower_[iSequence]) {
1964        numberOutside++;
1965        sumOutside-=solution_[iSequence]-lower_[iSequence];
1966        if (solution_[iSequence]<lower_[iSequence]-primalTolerance_)
1967          numberOutsideLarge++;
1968      } else if (solution_[iSequence]>lower_[iSequence]) {
1969        numberInside++;
1970        sumInside+=solution_[iSequence]-lower_[iSequence];
1971        if (solution_[iSequence]>lower_[iSequence]+primalTolerance_)
1972          numberInsideLarge++;
1973      }
1974    } else if (getStatus(iSequence) == atUpperBound) {
1975      assert (fabs(solution_[iSequence]-upper_[iSequence])<1000.0*primalTolerance_);
1976      if (solution_[iSequence]>upper_[iSequence]) {
1977        numberOutside++;
1978        sumOutside+=solution_[iSequence]-upper_[iSequence];
1979        if (solution_[iSequence]>upper_[iSequence]+primalTolerance_)
1980          numberOutsideLarge++;
1981      } else if (solution_[iSequence]<upper_[iSequence]) {
1982        numberInside++;
1983        sumInside-=solution_[iSequence]-upper_[iSequence];
1984        if (solution_[iSequence]<upper_[iSequence]-primalTolerance_)
1985          numberInsideLarge++;
1986      }
1987    } else if (getStatus(iSequence) == superBasic) {
1988      //assert (!valuesPass);
1989    }
1990  }
1991  if (numberInside+numberOutside)
1992    printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
1993           numberOutside,sumOutside,numberOutsideLarge,
1994           numberInside,sumInside,numberInsideLarge);
1995#endif
1996  int status = factorization_->factorize(this, solveType, valuesPass);
1997  if (status) {
1998    handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
1999      << status
2000      << CoinMessageEol;
2001#ifdef CLP_USEFUL_PRINTOUT
2002    printf("Basis singular - pivot tolerance %g\n",
2003      factorization_->pivotTolerance());
2004#endif
2005    return -1;
2006  } else if (!solveType) {
2007    // Initial basis - return number of singularities
2008    int numberSlacks = 0;
2009    for (iRow = 0; iRow < numberRows_; iRow++) {
2010      if (getRowStatus(iRow) == basic)
2011        numberSlacks++;
2012    }
2013    status = CoinMax(numberSlacks - totalSlacks, 0);
2014    // special case if all slack
2015    if (numberSlacks == numberRows_) {
2016      status = numberRows_ + 1;
2017    }
2018  }
2019
2020  // sparse methods
2021  //if (factorization_->sparseThreshold()) {
2022  // get default value
2023  factorization_->sparseThreshold(0);
2024  if (!(moreSpecialOptions_ & 1024))
2025    factorization_->goSparse();
2026  //}
2027
2028  return status;
2029}
2030/*
2031   This does basis housekeeping and does values for in/out variables.
2032   Can also decide to re-factorize
2033*/
2034int ClpSimplex::housekeeping(double objectiveChange)
2035{
2036  // save value of incoming and outgoing
2037  double oldIn = solution_[sequenceIn_];
2038  double oldOut = solution_[sequenceOut_];
2039  numberIterations_++;
2040  changeMade_++; // something has happened
2041  // incoming variable
2042  if (handler_->logLevel() > 7) {
2043    //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
2044    handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
2045      << directionOut_
2046      << directionIn_ << theta_
2047      << dualOut_ << dualIn_ << alpha_
2048      << CoinMessageEol;
2049    if (getStatus(sequenceIn_) == isFree) {
2050      handler_->message(CLP_SIMPLEX_FREEIN, messages_)
2051        << sequenceIn_
2052        << CoinMessageEol;
2053    }
2054  }
2055#if 0
2056     printf("h1 %d %d %g %g %g %g",
2057            directionOut_
2058            , directionIn_, theta_
2059            , dualOut_, dualIn_, alpha_);
2060#endif
2061  // change of incoming
2062  char rowcol[] = { 'R', 'C' };
2063  if (pivotRow_ >= 0)
2064    pivotVariable_[pivotRow_] = sequenceIn();
2065  if (upper_[sequenceIn_] > 1.0e20 && lower_[sequenceIn_] < -1.0e20)
2066    progressFlag_ |= 2; // making real progress
2067  solution_[sequenceIn_] = valueIn_;
2068  if (upper_[sequenceOut_] - lower_[sequenceOut_] < 1.0e-12)
2069    progressFlag_ |= 1; // making real progress
2070  if (sequenceIn_ != sequenceOut_) {
2071    if (alphaAccuracy_ > 0.0) {
2072      double value = fabs(alpha_);
2073      if (value > 1.0)
2074        alphaAccuracy_ *= value;
2075      else
2076        alphaAccuracy_ /= value;
2077    }
2078    //assert( getStatus(sequenceOut_)== basic);
2079    setStatus(sequenceIn_, basic);
2080    if (upper_[sequenceOut_] - lower_[sequenceOut_] > 0) {
2081      // As Nonlinear costs may have moved bounds (to more feasible)
2082      // Redo using value
2083      if (fabs(valueOut_ - lower_[sequenceOut_]) < fabs(valueOut_ - upper_[sequenceOut_])) {
2084        // going to lower
2085        setStatus(sequenceOut_, atLowerBound);
2086        oldOut = lower_[sequenceOut_];
2087      } else {
2088        // going to upper
2089        setStatus(sequenceOut_, atUpperBound);
2090        oldOut = upper_[sequenceOut_];
2091      }
2092    } else {
2093      // fixed
2094      setStatus(sequenceOut_, isFixed);
2095    }
2096    solution_[sequenceOut_] = valueOut_;
2097  } else {
2098    //if (objective_->type()<2)
2099    //assert (fabs(theta_)>1.0e-13);
2100    // flip from bound to bound
2101    // As Nonlinear costs may have moved bounds (to more feasible)
2102    // Redo using value
2103    if (fabs(valueIn_ - lower_[sequenceIn_]) < fabs(valueIn_ - upper_[sequenceIn_])) {
2104      // as if from upper bound
2105      setStatus(sequenceIn_, atLowerBound);
2106    } else {
2107      // as if from lower bound
2108      setStatus(sequenceIn_, atUpperBound);
2109    }
2110  }
2111
2112  // Update hidden stuff e.g. effective RHS and gub
2113  int invertNow = matrix_->updatePivot(this, oldIn, oldOut);
2114  objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
2115  if (handler_->logLevel() > 7) {
2116    //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
2117    handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
2118      << numberIterations_ << objectiveValue()
2119      << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
2120      << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
2121    handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
2122    handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
2123    handler_->message() << CoinMessageEol;
2124  }
2125#if 0
2126     if (numberIterations_ > 10000)
2127          printf(" it %d %g %c%d %c%d\n"
2128                 , numberIterations_, objectiveValue()
2129                 , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
2130                 , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
2131#endif
2132  if (trustedUserPointer_ && trustedUserPointer_->typeStruct == 1) {
2133    if (algorithm_ > 0 && integerType_ && !nonLinearCost_->numberInfeasibilities()) {
2134      if (fabs(theta_) > 1.0e-6 || !numberIterations_) {
2135        // For saving solutions
2136        typedef struct {
2137          int numberSolutions;
2138          int maximumSolutions;
2139          int numberColumns;
2140          double **solution;
2141          int *numberUnsatisfied;
2142        } clpSolution;
2143        clpSolution *solution = reinterpret_cast< clpSolution * >(trustedUserPointer_->data);
2144        if (solution->numberSolutions == solution->maximumSolutions) {
2145          int n = solution->maximumSolutions;
2146          int n2 = (n * 3) / 2 + 10;
2147          solution->maximumSolutions = n2;
2148          double **temp = new double *[n2];
2149          for (int i = 0; i < n; i++)
2150            temp[i] = solution->solution[i];
2151          delete[] solution->solution;
2152          solution->solution = temp;
2153          int *tempN = new int[n2];
2154          for (int i = 0; i < n; i++)
2155            tempN[i] = solution->numberUnsatisfied[i];
2156          delete[] solution->numberUnsatisfied;
2157          solution->numberUnsatisfied = tempN;
2158        }
2159        assert(numberColumns_ == solution->numberColumns);
2160        double *sol = new double[numberColumns_];
2161        solution->solution[solution->numberSolutions] = sol;
2162        int numberFixed = 0;
2163        int numberUnsat = 0;
2164        int numberSat = 0;
2165        double sumUnsat = 0.0;
2166        double tolerance = 10.0 * primalTolerance_;
2167        double mostAway = 0.0;
2168        for (int i = 0; i < numberColumns_; i++) {
2169          // Save anyway
2170          sol[i] = columnScale_ ? solution_[i] * columnScale_[i] : solution_[i];
2171          // rest is optional
2172          if (upper_[i] > lower_[i]) {
2173            double value = solution_[i];
2174            if (value > lower_[i] + tolerance && value < upper_[i] - tolerance && integerType_[i]) {
2175              // may have to modify value if scaled
2176              if (columnScale_)
2177                value *= columnScale_[i];
2178              double closest = floor(value + 0.5);
2179              // problem may be perturbed so relax test
2180              if (fabs(value - closest) > 1.0e-4) {
2181                numberUnsat++;
2182                sumUnsat += fabs(value - closest);
2183                if (mostAway < fabs(value - closest)) {
2184                  mostAway = fabs(value - closest);
2185                }
2186              } else {
2187                numberSat++;
2188              }
2189            } else {
2190              numberSat++;
2191            }
2192          } else {
2193            numberFixed++;
2194          }
2195        }
2196        solution->numberUnsatisfied[solution->numberSolutions++] = numberUnsat;
2197        COIN_DETAIL_PRINT(printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n",
2198          numberIterations_, numberUnsat, sumUnsat, mostAway, numberFixed, numberSat));
2199      }
2200    }
2201  }
2202  if (hitMaximumIterations())
2203    return 2;
2204#if 1
2205  //if (numberIterations_>14000)
2206  //handler_->setLogLevel(63);
2207  //if (numberIterations_>24000)
2208  //exit(77);
2209  // check for small cycles
2210  int in = sequenceIn_;
2211  int out = sequenceOut_;
2212  matrix_->correctSequence(this, in, out);
2213  int cycle = progress_.cycle(in, out,
2214    directionIn_, directionOut_);
2215  if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
2216    //if (cycle>0) {
2217    if (handler_->logLevel() >= 63)
2218      printf("Cycle of %d\n", cycle);
2219    // reset
2220    progress_.startCheck();
2221    double random = randomNumberGenerator_.randomDouble();
2222    int extra = static_cast< int >(9.999 * random);
2223    int off[] = { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 };
2224    if (factorization_->pivots() > cycle) {
2225      forceFactorization_ = CoinMax(1, cycle - off[extra]);
2226    } else {
2227      /* need to reject something
2228               should be better if don't reject incoming
2229               as it is in basis */
2230      int iSequence;
2231      //if (algorithm_ > 0)
2232      //   iSequence = sequenceIn_;
2233      //else
2234      iSequence = sequenceOut_;
2235      char x = isColumn(iSequence) ? 'C' : 'R';
2236      if (handler_->logLevel() >= 63)
2237        handler_->message(CLP_SIMPLEX_FLAG, messages_)
2238          << x << sequenceWithin(iSequence)
2239          << CoinMessageEol;
2240      setFlagged(iSequence);
2241      //printf("flagging %d\n",iSequence);
2242    }
2243    return 1;
2244  }
2245#endif
2246  // only time to re-factorize if one before real time
2247  // this is so user won't be surprised that maximumPivots has exact meaning
2248  int numberPivots = factorization_->pivots();
2249  int maximumPivots = factorization_->maximumPivots();
2250  int numberDense = factorization_->numberDense();
2251  bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 > 2 * maximumIterations());
2252  if (numberPivots == maximumPivots || maximumPivots < 2) {
2253    // If dense then increase
2254    if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots
2255      && false) {
2256      factorization_->maximumPivots(numberDense);
2257      dualRowPivot_->maximumPivotsChanged();
2258      primalColumnPivot_->maximumPivotsChanged();
2259      // and redo arrays
2260      for (int iRow = 0; iRow < 4; iRow++) {
2261        int length = rowArray_[iRow]->capacity() + numberDense - maximumPivots;
2262        rowArray_[iRow]->reserve(length);
2263      }
2264    }
2265#if CLP_FACTORIZATION_NEW_TIMING > 1
2266    factorization_->statsRefactor('M');
2267#endif
2268    return 1;
2269  } else if ((factorization_->timeToRefactorize() && !dontInvert)
2270    || invertNow) {
2271    //printf("ret after %d pivots\n",factorization_->pivots());
2272#if CLP_FACTORIZATION_NEW_TIMING > 1
2273    factorization_->statsRefactor('T');
2274#endif
2275    return 1;
2276  } else if (forceFactorization_ > 0 && factorization_->pivots() == forceFactorization_) {
2277    // relax
2278    forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
2279    if (forceFactorization_ > factorization_->maximumPivots())
2280      forceFactorization_ = -1; //off
2281#if CLP_FACTORIZATION_NEW_TIMING > 1
2282    factorization_->statsRefactor('F');
2283#endif
2284    return 1;
2285  } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type() < 15) {
2286    double random = randomNumberGenerator_.randomDouble();
2287    while (random < 0.45)
2288      random *= 2.0;
2289    int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
2290    if (factorization_->pivots() >= random * maxNumber) {
2291      return 1;
2292    } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && numberIterations_ < 1001000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
2293      return 1;
2294    } else {
2295      // carry on iterating
2296      return 0;
2297    }
2298  } else {
2299    // carry on iterating
2300    return 0;
2301  }
2302}
2303// Copy constructor.
2304ClpSimplex::ClpSimplex(const ClpSimplex &rhs, int scalingMode)
2305  : ClpModel(rhs, scalingMode)
2306  , bestPossibleImprovement_(0.0)
2307  , zeroTolerance_(1.0e-13)
2308  , columnPrimalSequence_(-2)
2309  , rowPrimalSequence_(-2)
2310  , bestObjectiveValue_(rhs.bestObjectiveValue_)
2311  , moreSpecialOptions_(2)
2312  , baseIteration_(0)
2313  , vectorMode_(0)
2314  , primalToleranceToGetOptimal_(-1.0)
2315  , largeValue_(1.0e15)
2316  , largestPrimalError_(0.0)
2317  , largestDualError_(0.0)
2318  , alphaAccuracy_(-1.0)
2319  , dualBound_(1.0e10)
2320  , alpha_(0.0)
2321  , theta_(0.0)
2322  , lowerIn_(0.0)
2323  , valueIn_(0.0)
2324  , upperIn_(-COIN_DBL_MAX)
2325  , dualIn_(0.0)
2326  , lowerOut_(-1)
2327  , valueOut_(-1)
2328  , upperOut_(-1)
2329  , dualOut_(-1)
2330  , dualTolerance_(1.0e-7)
2331  , primalTolerance_(1.0e-7)
2332  , sumDualInfeasibilities_(0.0)
2333  , sumPrimalInfeasibilities_(0.0)
2334  , infeasibilityCost_(1.0e10)
2335  , sumOfRelaxedDualInfeasibilities_(0.0)
2336  , sumOfRelaxedPrimalInfeasibilities_(0.0)
2337  , acceptablePivot_(1.0e-8)
2338  , lower_(NULL)
2339  , rowLowerWork_(NULL)
2340  , columnLowerWork_(NULL)
2341  , upper_(NULL)
2342  , rowUpperWork_(NULL)
2343  , columnUpperWork_(NULL)
2344  , cost_(NULL)
2345  , rowObjectiveWork_(NULL)
2346  , objectiveWork_(NULL)
2347  , sequenceIn_(-1)
2348  , directionIn_(-1)
2349  , sequenceOut_(-1)
2350  , directionOut_(-1)
2351  , pivotRow_(-1)
2352  , lastGoodIteration_(-100)
2353  , dj_(NULL)
2354  , rowReducedCost_(NULL)
2355  , reducedCostWork_(NULL)
2356  , solution_(NULL)
2357  , rowActivityWork_(NULL)
2358  , columnActivityWork_(NULL)
2359  , numberDualInfeasibilities_(0)
2360  , numberDualInfeasibilitiesWithoutFree_(0)
2361  , numberPrimalInfeasibilities_(100)
2362  , numberRefinements_(0)
2363  , pivotVariable_(NULL)
2364  , factorization_(NULL)
2365  , savedSolution_(NULL)
2366  , numberTimesOptimal_(0)
2367  , disasterArea_(NULL)
2368  , changeMade_(1)
2369  , algorithm_(0)
2370  , forceFactorization_(-1)
2371  , perturbation_(100)
2372  , nonLinearCost_(NULL)
2373  , lastBadIteration_(-999999)
2374  , lastFlaggedIteration_(-999999)
2375  , numberFake_(0)
2376  , numberChanged_(0)
2377  , progressFlag_(0)
2378  , firstFree_(-1)
2379  , numberExtraRows_(0)
2380  , maximumBasic_(0)
2381  , dontFactorizePivots_(0)
2382  , incomingInfeasibility_(1.0)
2383  , allowedInfeasibility_(10.0)
2384  , automaticScale_(0)
2385  , maximumPerturbationSize_(0)
2386  , perturbationArray_(NULL)
2387  , baseModel_(NULL)
2388#ifdef ABC_INHERIT
2389  , abcSimplex_(NULL)
2390  , abcState_(0)
2391#endif
2392{
2393  int i;
2394  for (i = 0; i < 6; i++) {
2395    rowArray_[i] = NULL;
2396    columnArray_[i] = NULL;
2397  }
2398  for (i = 0; i < 4; i++) {
2399    spareIntArray_[i] = 0;
2400    spareDoubleArray_[i] = 0.0;
2401  }
2402  saveStatus_ = NULL;
2403  factorization_ = NULL;
2404  dualRowPivot_ = NULL;
2405  primalColumnPivot_ = NULL;
2406  gutsOfDelete(0);
2407  delete nonLinearCost_;
2408  nonLinearCost_ = NULL;
2409  gutsOfCopy(rhs);
2410  solveType_ = 1; // say simplex based life form
2411}
2412// Copy constructor from model
2413ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode)
2414  : ClpModel(rhs, scalingMode)
2415  , bestPossibleImprovement_(0.0)
2416  , zeroTolerance_(1.0e-13)
2417  , columnPrimalSequence_(-2)
2418  , rowPrimalSequence_(-2)
2419  , bestObjectiveValue_(-COIN_DBL_MAX)
2420  , moreSpecialOptions_(2)
2421  , baseIteration_(0)
2422  , vectorMode_(0)
2423  , primalToleranceToGetOptimal_(-1.0)
2424  , largeValue_(1.0e15)
2425  , largestPrimalError_(0.0)
2426  , largestDualError_(0.0)
2427  , alphaAccuracy_(-1.0)
2428  , dualBound_(1.0e10)
2429  , alpha_(0.0)
2430  , theta_(0.0)
2431  , lowerIn_(0.0)
2432  , valueIn_(0.0)
2433  , upperIn_(-COIN_DBL_MAX)
2434  , dualIn_(0.0)
2435  , lowerOut_(-1)
2436  , valueOut_(-1)
2437  , upperOut_(-1)
2438  , dualOut_(-1)
2439  , dualTolerance_(1.0e-7)
2440  , primalTolerance_(1.0e-7)
2441  , sumDualInfeasibilities_(0.0)
2442  , sumPrimalInfeasibilities_(0.0)
2443  , infeasibilityCost_(1.0e10)
2444  , sumOfRelaxedDualInfeasibilities_(0.0)
2445  , sumOfRelaxedPrimalInfeasibilities_(0.0)
2446  , acceptablePivot_(1.0e-8)
2447  , lower_(NULL)
2448  , rowLowerWork_(NULL)
2449  , columnLowerWork_(NULL)
2450  , upper_(NULL)
2451  , rowUpperWork_(NULL)
2452  , columnUpperWork_(NULL)
2453  , cost_(NULL)
2454  , rowObjectiveWork_(NULL)
2455  , objectiveWork_(NULL)
2456  , sequenceIn_(-1)
2457  , directionIn_(-1)
2458  , sequenceOut_(-1)
2459  , directionOut_(-1)
2460  , pivotRow_(-1)
2461  , lastGoodIteration_(-100)
2462  , dj_(NULL)
2463  , rowReducedCost_(NULL)
2464  , reducedCostWork_(NULL)
2465  , solution_(NULL)
2466  , rowActivityWork_(NULL)
2467  , columnActivityWork_(NULL)
2468  , numberDualInfeasibilities_(0)
2469  , numberDualInfeasibilitiesWithoutFree_(0)
2470  , numberPrimalInfeasibilities_(100)
2471  , numberRefinements_(0)
2472  , pivotVariable_(NULL)
2473  , factorization_(NULL)
2474  , savedSolution_(NULL)
2475  , numberTimesOptimal_(0)
2476  , disasterArea_(NULL)
2477  , changeMade_(1)
2478  , algorithm_(0)
2479  , forceFactorization_(-1)
2480  , perturbation_(100)
2481  , nonLinearCost_(NULL)
2482  , lastBadIteration_(-999999)
2483  , lastFlaggedIteration_(-999999)
2484  , numberFake_(0)
2485  , numberChanged_(0)
2486  , progressFlag_(0)
2487  , firstFree_(-1)
2488  , numberExtraRows_(0)
2489  , maximumBasic_(0)
2490  , dontFactorizePivots_(0)
2491  , incomingInfeasibility_(1.0)
2492  , allowedInfeasibility_(10.0)
2493  , automaticScale_(0)
2494  , maximumPerturbationSize_(0)
2495  , perturbationArray_(NULL)
2496  , baseModel_(NULL)
2497#ifdef ABC_INHERIT
2498  , abcSimplex_(NULL)
2499  , abcState_(0)
2500#endif
2501{
2502  int i;
2503  for (i = 0; i < 6; i++) {
2504    rowArray_[i] = NULL;
2505    columnArray_[i] = NULL;
2506  }
2507  for (i = 0; i < 4; i++) {
2508    spareIntArray_[i] = 0;
2509    spareDoubleArray_[i] = 0.0;
2510  }
2511  saveStatus_ = NULL;
2512  // get an empty factorization so we can set tolerances etc
2513  getEmptyFactorization();
2514  // say Steepest pricing
2515  dualRowPivot_ = new ClpDualRowSteepest();
2516  // say Steepest pricing
2517  primalColumnPivot_ = new ClpPrimalColumnSteepest();
2518  solveType_ = 1; // say simplex based life form
2519}
2520// Assignment operator. This copies the data
2521ClpSimplex &
2522ClpSimplex::operator=(const ClpSimplex &rhs)
2523{
2524  if (this != &rhs) {
2525    gutsOfDelete(0);
2526    delete nonLinearCost_;
2527    nonLinearCost_ = NULL;
2528    ClpModel::operator=(rhs);
2529    gutsOfCopy(rhs);
2530  }
2531  return *this;
2532}
2533void ClpSimplex::gutsOfCopy(const ClpSimplex &rhs)
2534{
2535  assert(numberRows_ == rhs.numberRows_);
2536  assert(numberColumns_ == rhs.numberColumns_);
2537  numberExtraRows_ = rhs.numberExtraRows_;
2538  maximumBasic_ = rhs.maximumBasic_;
2539  dontFactorizePivots_ = rhs.dontFactorizePivots_;
2540  int numberRows2 = numberRows_ + numberExtraRows_;
2541  moreSpecialOptions_ = rhs.moreSpecialOptions_;
2542  if ((whatsChanged_ & 1) != 0) {
2543    int numberTotal = numberColumns_ + numberRows2;
2544    if ((specialOptions_ & 65536) != 0 && maximumRows_ >= 0) {
2545      assert(maximumInternalRows_ >= numberRows2);
2546      assert(maximumInternalColumns_ >= numberColumns_);
2547      numberTotal = 2 * (maximumInternalColumns_ + maximumInternalRows_);
2548    }
2549    lower_ = ClpCopyOfArray(rhs.lower_, numberTotal);
2550    rowLowerWork_ = lower_ + numberColumns_;
2551    columnLowerWork_ = lower_;
2552    upper_ = ClpCopyOfArray(rhs.upper_, numberTotal);
2553    rowUpperWork_ = upper_ + numberColumns_;
2554    columnUpperWork_ = upper_;
2555    cost_ = ClpCopyOfArray(rhs.cost_, numberTotal);
2556    objectiveWork_ = cost_;
2557    rowObjectiveWork_ = cost_ + numberColumns_;
2558    dj_ = ClpCopyOfArray(rhs.dj_, numberTotal);
2559    if (dj_) {
2560      reducedCostWork_ = dj_;
2561      rowReducedCost_ = dj_ + numberColumns_;
2562    }
2563    solution_ = ClpCopyOfArray(rhs.solution_, numberTotal);
2564    if (solution_) {
2565      columnActivityWork_ = solution_;
2566      rowActivityWork_ = solution_ + numberColumns_;
2567    }
2568    if (rhs.pivotVariable_) {
2569      pivotVariable_ = new int[numberRows2];
2570      CoinMemcpyN(rhs.pivotVariable_, numberRows2, pivotVariable_);
2571    } else {
2572      pivotVariable_ = NULL;
2573    }
2574    savedSolution_ = ClpCopyOfArray(rhs.savedSolution_, numberTotal);
2575    int i;
2576    for (i = 0; i < 6; i++) {
2577      rowArray_[i] = NULL;
2578      if (rhs.rowArray_[i])
2579        rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
2580      columnArray_[i] = NULL;
2581      if (rhs.columnArray_[i])
2582        columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
2583    }
2584    if (rhs.saveStatus_) {
2585      saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberTotal);
2586    }
2587  } else {
2588    lower_ = NULL;
2589    rowLowerWork_ = NULL;
2590    columnLowerWork_ = NULL;
2591    upper_ = NULL;
2592    rowUpperWork_ = NULL;
2593    columnUpperWork_ = NULL;
2594    cost_ = NULL;
2595    objectiveWork_ = NULL;
2596    rowObjectiveWork_ = NULL;
2597    dj_ = NULL;
2598    reducedCostWork_ = NULL;
2599    rowReducedCost_ = NULL;
2600    solution_ = NULL;
2601    columnActivityWork_ = NULL;
2602    rowActivityWork_ = NULL;
2603    pivotVariable_ = NULL;
2604    savedSolution_ = NULL;
2605    int i;
2606    for (i = 0; i < 6; i++) {
2607      rowArray_[i] = NULL;
2608      columnArray_[i] = NULL;
2609    }
2610    saveStatus_ = NULL;
2611  }
2612  if (rhs.factorization_) {
2613    setFactorization(*rhs.factorization_);
2614  } else {
2615    delete factorization_;
2616    factorization_ = NULL;
2617  }
2618  bestPossibleImprovement_ = rhs.bestPossibleImprovement_;
2619  columnPrimalSequence_ = rhs.columnPrimalSequence_;
2620  zeroTolerance_ = rhs.zeroTolerance_;
2621  rowPrimalSequence_ = rhs.rowPrimalSequence_;
2622  bestObjectiveValue_ = rhs.bestObjectiveValue_;
2623  baseIteration_ = rhs.baseIteration_;
2624  vectorMode_ = rhs.vectorMode_;
2625  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
2626  largeValue_ = rhs.largeValue_;
2627  largestPrimalError_ = rhs.largestPrimalError_;
2628  largestDualError_ = rhs.largestDualError_;
2629  alphaAccuracy_ = rhs.alphaAccuracy_;
2630  dualBound_ = rhs.dualBound_;
2631  alpha_ = rhs.alpha_;
2632  theta_ = rhs.theta_;
2633  lowerIn_ = rhs.lowerIn_;
2634  valueIn_ = rhs.valueIn_;
2635  upperIn_ = rhs.upperIn_;
2636  dualIn_ = rhs.dualIn_;
2637  sequenceIn_ = rhs.sequenceIn_;
2638  directionIn_ = rhs.directionIn_;
2639  lowerOut_ = rhs.lowerOut_;
2640  valueOut_ = rhs.valueOut_;
2641  upperOut_ = rhs.upperOut_;
2642  dualOut_ = rhs.dualOut_;
2643  sequenceOut_ = rhs.sequenceOut_;
2644  directionOut_ = rhs.directionOut_;
2645  pivotRow_ = rhs.pivotRow_;
2646  lastGoodIteration_ = rhs.lastGoodIteration_;
2647  numberRefinements_ = rhs.numberRefinements_;
2648  dualTolerance_ = rhs.dualTolerance_;
2649  primalTolerance_ = rhs.primalTolerance_;
2650  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
2651  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
2652  numberDualInfeasibilitiesWithoutFree_ = rhs.numberDualInfeasibilitiesWithoutFree_;
2653  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
2654  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
2655  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
2656  dualRowPivot_->setModel(this);
2657  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
2658  primalColumnPivot_->setModel(this);
2659  numberTimesOptimal_ = rhs.numberTimesOptimal_;
2660  disasterArea_ = NULL;
2661  changeMade_ = rhs.changeMade_;
2662  algorithm_ = rhs.algorithm_;
2663  forceFactorization_ = rhs.forceFactorization_;
2664  perturbation_ = rhs.perturbation_;
2665  infeasibilityCost_ = rhs.infeasibilityCost_;
2666  lastBadIteration_ = rhs.lastBadIteration_;
2667  lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
2668  numberFake_ = rhs.numberFake_;
2669  numberChanged_ = rhs.numberChanged_;
2670  progressFlag_ = rhs.progressFlag_;
2671  firstFree_ = rhs.firstFree_;
2672  incomingInfeasibility_ = rhs.incomingInfeasibility_;
2673  allowedInfeasibility_ = rhs.allowedInfeasibility_;
2674  automaticScale_ = rhs.automaticScale_;
2675#ifdef ABC_INHERIT
2676  abcSimplex_ = NULL;
2677  abcState_ = rhs.abcState_;
2678#endif
2679  maximumPerturbationSize_ = rhs.maximumPerturbationSize_;
2680  if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) {
2681    perturbationArray_ = CoinCopyOfArray(rhs.perturbationArray_,
2682      maximumPerturbationSize_);
2683  } else {
2684    maximumPerturbationSize_ = 0;
2685    perturbationArray_ = NULL;
2686  }
2687  if (rhs.baseModel_) {
2688    baseModel_ = new ClpSimplex(*rhs.baseModel_);
2689  } else {
2690    baseModel_ = NULL;
2691  }
2692  progress_ = rhs.progress_;
2693  for (int i = 0; i < 4; i++) {
2694    spareIntArray_[i] = rhs.spareIntArray_[i];
2695    spareDoubleArray_[i] = rhs.spareDoubleArray_[i];
2696  }
2697  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
2698  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
2699  acceptablePivot_ = rhs.acceptablePivot_;
2700  if (rhs.nonLinearCost_ != NULL)
2701    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
2702  else
2703    nonLinearCost_ = NULL;
2704  solveType_ = rhs.solveType_;
2705  eventHandler_->setSimplex(this);
2706}
2707// type == 0 do everything, most + pivot data, 2 factorization data as well
2708void ClpSimplex::gutsOfDelete(int type)
2709{
2710  if (!type || (specialOptions_ & 65536) == 0) {
2711    maximumInternalColumns_ = -1;
2712    maximumInternalRows_ = -1;
2713    delete[] lower_;
2714    lower_ = NULL;
2715    rowLowerWork_ = NULL;
2716    columnLowerWork_ = NULL;
2717    delete[] upper_;
2718    upper_ = NULL;
2719    rowUpperWork_ = NULL;
2720    columnUpperWork_ = NULL;
2721    delete[] cost_;
2722    cost_ = NULL;
2723    objectiveWork_ = NULL;
2724    rowObjectiveWork_ = NULL;
2725    delete[] dj_;
2726    dj_ = NULL;
2727    reducedCostWork_ = NULL;
2728    rowReducedCost_ = NULL;
2729    delete[] solution_;
2730    solution_ = NULL;
2731    rowActivityWork_ = NULL;
2732    columnActivityWork_ = NULL;
2733    delete[] savedSolution_;
2734    savedSolution_ = NULL;
2735  }
2736  if ((specialOptions_ & 2) == 0) {
2737    delete nonLinearCost_;
2738    nonLinearCost_ = NULL;
2739  }
2740  int i;
2741  if ((specialOptions_ & 65536) == 0) {
2742    for (i = 0; i < 6; i++) {
2743      delete rowArray_[i];
2744      rowArray_[i] = NULL;
2745      delete columnArray_[i];
2746      columnArray_[i] = NULL;
2747    }
2748  }
2749  delete[] saveStatus_;
2750  saveStatus_ = NULL;
2751  if (type != 1) {
2752    delete rowCopy_;
2753    rowCopy_ = NULL;
2754  }
2755  if (!type) {
2756    // delete everything
2757    setEmptyFactorization();
2758    delete[] pivotVariable_;
2759    pivotVariable_ = NULL;
2760    delete dualRowPivot_;
2761    dualRowPivot_ = NULL;
2762    delete primalColumnPivot_;
2763    primalColumnPivot_ = NULL;
2764    delete baseModel_;
2765    baseModel_ = NULL;
2766    delete[] perturbationArray_;
2767    perturbationArray_ = NULL;
2768    maximumPerturbationSize_ = 0;
2769  } else {
2770    // delete any size information in methods
2771    if (type > 1) {
2772      //assert (factorization_);
2773      if (factorization_)
2774        factorization_->clearArrays();
2775      delete[] pivotVariable_;
2776      pivotVariable_ = NULL;
2777    }
2778    dualRowPivot_->clearArrays();
2779    primalColumnPivot_->clearArrays();
2780  }
2781}
2782// This sets largest infeasibility and most infeasible
2783void ClpSimplex::checkPrimalSolution(const double *rowActivities,
2784  const double *columnActivities)
2785{
2786  double *solution;
2787  int iRow, iColumn;
2788
2789  objectiveValue_ = 0.0;
2790  // now look at primal solution
2791  solution = rowActivityWork_;
2792  sumPrimalInfeasibilities_ = 0.0;
2793  numberPrimalInfeasibilities_ = 0;
2794  double primalTolerance = primalTolerance_;
2795  double relaxedTolerance = primalTolerance_;
2796  // we can't really trust infeasibilities if there is primal error
2797  double error = CoinMin(1.0e-2, largestPrimalError_);
2798  // allow tolerance at least slightly bigger than standard
2799  relaxedTolerance = relaxedTolerance + error;
2800  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2801  for (iRow = 0; iRow < numberRows_; iRow++) {
2802    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
2803    double infeasibility = 0.0;
2804    objectiveValue_ += solution[iRow] * rowObjectiveWork_[iRow];
2805    if (solution[iRow] > rowUpperWork_[iRow]) {
2806      infeasibility = solution[iRow] - rowUpperWork_[iRow];
2807    } else if (solution[iRow] < rowLowerWork_[iRow]) {
2808      infeasibility = rowLowerWork_[iRow] - solution[iRow];
2809    }
2810    if (infeasibility > primalTolerance) {
2811      sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2812      if (infeasibility > relaxedTolerance)
2813        sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2814      numberPrimalInfeasibilities_++;
2815    }
2816    infeasibility = fabs(rowActivities[iRow] - solution[iRow]);
2817  }
2818  // Check any infeasibilities from dynamic rows
2819  matrix_->primalExpanded(this, 2);
2820  solution = columnActivityWork_;
2821  if (!matrix_->rhsOffset(this)) {
2822    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2823      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2824      double infeasibility = 0.0;
2825      objectiveValue_ += objectiveWork_[iColumn] * solution[iColumn];
2826      if (solution[iColumn] > columnUpperWork_[iColumn]) {
2827        infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2828      } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2829        infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2830      }
2831      if (infeasibility > primalTolerance) {
2832        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2833        if (infeasibility > relaxedTolerance)
2834          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2835        numberPrimalInfeasibilities_++;
2836      }
2837      infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2838    }
2839  } else {
2840    // as we are using effective rhs we only check basics
2841    // But we do need to get objective
2842    objectiveValue_ += innerProduct(objectiveWork_, numberColumns_, solution);
2843    for (int j = 0; j < numberRows_; j++) {
2844      int iColumn = pivotVariable_[j];
2845      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2846      double infeasibility = 0.0;
2847      if (solution[iColumn] > columnUpperWork_[iColumn]) {
2848        infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2849      } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2850        infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2851      }
2852      if (infeasibility > primalTolerance) {
2853        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2854        if (infeasibility > relaxedTolerance)
2855          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2856        numberPrimalInfeasibilities_++;
2857      }
2858      infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2859    }
2860  }
2861  objectiveValue_ += objective_->nonlinearOffset();
2862  objectiveValue_ /= (objectiveScale_ * rhsScale_);
2863}
2864void ClpSimplex::checkDualSolution()
2865{
2866
2867  int iRow, iColumn;
2868  sumDualInfeasibilities_ = 0.0;
2869  numberDualInfeasibilities_ = 0;
2870  numberDualInfeasibilitiesWithoutFree_ = 0;
2871  if (matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) {
2872    // pretend we found dual infeasibilities
2873    sumOfRelaxedDualInfeasibilities_ = 1.0;
2874    sumDualInfeasibilities_ = 1.0;
2875    numberDualInfeasibilities_ = 1;
2876    return;
2877  }
2878  int firstFreePrimal = -1;
2879  int firstFreeDual = -1;
2880  int numberSuperBasicWithDj = 0;
2881  bestPossibleImprovement_ = 0.0;
2882  // we can't really trust infeasibilities if there is dual error
2883  double error = CoinMin(1.0e-2, largestDualError_);
2884  // allow tolerance at least slightly bigger than standard
2885  double relaxedTolerance = dualTolerance_ + error;
2886  // allow bigger tolerance for possible improvement
2887  double possTolerance = 5.0 * relaxedTolerance;
2888  sumOfRelaxedDualInfeasibilities_ = 0.0;
2889
2890  // Check any djs from dynamic rows
2891  matrix_->dualExpanded(this, NULL, NULL, 3);
2892  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_;
2893  objectiveValue_ = 0.0;
2894  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2895    objectiveValue_ += objectiveWork_[iColumn] * columnActivityWork_[iColumn];
2896    if (getColumnStatus(iColumn) != basic && !flagged(iColumn)) {
2897      // not basic
2898      double distanceUp = columnUpperWork_[iColumn] - columnActivityWork_[iColumn];
2899      double distanceDown = columnActivityWork_[iColumn] - columnLowerWork_[iColumn];
2900      if (distanceUp > primalTolerance_) {
2901        double value = reducedCostWork_[iColumn];
2902        // Check if "free"
2903        if (distanceDown > primalTolerance_) {
2904          if (fabs(value) > 1.0e2 * relaxedTolerance) {
2905            numberSuperBasicWithDj++;
2906            if (firstFreeDual < 0)
2907              firstFreeDual = iColumn;
2908          }
2909          if (firstFreePrimal < 0)
2910            firstFreePrimal = iColumn;
2911        }
2912        // should not be negative
2913        if (value < 0.0) {
2914          value = -value;
2915          if (value > dualTolerance_) {
2916            if (getColumnStatus(iColumn) != isFree) {
2917              numberDualInfeasibilitiesWithoutFree_++;
2918              sumDualInfeasibilities_ += value - dualTolerance_;
2919              if (value > possTolerance)
2920                bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2921              if (value > relaxedTolerance)
2922                sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2923              numberDualInfeasibilities_++;
2924            } else {
2925              // free so relax a lot
2926              value *= 0.01;
2927              if (value > dualTolerance_) {
2928                sumDualInfeasibilities_ += value - dualTolerance_;
2929                if (value > possTolerance)
2930                  bestPossibleImprovement_ = 1.0e100;
2931                if (value > relaxedTolerance)
2932                  sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2933                numberDualInfeasibilities_++;
2934              }
2935            }
2936          }
2937        }
2938      }
2939      if (distanceDown > primalTolerance_) {
2940        double value = reducedCostWork_[iColumn];
2941        // should not be positive
2942        if (value > 0.0) {
2943          if (value > dualTolerance_) {
2944            sumDualInfeasibilities_ += value - dualTolerance_;
2945            if (value > possTolerance)
2946              bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2947            if (value > relaxedTolerance)
2948              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2949            numberDualInfeasibilities_++;
2950            if (getColumnStatus(iColumn) != isFree)
2951              numberDualInfeasibilitiesWithoutFree_++;
2952            // maybe we can make feasible by increasing tolerance
2953          }
2954        }
2955      }
2956    }
2957  }
2958  for (iRow = 0; iRow < numberRows_; iRow++) {
2959    objectiveValue_ += rowActivityWork_[iRow] * rowObjectiveWork_[iRow];
2960    if (getRowStatus(iRow) != basic && !flagged(iRow + numberColumns_)) {
2961      // not basic
2962      double distanceUp = rowUpperWork_[iRow] - rowActivityWork_[iRow];
2963      double distanceDown = rowActivityWork_[iRow] - rowLowerWork_[iRow];
2964      if (distanceUp > primalTolerance_) {
2965        double value = rowReducedCost_[iRow];
2966        // Check if "free"
2967        if (distanceDown > primalTolerance_) {
2968          if (fabs(value) > 1.0e2 * relaxedTolerance) {
2969            numberSuperBasicWithDj++;
2970            if (firstFreeDual < 0)
2971              firstFreeDual = iRow + numberColumns_;
2972          }
2973          if (firstFreePrimal < 0)
2974            firstFreePrimal = iRow + numberColumns_;
2975        }
2976        // should not be negative
2977        if (value < 0.0) {
2978          value = -value;
2979          if (value > dualTolerance_) {
2980            sumDualInfeasibilities_ += value - dualTolerance_;
2981            if (value > possTolerance)
2982              bestPossibleImprovement_ += value * CoinMin(distanceUp, 1.0e10);
2983            if (value > relaxedTolerance)
2984              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2985            numberDualInfeasibilities_++;
2986            if (getRowStatus(iRow) != isFree)
2987              numberDualInfeasibilitiesWithoutFree_++;
2988          }
2989        }
2990      }
2991      if (distanceDown > primalTolerance_) {
2992        double value = rowReducedCost_[iRow];
2993        // should not be positive
2994        if (value > 0.0) {
2995          if (value > dualTolerance_) {
2996            sumDualInfeasibilities_ += value - dualTolerance_;
2997            if (value > possTolerance)
2998              bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2999            if (value > relaxedTolerance)
3000              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
3001            numberDualInfeasibilities_++;
3002            if (getRowStatus(iRow) != isFree)
3003              numberDualInfeasibilitiesWithoutFree_++;
3004            // maybe we can make feasible by increasing tolerance
3005          }
3006        }
3007      }
3008    }
3009  }
3010  if (algorithm_ < 0 && firstFreeDual >= 0) {
3011    // dual
3012    firstFree_ = firstFreeDual;
3013  } else if (numberSuperBasicWithDj || (progress_.lastIterationNumber(0) <= 0)) {
3014    firstFree_ = firstFreePrimal;
3015  }
3016  objectiveValue_ += objective_->nonlinearOffset();
3017  objectiveValue_ /= (objectiveScale_ * rhsScale_);
3018}
3019/* This sets sum and number of infeasibilities (Dual and Primal) */
3020void ClpSimplex::checkBothSolutions()
3021{
3022  if ((matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) || matrix_->rhsOffset(this)) {
3023    // Say may be free or superbasic
3024    moreSpecialOptions_ &= ~8;
3025    // old way
3026    checkPrimalSolution(rowActivityWork_, columnActivityWork_);
3027    checkDualSolution();
3028    return;
3029  }
3030  int iSequence;
3031  assert(dualTolerance_ > 0.0 && dualTolerance_ < 1.0e10);
3032  assert(primalTolerance_ > 0.0 && primalTolerance_ < 1.0e10);
3033  objectiveValue_ = 0.0;
3034  sumPrimalInfeasibilities_ = 0.0;
3035  numberPrimalInfeasibilities_ = 0;
3036  double primalTolerance = primalTolerance_;
3037  double relaxedToleranceP = primalTolerance_;
3038  // we can't really trust infeasibilities if there is primal error
3039  double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_, 0.0 * primalTolerance_));
3040  // allow tolerance at least slightly bigger than standard
3041  relaxedToleranceP = relaxedToleranceP + error;
3042  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
3043  sumDualInfeasibilities_ = 0.0;
3044  numberDualInfeasibilities_ = 0;
3045  double dualTolerance = dualTolerance_;
3046  double relaxedToleranceD = dualTolerance;
3047  // we can't really trust infeasibilities if there is dual error
3048  error = CoinMin(1.0e-2, CoinMax(largestDualError_, 5.0 * dualTolerance_));
3049  // allow tolerance at least slightly bigger than standard
3050  relaxedToleranceD = relaxedToleranceD + error;
3051  // allow bigger tolerance for possible improvement
3052  double possTolerance = 5.0 * relaxedToleranceD;
3053  sumOfRelaxedDualInfeasibilities_ = 0.0;
3054  bestPossibleImprovement_ = 0.0;
3055
3056  // Check any infeasibilities from dynamic rows
3057  matrix_->primalExpanded(this, 2);
3058  // Check any djs from dynamic rows
3059  matrix_->dualExpanded(this, NULL, NULL, 3);
3060  int numberDualInfeasibilitiesFree = 0;
3061  int firstFreePrimal = -1;
3062  int firstFreeDual = -1;
3063  int numberSuperBasicWithDj = 0;
3064
3065  int numberTotal = numberRows_ + numberColumns_;
3066  // Say no free or superbasic
3067  moreSpecialOptions_ |= 8;
3068  //#define PRINT_INFEAS
3069#ifdef PRINT_INFEAS
3070  int seqInf[10];
3071#endif
3072  for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3073    //#define LIKELY_SUPERBASIC
3074#ifdef LIKELY_SUPERBASIC
3075    if (getStatus(iSequence) == isFree || getStatus(iSequence) == superBasic)
3076      moreSpecialOptions_ &= ~8; // Say superbasic variables exist
3077#endif
3078    double value = solution_[iSequence];
3079#ifdef COIN_DEBUG
3080    if (fabs(value) > 1.0e20)
3081      printf("%d values %g %g %g - status %d\n", iSequence, lower_[iSequence],
3082        solution_[iSequence], upper_[iSequence], status_[iSequence]);
3083#endif
3084    objectiveValue_ += value * cost_[iSequence];
3085    double distanceUp = upper_[iSequence] - value;
3086    double distanceDown = value - lower_[iSequence];
3087    if (distanceUp < -primalTolerance) {
3088      double infeasibility = -distanceUp;
3089#ifndef LIKELY_SUPERBASIC
3090      if (getStatus(iSequence) != basic)
3091        moreSpecialOptions_ &= ~8; // Say superbasic variables exist
3092#endif
3093      sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
3094      if (infeasibility > relaxedToleranceP)
3095        sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
3096#ifdef PRINT_INFEAS
3097      if (numberPrimalInfeasibilities_ < 10) {
3098        seqInf[numberPrimalInfeasibilities_] = iSequence;
3099      }
3100#endif
3101      numberPrimalInfeasibilities_++;
3102    } else if (distanceDown < -primalTolerance) {
3103      double infeasibility = -distanceDown;
3104#ifndef LIKELY_SUPERBASIC
3105      if (getStatus(iSequence) != basic)
3106        moreSpecialOptions_ &= ~8; // Say superbasic variables exist
3107#endif
3108      sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
3109      if (infeasibility > relaxedToleranceP)
3110        sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
3111#ifdef PRINT_INFEAS
3112      if (numberPrimalInfeasibilities_ < 10) {
3113        seqInf[numberPrimalInfeasibilities_] = iSequence;
3114      }
3115#endif
3116      numberPrimalInfeasibilities_++;
3117    } else {
3118      // feasible (so could be free)
3119      if (getStatus(iSequence) != basic && !flagged(iSequence)) {
3120        // not basic
3121        double djValue = dj_[iSequence];
3122        if (distanceDown < primalTolerance) {
3123          if (distanceUp > primalTolerance && djValue < -dualTolerance) {
3124            sumDualInfeasibilities_ -= djValue + dualTolerance;
3125            if (djValue < -possTolerance)
3126              bestPossibleImprovement_ -= distanceUp * djValue;
3127            if (djValue < -relaxedToleranceD)
3128              sumOfRelaxedDualInfeasibilities_ -= djValue + relaxedToleranceD;
3129            numberDualInfeasibilities_++;
3130          }
3131        } else if (distanceUp < primalTolerance) {
3132          if (djValue > dualTolerance) {
3133            sumDualInfeasibilities_ += djValue - dualTolerance;
3134            if (djValue > possTolerance)
3135              bestPossibleImprovement_ += distanceDown * djValue;
3136            if (djValue > relaxedToleranceD)
3137              sumOfRelaxedDualInfeasibilities_ += djValue - relaxedToleranceD;
3138            numberDualInfeasibilities_++;
3139          }
3140        } else {
3141          // may be free
3142          // Say free or superbasic
3143          moreSpecialOptions_ &= ~8;
3144          djValue *= 0.01;
3145          if (fabs(djValue) > dualTolerance) {
3146            if (getStatus(iSequence) == isFree)
3147              numberDualInfeasibilitiesFree++;
3148            sumDualInfeasibilities_ += fabs(djValue) - dualTolerance;
3149            bestPossibleImprovement_ = 1.0e100;
3150            numberDualInfeasibilities_++;
3151            if (fabs(djValue) > relaxedToleranceD) {
3152              sumOfRelaxedDualInfeasibilities_ += value - relaxedToleranceD;
3153              numberSuperBasicWithDj++;
3154              if (firstFreeDual < 0)
3155                firstFreeDual = iSequence;
3156            }
3157          }
3158          if (firstFreePrimal < 0)
3159            firstFreePrimal = iSequence;
3160        }
3161      }
3162    }
3163  }
3164  objectiveValue_ += objective_->nonlinearOffset();
3165  objectiveValue_ /= (objectiveScale_ * rhsScale_);
3166  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberDualInfeasibilitiesFree;
3167#ifdef PRINT_INFEAS
3168  if (numberPrimalInfeasibilities_ <= 10) {
3169    printf("---------------start-----------\n");
3170    if (!rowScale_) {
3171      for (int i = 0; i < numberPrimalInfeasibilities_; i++) {
3172        int iSeq = seqInf[i];
3173        double infeas;
3174        if (solution_[iSeq] < lower_[iSeq])
3175          infeas = lower_[iSeq] - solution_[iSeq];
3176        else
3177          infeas = solution_[iSeq] - upper_[iSeq];
3178        if (iSeq < numberColumns_) {
3179          printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g\n",
3180            iSeq, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas);
3181        } else {
3182          printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g\n",
3183            iSeq - numberColumns_, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas);
3184        }
3185      }
3186    } else {
3187      for (int i = 0; i < numberPrimalInfeasibilities_; i++) {
3188        int iSeq = seqInf[i];
3189        double infeas;
3190        if (solution_[iSeq] < lower_[iSeq])
3191          infeas = lower_[iSeq] - solution_[iSeq];
3192        else
3193          infeas = solution_[iSeq] - upper_[iSeq];
3194        double unscaled = infeas;
3195        if (iSeq < numberColumns_) {
3196          unscaled *= columnScale_[iSeq];
3197          printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3198            iSeq, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas, unscaled);
3199        } else {
3200          unscaled /= rowScale_[iSeq - numberColumns_];
3201          printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3202            iSeq - numberColumns_, lower_[iSeq], solution_[iSeq], upper_[iSeq], infeas, unscaled);
3203        }
3204      }
3205    }
3206  }
3207#endif
3208  if (algorithm_ < 0 && firstFreeDual >= 0) {
3209    // dual
3210    firstFree_ = firstFreeDual;
3211  } else if (numberSuperBasicWithDj || (progress_.lastIterationNumber(0) <= 0)) {
3212    firstFree_ = firstFreePrimal;
3213  }
3214}
3215/* Adds multiple of a column into an array */
3216void ClpSimplex::add(double *array,
3217  int sequence, double multiplier) const
3218{
3219  if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3220    //slack
3221    array[sequence - numberColumns_] -= multiplier;
3222  } else {
3223    // column
3224    matrix_->add(this, array, sequence, multiplier);
3225  }
3226}
3227/*
3228  Unpacks one column of the matrix into indexed array
3229*/
3230void ClpSimplex::unpack(CoinIndexedVector *rowArray) const
3231{
3232  rowArray->clear();
3233  if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3234    //slack
3235    rowArray->insert(sequenceIn_ - numberColumns_, -1.0);
3236  } else {
3237    // column
3238    matrix_->unpack(this, rowArray, sequenceIn_);
3239  }
3240}
3241void ClpSimplex::unpack(CoinIndexedVector *rowArray, int sequence) const
3242{
3243  rowArray->clear();
3244  if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3245    //slack
3246    rowArray->insert(sequence - numberColumns_, -1.0);
3247  } else {
3248    // column
3249    matrix_->unpack(this, rowArray, sequence);
3250  }
3251}
3252/*
3253  Unpacks one column of the matrix into indexed array
3254*/
3255void ClpSimplex::unpackPacked(CoinIndexedVector *rowArray)
3256{
3257  rowArray->clear();
3258  if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3259    //slack
3260    int *index = rowArray->getIndices();
3261    double *array = rowArray->denseVector();
3262    array[0] = -1.0;
3263    index[0] = sequenceIn_ - numberColumns_;
3264    rowArray->setNumElements(1);
3265    rowArray->setPackedMode(true);
3266  } else {
3267    // column
3268    matrix_->unpackPacked(this, rowArray, sequenceIn_);
3269  }
3270}
3271void ClpSimplex::unpackPacked(CoinIndexedVector *rowArray, int sequence)
3272{
3273  rowArray->clear();
3274  if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3275    //slack
3276    int *index = rowArray->getIndices();
3277    double *array = rowArray->denseVector();
3278    array[0] = -1.0;
3279    index[0] = sequence - numberColumns_;
3280    rowArray->setNumElements(1);
3281    rowArray->setPackedMode(true);
3282  } else {
3283    // column
3284    matrix_->unpackPacked(this, rowArray, sequence);
3285  }
3286}
3287//static int x_gaps[4]={0,0,0,0};
3288//static int scale_times[]={0,0,0,0};
3289bool ClpSimplex::createRim(int what, bool makeRowCopy, int startFinishOptions)
3290{
3291  bool goodMatrix = true;
3292  int saveLevel = handler_->logLevel();
3293  spareIntArray_[0] = 0;
3294  if (!matrix_->canGetRowCopy())
3295    makeRowCopy = false; // switch off row copy if can't produce
3296  // Arrays will be there and correct size unless what is 63
3297  bool newArrays = (what == 63);
3298  // We may be restarting with same size
3299  bool keepPivots = false;
3300  if (startFinishOptions == -1) {
3301    startFinishOptions = 0;
3302    keepPivots = true;
3303  }
3304  bool oldMatrix = ((startFinishOptions & 4) != 0 && (whatsChanged_ & 1) != 0);
3305  if (what == 63) {
3306    pivotRow_ = -1;
3307    if (!status_)
3308      createStatus();
3309    if (oldMatrix)
3310      newArrays = false;
3311    if (problemStatus_ == 10) {
3312      handler_->setLogLevel(0); // switch off messages
3313      if (rowArray_[0]) {
3314        // stuff is still there
3315        oldMatrix = true;
3316        newArrays = false;
3317        keepPivots = true;
3318        for (int iRow = 0; iRow < 4; iRow++) {
3319          rowArray_[iRow]->clear();
3320        }
3321        for (int iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
3322          columnArray_[iColumn]->clear();
3323        }
3324      }
3325    } else if (factorization_) {
3326      // match up factorization messages
3327      if (handler_->logLevel() < 3)
3328        factorization_->messageLevel(0);
3329      else
3330        factorization_->messageLevel(CoinMax(3, factorization_->messageLevel()));
3331      /* Faster to keep pivots rather than re-scan matrix.  Matrix may have changed
3332                  i.e. oldMatrix false but okay as long as same number rows and status array exists
3333               */
3334      if ((startFinishOptions & 2) != 0 && factorization_->numberRows() == numberRows_ && status_)
3335        keepPivots = true;
3336    }
3337    numberExtraRows_ = matrix_->generalExpanded(this, 2, maximumBasic_);
3338    if (numberExtraRows_ && newArrays) {
3339      // make sure status array large enough
3340      assert(status_);
3341      int numberOld = numberRows_ + numberColumns_;
3342      int numberNew = numberRows_ + numberColumns_ + numberExtraRows_;
3343      unsigned char *newStatus = new unsigned char[numberNew];
3344      memset(newStatus + numberOld, 0, numberExtraRows_);
3345      CoinMemcpyN(status_, numberOld, newStatus);
3346      delete[] status_;
3347      status_ = newStatus;
3348    }
3349  }
3350  int numberRows2 = numberRows_ + numberExtraRows_;
3351  int numberTotal = numberRows2 + numberColumns_;
3352  if ((specialOptions_ & 65536) != 0) {
3353    assert(!numberExtraRows_);
3354    if (!cost_ || numberRows2 > maximumInternalRows_ || numberColumns_ > maximumInternalColumns_) {
3355      newArrays = true;
3356      keepPivots = false;
3357      COIN_DETAIL_PRINT(printf("createrim a %d rows, %d maximum rows %d maxinternal\n",
3358        numberRows_, maximumRows_, maximumInternalRows_));
3359      int oldMaximumRows = maximumInternalRows_;
3360      int oldMaximumColumns = maximumInternalColumns_;
3361      if (cost_) {
3362        if (numberRows2 > maximumInternalRows_)
3363          maximumInternalRows_ = numberRows2;
3364        if (numberColumns_ > maximumInternalColumns_)
3365          maximumInternalColumns_ = numberColumns_;
3366      } else {
3367        maximumInternalRows_ = numberRows2;
3368        maximumInternalColumns_ = numberColumns_;
3369      }
3370      //maximumRows_=CoinMax(maximumInternalRows_,maximumRows_);
3371      //maximumColumns_=CoinMax(maximumInternalColumns_,maximumColumns_);
3372      assert(maximumInternalRows_ == maximumRows_);
3373      assert(maximumInternalColumns_ == maximumColumns_);
3374      COIN_DETAIL_PRINT(printf("createrim b %d rows, %d maximum rows, %d maxinternal\n",
3375        numberRows_, maximumRows_, maximumInternalRows_));
3376      int numberTotal2 = (maximumInternalRows_ + maximumInternalColumns_) * 2;
3377      delete[] cost_;
3378      cost_ = new double[numberTotal2];
3379      delete[] lower_;
3380      delete[] upper_;
3381      lower_ = new double[numberTotal2];
3382      upper_ = new double[numberTotal2];
3383      delete[] dj_;
3384      dj_ = new double[numberTotal2];
3385      delete[] solution_;
3386      solution_ = new double[numberTotal2];
3387      // ***** should be non NULL but seems to be too much
3388      //printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_);
3389      if (savedRowScale_) {
3390        assert(oldMaximumRows > 0);
3391        double *temp;
3392        temp = new double[4 * maximumRows_];
3393        CoinFillN(temp, 4 * maximumRows_, 1.0);
3394        CoinMemcpyN(savedRowScale_, numberRows_, temp);
3395        CoinMemcpyN(savedRowScale_ + oldMaximumRows, numberRows_, temp + maximumRows_);
3396        CoinMemcpyN(savedRowScale_ + 2 * oldMaximumRows, numberRows_, temp + 2 * maximumRows_);
3397        CoinMemcpyN(savedRowScale_ + 3 * oldMaximumRows, numberRows_, temp + 3 * maximumRows_);
3398        delete[] savedRowScale_;
3399        savedRowScale_ = temp;
3400        temp = new double[4 * maximumColumns_];
3401        CoinFillN(temp, 4 * maximumColumns_, 1.0);
3402        CoinMemcpyN(savedColumnScale_, numberColumns_, temp);
3403        CoinMemcpyN(savedColumnScale_ + oldMaximumColumns, numberColumns_, temp + maximumColumns_);
3404        CoinMemcpyN(savedColumnScale_ + 2 * oldMaximumColumns, numberColumns_, temp + 2 * maximumColumns_);
3405        CoinMemcpyN(savedColumnScale_ + 3 * oldMaximumColumns, numberColumns_, temp + 3 * maximumColumns_);
3406        delete[] savedColumnScale_;
3407        savedColumnScale_ = temp;
3408      }
3409    }
3410  }
3411  int i;
3412  bool doSanityCheck = true;
3413  if (what == 63) {
3414    // We may want to switch stuff off for speed
3415    if ((specialOptions_ & 256) != 0)
3416      makeRowCopy = false; // no row copy
3417    if ((specialOptions_ & 128) != 0)
3418      doSanityCheck = false; // no sanity check
3419    //check matrix
3420    if (!matrix_)
3421      matrix_ = new ClpPackedMatrix();
3422    int checkType = (doSanityCheck) ? 15 : 14;
3423    if (oldMatrix)
3424      checkType = 14;
3425    bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
3426    if (inCbcOrOther)
3427      checkType -= 4; // don't check for duplicates
3428    if (!matrix_->allElementsInRange(this, smallElement_, 1.0e20, checkType)) {
3429      problemStatus_ = 4;
3430      secondaryStatus_ = 8;
3431      //goodMatrix= false;
3432      return false;
3433    }
3434    bool rowCopyIsScaled;
3435    if (makeRowCopy) {
3436      if (!oldMatrix || !rowCopy_) {
3437        delete rowCopy_;
3438        // may return NULL if can't give row copy
3439        rowCopy_ = matrix_->reverseOrderedCopy();
3440        rowCopyIsScaled = false;
3441      } else {
3442        rowCopyIsScaled = true;
3443      }
3444    }
3445#if 0
3446          if (what == 63) {
3447               int k = rowScale_ ? 1 : 0;
3448               if (oldMatrix)
3449                    k += 2;
3450               scale_times[k]++;
3451               if ((scale_times[0] + scale_times[1] + scale_times[2] + scale_times[3]) % 1000 == 0)
3452                    printf("scale counts %d %d %d %d\n",
3453                           scale_times[0], scale_times[1], scale_times[2], scale_times[3]);
3454          }
3455#endif
3456    // do scaling if needed
3457    if (!oldMatrix && scalingFlag_ < 0) {
3458      if (scalingFlag_ < 0 && rowScale_) {
3459        //if (handler_->logLevel()>0)
3460        printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
3461          scalingFlag_);
3462        scalingFlag_ = 0;
3463      }
3464      delete[] rowScale_;
3465      delete[] columnScale_;
3466      rowScale_ = NULL;
3467      columnScale_ = NULL;
3468    }
3469    inverseRowScale_ = NULL;
3470    inverseColumnScale_ = NULL;
3471    if (scalingFlag_ > 0 && (specialOptions_ & 65536) != 0 && rowScale_ && rowScale_ == savedRowScale_)
3472      rowScale_ = NULL;
3473    if (scalingFlag_ > 0 && !rowScale_) {
3474      if ((specialOptions_ & 65536) != 0) {
3475        assert(!rowScale_);
3476        rowScale_ = savedRowScale_;
3477        columnScale_ = savedColumnScale_;
3478        // put back original
3479        if (savedRowScale_) {
3480          inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
3481          inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
3482          CoinMemcpyN(savedRowScale_ + 2 * maximumInternalRows_,
3483            numberRows2, savedRowScale_);
3484          CoinMemcpyN(savedRowScale_ + 3 * maximumInternalRows_,
3485            numberRows2, inverseRowScale_);
3486          CoinMemcpyN(savedColumnScale_ + 2 * maximumColumns_,
3487            numberColumns_, savedColumnScale_);
3488          CoinMemcpyN(savedColumnScale_ + 3 * maximumColumns_,
3489            numberColumns_, inverseColumnScale_);
3490        }
3491      }
3492      if (matrix_->scale(this, this))
3493        scalingFlag_ = -scalingFlag_; // not scaled after all
3494      if (rowScale_ && automaticScale_) {
3495        if (!savedRowScale_) {
3496          inverseRowScale_ = rowScale_ + numberRows2;
3497          inverseColumnScale_ = columnScale_ + numberColumns_;
3498        }
3499        // try automatic scaling
3500        double smallestObj = 1.0e100;
3501        double largestObj = 0.0;
3502        double largestRhs = 0.0;
3503        const double *obj = objective();
3504        for (i = 0; i < numberColumns_; i++) {
3505          double value = fabs(obj[i]);
3506          value *= columnScale_[i];
3507          if (value && columnLower_[i] != columnUpper_[i]) {
3508            smallestObj = CoinMin(smallestObj, value);
3509            largestObj = CoinMax(largestObj, value);
3510          }
3511          if (columnLower_[i] > 0.0 || columnUpper_[i] < 0.0) {
3512            double scale = 1.0 * inverseColumnScale_[i];
3513            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3514            if (columnLower_[i] > 0)
3515              largestRhs = CoinMax(largestRhs, columnLower_[i] * scale);
3516            if (columnUpper_[i] < 0.0)
3517              largestRhs = CoinMax(largestRhs, -columnUpper_[i] * scale);
3518          }
3519        }
3520        for (i = 0; i < numberRows_; i++) {
3521          if (rowLower_[i] > 0.0 || rowUpper_[i] < 0.0) {
3522            double scale = rowScale_[i];
3523            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3524            if (rowLower_[i] > 0)
3525              largestRhs = CoinMax(largestRhs, rowLower_[i] * scale);
3526            if (rowUpper_[i] < 0.0)
3527              largestRhs = CoinMax(largestRhs, -rowUpper_[i] * scale);
3528          }
3529        }
3530        COIN_DETAIL_PRINT(printf("small obj %g, large %g - rhs %g\n", smallestObj, largestObj, largestRhs));
3531        bool scalingDone = false;
3532        // look at element range
3533        double smallestNegative;
3534        double largestNegative;
3535        double smallestPositive;
3536        double largestPositive;
3537        matrix_->rangeOfElements(smallestNegative, largestNegative,
3538          smallestPositive, largestPositive);
3539        smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
3540        largestPositive = CoinMax(fabs(largestNegative), largestPositive);
3541        if (largestObj) {
3542          double ratio = largestObj / smallestObj;
3543          double scale = 1.0;
3544          if (ratio < 1.0e8) {
3545            // reasonable
3546            if (smallestObj < 1.0e-4) {
3547              // may as well scale up
3548              scalingDone = true;
3549              scale = 1.0e-3 / smallestObj;
3550            } else if (largestObj < 1.0e6 || (algorithm_ > 0 && largestObj < 1.0e-4 * infeasibilityCost_)) {
3551              //done=true;
3552            } else {
3553              scalingDone = true;
3554              if (algorithm_ < 0) {
3555                scale = 1.0e6 / largestObj;
3556              } else {
3557                scale = CoinMax(1.0e6, 1.0e-4 * infeasibilityCost_) / largestObj;
3558              }
3559            }
3560          } else if (ratio < 1.0e12) {
3561            // not so good
3562            if (smallestObj < 1.0e-7) {
3563              // may as well scale up
3564              scalingDone = true;
3565              scale = 1.0e-6 / smallestObj;
3566            } else if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3567              //done=true;
3568            } else {
3569              scalingDone = true;
3570              if (algorithm_ < 0) {
3571                scale = 1.0e7 / largestObj;
3572              } else {
3573                scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3574              }
3575            }
3576          } else {
3577            // Really nasty problem
3578            if (smallestObj < 1.0e-8) {
3579              // may as well scale up
3580              scalingDone = true;
3581              scale = 1.0e-7 / smallestObj;
3582              largestObj *= scale;
3583            }
3584            if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3585              //done=true;
3586            } else {
3587              scalingDone = true;
3588              if (algorithm_ < 0) {
3589                scale = 1.0e7 / largestObj;
3590              } else {
3591                scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3592              }
3593            }
3594          }
3595          objectiveScale_ = scale;
3596        }
3597        if (largestRhs > 1.0e12) {
3598          scalingDone = true;
3599          rhsScale_ = 1.0e9 / largestRhs;
3600        } else if (largestPositive > 1.0e-14 * smallestPositive && largestRhs > 1.0e6) {
3601          scalingDone = true;
3602          rhsScale_ = 1.0e6 / largestRhs;
3603        } else {
3604          rhsScale_ = 1.0;
3605        }
3606        if (scalingDone) {
3607          handler_->message(CLP_RIM_SCALE, messages_)
3608            << objectiveScale_ << rhsScale_
3609            << CoinMessageEol;
3610        }
3611      }
3612    } else if (makeRowCopy && scalingFlag_ > 0 && !rowCopyIsScaled) {
3613      matrix_->scaleRowCopy(this);
3614    }
3615    if (rowScale_ && !savedRowScale_) {
3616      inverseRowScale_ = rowScale_ + numberRows2;
3617      inverseColumnScale_ = columnScale_ + numberColumns_;
3618    }
3619    // See if we can try for faster row copy
3620    if (makeRowCopy && !oldMatrix) {
3621      ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix_);
3622      if (clpMatrix && numberThreads_)
3623        clpMatrix->specialRowCopy(this, rowCopy_);
3624      if (clpMatrix)
3625        clpMatrix->specialColumnCopy(this);
3626    }
3627  }
3628  if (what == 63) {
3629#if 0
3630          {
3631               x_gaps[0]++;
3632               ClpPackedMatrix* clpMatrix =
3633               dynamic_cast< ClpPackedMatrix*>(matrix_);
3634               if (clpMatrix) {
3635                    if (!clpMatrix->getPackedMatrix()->hasGaps())
3636                         x_gaps[1]++;
3637                    if ((clpMatrix->flags() & 2) == 0)
3638                         x_gaps[3]++;
3639               } else {
3640                    x_gaps[2]++;
3641               }
3642               if ((x_gaps[0] % 1000) == 0)
3643                    printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n",
3644                           x_gaps[0], x_gaps[1], x_gaps[2], x_gaps[3]);
3645          }
3646#endif
3647    if (newArrays && (specialOptions_ & 65536) == 0) {
3648      delete[] cost_;
3649      cost_ = new double[2 * numberTotal];
3650      delete[] lower_;
3651      delete[] upper_;
3652      lower_ = new double[numberTotal];
3653      upper_ = new double[numberTotal];
3654      delete[] dj_;
3655      dj_ = new double[numberTotal];
3656      delete[] solution_;
3657      solution_ = new double[numberTotal];
3658    }
3659    reducedCostWork_ = dj_;
3660    rowReducedCost_ = dj_ + numberColumns_;
3661    columnActivityWork_ = solution_;
3662    rowActivityWork_ = solution_ + numberColumns_;
3663    objectiveWork_ = cost_;
3664    rowObjectiveWork_ = cost_ + numberColumns_;
3665    rowLowerWork_ = lower_ + numberColumns_;
3666    columnLowerWork_ = lower_;
3667    rowUpperWork_ = upper_ + numberColumns_;
3668    columnUpperWork_ = upper_;
3669  }
3670  if ((what & 4) != 0) {
3671    double direction = optimizationDirection_ * objectiveScale_;
3672    const double *obj = objective();
3673    const double *rowScale = rowScale_;
3674    const double *columnScale = columnScale_;
3675    // and also scale by scale factors
3676    if (rowScale) {
3677      if (rowObjective_) {
3678        for (i = 0; i < numberRows_; i++)
3679          rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
3680      } else {
3681        memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3682      }
3683      // If scaled then do all columns later in one loop
3684      if (what != 63) {
3685        for (i = 0; i < numberColumns_; i++) {
3686          CoinAssert(fabs(obj[i]) < 1.0e25);
3687          objectiveWork_[i] = obj[i] * direction * columnScale[i];
3688        }
3689      }
3690    } else {
3691      if (rowObjective_) {
3692        for (i = 0; i < numberRows_; i++)
3693          rowObjectiveWork_[i] = rowObjective_[i] * direction;
3694      } else {
3695        memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3696      }
3697      for (i = 0; i < numberColumns_; i++) {
3698        CoinAssert(fabs(obj[i]) < 1.0e25);
3699        objectiveWork_[i] = obj[i] * direction;
3700      }
3701    }
3702  }
3703  if ((what & 1) != 0) {
3704    const double *rowScale = rowScale_;
3705    // clean up any mismatches on infinity
3706    // and fix any variables with tiny gaps
3707    double primalTolerance = dblParam_[ClpPrimalTolerance];
3708    if (rowScale) {
3709      // If scaled then do all columns later in one loop
3710      if (what != 63) {
3711        const double *inverseScale = inverseColumnScale_;
3712        for (i = 0; i < numberColumns_; i++) {
3713          double multiplier = rhsScale_ * inverseScale[i];
3714          double lowerValue = columnLower_[i];
3715          double upperValue = columnUpper_[i];
3716          if (lowerValue > -1.0e20) {
3717            columnLowerWork_[i] = lowerValue * multiplier;
3718            if (upperValue >= 1.0e20) {
3719              columnUpperWork_[i] = COIN_DBL_MAX;
3720            } else {
3721              columnUpperWork_[i] = upperValue * multiplier;
3722              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3723                if (columnLowerWork_[i] >= 0.0) {
3724                  columnUpperWork_[i] = columnLowerWork_[i];
3725                } else if (columnUpperWork_[i] <= 0.0) {
3726                  columnLowerWork_[i] = columnUpperWork_[i];
3727                } else {
3728                  columnUpperWork_[i] = 0.0;
3729                  columnLowerWork_[i] = 0.0;
3730                }
3731              }
3732            }
3733          } else if (upperValue < 1.0e20) {
3734            columnLowerWork_[i] = -COIN_DBL_MAX;
3735            columnUpperWork_[i] = upperValue * multiplier;
3736          } else {
3737            // free
3738            columnLowerWork_[i] = -COIN_DBL_MAX;
3739            columnUpperWork_[i] = COIN_DBL_MAX;
3740          }
3741        }
3742      }
3743      for (i = 0; i < numberRows_; i++) {
3744        double multiplier = rhsScale_ * rowScale[i];
3745        double lowerValue = rowLower_[i];
3746        double upperValue = rowUpper_[i];
3747        if (lowerValue > -1.0e20) {
3748          rowLowerWork_[i] = lowerValue * multiplier;
3749          if (upperValue >= 1.0e20) {
3750            rowUpperWork_[i] = COIN_DBL_MAX;
3751          } else {
3752            rowUpperWork_[i] = upperValue * multiplier;
3753            if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3754              if (rowLowerWork_[i] >= 0.0) {
3755                rowUpperWork_[i] = rowLowerWork_[i];
3756              } else if (rowUpperWork_[i] <= 0.0) {
3757                rowLowerWork_[i] = rowUpperWork_[i];
3758              } else {
3759                rowUpperWork_[i] = 0.0;
3760                rowLowerWork_[i] = 0.0;
3761              }
3762            }
3763          }
3764        } else if (upperValue < 1.0e20) {
3765          rowLowerWork_[i] = -COIN_DBL_MAX;
3766          rowUpperWork_[i] = upperValue * multiplier;
3767        } else {
3768          // free
3769          rowLowerWork_[i] = -COIN_DBL_MAX;
3770          rowUpperWork_[i] = COIN_DBL_MAX;
3771        }
3772      }
3773    } else if (rhsScale_ != 1.0) {
3774      for (i = 0; i < numberColumns_; i++) {
3775        double lowerValue = columnLower_[i];
3776        double upperValue = columnUpper_[i];
3777        if (lowerValue > -1.0e20) {
3778          columnLowerWork_[i] = lowerValue * rhsScale_;
3779          if (upperValue >= 1.0e20) {
3780            columnUpperWork_[i] = COIN_DBL_MAX;
3781          } else {
3782            columnUpperWork_[i] = upperValue * rhsScale_;
3783            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3784              if (columnLowerWork_[i] >= 0.0) {
3785                columnUpperWork_[i] = columnLowerWork_[i];
3786              } else if (columnUpperWork_[i] <= 0.0) {
3787                columnLowerWork_[i] = columnUpperWork_[i];
3788              } else {
3789                columnUpperWork_[i] = 0.0;
3790                columnLowerWork_[i] = 0.0;
3791              }
3792            }
3793          }
3794        } else if (upperValue < 1.0e20) {
3795          columnLowerWork_[i] = -COIN_DBL_MAX;
3796          columnUpperWork_[i] = upperValue * rhsScale_;
3797        } else {
3798          // free
3799          columnLowerWork_[i] = -COIN_DBL_MAX;
3800          columnUpperWork_[i] = COIN_DBL_MAX;
3801        }
3802      }
3803      for (i = 0; i < numberRows_; i++) {
3804        double lowerValue = rowLower_[i];
3805        double upperValue = rowUpper_[i];
3806        if (lowerValue > -1.0e20) {
3807          rowLowerWork_[i] = lowerValue * rhsScale_;
3808          if (upperValue >= 1.0e20) {
3809            rowUpperWork_[i] = COIN_DBL_MAX;
3810          } else {
3811            rowUpperWork_[i] = upperValue * rhsScale_;
3812            if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3813              if (rowLowerWork_[i] >= 0.0) {
3814                rowUpperWork_[i] = rowLowerWork_[i];
3815              } else if (rowUpperWork_[i] <= 0.0) {
3816                rowLowerWork_[i] = rowUpperWork_[i];
3817              } else {
3818                rowUpperWork_[i] = 0.0;
3819                rowLowerWork_[i] = 0.0;
3820              }
3821            }
3822          }
3823        } else if (upperValue < 1.0e20) {
3824          rowLowerWork_[i] = -COIN_DBL_MAX;
3825          rowUpperWork_[i] = upperValue * rhsScale_;
3826        } else {
3827          // free
3828          rowLowerWork_[i] = -COIN_DBL_MAX;
3829          rowUpperWork_[i] = COIN_DBL_MAX;
3830        }
3831      }
3832    } else {
3833      for (i = 0; i < numberColumns_; i++) {
3834        double lowerValue = columnLower_[i];
3835        double upperValue = columnUpper_[i];
3836        if (lowerValue > -1.0e20) {
3837          columnLowerWork_[i] = lowerValue;
3838          if (upperValue >= 1.0e20) {
3839            columnUpperWork_[i] = COIN_DBL_MAX;
3840          } else {
3841            columnUpperWork_[i] = upperValue;
3842            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3843              if (columnLowerWork_[i] >= 0.0) {
3844                columnUpperWork_[i] = columnLowerWork_[i];
3845              } else if (columnUpperWork_[i] <= 0.0) {
3846                columnLowerWork_[i] = columnUpperWork_[i];
3847              } else {
3848                columnUpperWork_[i] = 0.0;
3849                columnLowerWork_[i] = 0.0;
3850              }
3851            }
3852          }
3853        } else if (upperValue < 1.0e20) {
3854          columnLowerWork_[i] = -COIN_DBL_MAX;
3855          columnUpperWork_[i] = upperValue;
3856        } else {
3857          // free
3858          columnLowerWork_[i] = -COIN_DBL_MAX;
3859          columnUpperWork_[i] = COIN_DBL_MAX;
3860        }
3861      }
3862      for (i = 0; i < numberRows_; i++) {
3863        double lowerValue = rowLower_[i];
3864        double upperValue = rowUpper_[i];
3865        if (lowerValue > -1.0e20) {
3866          rowLowerWork_[i] = lowerValue;
3867          if (upperValue >= 1.0e20) {
3868            rowUpperWork_[i] = COIN_DBL_MAX;
3869          } else {
3870            rowUpperWork_[i] = upperValue;
3871            if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3872              if (rowLowerWork_[i] >= 0.0) {
3873                rowUpperWork_[i] = rowLowerWork_[i];
3874              } else if (rowUpperWork_[i] <= 0.0) {
3875                rowLowerWork_[i] = rowUpperWork_[i];
3876              } else {
3877                rowUpperWork_[i] = 0.0;
3878                rowLowerWork_[i] = 0.0;
3879              }
3880            }
3881          }
3882        } else if (upperValue < 1.0e20) {
3883          rowLowerWork_[i] = -COIN_DBL_MAX;
3884          rowUpperWork_[i] = upperValue;
3885        } else {
3886          // free
3887          rowLowerWork_[i] = -COIN_DBL_MAX;
3888          rowUpperWork_[i] = COIN_DBL_MAX;
3889        }
3890      }
3891    }
3892  }
3893  if (what == 63) {
3894    // move information to work arrays
3895    double direction = optimizationDirection_;
3896    // direction is actually scale out not scale in
3897    if (direction)
3898      direction = 1.0 / direction;
3899    if (direction != 1.0) {
3900      // reverse all dual signs
3901      for (i = 0; i < numberColumns_; i++)
3902        reducedCost_[i] *= direction;
3903      for (i = 0; i < numberRows_; i++)
3904        dual_[i] *= direction;
3905    }
3906    for (i = 0; i < numberRows_ + numberColumns_; i++) {
3907      setFakeBound(i, noFake);
3908    }
3909    if (rowScale_) {
3910      const double *obj = objective();
3911      double direction = optimizationDirection_ * objectiveScale_;
3912      // clean up any mismatches on infinity
3913      // and fix any variables with tiny gaps
3914      double primalTolerance = dblParam_[ClpPrimalTolerance];
3915      // on entry
3916      const double *inverseScale = inverseColumnScale_;
3917      for (i = 0; i < numberColumns_; i++) {
3918        CoinAssert(fabs(obj[i]) < 1.0e25);
3919        double scaleFactor = columnScale_[i];
3920        double multiplier = rhsScale_ * inverseScale[i];
3921        scaleFactor *= direction;
3922        objectiveWork_[i] = obj[i] * scaleFactor;
3923        reducedCostWork_[i] = reducedCost_[i] * scaleFactor;
3924        double lowerValue = columnLower_[i];
3925        double upperValue = columnUpper_[i];
3926        if (lowerValue > -1.0e20) {
3927          columnLowerWork_[i] = lowerValue * multiplier;
3928          if (upperValue >= 1.0e20) {
3929            columnUpperWork_[i] = COIN_DBL_MAX;
3930          } else {
3931            columnUpperWork_[i] = upperValue * multiplier;
3932            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3933              if (columnLowerWork_[i] >= 0.0) {
3934                columnUpperWork_[i] = columnLowerWork_[i];
3935              } else if (columnUpperWork_[i] <= 0.0) {
3936                columnLowerWork_[i] = columnUpperWork_[i];
3937              } else {
3938                columnUpperWork_[i] = 0.0;
3939                columnLowerWork_[i] = 0.0;
3940              }
3941            }
3942          }
3943        } else if (upperValue < 1.0e20) {
3944          columnLowerWork_[i] = -COIN_DBL_MAX;
3945          columnUpperWork_[i] = upperValue * multiplier;
3946        } else {
3947          // free
3948          columnLowerWork_[i] = -COIN_DBL_MAX;
3949          columnUpperWork_[i] = COIN_DBL_MAX;
3950        }
3951        double value = columnActivity_[i] * multiplier;
3952        if (fabs(value) > 1.0e20) {
3953          //printf("bad value of %g for column %d\n",value,i);
3954          setColumnStatus(i, superBasic);
3955          if (columnUpperWork_[i] < 0.0) {
3956            value = columnUpperWork_[i];
3957          } else if (columnLowerWork_[i] > 0.0) {
3958            value = columnLowerWork_[i];
3959          } else {
3960            value = 0.0;
3961          }
3962        }
3963        columnActivityWork_[i] = value;
3964      }
3965      inverseScale = inverseRowScale_;
3966      for (i = 0; i < numberRows_; i++) {
3967        dual_[i] *= inverseScale[i];
3968        dual_[i] *= objectiveScale_;
3969        rowReducedCost_[i] = dual_[i];
3970        double multiplier = rhsScale_ * rowScale_[i];
3971        double value = rowActivity_[i] * multiplier;
3972        if (fabs(value) > 1.0e20) {
3973          //printf("bad value of %g for row %d\n",value,i);
3974          setRowStatus(i, superBasic);
3975          if (rowUpperWork_[i] < 0.0) {
3976            value = rowUpperWork_[i];
3977          } else if (rowLowerWork_[i] > 0.0) {
3978            value = rowLowerWork_[i];
3979          } else {
3980            value = 0.0;
3981          }
3982        }
3983        rowActivityWork_[i] = value;
3984      }
3985    } else if (objectiveScale_ != 1.0 || rhsScale_ != 1.0) {
3986      // on entry
3987      for (i = 0; i < numberColumns_; i++) {
3988        double value = columnActivity_[i];
3989        value *= rhsScale_;
3990        if (fabs(value) > 1.0e20) {
3991          //printf("bad value of %g for column %d\n",value,i);
3992          setColumnStatus(i, superBasic);
3993          if (columnUpperWork_[i] < 0.0) {
3994            value = columnUpperWork_[i];
3995          } else if (columnLowerWork_[i] > 0.0) {
3996            value = columnLowerWork_[i];
3997          } else {
3998            value = 0.0;
3999          }
4000        }
4001        columnActivityWork_[i] = value;
4002        reducedCostWork_[i] = reducedCost_[i] * objectiveScale_;
4003      }
4004      for (i = 0; i < numberRows_; i++) {
4005        double value = rowActivity_[i];
4006        value *= rhsScale_;
4007        if (fabs(value) > 1.0e20) {
4008          //printf("bad value of %g for row %d\n",value,i);
4009          setRowStatus(i, superBasic);
4010          if (rowUpperWork_[i] < 0.0) {
4011            value = rowUpperWork_[i];
4012          } else if (rowLowerWork_[i] > 0.0) {
4013            value = rowLowerWork_[i];
4014          } else {
4015            value = 0.0;
4016          }
4017        }
4018        rowActivityWork_[i] = value;
4019        dual_[i] *= objectiveScale_;
4020        rowReducedCost_[i] = dual_[i];
4021      }
4022    } else {
4023      // on entry
4024      for (i = 0; i < numberColumns_; i++) {
4025        double value = columnActivity_[i];
4026        if (fabs(value) > 1.0e20) {
4027          //printf("bad value of %g for column %d\n",value,i);
4028          setColumnStatus(i, superBasic);
4029          if (columnUpperWork_[i] < 0.0) {
4030            value = columnUpperWork_[i];
4031          } else if (columnLowerWork_[i] > 0.0) {
4032            value = columnLowerWork_[i];
4033          } else {
4034            value = 0.0;
4035          }
4036        }
4037        columnActivityWork_[i] = value;
4038        reducedCostWork_[i] = reducedCost_[i];
4039      }
4040      for (i = 0; i < numberRows_; i++) {
4041        double value = rowActivity_[i];
4042        if (fabs(value) > 1.0e20) {
4043          //printf("bad value of %g for row %d\n",value,i);
4044          setRowStatus(i, superBasic);
4045          if (rowUpperWork_[i] < 0.0) {
4046            value = rowUpperWork_[i];
4047          } else if (rowLowerWork_[i] > 0.0) {
4048            value = rowLowerWork_[i];
4049          } else {
4050            value = 0.0;
4051          }
4052        }
4053        rowActivityWork_[i] = value;
4054        rowReducedCost_[i] = dual_[i];
4055      }
4056    }
4057  }
4058
4059  if (what == 63 && doSanityCheck) {
4060    // check rim of problem okay
4061    if (!sanityCheck())
4062      goodMatrix = false;
4063  }
4064  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
4065  // maybe we need to move scales to SimplexModel for factorization?
4066  if ((what == 63 && !pivotVariable_) || (newArrays && !keepPivots)) {
4067    delete[] pivotVariable_;
4068    pivotVariable_ = new int[numberRows2 + 1];
4069    for (int i = 0; i < numberRows2 + 1; i++)
4070      pivotVariable_[i] = -1;
4071  } else if (what == 63 && !keepPivots) {
4072    // just reset
4073    for (int i = 0; i < numberRows2; i++)
4074      pivotVariable_[i] = -1;
4075  } else if (what == 63) {
4076    // check pivots
4077    for (int i = 0; i < numberRows2; i++) {
4078      int iSequence = pivotVariable_[i];
4079      if (iSequence < 0 || (iSequence < numberRows_ + numberColumns_ && getStatus(iSequence) != basic)) {
4080        keepPivots = false;
4081        break;
4082      }
4083    }
4084    if (!keepPivots) {
4085      // reset
4086      for (int i = 0; i < numberRows2; i++)
4087        pivotVariable_[i] = -1;
4088    } else {
4089      // clean
4090      for (int i = 0; i < numberColumns_ + numberRows_; i++) {
4091        Status status = getStatus(i);
4092        if (status != basic) {
4093          if (upper_[i] == lower_[i]) {
4094            setStatus(i, isFixed);
4095            solution_[i] = lower_[i];
4096          } else if (status == atLowerBound) {
4097            if (lower_[i] > -1.0e20) {
4098              solution_[i] = lower_[i];
4099            } else {
4100              //printf("seq %d at lower of %g\n",i,lower_[i]);
4101              if (upper_[i] < 1.0e20) {
4102                solution_[i] = upper_[i];
4103                setStatus(i, atUpperBound);
4104              } else {
4105                setStatus(i, isFree);
4106              }
4107            }
4108          } else if (status == atUpperBound) {
4109            if (upper_[i] < 1.0e20) {
4110              solution_[i] = upper_[i];
4111            } else {
4112              //printf("seq %d at upper of %g\n",i,upper_[i]);
4113              if (lower_[i] > -1.0e20) {
4114                solution_[i] = lower_[i];
4115                setStatus(i, atLowerBound);
4116              } else {
4117                setStatus(i, isFree);
4118              }
4119            }
4120          } else if (status == isFixed && upper_[i] > lower_[i]) {
4121            // was fixed - not now
4122            if (solution_[i] <= lower_[i] + primalTolerance_) {
4123              setStatus(i, atLowerBound);
4124            } else if (solution_[i] >= upper_[i] - primalTolerance_) {
4125              setStatus(i, atUpperBound);
4126            } else {
4127              setStatus(i, superBasic);
4128            }
4129          }
4130        }
4131      }
4132    }
4133  }
4134
4135  if (what == 63) {
4136    // Safer to get new arrays anyway so test on newArrays removed
4137    // get some arrays
4138    int iRow, iColumn;
4139    // these are "indexed" arrays so we always know where nonzeros are
4140    /**********************************************************
4141               rowArray_[3] is long enough for rows+columns (2 also maybe)
4142               rowArray_[1] is long enough for max(rows,columns)
4143          *********************************************************/
4144    for (iRow = 0; iRow < 4; iRow++) {
4145      int length = numberRows2 + factorization_->maximumPivots();
4146      if (iRow > SHORT_REGION || objective_->type() > 1)
4147        length += numberColumns_;
4148      else if (iRow == 1)
4149        length = CoinMax(length, numberColumns_);
4150      if ((specialOptions_ & 65536) == 0 || !rowArray_[iRow]) {
4151        delete rowArray_[iRow];
4152        rowArray_[iRow] = new CoinIndexedVector();
4153      }
4154      rowArray_[iRow]->reserve(length);
4155    }
4156
4157    for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
4158      if ((specialOptions_ & 65536) == 0 || !columnArray_[iColumn]) {
4159        delete columnArray_[iColumn];
4160        columnArray_[iColumn] = new CoinIndexedVector();
4161      }
4162      columnArray_[iColumn]->reserve(numberColumns_ + numberRows2);
4163    }
4164  }
4165  if (problemStatus_ == 10) {
4166    problemStatus_ = -1;
4167    handler_->setLogLevel(saveLevel); // switch back messages
4168  }
4169  if ((what & 5) != 0)
4170    matrix_->generalExpanded(this, 9, what); // update costs and bounds if necessary
4171  if (goodMatrix && (specialOptions_ & 65536) != 0) {
4172    int save = maximumColumns_ + maximumRows_;
4173    CoinMemcpyN(cost_, numberTotal, cost_ + save);
4174    CoinMemcpyN(lower_, numberTotal, lower_ + save);
4175    CoinMemcpyN(upper_, numberTotal, upper_ + save);
4176    CoinMemcpyN(dj_, numberTotal, dj_ + save);
4177    CoinMemcpyN(solution_, numberTotal, solution_ + save);
4178    if (rowScale_ && !savedRowScale_) {
4179      double *temp;
4180      temp = new double[4 * maximumRows_];
4181      CoinFillN(temp, 4 * maximumRows_, 1.0);
4182      CoinMemcpyN(rowScale_, numberRows2, temp);
4183      CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + maximumRows_);
4184      CoinMemcpyN(rowScale_, numberRows2, temp + 2 * maximumRows_);
4185      CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + 3 * maximumRows_);
4186      delete[] rowScale_;
4187      savedRowScale_ = temp;
4188      rowScale_ = savedRowScale_;
4189      inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
4190      temp = new double[4 * maximumColumns_];
4191      CoinFillN(temp, 4 * maximumColumns_, 1.0);
4192      CoinMemcpyN(columnScale_, numberColumns_, temp);
4193      CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + maximumColumns_);
4194      CoinMemcpyN(columnScale_, numberColumns_, temp + 2 * maximumColumns_);
4195      CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + 3 * maximumColumns_);
4196      delete[] columnScale_;
4197      savedColumnScale_ = temp;
4198      columnScale_ = savedColumnScale_;
4199      inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
4200    }
4201  }
4202#ifdef CLP_USER_DRIVEN
4203  eventHandler_->event(ClpEventHandler::endOfCreateRim);
4204#endif
4205  return goodMatrix;
4206}
4207// Does rows and columns
4208void ClpSimplex::createRim1(bool initial)
4209{
4210  int i;
4211  int numberRows2 = numberRows_ + numberExtraRows_;
4212  int numberTotal = numberRows2 + numberColumns_;
4213  if ((specialOptions_ & 65536) != 0 && true) {
4214    assert(!initial);
4215    int save = maximumColumns_ + maximumRows_;
4216    CoinMemcpyN(lower_ + save, numberTotal, lower_);
4217    CoinMemcpyN(upper_ + save, numberTotal, upper_);
4218    return;
4219  }
4220  const double *rowScale = rowScale_;
4221  // clean up any mismatches on infinity
4222  // and fix any variables with tiny gaps
4223  double primalTolerance = dblParam_[ClpPrimalTolerance];
4224  if (rowScale) {
4225    // If scaled then do all columns later in one loop
4226    if (!initial) {
4227      const double *inverseScale = inverseColumnScale_;
4228      for (i = 0; i < numberColumns_; i++) {
4229        double multiplier = rhsScale_ * inverseScale[i];
4230        double lowerValue = columnLower_[i];
4231        double upperValue = columnUpper_[i];
4232        if (lowerValue > -1.0e20) {
4233          columnLowerWork_[i] = lowerValue * multiplier;
4234          if (upperValue >= 1.0e20) {
4235            columnUpperWork_[i] = COIN_DBL_MAX;
4236          } else {
4237            columnUpperWork_[i] = upperValue * multiplier;
4238            if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4239              if (columnLowerWork_[i] >= 0.0) {
4240                columnUpperWork_[i] = columnLowerWork_[i];
4241              } else if (columnUpperWork_[i] <= 0.0) {
4242                columnLowerWork_[i] = columnUpperWork_[i];
4243              } else {
4244                columnUpperWork_[i] = 0.0;
4245                columnLowerWork_[i] = 0.0;
4246              }
4247            }
4248          }
4249        } else if (upperValue < 1.0e20) {
4250          columnLowerWork_[i] = -COIN_DBL_MAX;
4251          columnUpperWork_[i] = upperValue * multiplier;
4252        } else {
4253          // free
4254          columnLowerWork_[i] = -COIN_DBL_MAX;
4255          columnUpperWork_[i] = COIN_DBL_MAX;
4256        }
4257      }
4258    }
4259    for (i = 0; i < numberRows_; i++) {
4260      double multiplier = rhsScale_ * rowScale[i];
4261      double lowerValue = rowLower_[i];
4262      double upperValue = rowUpper_[i];
4263      if (lowerValue > -1.0e20) {
4264        rowLowerWork_[i] = lowerValue * multiplier;
4265        if (upperValue >= 1.0e20) {
4266          rowUpperWork_[i] = COIN_DBL_MAX;
4267        } else {
4268          rowUpperWork_[i] = upperValue * multiplier;
4269          if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4270            if (rowLowerWork_[i] >= 0.0) {
4271              rowUpperWork_[i] = rowLowerWork_[i];
4272            } else if (rowUpperWork_[i] <= 0.0) {
4273              rowLowerWork_[i] = rowUpperWork_[i];
4274            } else {
4275              rowUpperWork_[i] = 0.0;
4276              rowLowerWork_[i] = 0.0;
4277            }
4278          }
4279        }
4280      } else if (upperValue < 1.0e20) {
4281        rowLowerWork_[i] = -COIN_DBL_MAX;
4282        rowUpperWork_[i] = upperValue * multiplier;
4283      } else {
4284        // free
4285        rowLowerWork_[i] = -COIN_DBL_MAX;
4286        rowUpperWork_[i] = COIN_DBL_MAX;
4287      }
4288    }
4289  } else if (rhsScale_ != 1.0) {
4290    for (i = 0; i < numberColumns_; i++) {
4291      double lowerValue = columnLower_[i];
4292      double upperValue = columnUpper_[i];
4293      if (lowerValue > -1.0e20) {
4294        columnLowerWork_[i] = lowerValue * rhsScale_;
4295        if (upperValue >= 1.0e20) {
4296          columnUpperWork_[i] = COIN_DBL_MAX;
4297        } else {
4298          columnUpperWork_[i] = upperValue * rhsScale_;
4299          if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4300            if (columnLowerWork_[i] >= 0.0) {
4301              columnUpperWork_[i] = columnLowerWork_[i];
4302            } else if (columnUpperWork_[i] <= 0.0) {
4303              columnLowerWork_[i] = columnUpperWork_[i];
4304            } else {
4305              columnUpperWork_[i] = 0.0;
4306              columnLowerWork_[i] = 0.0;
4307            }
4308          }
4309        }
4310      } else if (upperValue < 1.0e20) {
4311        columnLowerWork_[i] = -COIN_DBL_MAX;
4312        columnUpperWork_[i] = upperValue * rhsScale_;
4313      } else {
4314        // free
4315        columnLowerWork_[i] = -COIN_DBL_MAX;
4316        columnUpperWork_[i] = COIN_DBL_MAX;
4317      }
4318    }
4319    for (i = 0; i < numberRows_; i++) {
4320      double lowerValue = rowLower_[i];
4321      double upperValue = rowUpper_[i];
4322      if (lowerValue > -1.0e20) {
4323        rowLowerWork_[i] = lowerValue * rhsScale_;
4324        if (upperValue >= 1.0e20) {
4325          rowUpperWork_[i] = COIN_DBL_MAX;
4326        } else {
4327          rowUpperWork_[i] = upperValue * rhsScale_;
4328          if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4329            if (rowLowerWork_[i] >= 0.0) {
4330              rowUpperWork_[i] = rowLowerWork_[i];
4331            } else if (rowUpperWork_[i] <= 0.0) {
4332              rowLowerWork_[i] = rowUpperWork_[i];
4333            } else {
4334              rowUpperWork_[i] = 0.0;
4335              rowLowerWork_[i] = 0.0;
4336            }
4337          }
4338        }
4339      } else if (upperValue < 1.0e20) {
4340        rowLowerWork_[i] = -COIN_DBL_MAX;
4341        rowUpperWork_[i] = upperValue * rhsScale_;
4342      } else {
4343        // free
4344        rowLowerWork_[i] = -COIN_DBL_MAX;
4345        rowUpperWork_[i] = COIN_DBL_MAX;
4346      }
4347    }
4348  } else {
4349    for (i = 0; i < numberColumns_; i++) {
4350      double lowerValue = columnLower_[i];
4351      double upperValue = columnUpper_[i];
4352      if (lowerValue > -1.0e20) {
4353        columnLowerWork_[i] = lowerValue;
4354        if (upperValue >= 1.0e20) {
4355          columnUpperWork_[i] = COIN_DBL_MAX;
4356        } else {
4357          columnUpperWork_[i] = upperValue;
4358          if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4359            if (columnLowerWork_[i] >= 0.0) {
4360              columnUpperWork_[i] = columnLowerWork_[i];
4361            } else if (columnUpperWork_[i] <= 0.0) {
4362              columnLowerWork_[i] = columnUpperWork_[i];
4363            } else {
4364              columnUpperWork_[i] = 0.0;
4365              columnLowerWork_[i] = 0.0;
4366            }
4367          }
4368        }
4369      } else if (upperValue < 1.0e20) {
4370        columnLowerWork_[i] = -COIN_DBL_MAX;
4371        columnUpperWork_[i] = upperValue;
4372      } else {
4373        // free
4374        columnLowerWork_[i] = -COIN_DBL_MAX;
4375        columnUpperWork_[i] = COIN_DBL_MAX;
4376      }
4377    }
4378    for (i = 0; i < numberRows_; i++) {
4379      double lowerValue = rowLower_[i];
4380      double upperValue = rowUpper_[i];
4381      if (lowerValue > -1.0e20) {
4382        rowLowerWork_[i] = lowerValue;
4383        if (upperValue >= 1.0e20) {
4384          rowUpperWork_[i] = COIN_DBL_MAX;
4385        } else {
4386          rowUpperWork_[i] = upperValue;
4387          if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4388            if (rowLowerWork_[i] >= 0.0) {
4389              rowUpperWork_[i] = rowLowerWork_[i];
4390            } else if (rowUpperWork_[i] <= 0.0) {
4391              rowLowerWork_[i] = rowUpperWork_[i];
4392            } else {
4393              rowUpperWork_[i] = 0.0;
4394              rowLowerWork_[i] = 0.0;
4395            }
4396          }
4397        }
4398      } else if (upperValue < 1.0e20) {
4399        rowLowerWork_[i] = -COIN_DBL_MAX;
4400        rowUpperWork_[i] = upperValue;
4401      } else {
4402        // free
4403        rowLowerWork_[i] = -COIN_DBL_MAX;
4404        rowUpperWork_[i] = COIN_DBL_MAX;
4405      }
4406    }
4407  }
4408#ifndef NDEBUG
4409  if ((specialOptions_ & 65536) != 0 && false) {
4410    assert(!initial);
4411    int save = maximumColumns_ + maximumRows_;
4412    for (int i = 0; i < numberTotal; i++) {
4413      assert(fabs(lower_[i] - lower_[i + save]) < 1.0e-5);
4414      assert(fabs(upper_[i] - upper_[i + save]) < 1.0e-5);
4415    }
4416  }
4417#endif
4418}
4419// Does objective
4420void ClpSimplex::createRim4(bool initial)
4421{
4422  int i;
4423  int numberRows2 = numberRows_ + numberExtraRows_;
4424  int numberTotal = numberRows2 + numberColumns_;
4425  if ((specialOptions_ & 65536) != 0 && true) {
4426    assert(!initial);
4427    int save = maximumColumns_ + maximumRows_;
4428    CoinMemcpyN(cost_ + save, numberTotal, cost_);
4429    return;
4430  }
4431  double direction = optimizationDirection_ * objectiveScale_;
4432  const double *obj = objective();
4433  const double *rowScale = rowScale_;
4434  const double *columnScale = columnScale_;
4435  // and also scale by scale factors
4436  if (rowScale) {
4437    if (rowObjective_) {
4438      for (i = 0; i < numberRows_; i++)
4439        rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
4440    } else {
4441      memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4442    }
4443    // If scaled then do all columns later in one loop
4444    if (!initial) {
4445      for (i = 0; i < numberColumns_; i++) {
4446        CoinAssert(fabs(obj[i]) < 1.0e25);
4447        objectiveWork_[i] = obj[i] * direction * columnScale[i];
4448      }
4449    }
4450  } else {
4451    if (rowObjective_) {
4452      for (i = 0; i < numberRows_; i++)
4453        rowObjectiveWork_[i] = rowObjective_[i] * direction;
4454    } else {
4455      memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4456    }
4457    for (i = 0; i < numberColumns_; i++) {
4458      CoinAssert(fabs(obj[i]) < 1.0e25);
4459      objectiveWork_[i] = obj[i] * direction;
4460    }
4461  }
4462}
4463// Does rows and columns and objective
4464void ClpSimplex::createRim5(bool initial)
4465{
4466  createRim4(initial);
4467  createRim1(initial);
4468}
4469void ClpSimplex::deleteRim(int getRidOfFactorizationData)
4470{
4471#ifdef CLP_USER_DRIVEN
4472  eventHandler_->event(ClpEventHandler::beforeDeleteRim);
4473#endif
4474  // Just possible empty problem
4475  int numberRows = numberRows_;
4476  int numberColumns = numberColumns_;
4477  if (!numberRows || !numberColumns) {
4478    numberRows = 0;
4479    if (objective_->type() < 2)
4480      numberColumns = 0;
4481  }
4482  int i;
4483  if (problemStatus_ != 1 && problemStatus_ != 2) {
4484    delete[] ray_;
4485    ray_ = NULL;
4486  }
4487  // set upperOut_ to furthest away from bound so can use in dual for dualBound_
4488  upperOut_ = 1.0;
4489#if 0
4490     {
4491          int nBad = 0;
4492          for (i = 0; i < numberColumns; i++) {
4493               if (lower_[i] == upper_[i] && getColumnStatus(i) == basic)
4494                    nBad++;
4495          }
4496          if (nBad)
4497               printf("yy %d basic fixed\n", nBad);
4498     }
4499#endif
4500  if ((moreSpecialOptions_ & 4194304) != 0) {
4501    // preset tolerances were changed
4502    moreSpecialOptions_ &= ~4194304;
4503    primalTolerance_ = 1.0e-7;
4504    dblParam_[ClpPrimalTolerance] = primalTolerance_;
4505    dualTolerance_ = 1.0e-7;
4506    dblParam_[ClpDualTolerance] = dualTolerance_;
4507  }
4508  // ray may be null if in branch and bound
4509  if (rowScale_ && solution_) {
4510    // Collect infeasibilities
4511    int numberPrimalScaled = 0;
4512    int numberPrimalUnscaled = 0;
4513    int numberDualScaled = 0;
4514    int numberDualUnscaled = 0;
4515    double scaleC = 1.0 / objectiveScale_;
4516    double scaleR = 1.0 / rhsScale_;
4517    const double *inverseScale = inverseColumnScale_;
4518    for (i = 0; i < numberColumns; i++) {
4519      double scaleFactor = columnScale_[i];
4520      double valueScaled = columnActivityWork_[i];
4521      double lowerScaled = columnLowerWork_[i];
4522      double upperScaled = columnUpperWork_[i];
4523      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4524        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4525          numberPrimalScaled++;
4526        else
4527          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4528      }
4529      columnActivity_[i] = valueScaled * scaleFactor * scaleR;
4530      double value = columnActivity_[i];
4531      if (value < columnLower_[i] - primalTolerance_)
4532        numberPrimalUnscaled++;
4533      else if (value > columnUpper_[i] + primalTolerance_)
4534        numberPrimalUnscaled++;
4535      double valueScaledDual = reducedCostWork_[i];
4536      if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4537        numberDualScaled++;
4538      if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4539        numberDualScaled++;
4540      reducedCost_[i] = (valueScaledDual * scaleC) * inverseScale[i];
4541      double valueDual = reducedCost_[i];
4542      if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4543        numberDualUnscaled++;
4544      if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4545        numberDualUnscaled++;
4546    }
4547    inverseScale = inverseRowScale_;
4548    for (i = 0; i < numberRows; i++) {
4549      double scaleFactor = rowScale_[i];
4550      double valueScaled = rowActivityWork_[i];
4551      double lowerScaled = rowLowerWork_[i];
4552      double upperScaled = rowUpperWork_[i];
4553      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4554        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4555          numberPrimalScaled++;
4556        else
4557          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4558      }
4559      rowActivity_[i] = (valueScaled * scaleR) * inverseScale[i];
4560      double value = rowActivity_[i];
4561      if (value < rowLower_[i] - primalTolerance_)
4562        numberPrimalUnscaled++;
4563      else if (value > rowUpper_[i] + primalTolerance_)
4564        numberPrimalUnscaled++;
4565      double valueScaledDual = dual_[i] + rowObjectiveWork_[i];
4566      ;
4567      if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4568        numberDualScaled++;
4569      if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4570        numberDualScaled++;
4571      dual_[i] *= scaleFactor * scaleC;
4572      double valueDual = dual_[i];
4573      if (rowObjective_)
4574        valueDual += rowObjective_[i];
4575      if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4576        numberDualUnscaled++;
4577      if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4578        numberDualUnscaled++;
4579    }
4580    if (!problemStatus_ && !secondaryStatus_) {
4581      // See if we need to set secondary status
4582      if (numberPrimalUnscaled) {
4583        if (numberDualUnscaled)
4584          secondaryStatus_ = 4;
4585        else
4586          secondaryStatus_ = 2;
4587      } else {
4588        if (numberDualUnscaled)
4589          secondaryStatus_ = 3;
4590      }
4591    }
4592    if (problemStatus_ == 2 && ray_) {
4593      for (i = 0; i < numberColumns; i++) {
4594        ray_[i] *= columnScale_[i];
4595      }
4596    } else if (problemStatus_ == 1 && ray_) {
4597      for (i = 0; i < numberRows; i++) {
4598        ray_[i] *= rowScale_[i];
4599      }
4600    }
4601  } else if (rhsScale_ != 1.0 || objectiveScale_ != 1.0) {
4602    // Collect infeasibilities
4603    int numberPrimalScaled = 0;
4604    int numberPrimalUnscaled = 0;
4605    int numberDualScaled = 0;
4606    int numberDualUnscaled = 0;
4607    double scaleC = 1.0 / objectiveScale_;
4608    double scaleR = 1.0 / rhsScale_;
4609    for (i = 0; i < numberColumns; i++) {
4610      double valueScaled = columnActivityWork_[i];
4611      double lowerScaled = columnLowerWork_[i];
4612      double upperScaled = columnUpperWork_[i];
4613      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4614        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4615          numberPrimalScaled++;
4616        else
4617          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4618      }
4619      columnActivity_[i] = valueScaled * scaleR;
4620      double value = columnActivity_[i];
4621      if (value < columnLower_[i] - primalTolerance_)
4622        numberPrimalUnscaled++;
4623      else if (value > columnUpper_[i] + primalTolerance_)
4624        numberPrimalUnscaled++;
4625      double valueScaledDual = reducedCostWork_[i];
4626      if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4627        numberDualScaled++;
4628      if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4629        numberDualScaled++;
4630      reducedCost_[i] = valueScaledDual * scaleC;
4631      double valueDual = reducedCost_[i];
4632      if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4633        numberDualUnscaled++;
4634      if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4635        numberDualUnscaled++;
4636    }
4637    for (i = 0; i < numberRows; i++) {
4638      double valueScaled = rowActivityWork_[i];
4639      double lowerScaled = rowLowerWork_[i];
4640      double upperScaled = rowUpperWork_[i];
4641      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4642        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
4643          numberPrimalScaled++;
4644        else
4645          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4646      }
4647      rowActivity_[i] = valueScaled * scaleR;
4648      double value = rowActivity_[i];
4649      if (value < rowLower_[i] - primalTolerance_)
4650        numberPrimalUnscaled++;
4651      else if (value > rowUpper_[i] + primalTolerance_)
4652        numberPrimalUnscaled++;
4653      double valueScaledDual = dual_[i] + rowObjectiveWork_[i];
4654      ;
4655      if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4656        numberDualScaled++;
4657      if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4658        numberDualScaled++;
4659      dual_[i] *= scaleC;
4660      double valueDual = dual_[i];
4661      if (rowObjective_)
4662        valueDual += rowObjective_[i];
4663      if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4664        numberDualUnscaled++;
4665      if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4666        numberDualUnscaled++;
4667    }
4668    if (!problemStatus_ && !secondaryStatus_) {
4669      // See if we need to set secondary status
4670      if (numberPrimalUnscaled) {
4671        if (numberDualUnscaled)
4672          secondaryStatus_ = 4;
4673        else
4674          secondaryStatus_ = 2;
4675      } else {
4676        if (numberDualUnscaled)
4677          secondaryStatus_ = 3;
4678      }
4679    }
4680  } else {
4681    if (columnActivityWork_) {
4682      for (i = 0; i < numberColumns; i++) {
4683        double value = columnActivityWork_[i];
4684        double lower = columnLowerWork_[i];
4685        double upper = columnUpperWork_[i];
4686        if (lower > -1.0e20 || upper < 1.0e20) {
4687          if (value > lower && value < upper)
4688            upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4689        }
4690        columnActivity_[i] = columnActivityWork_[i];
4691        reducedCost_[i] = reducedCostWork_[i];
4692      }
4693      for (i = 0; i < numberRows; i++) {
4694        double value = rowActivityWork_[i];
4695        double lower = rowLowerWork_[i];
4696        double upper = rowUpperWork_[i];
4697        if (lower > -1.0e20 || upper < 1.0e20) {
4698          if (value > lower && value < upper)
4699            upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4700        }
4701        rowActivity_[i] = rowActivityWork_[i];
4702      }
4703    }
4704  }
4705  // switch off scalefactor if auto
4706  if (automaticScale_) {
4707    rhsScale_ = 1.0;
4708    objectiveScale_ = 1.0;
4709  }
4710  if (optimizationDirection_ != 1.0) {
4711    // and modify all dual signs
4712    for (i = 0; i < numberColumns; i++)
4713      reducedCost_[i] *= optimizationDirection_;
4714    for (i = 0; i < numberRows; i++)
4715      dual_[i] *= optimizationDirection_;
4716  }
4717  // scaling may have been turned off
4718  scalingFlag_ = abs(scalingFlag_);
4719  if (getRidOfFactorizationData > 0) {
4720    gutsOfDelete(getRidOfFactorizationData + 1);
4721  } else {
4722    // at least get rid of nonLinearCost_
4723    delete nonLinearCost_;
4724    nonLinearCost_ = NULL;
4725  }
4726  if (!rowObjective_ && problemStatus_ == 0 && objective_->type() == 1 && numberRows && numberColumns) {
4727    // Redo objective value
4728    double objectiveValue = 0.0;
4729    const double *cost = objective();
4730    for (int i = 0; i < numberColumns; i++) {
4731      double value = columnActivity_[i];
4732      objectiveValue += value * cost[i];
4733    }
4734    //if (fabs(objectiveValue_ -objectiveValue*optimizationDirection())>1.0e-5)
4735    //printf("old obj %g new %g\n",objectiveValue_, objectiveValue*optimizationDirection());
4736    objectiveValue_ = objectiveValue * optimizationDirection();
4737  }
4738  // get rid of data
4739  matrix_->generalExpanded(this, 13, scalingFlag_);
4740}
4741void ClpSimplex::setDualBound(double value)
4742{
4743  if (value > 0.0)
4744    dualBound_ = value;
4745}
4746void ClpSimplex::setInfeasibilityCost(double value)
4747{
4748  if (value > 0.0)
4749    infeasibilityCost_ = value;
4750}
4751void ClpSimplex::setNumberRefinements(int value)
4752{
4753  if (value >= 0 && value < 10)
4754    numberRefinements_ = value;
4755}
4756// Sets row pivot choice algorithm in dual
4757void ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot &choice)
4758{
4759  delete dualRowPivot_;
4760  dualRowPivot_ = choice.clone(true);
4761  dualRowPivot_->setModel(this);
4762}
4763// Sets column pivot choice algorithm in primal
4764void ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot &choice)
4765{
4766  delete primalColumnPivot_;
4767  primalColumnPivot_ = choice.clone(true);
4768  primalColumnPivot_->setModel(this);
4769}
4770void ClpSimplex::setFactorization(ClpFactorization &factorization)
4771{
4772  if (factorization_)
4773    factorization_->setFactorization(factorization);
4774  else
4775    factorization_ = new ClpFactorization(factorization,
4776      numberRows_);
4777}
4778
4779// Swaps factorization
4780ClpFactorization *
4781ClpSimplex::swapFactorization(ClpFactorization *factorization)
4782{
4783  ClpFactorization *swap = factorization_;
4784  factorization_ = factorization;
4785  return swap;
4786}
4787// Copies in factorization to existing one
4788void ClpSimplex::copyFactorization(ClpFactorization &factorization)
4789{
4790  *factorization_ = factorization;
4791}
4792/* Perturbation:
4793   -50 to +50 - perturb by this power of ten (-6 sounds good)
4794   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
4795   101 - we are perturbed
4796   102 - don't try perturbing again
4797   default is 100
4798*/
4799void ClpSimplex::setPerturbation(int value)
4800{
4801  if (value <= 102 && value >= -1000) {
4802    perturbation_ = value;
4803  }
4804}
4805// Sparsity on or off
4806bool ClpSimplex::sparseFactorization() const
4807{
4808  return factorization_->sparseThreshold() != 0;
4809}
4810void ClpSimplex::setSparseFactorization(bool value)
4811{
4812  if (value) {
4813    if (!factorization_->sparseThreshold())
4814      factorization_->goSparse();
4815  } else {
4816    factorization_->sparseThreshold(0);
4817  }
4818}
4819void checkCorrect(ClpSimplex * /*model*/, int iRow,
4820  const double *element, const int *rowStart, const int *rowLength,
4821  const int *column,
4822  const double *columnLower_, const double *columnUpper_,
4823  int /*infiniteUpperC*/,
4824  int /*infiniteLowerC*/,
4825  double &maximumUpC,
4826  double &maximumDownC)
4827{
4828  int infiniteUpper = 0;
4829  int infiniteLower = 0;
4830  double maximumUp = 0.0;
4831  double maximumDown = 0.0;
4832  CoinBigIndex rStart = rowStart[iRow];
4833  CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4834  CoinBigIndex j;
4835  double large = 1.0e15;
4836  int iColumn;
4837  // Compute possible lower and upper ranges
4838
4839  for (j = rStart; j < rEnd; ++j) {
4840    double value = element[j];
4841    iColumn = column[j];
4842    if (value > 0.0) {
4843      if (columnUpper_[iColumn] >= large) {
4844        ++infiniteUpper;
4845      } else {
4846        maximumUp += columnUpper_[iColumn] * value;
4847      }
4848      if (columnLower_[iColumn] <= -large) {
4849        ++infiniteLower;
4850      } else {
4851        maximumDown += columnLower_[iColumn] * value;
4852      }
4853    } else if (value < 0.0) {
4854      if (columnUpper_[iColumn] >= large) {
4855        ++infiniteLower;
4856      } else {
4857        maximumDown += columnUpper_[iColumn] * value;
4858      }
4859      if (columnLower_[iColumn] <= -large) {
4860        ++infiniteUpper;
4861      } else {
4862        maximumUp += columnLower_[iColumn] * value;
4863      }
4864    }
4865  }
4866  //assert (infiniteLowerC==infiniteLower);
4867  //assert (infiniteUpperC==infiniteUpper);
4868  if (fabs(maximumUp - maximumUpC) > 1.0e-12 * CoinMax(fabs(maximumUp), fabs(maximumUpC)))
4869    COIN_DETAIL_PRINT(printf("row %d comp up %g, true up %g\n", iRow,
4870      maximumUpC, maximumUp));
4871  if (fabs(maximumDown - maximumDownC) > 1.0e-12 * CoinMax(fabs(maximumDown), fabs(maximumDownC)))
4872    COIN_DETAIL_PRINT(printf("row %d comp down %g, true down %g\n", iRow,
4873      maximumDownC, maximumDown));
4874  maximumUpC = maximumUp;
4875  maximumDownC = maximumDown;
4876}
4877
4878/* Tightens primal bounds to make dual faster.  Unless
4879   fixed, bounds are slightly looser than they could be.
4880   This is to make dual go faster and is probably not needed
4881   with a presolve.  Returns non-zero if problem infeasible
4882
4883   Fudge for branch and bound - put bounds on columns of factor *
4884   largest value (at continuous) - should improve stability
4885   in branch and bound on infeasible branches (0.0 is off)
4886*/
4887int ClpSimplex::tightenPrimalBounds(double factor, int doTight, bool tightIntegers)
4888{
4889
4890  // Get a row copy in standard format
4891  CoinPackedMatrix copy;
4892  copy.setExtraGap(0.0);
4893  copy.setExtraMajor(0.0);
4894  copy.reverseOrderedCopyOf(*matrix());
4895  // Matrix may have been created so get rid of it
4896  matrix_->releasePackedMatrix();
4897  // get matrix data pointers
4898  const int *column = copy.getIndices();
4899  const CoinBigIndex *rowStart = copy.getVectorStarts();
4900  const int *rowLength = copy.getVectorLengths();
4901  const double *element = copy.getElements();
4902  int numberChanged = 1, iPass = 0;
4903  double large = largeValue(); // treat bounds > this as infinite
4904#ifndef NDEBUG
4905  double large2 = 1.0e10 * large;
4906#endif
4907  int numberInfeasible = 0;
4908  int totalTightened = 0;
4909
4910  double tolerance = primalTolerance();
4911
4912  // Save column bounds
4913  double *saveLower = new double[numberColumns_];
4914  CoinMemcpyN(columnLower_, numberColumns_, saveLower);
4915  double *saveUpper = new double[numberColumns_];
4916  CoinMemcpyN(columnUpper_, numberColumns_, saveUpper);
4917
4918  int iRow, iColumn;
4919  // If wanted compute a reasonable dualBound_
4920  if (factor == COIN_DBL_MAX) {
4921    factor = 0.0;
4922    if (dualBound_ == 1.0e10) {
4923      // get largest scaled away from bound
4924      double largest = 1.0e-12;
4925      double largestScaled = 1.0e-12;
4926      int iRow;
4927      for (iRow = 0; iRow < numberRows_; iRow++) {
4928        double value = rowActivity_[iRow];
4929        double above = value - rowLower_[iRow];
4930        double below = rowUpper_[iRow] - value;
4931        if (above < 1.0e12) {
4932          largest = CoinMax(largest, above);
4933        }
4934        if (below < 1.0e12) {
4935          largest = CoinMax(largest, below);
4936        }
4937        if (rowScale_) {
4938          double multiplier = rowScale_[iRow];
4939          above *= multiplier;
4940          below *= multiplier;
4941        }
4942        if (above < 1.0e12) {
4943          largestScaled = CoinMax(largestScaled, above);
4944        }
4945        if (below < 1.0e12) {
4946          largestScaled = CoinMax(largestScaled, below);
4947        }
4948      }
4949
4950      int iColumn;
4951      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4952        double value = columnActivity_[iColumn];
4953        double above = value - columnLower_[iColumn];
4954        double below = columnUpper_[iColumn] - value;
4955        if (above < 1.0e12) {
4956          largest = CoinMax(largest, above);
4957        }
4958        if (below < 1.0e12) {
4959          largest = CoinMax(largest, below);
4960        }
4961        if (columnScale_) {
4962          double multiplier = 1.0 / columnScale_[iColumn];
4963          above *= multiplier;
4964          below *= multiplier;
4965        }
4966        if (above < 1.0e12) {
4967          largestScaled = CoinMax(largestScaled, above);
4968        }
4969        if (below < 1.0e12) {
4970          largestScaled = CoinMax(largestScaled, below);
4971        }
4972      }
4973      std::cout << "Largest (scaled) away from bound " << largestScaled
4974                << " unscaled " << largest << std::endl;
4975      dualBound_ = CoinMax(1.0001e7, CoinMin(100.0 * largest, 1.00001e10));
4976    }
4977  }
4978
4979  // If wanted - tighten column bounds using solution
4980  if (factor) {
4981    double largest = 0.0;
4982    if (factor > 0.0) {
4983      assert(factor > 1.0);
4984      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4985        if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
4986          largest = CoinMax(largest, fabs(columnActivity_[iColumn]));
4987        }
4988      }
4989      largest *= factor;
4990    } else {
4991      // absolute
4992      largest = -factor;
4993    }
4994    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4995      if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
4996        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn], largest);
4997        columnLower_[iColumn] = CoinMax(columnLower_[iColumn], -largest);
4998      }
4999    }
5000  }
5001#define MAXPASS 10
5002
5003  // Loop round seeing if we can tighten bounds
5004  // Would be faster to have a stack of possible rows
5005  // and we put altered rows back on stack
5006  int numberCheck = -1;
5007  while (numberChanged > numberCheck) {
5008
5009    numberChanged = 0; // Bounds tightened this pass
5010
5011    if (iPass == MAXPASS)
5012      break;
5013    iPass++;
5014
5015    for (iRow = 0; iRow < numberRows_; iRow++) {
5016
5017      if (rowLower_[iRow] > -large || rowUpper_[iRow] < large) {
5018
5019        // possible row
5020        int infiniteUpper = 0;
5021        int infiniteLower = 0;
5022        double maximumUp = 0.0;
5023        double maximumDown = 0.0;
5024        double newBound;
5025        CoinBigIndex rStart = rowStart[iRow];
5026        CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
5027        CoinBigIndex j;
5028        // Compute possible lower and upper ranges
5029
5030        for (j = rStart; j < rEnd; ++j) {
5031          double value = element[j];
5032          iColumn = column[j];
5033          if (value > 0.0) {
5034            if (columnUpper_[iColumn] >= large) {
5035              ++infiniteUpper;
5036            } else {
5037              maximumUp += columnUpper_[iColumn] * value;
5038            }
5039            if (columnLower_[iColumn] <= -large) {
5040              ++infiniteLower;
5041            } else {
5042              maximumDown += columnLower_[iColumn] * value;
5043            }
5044          } else if (value < 0.0) {
5045            if (columnUpper_[iColumn] >= large) {
5046              ++infiniteLower;
5047            } else {
5048              maximumDown += columnUpper_[iColumn] * value;
5049            }
5050            if (columnLower_[iColumn] <= -large) {
5051              ++infiniteUpper;
5052            } else {
5053              maximumUp += columnLower_[iColumn] * value;
5054            }
5055          }
5056        }
5057        // Build in a margin of error
5058        maximumUp += 1.0e-8 * fabs(maximumUp);
5059        maximumDown -= 1.0e-8 * fabs(maximumDown);
5060        double maxUp = maximumUp + infiniteUpper * 1.0e31;
5061        double maxDown = maximumDown - infiniteLower * 1.0e31;
5062        if (maxUp <= rowUpper_[iRow] + tolerance && maxDown >= rowLower_[iRow] - tolerance) {
5063
5064          // Row is redundant - make totally free
5065          // NO - can't do this for postsolve
5066          // rowLower_[iRow]=-COIN_DBL_MAX;
5067          // rowUpper_[iRow]=COIN_DBL_MAX;
5068          //printf("Redundant row in presolveX %d\n",iRow);
5069
5070        } else {
5071          if (maxUp < rowLower_[iRow] - 100.0 * tolerance || maxDown > rowUpper_[iRow] + 100.0 * tolerance) {
5072            // problem is infeasible - exit at once
5073            numberInfeasible++;
5074            break;
5075          }
5076          double lower = rowLower_[iRow];
5077          double upper = rowUpper_[iRow];
5078          for (j = rStart; j < rEnd; ++j) {
5079            double value = element[j];
5080            iColumn = column[j];
5081            double nowLower = columnLower_[iColumn];
5082            double nowUpper = columnUpper_[iColumn];
5083            if (value > 0.0) {
5084              // positive value
5085              if (lower > -large) {
5086                if (!infiniteUpper) {
5087                  assert(nowUpper < large2);
5088                  newBound = nowUpper + (lower - maximumUp) / value;
5089                  // relax if original was large
5090                  if (fabs(maximumUp) > 1.0e8)
5091                    newBound -= 1.0e-12 * fabs(maximumUp);
5092                } else if (infiniteUpper == 1 && nowUpper > large) {
5093                  newBound = (lower - maximumUp) / value;
5094                  // relax if original was large
5095                  if (fabs(maximumUp) > 1.0e8)
5096                    newBound -= 1.0e-12 * fabs(maximumUp);
5097                } else {
5098                  newBound = -COIN_DBL_MAX;
5099                }
5100                if (newBound > nowLower + 1.0e-12 && newBound > -large) {
5101                  // Tighten the lower bound
5102                  numberChanged++;
5103                  // check infeasible (relaxed)
5104                  if (nowUpper < newBound) {
5105                    if (nowUpper - newBound < -100.0 * tolerance)
5106                      numberInfeasible++;
5107                    else
5108                      newBound = nowUpper;
5109                  }
5110                  columnLower_[iColumn] = newBound;
5111                  // adjust
5112                  double now;
5113                  if (nowLower < -large) {
5114                    now = 0.0;
5115                    infiniteLower--;
5116                  } else {
5117                    now = nowLower;
5118                  }
5119                  maximumDown += (newBound - now) * value;
5120                  nowLower = newBound;
5121#ifdef DEBUG
5122                  checkCorrect(this, iRow,
5123                    element, rowStart, rowLength,
5124                    column,
5125                    columnLower_, columnUpper_,
5126                    infiniteUpper,
5127                    infiniteLower,
5128                    maximumUp,
5129                    maximumDown);
5130#endif
5131                }
5132              }
5133              if (upper < large) {
5134                if (!infiniteLower) {
5135                  assert(nowLower > -large2);
5136                  newBound = nowLower + (upper - maximumDown) / value;
5137                  // relax if original was large
5138                  if (fabs(maximumDown) > 1.0e8)
5139                    newBound += 1.0e-12 * fabs(maximumDown);
5140                } else if (infiniteLower == 1 && nowLower < -large) {
5141                  newBound = (upper - maximumDown) / value;
5142                  // relax if original was large
5143                  if (fabs(maximumDown) > 1.0e8)
5144                    newBound += 1.0e-12 * fabs(maximumDown);
5145                } else {
5146                  newBound = COIN_DBL_MAX;
5147                }
5148                if (newBound < nowUpper - 1.0e-12 && newBound < large) {
5149                  // Tighten the upper bound
5150                  numberChanged++;
5151                  // check infeasible (relaxed)
5152                  if (nowLower > newBound) {
5153                    if (newBound - nowLower < -100.0 * tolerance)
5154                      numberInfeasible++;
5155                    else
5156                      newBound = nowLower;
5157                  }
5158                  columnUpper_[iColumn] = newBound;
5159                  // adjust
5160                  double now;
5161                  if (nowUpper > large) {
5162                    now = 0.0;
5163                    infiniteUpper--;
5164                  } else {
5165                    now = nowUpper;
5166                  }
5167                  maximumUp += (newBound - now) * value;
5168                  nowUpper = newBound;
5169#ifdef DEBUG
5170                  checkCorrect(this, iRow,
5171                    element, rowStart, rowLength,
5172                    column,
5173                    columnLower_, columnUpper_,
5174                    infiniteUpper,
5175                    infiniteLower,
5176                    maximumUp,
5177                    maximumDown);
5178#endif
5179                }
5180              }
5181            } else {
5182              // negative value
5183              if (lower > -large) {
5184                if (!infiniteUpper) {
5185                  assert(nowLower < large2);
5186                  newBound = nowLower + (lower - maximumUp) / value;
5187                  // relax if original was large
5188                  if (fabs(maximumUp) > 1.0e8)
5189                    newBound += 1.0e-12 * fabs(maximumUp);
5190                } else if (infiniteUpper == 1 && nowLower < -large) {
5191                  newBound = (lower - maximumUp) / value;
5192                  // relax if original was large
5193                  if (fabs(maximumUp) > 1.0e8)
5194                    newBound += 1.0e-12 * fabs(maximumUp);
5195                } else {
5196                  newBound = COIN_DBL_MAX;
5197                }
5198                if (newBound < nowUpper - 1.0e-12 && newBound < large) {
5199                  // Tighten the upper bound
5200                  numberChanged++;
5201                  // check infeasible (relaxed)
5202                  if (nowLower > newBound) {
5203                    if (newBound - nowLower < -100.0 * tolerance)
5204                      numberInfeasible++;
5205                    else
5206                      newBound = nowLower;
5207                  }
5208                  columnUpper_[iColumn] = newBound;
5209                  // adjust
5210                  double now;
5211                  if (nowUpper > large) {
5212                    now = 0.0;
5213                    infiniteLower--;
5214                  } else {
5215                    now = nowUpper;
5216                  }
5217                  maximumDown += (newBound - now) * value;
5218                  nowUpper = newBound;
5219#ifdef DEBUG
5220                  checkCorrect(this, iRow,
5221                    element, rowStart, rowLength,
5222                    column,
5223                    columnLower_, columnUpper_,
5224                    infiniteUpper,
5225                    infiniteLower,
5226                    maximumUp,
5227                    maximumDown);
5228#endif
5229                }
5230              }
5231              if (upper < large) {
5232                if (!infiniteLower) {
5233                  assert(nowUpper < large2);
5234                  newBound = nowUpper + (upper - maximumDown) / value;
5235                  // relax if original was large
5236                  if (fabs(maximumDown) > 1.0e8)
5237                    newBound -= 1.0e-12 * fabs(maximumDown);
5238                } else if (infiniteLower == 1 && nowUpper > large) {
5239                  newBound = (upper - maximumDown) / value;
5240                  // relax if original was large
5241                  if (fabs(maximumDown) > 1.0e8)
5242                    newBound -= 1.0e-12 * fabs(maximumDown);
5243                } else {
5244                  newBound = -COIN_DBL_MAX;
5245                }
5246                if (newBound > nowLower + 1.0e-12 && newBound > -large) {
5247                  // Tighten the lower bound
5248                  numberChanged++;
5249                  // check infeasible (relaxed)
5250                  if (nowUpper < newBound) {
5251                    if (nowUpper - newBound < -100.0 * tolerance)
5252                      numberInfeasible++;
5253                    else
5254                      newBound = nowUpper;
5255                  }
5256                  columnLower_[iColumn] = newBound;
5257                  // adjust
5258                  double now;
5259                  if (nowLower < -large) {
5260                    now = 0.0;
5261                    infiniteUpper--;
5262                  } else {
5263                    now = nowLower;
5264                  }
5265                  maximumUp += (newBound - now) * value;
5266                  nowLower = newBound;
5267#ifdef DEBUG
5268                  checkCorrect(this, iRow,
5269                    element, rowStart, rowLength,
5270                    column,
5271                    columnLower_, columnUpper_,
5272                    infiniteUpper,
5273                    infiniteLower,
5274                    maximumUp,
5275                    maximumDown);
5276#endif
5277                }
5278              }
5279            }
5280          }
5281        }
5282      }
5283    }
5284    totalTightened += numberChanged;
5285    if (iPass == 1)
5286      numberCheck = numberChanged >> 4;
5287    if (numberInfeasible)
5288      break;
5289  }
5290  if (!numberInfeasible) {
5291    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN, messages_)
5292      << totalTightened
5293      << CoinMessageEol;
5294    // Set bounds slightly loose
5295    double useTolerance = 1.0e-3;
5296    if (doTight > 0) {
5297      if (doTight > 10) {
5298        useTolerance = 0.0;
5299      } else {
5300        while (doTight) {
5301          useTolerance *= 0.1;
5302          doTight--;
5303        }
5304      }
5305    }
5306    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5307      if (saveUpper[iColumn] > saveLower[iColumn] + tolerance) {
5308        // Make large bounds stay infinite
5309        if (saveUpper[iColumn] > 1.0e30 && columnUpper_[iColumn] > 1.0e10) {
5310          columnUpper_[iColumn] = COIN_DBL_MAX;
5311        }
5312        if (saveLower[iColumn] < -1.0e30 && columnLower_[iColumn] < -1.0e10) {
5313          columnLower_[iColumn] = -COIN_DBL_MAX;
5314        }
5315#ifdef KEEP_GOING_IF_FIXED
5316        double multiplier = 5.0e-3 * floor(100.0 * randomNumberGenerator_.randomDouble()) + 1.0;
5317        multiplier *= 100.0;
5318#else
5319        double multiplier = 100.0;
5320#endif
5321        if (columnUpper_[iColumn] - columnLower_[iColumn] < useTolerance + 1.0e-8) {
5322          // relax enough so will have correct dj
5323#if 1
5324          columnLower_[iColumn] = CoinMax(saveLower[iColumn],
5325            columnLower_[iColumn] - multiplier * useTolerance);
5326          columnUpper_[iColumn] = CoinMin(saveUpper[iColumn],
5327            columnUpper_[iColumn] + multiplier * useTolerance);
5328#else
5329          if (fabs(columnUpper_[iColumn]) < fabs(columnLower_[iColumn])) {
5330            if (columnUpper_[iColumn] - multiplier * useTolerance > saveLower[iColumn]) {
5331              columnLower_[iColumn] = columnUpper_[iColumn] - multiplier * useTolerance;
5332            } else {
5333              columnLower_[iColumn] = saveLower[iColumn];
5334              columnUpper_[iColumn] = CoinMin(saveUpper[iColumn],
5335                saveLower[iColumn] + multiplier * useTolerance);
5336            }
5337          } else {
5338            if (columnLower_[iColumn] + multiplier * useTolerance < saveUpper[iColumn]) {
5339              columnUpper_[iColumn] = columnLower_[iColumn] + multiplier * useTolerance;
5340            } else {
5341              columnUpper_[iColumn] = saveUpper[iColumn];
5342              columnLower_[iColumn] = CoinMax(saveLower[iColumn],
5343                saveUpper[iColumn] - multiplier * useTolerance);
5344            }
5345          }
5346#endif
5347        } else {
5348          if (columnUpper_[iColumn] < saveUpper[iColumn]) {
5349            // relax a bit
5350            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn] + multiplier * useTolerance,
5351              saveUpper[iColumn]);
5352          }
5353          if (columnLower_[iColumn] > saveLower[iColumn]) {
5354            // relax a bit
5355            columnLower_[iColumn] = CoinMax(columnLower_[iColumn] - multiplier * useTolerance,
5356              saveLower[iColumn]);
5357          }
5358        }
5359      }
5360    }
5361    if (tightIntegers && integerType_) {
5362      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5363        if (integerType_[iColumn]) {
5364          double value;
5365          value = floor(columnLower_[iColumn] + 0.5);
5366          if (fabs(value - columnLower_[iColumn]) > primalTolerance_)
5367            value = ceil(columnLower_[iColumn]);
5368          columnLower_[iColumn] = value;
5369          value = floor(columnUpper_[iColumn] + 0.5);
5370          if (fabs(value - columnUpper_[iColumn]) > primalTolerance_)
5371            value = floor(columnUpper_[iColumn]);
5372          columnUpper_[iColumn] = value;
5373          if (columnLower_[iColumn] > columnUpper_[iColumn])
5374            numberInfeasible++;
5375        }
5376      }
5377      if (numberInfeasible) {
5378        handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
5379          << numberInfeasible
5380          << CoinMessageEol;
5381        // restore column bounds
5382        CoinMemcpyN(saveLower, numberColumns_, columnLower_);
5383        CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
5384      }
5385    }
5386  } else {
5387    handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
5388      << numberInfeasible
5389      << CoinMessageEol;
5390    // restore column bounds
5391    CoinMemcpyN(saveLower, numberColumns_, columnLower_);
5392    CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
5393  }
5394  delete[] saveLower;
5395  delete[] saveUpper;
5396  return (numberInfeasible);
5397}
5398//#define SAVE_AND_RESTORE
5399// dual
5400#include "ClpSimplexDual.hpp"
5401#include "ClpSimplexPrimal.hpp"
5402#ifndef SAVE_AND_RESTORE
5403int ClpSimplex::dual(int ifValuesPass, int startFinishOptions)
5404#else
5405int ClpSimplex::dual(int ifValuesPass, int startFinishOptions)
5406{
5407  // May be empty problem
5408  if (numberRow