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

Last change on this file since 1928 was 1928, checked in by stefan, 7 years ago

revert r1925 (so we can merge r1926) and sync with rev 1927

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 506.9 KB
Line 
1/* $Id: ClpSimplex.cpp 1928 2013-04-06 12:54:16Z 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                    // try automatic scaling
3316                    double smallestObj = 1.0e100;
3317                    double largestObj = 0.0;
3318                    double largestRhs = 0.0;
3319                    const double * obj = objective();
3320                    for (i = 0; i < numberColumns_; i++) {
3321                         double value = fabs(obj[i]);
3322                         value *= columnScale_[i];
3323                         if (value && columnLower_[i] != columnUpper_[i]) {
3324                              smallestObj = CoinMin(smallestObj, value);
3325                              largestObj = CoinMax(largestObj, value);
3326                         }
3327                         if (columnLower_[i] > 0.0 || columnUpper_[i] < 0.0) {
3328                              double scale = 1.0 * inverseColumnScale_[i];
3329                              //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3330                              if (columnLower_[i] > 0)
3331                                   largestRhs = CoinMax(largestRhs, columnLower_[i] * scale);
3332                              if (columnUpper_[i] < 0.0)
3333                                   largestRhs = CoinMax(largestRhs, -columnUpper_[i] * scale);
3334                         }
3335                    }
3336                    for (i = 0; i < numberRows_; i++) {
3337                         if (rowLower_[i] > 0.0 || rowUpper_[i] < 0.0) {
3338                              double scale = rowScale_[i];
3339                              //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3340                              if (rowLower_[i] > 0)
3341                                   largestRhs = CoinMax(largestRhs, rowLower_[i] * scale);
3342                              if (rowUpper_[i] < 0.0)
3343                                   largestRhs = CoinMax(largestRhs, -rowUpper_[i] * scale);
3344                         }
3345                    }
3346                    COIN_DETAIL_PRINT(printf("small obj %g, large %g - rhs %g\n", smallestObj, largestObj, largestRhs));
3347                    bool scalingDone = false;
3348                    // look at element range
3349                    double smallestNegative;
3350                    double largestNegative;
3351                    double smallestPositive;
3352                    double largestPositive;
3353                    matrix_->rangeOfElements(smallestNegative, largestNegative,
3354                                             smallestPositive, largestPositive);
3355                    smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
3356                    largestPositive = CoinMax(fabs(largestNegative), largestPositive);
3357                    if (largestObj) {
3358                         double ratio = largestObj / smallestObj;
3359                         double scale = 1.0;
3360                         if (ratio < 1.0e8) {
3361                              // reasonable
3362                              if (smallestObj < 1.0e-4) {
3363                                   // may as well scale up
3364                                   scalingDone = true;
3365                                   scale = 1.0e-3 / smallestObj;
3366                              } else if (largestObj < 1.0e6 || (algorithm_ > 0 && largestObj < 1.0e-4 * infeasibilityCost_)) {
3367                                   //done=true;
3368                              } else {
3369                                   scalingDone = true;
3370                                   if (algorithm_ < 0) {
3371                                        scale = 1.0e6 / largestObj;
3372                                   } else {
3373                                        scale = CoinMax(1.0e6, 1.0e-4 * infeasibilityCost_) / largestObj;
3374                                   }
3375                              }
3376                         } else if (ratio < 1.0e12) {
3377                              // not so good
3378                              if (smallestObj < 1.0e-7) {
3379                                   // may as well scale up
3380                                   scalingDone = true;
3381                                   scale = 1.0e-6 / smallestObj;
3382                              } else if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3383                                   //done=true;
3384                              } else {
3385                                   scalingDone = true;
3386                                   if (algorithm_ < 0) {
3387                                        scale = 1.0e7 / largestObj;
3388                                   } else {
3389                                        scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3390                                   }
3391                              }
3392                         } else {
3393                              // Really nasty problem
3394                              if (smallestObj < 1.0e-8) {
3395                                   // may as well scale up
3396                                   scalingDone = true;
3397                                   scale = 1.0e-7 / smallestObj;
3398                                   largestObj *= scale;
3399                              }
3400                              if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3401                                   //done=true;
3402                              } else {
3403                                   scalingDone = true;
3404                                   if (algorithm_ < 0) {
3405                                        scale = 1.0e7 / largestObj;
3406                                   } else {
3407                                        scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3408                                   }
3409                              }
3410                         }
3411                         objectiveScale_ = scale;
3412                    }
3413                    if (largestRhs > 1.0e12) {
3414                         scalingDone = true;
3415                         rhsScale_ = 1.0e9 / largestRhs;
3416                    } else if (largestPositive > 1.0e-14 * smallestPositive && largestRhs > 1.0e6) {
3417                         scalingDone = true;
3418                         rhsScale_ = 1.0e6 / largestRhs;
3419                    } else {
3420                         rhsScale_ = 1.0;
3421                    }
3422                    if (scalingDone) {
3423                         handler_->message(CLP_RIM_SCALE, messages_)
3424                                   << objectiveScale_ << rhsScale_
3425                                   << CoinMessageEol;
3426                    }
3427               }
3428          } else if (makeRowCopy && scalingFlag_ > 0 && !rowCopyIsScaled) {
3429               matrix_->scaleRowCopy(this);
3430          }
3431          if (rowScale_ && !savedRowScale_) {
3432               inverseRowScale_ = rowScale_ + numberRows2;
3433               inverseColumnScale_ = columnScale_ + numberColumns_;
3434          }
3435          // See if we can try for faster row copy
3436          if (makeRowCopy && !oldMatrix) {
3437               ClpPackedMatrix* clpMatrix =
3438                    dynamic_cast< ClpPackedMatrix*>(matrix_);
3439               if (clpMatrix && numberThreads_)
3440                    clpMatrix->specialRowCopy(this, rowCopy_);
3441               if (clpMatrix)
3442                    clpMatrix->specialColumnCopy(this);
3443          }
3444     }
3445     if (what == 63) {
3446#if 0
3447          {
3448               x_gaps[0]++;
3449               ClpPackedMatrix* clpMatrix =
3450               dynamic_cast< ClpPackedMatrix*>(matrix_);
3451               if (clpMatrix) {
3452                    if (!clpMatrix->getPackedMatrix()->hasGaps())
3453                         x_gaps[1]++;
3454                    if ((clpMatrix->flags() & 2) == 0)
3455                         x_gaps[3]++;
3456               } else {
3457                    x_gaps[2]++;
3458               }
3459               if ((x_gaps[0] % 1000) == 0)
3460                    printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n",
3461                           x_gaps[0], x_gaps[1], x_gaps[2], x_gaps[3]);
3462          }
3463#endif
3464          if (newArrays && (specialOptions_ & 65536) == 0) {
3465               delete [] cost_;
3466               cost_ = new double[2*numberTotal];
3467               delete [] lower_;
3468               delete [] upper_;
3469               lower_ = new double[numberTotal];
3470               upper_ = new double[numberTotal];
3471               delete [] dj_;
3472               dj_ = new double[numberTotal];
3473               delete [] solution_;
3474               solution_ = new double[numberTotal];
3475          }
3476          reducedCostWork_ = dj_;
3477          rowReducedCost_ = dj_ + numberColumns_;
3478          columnActivityWork_ = solution_;
3479          rowActivityWork_ = solution_ + numberColumns_;
3480          objectiveWork_ = cost_;
3481          rowObjectiveWork_ = cost_ + numberColumns_;
3482          rowLowerWork_ = lower_ + numberColumns_;
3483          columnLowerWork_ = lower_;
3484          rowUpperWork_ = upper_ + numberColumns_;
3485          columnUpperWork_ = upper_;
3486     }
3487     if ((what & 4) != 0) {
3488          double direction = optimizationDirection_ * objectiveScale_;
3489          const double * obj = objective();
3490          const double * rowScale = rowScale_;
3491          const double * columnScale = columnScale_;
3492          // and also scale by scale factors
3493          if (rowScale) {
3494               if (rowObjective_) {
3495                    for (i = 0; i < numberRows_; i++)
3496                         rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
3497               } else {
3498                    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3499               }
3500               // If scaled then do all columns later in one loop
3501               if (what != 63) {
3502                    for (i = 0; i < numberColumns_; i++) {
3503                         CoinAssert(fabs(obj[i]) < 1.0e25);
3504                         objectiveWork_[i] = obj[i] * direction * columnScale[i];
3505                    }
3506               }
3507          } else {
3508               if (rowObjective_) {
3509                    for (i = 0; i < numberRows_; i++)
3510                         rowObjectiveWork_[i] = rowObjective_[i] * direction;
3511               } else {
3512                    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3513               }
3514               for (i = 0; i < numberColumns_; i++) {
3515                    CoinAssert(fabs(obj[i]) < 1.0e25);
3516                    objectiveWork_[i] = obj[i] * direction;
3517               }
3518          }
3519     }
3520     if ((what & 1) != 0) {
3521          const double * rowScale = rowScale_;
3522          // clean up any mismatches on infinity
3523          // and fix any variables with tiny gaps
3524          double primalTolerance = dblParam_[ClpPrimalTolerance];
3525          if(rowScale) {
3526               // If scaled then do all columns later in one loop
3527               if (what != 63) {
3528                    const double * inverseScale = inverseColumnScale_;
3529                    for (i = 0; i < numberColumns_; i++) {
3530                         double multiplier = rhsScale_ * inverseScale[i];
3531                         double lowerValue = columnLower_[i];
3532                         double upperValue = columnUpper_[i];
3533                         if (lowerValue > -1.0e20) {
3534                              columnLowerWork_[i] = lowerValue * multiplier;
3535                              if (upperValue >= 1.0e20) {
3536                                   columnUpperWork_[i] = COIN_DBL_MAX;
3537                              } else {
3538                                   columnUpperWork_[i] = upperValue * multiplier;
3539                                   if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3540                                        if (columnLowerWork_[i] >= 0.0) {
3541                                             columnUpperWork_[i] = columnLowerWork_[i];
3542                                        } else if (columnUpperWork_[i] <= 0.0) {
3543                                             columnLowerWork_[i] = columnUpperWork_[i];
3544                                        } else {
3545                                             columnUpperWork_[i] = 0.0;
3546                                             columnLowerWork_[i] = 0.0;
3547                                        }
3548                                   }
3549                              }
3550                         } else if (upperValue < 1.0e20) {
3551                              columnLowerWork_[i] = -COIN_DBL_MAX;
3552                              columnUpperWork_[i] = upperValue * multiplier;
3553                         } else {
3554                              // free
3555                              columnLowerWork_[i] = -COIN_DBL_MAX;
3556                              columnUpperWork_[i] = COIN_DBL_MAX;
3557                         }
3558                    }
3559               }
3560               for (i = 0; i < numberRows_; i++) {
3561                    double multiplier = rhsScale_ * rowScale[i];
3562                    double lowerValue = rowLower_[i];
3563                    double upperValue = rowUpper_[i];
3564                    if (lowerValue > -1.0e20) {
3565                         rowLowerWork_[i] = lowerValue * multiplier;
3566                         if (upperValue >= 1.0e20) {
3567                              rowUpperWork_[i] = COIN_DBL_MAX;
3568                         } else {
3569                              rowUpperWork_[i] = upperValue * multiplier;
3570                              if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3571                                   if (rowLowerWork_[i] >= 0.0) {
3572                                        rowUpperWork_[i] = rowLowerWork_[i];
3573                                   } else if (rowUpperWork_[i] <= 0.0) {
3574                                        rowLowerWork_[i] = rowUpperWork_[i];
3575                                   } else {
3576                                        rowUpperWork_[i] = 0.0;
3577                                        rowLowerWork_[i] = 0.0;
3578                                   }
3579                              }
3580                         }
3581                    } else if (upperValue < 1.0e20) {
3582                         rowLowerWork_[i] = -COIN_DBL_MAX;
3583                         rowUpperWork_[i] = upperValue * multiplier;
3584                    } else {
3585                         // free
3586                         rowLowerWork_[i] = -COIN_DBL_MAX;
3587                         rowUpperWork_[i] = COIN_DBL_MAX;
3588                    }
3589               }
3590          } else if (rhsScale_ != 1.0) {
3591               for (i = 0; i < numberColumns_; i++) {
3592                    double lowerValue = columnLower_[i];
3593                    double upperValue = columnUpper_[i];
3594                    if (lowerValue > -1.0e20) {
3595                         columnLowerWork_[i] = lowerValue * rhsScale_;
3596                         if (upperValue >= 1.0e20) {
3597                              columnUpperWork_[i] = COIN_DBL_MAX;
3598                         } else {
3599                              columnUpperWork_[i] = upperValue * rhsScale_;
3600                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3601                                   if (columnLowerWork_[i] >= 0.0) {
3602                                        columnUpperWork_[i] = columnLowerWork_[i];
3603                                   } else if (columnUpperWork_[i] <= 0.0) {
3604                                        columnLowerWork_[i] = columnUpperWork_[i];
3605                                   } else {
3606                                        columnUpperWork_[i] = 0.0;
3607                                        columnLowerWork_[i] = 0.0;
3608                                   }
3609                              }
3610                         }
3611                    } else if (upperValue < 1.0e20) {
3612                         columnLowerWork_[i] = -COIN_DBL_MAX;
3613                         columnUpperWork_[i] = upperValue * rhsScale_;
3614                    } else {
3615                         // free
3616                         columnLowerWork_[i] = -COIN_DBL_MAX;
3617                         columnUpperWork_[i] = COIN_DBL_MAX;
3618                    }
3619               }
3620               for (i = 0; i < numberRows_; i++) {
3621                    double lowerValue = rowLower_[i];
3622                    double upperValue = rowUpper_[i];
3623                    if (lowerValue > -1.0e20) {
3624                         rowLowerWork_[i] = lowerValue * rhsScale_;
3625                         if (upperValue >= 1.0e20) {
3626                              rowUpperWork_[i] = COIN_DBL_MAX;
3627                         } else {
3628                              rowUpperWork_[i] = upperValue * rhsScale_;
3629                              if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3630                                   if (rowLowerWork_[i] >= 0.0) {
3631                                        rowUpperWork_[i] = rowLowerWork_[i];
3632                                   } else if (rowUpperWork_[i] <= 0.0) {
3633                                        rowLowerWork_[i] = rowUpperWork_[i];
3634                                   } else {
3635                                        rowUpperWork_[i] = 0.0;
3636                                        rowLowerWork_[i] = 0.0;
3637                                   }
3638                              }
3639                         }
3640                    } else if (upperValue < 1.0e20) {
3641                         rowLowerWork_[i] = -COIN_DBL_MAX;
3642                         rowUpperWork_[i] = upperValue * rhsScale_;
3643                    } else {
3644                         // free
3645                         rowLowerWork_[i] = -COIN_DBL_MAX;
3646                         rowUpperWork_[i] = COIN_DBL_MAX;
3647                    }
3648               }
3649          } else {
3650               for (i = 0; i < numberColumns_; i++) {
3651                    double lowerValue = columnLower_[i];
3652                    double upperValue = columnUpper_[i];
3653                    if (lowerValue > -1.0e20) {
3654                         columnLowerWork_[i] = lowerValue;
3655                         if (upperValue >= 1.0e20) {
3656                              columnUpperWork_[i] = COIN_DBL_MAX;
3657                         } else {
3658                              columnUpperWork_[i] = upperValue;
3659                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3660                                   if (columnLowerWork_[i] >= 0.0) {
3661                                        columnUpperWork_[i] = columnLowerWork_[i];
3662                                   } else if (columnUpperWork_[i] <= 0.0) {
3663                                        columnLowerWork_[i] = columnUpperWork_[i];
3664                                   } else {
3665                                        columnUpperWork_[i] = 0.0;
3666                                        columnLowerWork_[i] = 0.0;
3667                                   }
3668                              }
3669                         }
3670                    } else if (upperValue < 1.0e20) {
3671                         columnLowerWork_[i] = -COIN_DBL_MAX;
3672                         columnUpperWork_[i] = upperValue;
3673                    } else {
3674                         // free
3675                         columnLowerWork_[i] = -COIN_DBL_MAX;
3676                         columnUpperWork_[i] = COIN_DBL_MAX;
3677                    }
3678               }
3679               for (i = 0; i < numberRows_; i++) {
3680                    double lowerValue = rowLower_[i];
3681                    double upperValue = rowUpper_[i];
3682                    if (lowerValue > -1.0e20) {
3683                         rowLowerWork_[i] = lowerValue;
3684                         if (upperValue >= 1.0e20) {
3685                              rowUpperWork_[i] = COIN_DBL_MAX;
3686                         } else {
3687                              rowUpperWork_[i] = upperValue;
3688                              if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3689                                   if (rowLowerWork_[i] >= 0.0) {
3690                                        rowUpperWork_[i] = rowLowerWork_[i];
3691                                   } else if (rowUpperWork_[i] <= 0.0) {
3692                                        rowLowerWork_[i] = rowUpperWork_[i];
3693                                   } else {
3694                                        rowUpperWork_[i] = 0.0;
3695                                        rowLowerWork_[i] = 0.0;
3696                                   }
3697                              }
3698                         }
3699                    } else if (upperValue < 1.0e20) {
3700                         rowLowerWork_[i] = -COIN_DBL_MAX;
3701                         rowUpperWork_[i] = upperValue;
3702                    } else {
3703                         // free
3704                         rowLowerWork_[i] = -COIN_DBL_MAX;
3705                         rowUpperWork_[i] = COIN_DBL_MAX;
3706                    }
3707               }
3708          }
3709     }
3710     if (what == 63) {
3711          // move information to work arrays
3712          double direction = optimizationDirection_;
3713          // direction is actually scale out not scale in
3714          if (direction)
3715               direction = 1.0 / direction;
3716          if (direction != 1.0) {
3717               // reverse all dual signs
3718               for (i = 0; i < numberColumns_; i++)
3719                    reducedCost_[i] *= direction;
3720               for (i = 0; i < numberRows_; i++)
3721                    dual_[i] *= direction;
3722          }
3723          for (i = 0; i < numberRows_ + numberColumns_; i++) {
3724               setFakeBound(i, noFake);
3725          }
3726          if (rowScale_) {
3727               const double * obj = objective();
3728               double direction = optimizationDirection_ * objectiveScale_;
3729               // clean up any mismatches on infinity
3730               // and fix any variables with tiny gaps
3731               double primalTolerance = dblParam_[ClpPrimalTolerance];
3732               // on entry
3733               const double * inverseScale = inverseColumnScale_;
3734               for (i = 0; i < numberColumns_; i++) {
3735                    CoinAssert(fabs(obj[i]) < 1.0e25);
3736                    double scaleFactor = columnScale_[i];
3737                    double multiplier = rhsScale_ * inverseScale[i];
3738                    scaleFactor *= direction;
3739                    objectiveWork_[i] = obj[i] * scaleFactor;
3740                    reducedCostWork_[i] = reducedCost_[i] * scaleFactor;
3741                    double lowerValue = columnLower_[i];
3742                    double upperValue = columnUpper_[i];
3743                    if (lowerValue > -1.0e20) {
3744                         columnLowerWork_[i] = lowerValue * multiplier;
3745                         if (upperValue >= 1.0e20) {
3746                              columnUpperWork_[i] = COIN_DBL_MAX;
3747                         } else {
3748                              columnUpperWork_[i] = upperValue * multiplier;
3749                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3750                                   if (columnLowerWork_[i] >= 0.0) {
3751                                        columnUpperWork_[i] = columnLowerWork_[i];
3752                                   } else if (columnUpperWork_[i] <= 0.0) {
3753                                        columnLowerWork_[i] = columnUpperWork_[i];
3754                                   } else {
3755                                        columnUpperWork_[i] = 0.0;
3756                                        columnLowerWork_[i] = 0.0;
3757                                   }
3758                              }
3759                         }
3760                    } else if (upperValue < 1.0e20) {
3761                         columnLowerWork_[i] = -COIN_DBL_MAX;
3762                         columnUpperWork_[i] = upperValue * multiplier;
3763                    } else {
3764                         // free
3765                         columnLowerWork_[i] = -COIN_DBL_MAX;
3766                         columnUpperWork_[i] = COIN_DBL_MAX;
3767                    }
3768                    double value = columnActivity_[i] * multiplier;
3769                    if (fabs(value) > 1.0e20) {
3770                         //printf("bad value of %g for column %d\n",value,i);
3771                         setColumnStatus(i, superBasic);
3772                         if (columnUpperWork_[i] < 0.0) {
3773                              value = columnUpperWork_[i];
3774                         } else if (columnLowerWork_[i] > 0.0) {
3775                              value = columnLowerWork_[i];
3776                         } else {
3777                              value = 0.0;
3778                         }
3779                    }
3780                    columnActivityWork_[i] = value;
3781               }
3782               inverseScale = inverseRowScale_;
3783               for (i = 0; i < numberRows_; i++) {
3784                    dual_[i] *= inverseScale[i];
3785                    dual_[i] *= objectiveScale_;
3786                    rowReducedCost_[i] = dual_[i];
3787                    double multiplier = rhsScale_ * rowScale_[i];
3788                    double value = rowActivity_[i] * multiplier;
3789                    if (fabs(value) > 1.0e20) {
3790                         //printf("bad value of %g for row %d\n",value,i);
3791                         setRowStatus(i, superBasic);
3792                         if (rowUpperWork_[i] < 0.0) {
3793                              value = rowUpperWork_[i];
3794                         } else if (rowLowerWork_[i] > 0.0) {
3795                              value = rowLowerWork_[i];
3796                         } else {
3797                              value = 0.0;
3798                         }
3799                    }
3800                    rowActivityWork_[i] = value;
3801               }
3802          } else if (objectiveScale_ != 1.0 || rhsScale_ != 1.0) {
3803               // on entry
3804               for (i = 0; i < numberColumns_; i++) {
3805                    double value = columnActivity_[i];
3806                    value *= rhsScale_;
3807                    if (fabs(value) > 1.0e20) {
3808                         //printf("bad value of %g for column %d\n",value,i);
3809                         setColumnStatus(i, superBasic);
3810                         if (columnUpperWork_[i] < 0.0) {
3811                              value = columnUpperWork_[i];
3812                         } else if (columnLowerWork_[i] > 0.0) {
3813                              value = columnLowerWork_[i];
3814                         } else {
3815                              value = 0.0;
3816                         }
3817                    }
3818                    columnActivityWork_[i] = value;
3819                    reducedCostWork_[i] = reducedCost_[i] * objectiveScale_;
3820               }
3821               for (i = 0; i < numberRows_; i++) {
3822                    double value = rowActivity_[i];
3823                    value *= rhsScale_;
3824                    if (fabs(value) > 1.0e20) {
3825                         //printf("bad value of %g for row %d\n",value,i);
3826                         setRowStatus(i, superBasic);
3827                         if (rowUpperWork_[i] < 0.0) {
3828                              value = rowUpperWork_[i];
3829                         } else if (rowLowerWork_[i] > 0.0) {
3830                              value = rowLowerWork_[i];
3831                         } else {
3832                              value = 0.0;
3833                         }
3834                    }
3835                    rowActivityWork_[i] = value;
3836                    dual_[i] *= objectiveScale_;
3837                    rowReducedCost_[i] = dual_[i];
3838               }
3839          } else {
3840               // on entry
3841               for (i = 0; i < numberColumns_; i++) {
3842                    double value = columnActivity_[i];
3843                    if (fabs(value) > 1.0e20) {
3844                         //printf("bad value of %g for column %d\n",value,i);
3845                         setColumnStatus(i, superBasic);
3846                         if (columnUpperWork_[i] < 0.0) {
3847                              value = columnUpperWork_[i];
3848                         } else if (columnLowerWork_[i] > 0.0) {
3849                              value = columnLowerWork_[i];
3850                         } else {
3851                              value = 0.0;
3852                         }
3853                    }
3854                    columnActivityWork_[i] = value;
3855                    reducedCostWork_[i] = reducedCost_[i];
3856               }
3857               for (i = 0; i < numberRows_; i++) {
3858                    double value = rowActivity_[i];
3859                    if (fabs(value) > 1.0e20) {
3860                         //printf("bad value of %g for row %d\n",value,i);
3861                         setRowStatus(i, superBasic);
3862                         if (rowUpperWork_[i] < 0.0) {
3863                              value = rowUpperWork_[i];
3864                         } else if (rowLowerWork_[i] > 0.0) {
3865                              value = rowLowerWork_[i];
3866                         } else {
3867                              value = 0.0;
3868                         }
3869                    }
3870                    rowActivityWork_[i] = value;
3871                    rowReducedCost_[i] = dual_[i];
3872               }
3873          }
3874     }
3875
3876     if (what == 63 && doSanityCheck) {
3877          // check rim of problem okay
3878          if (!sanityCheck())
3879               goodMatrix = false;
3880     }
3881     // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
3882     // maybe we need to move scales to SimplexModel for factorization?
3883     if ((what == 63 && !pivotVariable_) || (newArrays && !keepPivots)) {
3884          delete [] pivotVariable_;
3885          pivotVariable_ = new int[numberRows2];
3886          for (int i = 0; i < numberRows2; i++)
3887               pivotVariable_[i] = -1;
3888     } else if (what == 63 && !keepPivots) {
3889          // just reset
3890          for (int i = 0; i < numberRows2; i++)
3891               pivotVariable_[i] = -1;
3892     } else if (what == 63) {
3893          // check pivots
3894          for (int i = 0; i < numberRows2; i++) {
3895               int iSequence = pivotVariable_[i];
3896               if (iSequence<0||(iSequence < numberRows_ + numberColumns_ &&
3897                                 getStatus(iSequence) != basic)) {
3898                    keepPivots = false;
3899                    break;
3900               }
3901          }
3902          if (!keepPivots) {
3903               // reset
3904               for (int i = 0; i < numberRows2; i++)
3905                    pivotVariable_[i] = -1;
3906          } else {
3907               // clean
3908               for (int i = 0; i < numberColumns_ + numberRows_; i++) {
3909                    Status status = getStatus(i);
3910                    if (status != basic) {
3911                         if (upper_[i] == lower_[i]) {
3912                              setStatus(i, isFixed);
3913                              solution_[i] = lower_[i];
3914                         } else if (status == atLowerBound) {
3915                              if (lower_[i] > -1.0e20) {
3916                                   solution_[i] = lower_[i];
3917                              } else {
3918                                   //printf("seq %d at lower of %g\n",i,lower_[i]);
3919                                   if (upper_[i] < 1.0e20) {
3920                                        solution_[i] = upper_[i];
3921                                        setStatus(i, atUpperBound);
3922                                   } else {
3923                                        setStatus(i, isFree);
3924                                   }
3925                              }
3926                         } else if (status == atUpperBound) {
3927                              if (upper_[i] < 1.0e20) {
3928                                   solution_[i] = upper_[i];
3929                              } else {
3930                                   //printf("seq %d at upper of %g\n",i,upper_[i]);
3931                                   if (lower_[i] > -1.0e20) {
3932                                        solution_[i] = lower_[i];
3933                                        setStatus(i, atLowerBound);
3934                                   } else {
3935                                        setStatus(i, isFree);
3936                                   }
3937                              }
3938                         } else if (status == isFixed && upper_[i] > lower_[i]) {
3939                              // was fixed - not now
3940                              if (solution_[i] <= lower_[i]) {
3941                                   setStatus(i, atLowerBound);
3942                              } else if (solution_[i] >= upper_[i]) {
3943                                   setStatus(i, atUpperBound);
3944                              } else {
3945                                   setStatus(i, superBasic);
3946                              }
3947                         }
3948                    }
3949               }
3950          }
3951     }
3952
3953     if (what == 63) {
3954          if (newArrays) {
3955               // get some arrays
3956               int iRow, iColumn;
3957               // these are "indexed" arrays so we always know where nonzeros are
3958               /**********************************************************
3959               rowArray_[3] is long enough for rows+columns
3960               rowArray_[1] is long enough for max(rows,columns)
3961               *********************************************************/
3962               for (iRow = 0; iRow < 4; iRow++) {
3963                    int length = numberRows2 + factorization_->maximumPivots();
3964                    if (iRow == 3 || objective_->type() > 1)
3965                         length += numberColumns_;
3966                    else if (iRow == 1)
3967                         length = CoinMax(length, numberColumns_);
3968                    if ((specialOptions_ & 65536) == 0 || !rowArray_[iRow]) {
3969                         delete rowArray_[iRow];
3970                         rowArray_[iRow] = new CoinIndexedVector();
3971                    }
3972                    rowArray_[iRow]->reserve(length);
3973               }
3974
3975               for (iColumn = 0; iColumn < 2; iColumn++) {
3976                    if ((specialOptions_ & 65536) == 0 || !columnArray_[iColumn]) {
3977                         delete columnArray_[iColumn];
3978                         columnArray_[iColumn] = new CoinIndexedVector();
3979                    }
3980                    columnArray_[iColumn]->reserve(numberColumns_+numberRows2);
3981               }
3982          } else {
3983               int iRow, iColumn;
3984               for (iRow = 0; iRow < 4; iRow++) {
3985                    int length = numberRows2 + factorization_->maximumPivots();
3986                    if (iRow == 3 || objective_->type() > 1)
3987                         length += numberColumns_;
3988                    if(rowArray_[iRow]->capacity() >= length) {
3989                         rowArray_[iRow]->clear();
3990                    } else {
3991                         // model size or maxinv changed
3992                         rowArray_[iRow]->reserve(length);
3993                    }
3994#ifndef NDEBUG
3995                    rowArray_[iRow]->checkClear();
3996#endif
3997               }
3998
3999               for (iColumn = 0; iColumn < 2; iColumn++) {
4000                    int length = numberColumns_;
4001                    if (iColumn)
4002                         length = CoinMax(numberRows2, numberColumns_);
4003                    if(columnArray_[iColumn]->capacity() >= length) {
4004                         columnArray_[iColumn]->clear();
4005                    } else {
4006                         // model size or maxinv changed
4007                         columnArray_[iColumn]->reserve(length);
4008                    }
4009#ifndef NDEBUG
4010                    columnArray_[iColumn]->checkClear();
4011#endif
4012               }
4013          }
4014     }
4015     if (problemStatus_ == 10) {
4016          problemStatus_ = -1;
4017          handler_->setLogLevel(saveLevel); // switch back messages
4018     }
4019     if ((what & 5) != 0)
4020          matrix_->generalExpanded(this, 9, what); // update costs and bounds if necessary
4021     if (goodMatrix && (specialOptions_ & 65536) != 0) {
4022          int save = maximumColumns_ + maximumRows_;
4023          CoinMemcpyN(cost_, numberTotal, cost_ + save);
4024          CoinMemcpyN(lower_, numberTotal, lower_ + save);
4025          CoinMemcpyN(upper_, numberTotal, upper_ + save);
4026          CoinMemcpyN(dj_, numberTotal, dj_ + save);
4027          CoinMemcpyN(solution_, numberTotal, solution_ + save);
4028          if (rowScale_ && !savedRowScale_) {
4029               double * temp;
4030               temp = new double [4*maximumRows_];
4031               CoinFillN(temp, 4 * maximumRows_, 1.0);
4032               CoinMemcpyN(rowScale_, numberRows2, temp);
4033               CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + maximumRows_);
4034               CoinMemcpyN(rowScale_, numberRows2, temp + 2 * maximumRows_);
4035               CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + 3 * maximumRows_);
4036               delete [] rowScale_;
4037               savedRowScale_ = temp;
4038               rowScale_ = savedRowScale_;
4039               inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
4040               temp = new double [4*maximumColumns_];
4041               CoinFillN(temp, 4 * maximumColumns_, 1.0);
4042               CoinMemcpyN(columnScale_, numberColumns_, temp);
4043               CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + maximumColumns_);
4044               CoinMemcpyN(columnScale_, numberColumns_, temp + 2 * maximumColumns_);
4045               CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + 3 * maximumColumns_);
4046               delete [] columnScale_;
4047               savedColumnScale_ = temp;
4048               columnScale_ = savedColumnScale_;
4049               inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
4050          }
4051     }
4052#ifdef CLP_USER_DRIVEN
4053     eventHandler_->event(ClpEventHandler::endOfCreateRim);
4054#endif
4055     return goodMatrix;
4056}
4057// Does rows and columns
4058void
4059ClpSimplex::createRim1(bool initial)
4060{
4061     int i;
4062     int numberRows2 = numberRows_ + numberExtraRows_;
4063     int numberTotal = numberRows2 + numberColumns_;
4064     if ((specialOptions_ & 65536) != 0 && true) {
4065          assert (!initial);
4066          int save = maximumColumns_ + maximumRows_;
4067          CoinMemcpyN(lower_ + save, numberTotal, lower_);
4068          CoinMemcpyN(upper_ + save, numberTotal, upper_);
4069          return;
4070     }
4071     const double * rowScale = rowScale_;
4072     // clean up any mismatches on infinity
4073     // and fix any variables with tiny gaps
4074     double primalTolerance = dblParam_[ClpPrimalTolerance];
4075     if(rowScale) {
4076          // If scaled then do all columns later in one loop
4077          if (!initial) {
4078               const double * inverseScale = inverseColumnScale_;
4079               for (i = 0; i < numberColumns_; i++) {
4080                    double multiplier = rhsScale_ * inverseScale[i];
4081                    double lowerValue = columnLower_[i];
4082                    double upperValue = columnUpper_[i];
4083                    if (lowerValue > -1.0e20) {
4084                         columnLowerWork_[i] = lowerValue * multiplier;
4085                         if (upperValue >= 1.0e20) {
4086                              columnUpperWork_[i] = COIN_DBL_MAX;
4087                         } else {
4088                              columnUpperWork_[i] = upperValue * multiplier;
4089                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4090                                   if (columnLowerWork_[i] >= 0.0) {
4091                                        columnUpperWork_[i] = columnLowerWork_[i];
4092                                   } else if (columnUpperWork_[i] <= 0.0) {
4093                                        columnLowerWork_[i] = columnUpperWork_[i];
4094                                   } else {
4095                                        columnUpperWork_[i] = 0.0;
4096                                        columnLowerWork_[i] = 0.0;
4097                                   }
4098                              }
4099                         }
4100                    } else if (upperValue < 1.0e20) {
4101                         columnLowerWork_[i] = -COIN_DBL_MAX;
4102                         columnUpperWork_[i] = upperValue * multiplier;
4103                    } else {
4104                         // free
4105                         columnLowerWork_[i] = -COIN_DBL_MAX;
4106                         columnUpperWork_[i] = COIN_DBL_MAX;
4107                    }
4108               }
4109          }
4110          for (i = 0; i < numberRows_; i++) {
4111               double multiplier = rhsScale_ * rowScale[i];
4112               double lowerValue = rowLower_[i];
4113               double upperValue = rowUpper_[i];
4114               if (lowerValue > -1.0e20) {
4115                    rowLowerWork_[i] = lowerValue * multiplier;
4116                    if (upperValue >= 1.0e20) {
4117                         rowUpperWork_[i] = COIN_DBL_MAX;
4118                    } else {
4119                         rowUpperWork_[i] = upperValue * multiplier;
4120                         if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4121                              if (rowLowerWork_[i] >= 0.0) {
4122                                   rowUpperWork_[i] = rowLowerWork_[i];
4123                              } else if (rowUpperWork_[i] <= 0.0) {
4124                                   rowLowerWork_[i] = rowUpperWork_[i];
4125                              } else {
4126                                   rowUpperWork_[i] = 0.0;
4127                                   rowLowerWork_[i] = 0.0;
4128                              }
4129                         }
4130                    }
4131               } else if (upperValue < 1.0e20) {
4132                    rowLowerWork_[i] = -COIN_DBL_MAX;
4133                    rowUpperWork_[i] = upperValue * multiplier;
4134               } else {
4135                    // free
4136                    rowLowerWork_[i] = -COIN_DBL_MAX;
4137                    rowUpperWork_[i] = COIN_DBL_MAX;
4138               }
4139          }
4140     } else if (rhsScale_ != 1.0) {
4141          for (i = 0; i < numberColumns_; i++) {
4142               double lowerValue = columnLower_[i];
4143               double upperValue = columnUpper_[i];
4144               if (lowerValue > -1.0e20) {
4145                    columnLowerWork_[i] = lowerValue * rhsScale_;
4146                    if (upperValue >= 1.0e20) {
4147                         columnUpperWork_[i] = COIN_DBL_MAX;
4148                    } else {
4149                         columnUpperWork_[i] = upperValue * rhsScale_;
4150                         if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4151                              if (columnLowerWork_[i] >= 0.0) {
4152                                   columnUpperWork_[i] = columnLowerWork_[i];
4153                              } else if (columnUpperWork_[i] <= 0.0) {
4154                                   columnLowerWork_[i] = columnUpperWork_[i];
4155                              } else {
4156                                   columnUpperWork_[i] = 0.0;
4157                                   columnLowerWork_[i] = 0.0;
4158                              }
4159                         }
4160                    }
4161               } else if (upperValue < 1.0e20) {
4162                    columnLowerWork_[i] = -COIN_DBL_MAX;
4163                    columnUpperWork_[i] = upperValue * rhsScale_;
4164               } else {
4165                    // free
4166                    columnLowerWork_[i] = -COIN_DBL_MAX;
4167                    columnUpperWork_[i] = COIN_DBL_MAX;
4168               }
4169          }
4170          for (i = 0; i < numberRows_; i++) {
4171               double lowerValue = rowLower_[i];
4172               double upperValue = rowUpper_[i];
4173               if (lowerValue > -1.0e20) {
4174                    rowLowerWork_[i] = lowerValue * rhsScale_;
4175                    if (upperValue >= 1.0e20) {
4176                         rowUpperWork_[i] = COIN_DBL_MAX;
4177                    } else {
4178                         rowUpperWork_[i] = upperValue * rhsScale_;
4179                         if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4180                              if (rowLowerWork_[i] >= 0.0) {
4181                                   rowUpperWork_[i] = rowLowerWork_[i];
4182                              } else if (rowUpperWork_[i] <= 0.0) {
4183                                   rowLowerWork_[i] = rowUpperWork_[i];
4184                              } else {
4185                                   rowUpperWork_[i] = 0.0;
4186                                   rowLowerWork_[i] = 0.0;
4187                              }
4188                         }
4189                    }
4190               } else if (upperValue < 1.0e20) {
4191                    rowLowerWork_[i] = -COIN_DBL_MAX;
4192                    rowUpperWork_[i] = upperValue * rhsScale_;
4193               } else {
4194                    // free
4195                    rowLowerWork_[i] = -COIN_DBL_MAX;
4196                    rowUpperWork_[i] = COIN_DBL_MAX;
4197               }
4198          }
4199     } else {
4200          for (i = 0; i < numberColumns_; i++) {
4201               double lowerValue = columnLower_[i];
4202               double upperValue = columnUpper_[i];
4203               if (lowerValue > -1.0e20) {
4204                    columnLowerWork_[i] = lowerValue;
4205                    if (upperValue >= 1.0e20) {
4206                         columnUpperWork_[i] = COIN_DBL_MAX;
4207                    } else {
4208                         columnUpperWork_[i] = upperValue;
4209                         if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4210                              if (columnLowerWork_[i] >= 0.0) {
4211                                   columnUpperWork_[i] = columnLowerWork_[i];
4212                              } else if (columnUpperWork_[i] <= 0.0) {
4213                                   columnLowerWork_[i] = columnUpperWork_[i];
4214                              } else {
4215                                   columnUpperWork_[i] = 0.0;
4216                                   columnLowerWork_[i] = 0.0;
4217                              }
4218                         }
4219                    }
4220               } else if (upperValue < 1.0e20) {
4221                    columnLowerWork_[i] = -COIN_DBL_MAX;
4222                    columnUpperWork_[i] = upperValue;
4223               } else {
4224                    // free
4225                    columnLowerWork_[i] = -COIN_DBL_MAX;
4226                    columnUpperWork_[i] = COIN_DBL_MAX;
4227               }
4228          }
4229          for (i = 0; i < numberRows_; i++) {
4230               double lowerValue = rowLower_[i];
4231               double upperValue = rowUpper_[i];
4232               if (lowerValue > -1.0e20) {
4233                    rowLowerWork_[i] = lowerValue;
4234                    if (upperValue >= 1.0e20) {
4235                         rowUpperWork_[i] = COIN_DBL_MAX;
4236                    } else {
4237                         rowUpperWork_[i] = upperValue;
4238                         if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4239                              if (rowLowerWork_[i] >= 0.0) {
4240                                   rowUpperWork_[i] = rowLowerWork_[i];
4241                              } else if (rowUpperWork_[i] <= 0.0) {
4242                                   rowLowerWork_[i] = rowUpperWork_[i];
4243                              } else {
4244                                   rowUpperWork_[i] = 0.0;
4245                                   rowLowerWork_[i] = 0.0;
4246                              }
4247                         }
4248                    }
4249               } else if (upperValue < 1.0e20) {
4250                    rowLowerWork_[i] = -COIN_DBL_MAX;
4251                    rowUpperWork_[i] = upperValue;
4252               } else {
4253                    // free
4254                    rowLowerWork_[i] = -COIN_DBL_MAX;
4255                    rowUpperWork_[i] = COIN_DBL_MAX;
4256               }
4257          }
4258     }
4259#ifndef NDEBUG
4260     if ((specialOptions_ & 65536) != 0 && false) {
4261          assert (!initial);
4262          int save = maximumColumns_ + maximumRows_;
4263          for (int i = 0; i < numberTotal; i++) {
4264               assert (fabs(lower_[i] - lower_[i+save]) < 1.0e-5);
4265               assert (fabs(upper_[i] - upper_[i+save]) < 1.0e-5);
4266          }
4267     }
4268#endif
4269}
4270// Does objective
4271void
4272ClpSimplex::createRim4(bool initial)
4273{
4274     int i;
4275     int numberRows2 = numberRows_ + numberExtraRows_;
4276     int numberTotal = numberRows2 + numberColumns_;
4277     if ((specialOptions_ & 65536) != 0 && true) {
4278          assert (!initial);
4279          int save = maximumColumns_ + maximumRows_;
4280          CoinMemcpyN(cost_ + save, numberTotal, cost_);
4281          return;
4282     }
4283     double direction = optimizationDirection_ * objectiveScale_;
4284     const double * obj = objective();
4285     const double * rowScale = rowScale_;
4286     const double * columnScale = columnScale_;
4287     // and also scale by scale factors
4288     if (rowScale) {
4289          if (rowObjective_) {
4290               for (i = 0; i < numberRows_; i++)
4291                    rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
4292          } else {
4293               memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4294          }
4295          // If scaled then do all columns later in one loop
4296          if (!initial) {
4297               for (i = 0; i < numberColumns_; i++) {
4298                    CoinAssert(fabs(obj[i]) < 1.0e25);
4299                    objectiveWork_[i] = obj[i] * direction * columnScale[i];
4300               }
4301          }
4302     } else {
4303          if (rowObjective_) {
4304               for (i = 0; i < numberRows_; i++)
4305                    rowObjectiveWork_[i] = rowObjective_[i] * direction;
4306          } else {
4307               memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4308          }
4309          for (i = 0; i < numberColumns_; i++) {
4310               CoinAssert(fabs(obj[i]) < 1.0e25);
4311               objectiveWork_[i] = obj[i] * direction;
4312          }
4313     }
4314}
4315// Does rows and columns and objective
4316void
4317ClpSimplex::createRim5(bool initial)
4318{
4319     createRim4(initial);
4320     createRim1(initial);
4321}
4322void
4323ClpSimplex::deleteRim(int getRidOfFactorizationData)
4324{
4325     // Just possible empty problem
4326     int numberRows = numberRows_;
4327     int numberColumns = numberColumns_;
4328     if (!numberRows || !numberColumns) {
4329          numberRows = 0;
4330          if (objective_->type() < 2)
4331               numberColumns = 0;
4332     }
4333     int i;
4334     if (problemStatus_ != 1 && problemStatus_ != 2) {
4335          delete [] ray_;
4336          ray_ = NULL;
4337     }
4338     // set upperOut_ to furthest away from bound so can use in dual for dualBound_
4339     upperOut_ = 1.0;
4340#if 0
4341     {
4342          int nBad = 0;
4343          for (i = 0; i < numberColumns; i++) {
4344               if (lower_[i] == upper_[i] && getColumnStatus(i) == basic)
4345                    nBad++;
4346          }
4347          if (nBad)
4348               printf("yy %d basic fixed\n", nBad);
4349     }
4350#endif
4351     // ray may be null if in branch and bound
4352     if (rowScale_) {
4353          // Collect infeasibilities
4354          int numberPrimalScaled = 0;
4355          int numberPrimalUnscaled = 0;
4356          int numberDualScaled = 0;
4357          int numberDualUnscaled = 0;
4358          double scaleC = 1.0 / objectiveScale_;
4359          double scaleR = 1.0 / rhsScale_;
4360          const double * inverseScale = inverseColumnScale_;
4361          for (i = 0; i < numberColumns; i++) {
4362               double scaleFactor = columnScale_[i];
4363               double valueScaled = columnActivityWork_[i];
4364               double lowerScaled = columnLowerWork_[i];
4365               double upperScaled = columnUpperWork_[i];
4366               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4367                    if (valueScaled < lowerScaled - primalTolerance_ ||
4368                              valueScaled > upperScaled + primalTolerance_)
4369                         numberPrimalScaled++;
4370                    else
4371                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4372               }
4373               columnActivity_[i] = valueScaled * scaleFactor * scaleR;
4374               double value = columnActivity_[i];
4375               if (value < columnLower_[i] - primalTolerance_)
4376                    numberPrimalUnscaled++;
4377               else if (value > columnUpper_[i] + primalTolerance_)
4378                    numberPrimalUnscaled++;
4379               double valueScaledDual = reducedCostWork_[i];
4380               if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4381                    numberDualScaled++;
4382               if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4383                    numberDualScaled++;
4384               reducedCost_[i] = (valueScaledDual * scaleC) * inverseScale[i];
4385               double valueDual = reducedCost_[i];
4386               if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4387                    numberDualUnscaled++;
4388               if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4389                    numberDualUnscaled++;
4390          }
4391          inverseScale = inverseRowScale_;
4392          for (i = 0; i < numberRows; i++) {
4393               double scaleFactor = rowScale_[i];
4394               double valueScaled = rowActivityWork_[i];
4395               double lowerScaled = rowLowerWork_[i];
4396               double upperScaled = rowUpperWork_[i];
4397               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4398                    if (valueScaled < lowerScaled - primalTolerance_ ||
4399                              valueScaled > upperScaled + primalTolerance_)
4400                         numberPrimalScaled++;
4401                    else
4402                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4403               }
4404               rowActivity_[i] = (valueScaled * scaleR) * inverseScale[i];
4405               double value = rowActivity_[i];
4406               if (value < rowLower_[i] - primalTolerance_)
4407                    numberPrimalUnscaled++;
4408               else if (value > rowUpper_[i] + primalTolerance_)
4409                    numberPrimalUnscaled++;
4410               double valueScaledDual = dual_[i] + rowObjectiveWork_[i];;
4411               if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4412                    numberDualScaled++;
4413               if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4414                    numberDualScaled++;
4415               dual_[i] *= scaleFactor * scaleC;
4416               double valueDual = dual_[i];
4417               if (rowObjective_)
4418                    valueDual += rowObjective_[i];
4419               if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4420                    numberDualUnscaled++;
4421               if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4422                    numberDualUnscaled++;
4423          }
4424          if (!problemStatus_ && !secondaryStatus_) {
4425               // See if we need to set secondary status
4426               if (numberPrimalUnscaled) {
4427                    if (numberDualUnscaled)
4428                         secondaryStatus_ = 4;
4429                    else
4430                         secondaryStatus_ = 2;
4431               } else {
4432                    if (numberDualUnscaled)
4433                         secondaryStatus_ = 3;
4434               }
4435          }
4436          if (problemStatus_ == 2 && ray_) {
4437               for (i = 0; i < numberColumns; i++) {
4438                    ray_[i] *= columnScale_[i];
4439               }
4440          } else if (problemStatus_ == 1 && ray_) {
4441               for (i = 0; i < numberRows; i++) {
4442                    ray_[i] *= rowScale_[i];
4443               }
4444          }
4445     } else if (rhsScale_ != 1.0 || objectiveScale_ != 1.0) {
4446          // Collect infeasibilities
4447          int numberPrimalScaled = 0;
4448          int numberPrimalUnscaled = 0;
4449          int numberDualScaled = 0;
4450          int numberDualUnscaled = 0;
4451          double scaleC = 1.0 / objectiveScale_;
4452          double scaleR = 1.0 / rhsScale_;
4453          for (i = 0; i < numberColumns; i++) {
4454               double valueScaled = columnActivityWork_[i];
4455               double lowerScaled = columnLowerWork_[i];
4456               double upperScaled = columnUpperWork_[i];
4457               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4458                    if (valueScaled < lowerScaled - primalTolerance_ ||
4459                              valueScaled > upperScaled + primalTolerance_)
4460                         numberPrimalScaled++;
4461                    else
4462                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4463               }
4464               columnActivity_[i] = valueScaled * scaleR;
4465               double value = columnActivity_[i];
4466               if (value < columnLower_[i] - primalTolerance_)
4467                    numberPrimalUnscaled++;
4468               else if (value > columnUpper_[i] + primalTolerance_)
4469                    numberPrimalUnscaled++;
4470               double valueScaledDual = reducedCostWork_[i];
4471               if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4472                    numberDualScaled++;
4473               if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4474                    numberDualScaled++;
4475               reducedCost_[i] = valueScaledDual * scaleC;
4476               double valueDual = reducedCost_[i];
4477               if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4478                    numberDualUnscaled++;
4479               if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4480                    numberDualUnscaled++;
4481          }
4482          for (i = 0; i < numberRows; i++) {
4483               double valueScaled = rowActivityWork_[i];
4484               double lowerScaled = rowLowerWork_[i];
4485               double upperScaled = rowUpperWork_[i];
4486               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4487                    if (valueScaled < lowerScaled - primalTolerance_ ||
4488                              valueScaled > upperScaled + primalTolerance_)
4489                         numberPrimalScaled++;
4490                    else
4491                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4492               }
4493               rowActivity_[i] = valueScaled * scaleR;
4494               double value = rowActivity_[i];
4495               if (value < rowLower_[i] - primalTolerance_)
4496                    numberPrimalUnscaled++;
4497               else if (value > rowUpper_[i] + primalTolerance_)
4498                    numberPrimalUnscaled++;
4499               double valueScaledDual = dual_[i] + rowObjectiveWork_[i];;
4500               if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4501                    numberDualScaled++;
4502               if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4503                    numberDualScaled++;
4504               dual_[i] *= scaleC;
4505               double valueDual = dual_[i];
4506               if (rowObjective_)
4507                    valueDual += rowObjective_[i];
4508               if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4509                    numberDualUnscaled++;
4510               if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4511                    numberDualUnscaled++;
4512          }
4513          if (!problemStatus_ && !secondaryStatus_) {
4514               // See if we need to set secondary status
4515               if (numberPrimalUnscaled) {
4516                    if (numberDualUnscaled)
4517                         secondaryStatus_ = 4;
4518                    else
4519                         secondaryStatus_ = 2;
4520               } else {
4521                    if (numberDualUnscaled)
4522                         secondaryStatus_ = 3;
4523               }
4524          }
4525     } else {
4526          if (columnActivityWork_) {
4527               for (i = 0; i < numberColumns; i++) {
4528                    double value = columnActivityWork_[i];
4529                    double lower = columnLowerWork_[i];
4530                    double upper = columnUpperWork_[i];
4531                    if (lower > -1.0e20 || upper < 1.0e20) {
4532                         if (value > lower && value < upper)
4533                              upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4534                    }
4535                    columnActivity_[i] = columnActivityWork_[i];
4536                    reducedCost_[i] = reducedCostWork_[i];
4537               }
4538               for (i = 0; i < numberRows; i++) {
4539                    double value = rowActivityWork_[i];
4540                    double lower = rowLowerWork_[i];
4541                    double upper = rowUpperWork_[i];
4542                    if (lower > -1.0e20 || upper < 1.0e20) {
4543                         if (value > lower && value < upper)
4544                              upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4545                    }
4546                    rowActivity_[i] = rowActivityWork_[i];
4547               }
4548          }
4549     }
4550     // switch off scalefactor if auto
4551     if (automaticScale_) {
4552          rhsScale_ = 1.0;
4553          objectiveScale_ = 1.0;
4554     }
4555     if (optimizationDirection_ != 1.0) {
4556          // and modify all dual signs
4557          for (i = 0; i < numberColumns; i++)
4558               reducedCost_[i] *= optimizationDirection_;
4559          for (i = 0; i < numberRows; i++)
4560               dual_[i] *= optimizationDirection_;
4561     }
4562     // scaling may have been turned off
4563     scalingFlag_ = abs(scalingFlag_);
4564     if(getRidOfFactorizationData > 0) {
4565          gutsOfDelete(getRidOfFactorizationData + 1);
4566     } else {
4567          // at least get rid of nonLinearCost_
4568          delete nonLinearCost_;
4569          nonLinearCost_ = NULL;
4570     }
4571     if (!rowObjective_ && problemStatus_ == 0 && objective_->type() == 1 &&
4572               numberRows && numberColumns) {
4573          // Redo objective value
4574          double objectiveValue = 0.0;
4575          const double * cost = objective();
4576          for (int i = 0; i < numberColumns; i++) {
4577               double value = columnActivity_[i];
4578               objectiveValue += value * cost[i];
4579          }
4580          //if (fabs(objectiveValue_ -objectiveValue*optimizationDirection())>1.0e-5)
4581          //printf("old obj %g new %g\n",objectiveValue_, objectiveValue*optimizationDirection());
4582          objectiveValue_ = objectiveValue * optimizationDirection();
4583     }
4584     // get rid of data
4585     matrix_->generalExpanded(this, 13, scalingFlag_);
4586}
4587void
4588ClpSimplex::setDualBound(double value)
4589{
4590     if (value > 0.0)
4591          dualBound_ = value;
4592}
4593void
4594ClpSimplex::setInfeasibilityCost(double value)
4595{
4596     if (value > 0.0)
4597          infeasibilityCost_ = value;
4598}
4599void ClpSimplex::setNumberRefinements( int value)
4600{
4601     if (value >= 0 && value < 10)
4602          numberRefinements_ = value;
4603}
4604// Sets row pivot choice algorithm in dual
4605void
4606ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
4607{
4608     delete dualRowPivot_;
4609     dualRowPivot_ = choice.clone(true);
4610     dualRowPivot_->setModel(this);
4611}
4612// Sets row pivot choice algorithm in dual
4613void
4614ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
4615{
4616     delete primalColumnPivot_;
4617     primalColumnPivot_ = choice.clone(true);
4618     primalColumnPivot_->setModel(this);
4619}
4620void
4621ClpSimplex::setFactorization( ClpFactorization & factorization)
4622{
4623     if (factorization_)
4624          factorization_->setFactorization(factorization);
4625     else
4626          factorization_ = new ClpFactorization(factorization,
4627                                                numberRows_);
4628}
4629
4630// Swaps factorization
4631ClpFactorization *
4632ClpSimplex::swapFactorization( ClpFactorization * factorization)
4633{
4634     ClpFactorization * swap = factorization_;
4635     factorization_ = factorization;
4636     return swap;
4637}
4638// Copies in factorization to existing one
4639void
4640ClpSimplex::copyFactorization( ClpFactorization & factorization)
4641{
4642     *factorization_ = factorization;
4643}
4644/* Perturbation:
4645   -50 to +50 - perturb by this power of ten (-6 sounds good)
4646   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
4647   101 - we are perturbed
4648   102 - don't try perturbing again
4649   default is 100
4650*/
4651void
4652ClpSimplex::setPerturbation(int value)
4653{
4654     if(value <= 100 && value >= -1000) {
4655          perturbation_ = value;
4656     }
4657}
4658// Sparsity on or off
4659bool
4660ClpSimplex::sparseFactorization() const
4661{
4662     return factorization_->sparseThreshold() != 0;
4663}
4664void
4665ClpSimplex::setSparseFactorization(bool value)
4666{
4667     if (value) {
4668          if (!factorization_->sparseThreshold())
4669               factorization_->goSparse();
4670     } else {
4671          factorization_->sparseThreshold(0);
4672     }
4673}
4674void checkCorrect(ClpSimplex * /*model*/, int iRow,
4675                  const double * element, const int * rowStart, const int * rowLength,
4676                  const int * column,
4677                  const double * columnLower_, const double * columnUpper_,
4678                  int /*infiniteUpperC*/,
4679                  int /*infiniteLowerC*/,
4680                  double &maximumUpC,
4681                  double &maximumDownC)
4682{
4683     int infiniteUpper = 0;
4684     int infiniteLower = 0;
4685     double maximumUp = 0.0;
4686     double maximumDown = 0.0;
4687     CoinBigIndex rStart = rowStart[iRow];
4688     CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4689     CoinBigIndex j;
4690     double large = 1.0e15;
4691     int iColumn;
4692     // Compute possible lower and upper ranges
4693
4694     for (j = rStart; j < rEnd; ++j) {
4695          double value = element[j];
4696          iColumn = column[j];
4697          if (value > 0.0) {
4698               if (columnUpper_[iColumn] >= large) {
4699                    ++infiniteUpper;
4700               } else {
4701                    maximumUp += columnUpper_[iColumn] * value;
4702               }
4703               if (columnLower_[iColumn] <= -large) {
4704                    ++infiniteLower;
4705               } else {
4706                    maximumDown += columnLower_[iColumn] * value;
4707               }
4708          } else if (value < 0.0) {
4709               if (columnUpper_[iColumn] >= large) {
4710                    ++infiniteLower;
4711               } else {
4712                    maximumDown += columnUpper_[iColumn] * value;
4713               }
4714               if (columnLower_[iColumn] <= -large) {
4715                    ++infiniteUpper;
4716               } else {
4717                    maximumUp += columnLower_[iColumn] * value;
4718               }
4719          }
4720     }
4721     //assert (infiniteLowerC==infiniteLower);
4722     //assert (infiniteUpperC==infiniteUpper);
4723     if (fabs(maximumUp - maximumUpC) > 1.0e-12 * CoinMax(fabs(maximumUp), fabs(maximumUpC)))
4724          COIN_DETAIL_PRINT(printf("row %d comp up %g, true up %g\n", iRow,
4725                                   maximumUpC, maximumUp));
4726     if (fabs(maximumDown - maximumDownC) > 1.0e-12 * CoinMax(fabs(maximumDown), fabs(maximumDownC)))
4727          COIN_DETAIL_PRINT(printf("row %d comp down %g, true down %g\n", iRow,
4728                 maximumDownC, maximumDown));
4729     maximumUpC = maximumUp;
4730     maximumDownC = maximumDown;
4731}
4732
4733/* Tightens primal bounds to make dual faster.  Unless
4734   fixed, bounds are slightly looser than they could be.
4735   This is to make dual go faster and is probably not needed
4736   with a presolve.  Returns non-zero if problem infeasible
4737
4738   Fudge for branch and bound - put bounds on columns of factor *
4739   largest value (at continuous) - should improve stability
4740   in branch and bound on infeasible branches (0.0 is off)
4741*/
4742int
4743ClpSimplex::tightenPrimalBounds(double factor, int doTight, bool tightIntegers)
4744{
4745
4746     // Get a row copy in standard format
4747     CoinPackedMatrix copy;
4748     copy.setExtraGap(0.0);
4749     copy.setExtraMajor(0.0);
4750     copy.reverseOrderedCopyOf(*matrix());
4751     // Matrix may have been created so get rid of it
4752     matrix_->releasePackedMatrix();
4753     // get matrix data pointers
4754     const int * column = copy.getIndices();
4755     const CoinBigIndex * rowStart = copy.getVectorStarts();
4756     const int * rowLength = copy.getVectorLengths();
4757     const double * element = copy.getElements();
4758     int numberChanged = 1, iPass = 0;
4759     double large = largeValue(); // treat bounds > this as infinite
4760#ifndef NDEBUG
4761     double large2 = 1.0e10 * large;
4762#endif
4763     int numberInfeasible = 0;
4764     int totalTightened = 0;
4765
4766     double tolerance = primalTolerance();
4767
4768
4769     // Save column bounds
4770     double * saveLower = new double [numberColumns_];
4771     CoinMemcpyN(columnLower_, numberColumns_, saveLower);
4772     double * saveUpper = new double [numberColumns_];
4773     CoinMemcpyN(columnUpper_, numberColumns_, saveUpper);
4774
4775     int iRow, iColumn;
4776     // If wanted compute a reasonable dualBound_
4777     if (factor == COIN_DBL_MAX) {
4778          factor = 0.0;
4779          if (dualBound_ == 1.0e10) {
4780               // get largest scaled away from bound
4781               double largest = 1.0e-12;
4782               double largestScaled = 1.0e-12;
4783               int iRow;
4784               for (iRow = 0; iRow < numberRows_; iRow++) {
4785                    double value = rowActivity_[iRow];
4786                    double above = value - rowLower_[iRow];
4787                    double below = rowUpper_[iRow] - value;
4788                    if (above < 1.0e12) {
4789                         largest = CoinMax(largest, above);
4790                    }
4791                    if (below < 1.0e12) {
4792                         largest = CoinMax(largest, below);
4793                    }
4794                    if (rowScale_) {
4795                         double multiplier = rowScale_[iRow];
4796                         above *= multiplier;
4797                         below *= multiplier;
4798                    }
4799                    if (above < 1.0e12) {
4800                         largestScaled = CoinMax(largestScaled, above);
4801                    }
4802                    if (below < 1.0e12) {
4803                         largestScaled = CoinMax(largestScaled, below);
4804                    }
4805               }
4806
4807               int iColumn;
4808               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4809                    double value = columnActivity_[iColumn];
4810                    double above = value - columnLower_[iColumn];
4811                    double below = columnUpper_[iColumn] - value;
4812                    if (above < 1.0e12) {
4813                         largest = CoinMax(largest, above);
4814                    }
4815                    if (below < 1.0e12) {
4816                         largest = CoinMax(largest, below);
4817                    }
4818                    if (columnScale_) {
4819                         double multiplier = 1.0 / columnScale_[iColumn];
4820                         above *= multiplier;
4821                         below *= multiplier;
4822                    }
4823                    if (above < 1.0e12) {
4824                         largestScaled = CoinMax(largestScaled, above);
4825                    }
4826                    if (below < 1.0e12) {
4827                         largestScaled = CoinMax(largestScaled, below);
4828                    }
4829               }
4830               std::cout << "Largest (scaled) away from bound " << largestScaled
4831                         << " unscaled " << largest << std::endl;
4832               dualBound_ = CoinMax(1.0001e7, CoinMin(100.0 * largest, 1.00001e10));
4833          }
4834     }
4835
4836     // If wanted - tighten column bounds using solution
4837     if (factor) {
4838          double largest = 0.0;
4839          if (factor > 0.0) {
4840               assert (factor > 1.0);
4841               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4842                    if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
4843                         largest = CoinMax(largest, fabs(columnActivity_[iColumn]));
4844                    }
4845               }
4846               largest *= factor;
4847          } else {
4848               // absolute
4849               largest = - factor;
4850          }
4851          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4852               if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
4853                    columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn], largest);
4854                    columnLower_[iColumn] = CoinMax(columnLower_[iColumn], -largest);
4855               }
4856          }
4857     }
4858#define MAXPASS 10
4859
4860     // Loop round seeing if we can tighten bounds
4861     // Would be faster to have a stack of possible rows
4862     // and we put altered rows back on stack
4863     int numberCheck = -1;
4864     while(numberChanged > numberCheck) {
4865
4866          numberChanged = 0; // Bounds tightened this pass
4867
4868          if (iPass == MAXPASS) break;
4869          iPass++;
4870
4871          for (iRow = 0; iRow < numberRows_; iRow++) {
4872
4873               if (rowLower_[iRow] > -large || rowUpper_[iRow] < large) {
4874
4875                    // possible row
4876                    int infiniteUpper = 0;
4877                    int infiniteLower = 0;
4878                    double maximumUp = 0.0;
4879                    double maximumDown = 0.0;
4880                    double newBound;
4881                    CoinBigIndex rStart = rowStart[iRow];
4882                    CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4883                    CoinBigIndex j;
4884                    // Compute possible lower and upper ranges
4885
4886                    for (j = rStart; j < rEnd; ++j) {
4887                         double value = element[j];
4888                         iColumn = column[j];
4889                         if (value > 0.0) {
4890                              if (columnUpper_[iColumn] >= large) {
4891                                   ++infiniteUpper;
4892                              } else {
4893                                   maximumUp += columnUpper_[iColumn] * value;
4894                              }
4895                              if (columnLower_[iColumn] <= -large) {
4896                                   ++infiniteLower;
4897                              } else {
4898                                   maximumDown += columnLower_[iColumn] * value;
4899                              }
4900                         } else if (value < 0.0) {
4901                              if (columnUpper_[iColumn] >= large) {
4902                                   ++infiniteLower;
4903                              } else {
4904                                   maximumDown += columnUpper_[iColumn] * value;
4905                              }
4906                              if (columnLower_[iColumn] <= -large) {
4907                                   ++infiniteUpper;
4908                              } else {
4909                                   maximumUp += columnLower_[iColumn] * value;
4910                              }
4911                         }
4912                    }
4913                    // Build in a margin of error
4914                    maximumUp += 1.0e-8 * fabs(maximumUp);
4915                    maximumDown -= 1.0e-8 * fabs(maximumDown);
4916                    double maxUp = maximumUp + infiniteUpper * 1.0e31;
4917                    double maxDown = maximumDown - infiniteLower * 1.0e31;
4918                    if (maxUp <= rowUpper_[iRow] + tolerance &&
4919                              maxDown >= rowLower_[iRow] - tolerance) {
4920
4921                         // Row is redundant - make totally free
4922                         // NO - can't do this for postsolve
4923                         // rowLower_[iRow]=-COIN_DBL_MAX;
4924                         // rowUpper_[iRow]=COIN_DBL_MAX;
4925                         //printf("Redundant row in presolveX %d\n",iRow);
4926
4927                    } else {
4928                         if (maxUp < rowLower_[iRow] - 100.0 * tolerance ||
4929                                   maxDown > rowUpper_[iRow] + 100.0 * tolerance) {
4930                              // problem is infeasible - exit at once
4931                              numberInfeasible++;
4932                              break;
4933                         }
4934                         double lower = rowLower_[iRow];
4935                         double upper = rowUpper_[iRow];
4936                         for (j = rStart; j < rEnd; ++j) {
4937                              double value = element[j];
4938                              iColumn = column[j];
4939                              double nowLower = columnLower_[iColumn];
4940                              double nowUpper = columnUpper_[iColumn];
4941                              if (value > 0.0) {
4942                                   // positive value
4943                                   if (lower > -large) {
4944                                        if (!infiniteUpper) {
4945                                             assert(nowUpper < large2);
4946                                             newBound = nowUpper +
4947                                                        (lower - maximumUp) / value;
4948                                             // relax if original was large
4949                                             if (fabs(maximumUp) > 1.0e8)
4950                                                  newBound -= 1.0e-12 * fabs(maximumUp);
4951                                        } else if (infiniteUpper == 1 && nowUpper > large) {
4952                                             newBound = (lower - maximumUp) / value;
4953                                             // relax if original was large
4954                                             if (fabs(maximumUp) > 1.0e8)
4955                                                  newBound -= 1.0e-12 * fabs(maximumUp);
4956                                        } else {
4957                                             newBound = -COIN_DBL_MAX;
4958                                        }
4959                                        if (newBound > nowLower + 1.0e-12 && newBound > -large) {
4960                                             // Tighten the lower bound
4961                                             numberChanged++;
4962                                             // check infeasible (relaxed)
4963                                             if (nowUpper < newBound) {
4964                                                  if (nowUpper - newBound <
4965                                                            -100.0 * tolerance)
4966                                                       numberInfeasible++;
4967                                                  else
4968                                                       newBound = nowUpper;
4969                                             }
4970                                             columnLower_[iColumn] = newBound;
4971                                             // adjust
4972                                             double now;
4973                                             if (nowLower < -large) {
4974                                                  now = 0.0;
4975                                                  infiniteLower--;
4976                                             } else {
4977                                                  now = nowLower;
4978                                             }
4979                                             maximumDown += (newBound - now) * value;
4980                                             nowLower = newBound;
4981#ifdef DEBUG
4982                                             checkCorrect(this, iRow,
4983                                                          element, rowStart, rowLength,
4984                                                          column,
4985                                                          columnLower_,  columnUpper_,
4986                                                          infiniteUpper,
4987                                                          infiniteLower,
4988                                                          maximumUp,
4989                                                          maximumDown);
4990#endif
4991                                        }
4992                                   }
4993                                   if (upper < large) {
4994                                        if (!infiniteLower) {
4995                                             assert(nowLower > - large2);
4996                                             newBound = nowLower +
4997                                                        (upper - maximumDown) / value;
4998                                             // relax if original was large
4999                                             if (fabs(maximumDown) > 1.0e8)
5000                                                  newBound += 1.0e-12 * fabs(maximumDown);
5001                                        } else if (infiniteLower == 1 && nowLower < -large) {
5002                                             newBound =   (upper - maximumDown) / value;
5003                                             // relax if original was large
5004                                             if (fabs(maximumDown) > 1.0e8)
5005                                                  newBound += 1.0e-12 * fabs(maximumDown);
5006                                        } else {
5007                                             newBound = COIN_DBL_MAX;
5008                                        }
5009                                        if (newBound < nowUpper - 1.0e-12 && newBound < large) {
5010                                             // Tighten the upper bound
5011                                             numberChanged++;
5012                                             // check infeasible (relaxed)
5013                                             if (nowLower > newBound) {
5014                                                  if (newBound - nowLower <
5015                                                            -100.0 * tolerance)
5016                                                       numberInfeasible++;
5017                                                  else
5018                                                       newBound = nowLower;
5019                                             }
5020                                             columnUpper_[iColumn] = newBound;
5021                                             // adjust
5022                                             double now;
5023                                             if (nowUpper > large) {
5024                                                  now = 0.0;
5025                                                  infiniteUpper--;
5026                                             } else {
5027                                                  now = nowUpper;
5028                                             }
5029                                             maximumUp += (newBound - now) * value;
5030                                             nowUpper = newBound;
5031#ifdef DEBUG
5032                                             checkCorrect(this, iRow,
5033                                                          element, rowStart, rowLength,
5034                                                          column,
5035                                                          columnLower_,  columnUpper_,
5036                                                          infiniteUpper,
5037                                                          infiniteLower,
5038                                                          maximumUp,
5039                                                          maximumDown);
5040#endif
5041                                        }
5042                                   }
5043                              } else {
5044                                   // negative value
5045                                   if (lower > -large) {
5046                                        if (!infiniteUpper) {
5047                                             assert(nowLower < large2);
5048                                             newBound = nowLower +
5049                                                        (lower - maximumUp) / value;
5050                                             // relax if original was large
5051                                             if (fabs(maximumUp) > 1.0e8)
5052                                                  newBound += 1.0e-12 * fabs(maximumUp);
5053                                        } else if (infiniteUpper == 1 && nowLower < -large) {
5054                                             newBound = (lower - maximumUp) / value;
5055                                             // relax if original was large
5056                                             if (fabs(maximumUp) > 1.0e8)
5057                                                  newBound += 1.0e-12 * fabs(maximumUp);
5058                                        } else {
5059                                             newBound = COIN_DBL_MAX;
5060                                        }
5061                                        if (newBound < nowUpper - 1.0e-12 && newBound < large) {
5062                                             // Tighten the upper bound
5063                                             numberChanged++;
5064                                             // check infeasible (relaxed)
5065                                             if (nowLower > newBound) {
5066                                                  if (newBound - nowLower <
5067                                                            -100.0 * tolerance)
5068                                                       numberInfeasible++;
5069                                                  else
5070                                                       newBound = nowLower;
5071                                             }
5072                                             columnUpper_[iColumn] = newBound;
5073                                             // adjust
5074                                             double now;
5075                                             if (nowUpper > large) {
5076                                                  now = 0.0;
5077                                                  infiniteLower--;
5078                                             } else {
5079                                                  now = nowUpper;
5080                                             }
5081                                             maximumDown += (newBound - now) * value;
5082                                             nowUpper = newBound;
5083#ifdef DEBUG
5084                                             checkCorrect(this, iRow,
5085                                                          element, rowStart, rowLength,
5086                                                          column,
5087                                                          columnLower_,  columnUpper_,
5088                                                          infiniteUpper,
5089                                                          infiniteLower,
5090                                                          maximumUp,
5091                                                          maximumDown);
5092#endif
5093                                        }
5094                                   }
5095                                   if (upper < large) {
5096                                        if (!infiniteLower) {
5097                                             assert(nowUpper < large2);
5098                                             newBound = nowUpper +
5099                                                        (upper - maximumDown) / value;
5100                                             // relax if original was large
5101                                             if (fabs(maximumDown) > 1.0e8)
5102                                                  newBound -= 1.0e-12 * fabs(maximumDown);
5103                                        } else if (infiniteLower == 1 && nowUpper > large) {
5104                                             newBound =   (upper - maximumDown) / value;
5105                                             // relax if original was large
5106                                             if (fabs(maximumDown) > 1.0e8)
5107                                                  newBound -= 1.0e-12 * fabs(maximumDown);
5108                                        } else {
5109                                             newBound = -COIN_DBL_MAX;
5110                                        }
5111                                        if (newBound > nowLower + 1.0e-12 && newBound > -large) {
5112                                             // Tighten the lower bound
5113                                             numberChanged++;
5114                                             // check infeasible (relaxed)
5115                                             if (nowUpper < newBound) {
5116                                                  if (nowUpper - newBound <
5117                                                            -100.0 * tolerance)
5118                                                       numberInfeasible++;
5119