source: stable/1.13/Clp/src/ClpSimplex.cpp @ 1649

Last change on this file since 1649 was 1649, checked in by stefan, 8 years ago

sync with trunk rev1648, except for half of chgset1635 and except for chgset1643

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