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

Last change on this file since 2453 was 2453, checked in by stefan, 2 months ago

sync with trunk r2451

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