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

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

fix some ampl stuff, add ClpSolver? and a few fixes

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