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

Last change on this file since 2082 was 2082, checked in by forrest, 5 years ago

out extraneous printfs

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