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

Last change on this file since 2383 was 2383, checked in by forrest, 4 months ago

allow compressed .lp files

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