source: stable/1.15/Clp/src/ClpSimplex.cpp @ 1959

Last change on this file since 1959 was 1959, checked in by stefan, 6 years ago

merge r1945, r1947, and r1950..r1958 from trunk

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