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

Last change on this file since 2018 was 1989, checked in by forrest, 6 years ago

Minor stability fixes and an option to allow perturbation after presolve

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