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

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

make more threadsafe

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