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

Last change on this file since 2292 was 2292, checked in by forrest, 2 years ago

to stop seg fault if disableFactorization out of place

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 529.4 KB
Line 
1/* $Id: ClpSimplex.cpp 2292 2017-12-04 16:44:17Z 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_&2097152)==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 (fabs(lower-value) <= fabs(upper-value)) {
1861                                        setColumnStatus(iColumn, atLowerBound);
1862                                        columnActivityWork_[iColumn] = lower;
1863                                   } else {
1864                                        setColumnStatus(iColumn, atUpperBound);
1865                                        columnActivityWork_[iColumn] = upper;
1866                                   }
1867                              }
1868                         } else {
1869                              setColumnStatus(iColumn, isFree);
1870                              columnActivityWork_[iColumn] = 0.0;
1871                         }
1872                    }
1873               }
1874          } else {
1875               // values pass has less coding
1876               // make row activities correct and clean basis a bit
1877               cleanStatus();
1878               if (status_) {
1879                    int numberBasic = 0;
1880                    for (iRow = 0; iRow < numberRows_; iRow++) {
1881                         if (getRowStatus(iRow) == basic)
1882                              numberBasic++;
1883                    }
1884                    totalSlacks = numberBasic;
1885#if 0
1886                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1887                         if (getColumnStatus(iColumn) == basic)
1888                              numberBasic++;
1889                    }
1890#endif
1891               } else {
1892                    // all slack basis
1893                    int numberBasic = 0;
1894                    if (!status_) {
1895                         createStatus();
1896                    }
1897                    for (iRow = 0; iRow < numberRows_; iRow++) {
1898                         setRowStatus(iRow, basic);
1899                         numberBasic++;
1900                    }
1901                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1902                         setColumnStatus(iColumn, superBasic);
1903                         // but put to bound if close
1904                         if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1905                                   <= primalTolerance_) {
1906                              columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1907                              setColumnStatus(iColumn, atLowerBound);
1908                         } else if (fabs(columnActivityWork_[iColumn]
1909                                         - columnUpperWork_[iColumn])
1910                                    <= primalTolerance_) {
1911                              columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1912                              setColumnStatus(iColumn, atUpperBound);
1913                         }
1914                    }
1915               }
1916          }
1917          numberRefinements_ = 1;
1918          // set fixed if they are
1919          for (iRow = 0; iRow < numberRows_; iRow++) {
1920               if (getRowStatus(iRow) != basic ) {
1921                    if (rowLowerWork_[iRow] == rowUpperWork_[iRow]) {
1922                         rowActivityWork_[iRow] = rowLowerWork_[iRow];
1923                         setRowStatus(iRow, isFixed);
1924                    }
1925               }
1926          }
1927          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1928               if (getColumnStatus(iColumn) != basic ) {
1929                    if (columnLowerWork_[iColumn] == columnUpperWork_[iColumn]) {
1930                         columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1931                         setColumnStatus(iColumn, isFixed);
1932                    }
1933               }
1934          }
1935     }
1936     //for (iRow=0;iRow<numberRows_+numberColumns_;iRow++) {
1937     //if (fabs(solution_[iRow])>1.0e10) {
1938     //  printf("large %g at %d - status %d\n",
1939     //         solution_[iRow],iRow,status_[iRow]);
1940     //}
1941     //}
1942#    if 0 //ndef _MSC_VER
1943         // The local static var k is a problem when trying to build a DLL. Since this is
1944         // just for debugging (likely done on *nix), just hide it from Windows
1945         // -- lh, 101016 --
1946     if (0)  {
1947          static int k = 0;
1948          printf("start basis\n");
1949          int i;
1950          for (i = 0; i < numberRows_; i++)
1951               printf ("xx %d %d\n", i, pivotVariable_[i]);
1952          for (i = 0; i < numberRows_ + numberColumns_; i++)
1953               if (getColumnStatus(i) == basic)
1954                    printf ("yy %d basic\n", i);
1955          if (k > 20)
1956               exit(0);
1957          k++;
1958     }
1959#    endif
1960#if 0 //ndef NDEBUG
1961  // Make sure everything is clean
1962  double sumOutside=0.0;
1963  int numberOutside=0;
1964  //double sumOutsideLarge=0.0;
1965  int numberOutsideLarge=0;
1966  double sumInside=0.0;
1967  int numberInside=0;
1968  //double sumInsideLarge=0.0;
1969  int numberInsideLarge=0;
1970  int numberTotal=numberRows_+numberColumns_;
1971  for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
1972    if(getStatus(iSequence) == isFixed) {
1973      // double check fixed
1974      assert (upper_[iSequence] == lower_[iSequence]);
1975      assert (fabs(solution_[iSequence]-lower_[iSequence])<primalTolerance_);
1976    } else if (getStatus(iSequence) == isFree) {
1977      assert (upper_[iSequence] == COIN_DBL_MAX && lower_[iSequence]==-COIN_DBL_MAX);
1978    } else if (getStatus(iSequence) == atLowerBound) {
1979      assert (fabs(solution_[iSequence]-lower_[iSequence])<1000.0*primalTolerance_);
1980      if (solution_[iSequence]<lower_[iSequence]) {
1981        numberOutside++;
1982        sumOutside-=solution_[iSequence]-lower_[iSequence];
1983        if (solution_[iSequence]<lower_[iSequence]-primalTolerance_)
1984          numberOutsideLarge++;
1985      } else if (solution_[iSequence]>lower_[iSequence]) {
1986        numberInside++;
1987        sumInside+=solution_[iSequence]-lower_[iSequence];
1988        if (solution_[iSequence]>lower_[iSequence]+primalTolerance_)
1989          numberInsideLarge++;
1990      }
1991    } else if (getStatus(iSequence) == atUpperBound) {
1992      assert (fabs(solution_[iSequence]-upper_[iSequence])<1000.0*primalTolerance_);
1993      if (solution_[iSequence]>upper_[iSequence]) {
1994        numberOutside++;
1995        sumOutside+=solution_[iSequence]-upper_[iSequence];
1996        if (solution_[iSequence]>upper_[iSequence]+primalTolerance_)
1997          numberOutsideLarge++;
1998      } else if (solution_[iSequence]<upper_[iSequence]) {
1999        numberInside++;
2000        sumInside-=solution_[iSequence]-upper_[iSequence];
2001        if (solution_[iSequence]<upper_[iSequence]-primalTolerance_)
2002          numberInsideLarge++;
2003      }
2004    } else if (getStatus(iSequence) == superBasic) {
2005      //assert (!valuesPass);
2006    }
2007  }
2008  if (numberInside+numberOutside)
2009    printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
2010           numberOutside,sumOutside,numberOutsideLarge,
2011           numberInside,sumInside,numberInsideLarge);
2012#endif
2013     int status = factorization_->factorize(this, solveType, valuesPass);
2014     if (status) {
2015          handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
2016                    << status
2017                    << CoinMessageEol;
2018#ifdef CLP_USEFUL_PRINTOUT
2019          printf("Basis singular - pivot tolerance %g\n",
2020                 factorization_->pivotTolerance());
2021#endif
2022          return -1;
2023     } else if (!solveType) {
2024          // Initial basis - return number of singularities
2025          int numberSlacks = 0;
2026          for (iRow = 0; iRow < numberRows_; iRow++) {
2027               if (getRowStatus(iRow) == basic)
2028                    numberSlacks++;
2029          }
2030          status = CoinMax(numberSlacks - totalSlacks, 0);
2031          // special case if all slack
2032          if (numberSlacks == numberRows_) {
2033               status = numberRows_ + 1;
2034          }
2035     }
2036
2037     // sparse methods
2038     //if (factorization_->sparseThreshold()) {
2039     // get default value
2040     factorization_->sparseThreshold(0);
2041     if (!(moreSpecialOptions_&1024))
2042       factorization_->goSparse();
2043     //}
2044
2045     return status;
2046}
2047/*
2048   This does basis housekeeping and does values for in/out variables.
2049   Can also decide to re-factorize
2050*/
2051int
2052ClpSimplex::housekeeping(double objectiveChange)
2053{
2054     // save value of incoming and outgoing
2055     double oldIn = solution_[sequenceIn_];
2056     double oldOut = solution_[sequenceOut_];
2057     numberIterations_++;
2058     changeMade_++; // something has happened
2059     // incoming variable
2060     if (handler_->logLevel() > 7) {
2061          //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
2062          handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
2063                    << directionOut_
2064                    << directionIn_ << theta_
2065                    << dualOut_ << dualIn_ << alpha_
2066                    << CoinMessageEol;
2067          if (getStatus(sequenceIn_) == isFree) {
2068               handler_->message(CLP_SIMPLEX_FREEIN, messages_)
2069                         << sequenceIn_
2070                         << CoinMessageEol;
2071          }
2072     }
2073#if 0
2074     printf("h1 %d %d %g %g %g %g",
2075            directionOut_
2076            , directionIn_, theta_
2077            , dualOut_, dualIn_, alpha_);
2078#endif
2079     // change of incoming
2080     char rowcol[] = {'R', 'C'};
2081     if (pivotRow_ >= 0)
2082          pivotVariable_[pivotRow_] = sequenceIn();
2083     if (upper_[sequenceIn_] > 1.0e20 && lower_[sequenceIn_] < -1.0e20)
2084          progressFlag_ |= 2; // making real progress
2085     solution_[sequenceIn_] = valueIn_;
2086     if (upper_[sequenceOut_] - lower_[sequenceOut_] < 1.0e-12)
2087          progressFlag_ |= 1; // making real progress
2088     if (sequenceIn_ != sequenceOut_) {
2089          if (alphaAccuracy_ > 0.0) {
2090               double value = fabs(alpha_);
2091               if (value > 1.0)
2092                    alphaAccuracy_ *= value;
2093               else
2094                    alphaAccuracy_ /= value;
2095          }
2096          //assert( getStatus(sequenceOut_)== basic);
2097          setStatus(sequenceIn_, basic);
2098          if (upper_[sequenceOut_] - lower_[sequenceOut_] > 0) {
2099               // As Nonlinear costs may have moved bounds (to more feasible)
2100               // Redo using value
2101               if (fabs(valueOut_ - lower_[sequenceOut_]) < fabs(valueOut_ - upper_[sequenceOut_])) {
2102                    // going to lower
2103                    setStatus(sequenceOut_, atLowerBound);
2104                    oldOut = lower_[sequenceOut_];
2105               } else {
2106                    // going to upper
2107                    setStatus(sequenceOut_, atUpperBound);
2108                    oldOut = upper_[sequenceOut_];
2109               }
2110          } else {
2111               // fixed
2112               setStatus(sequenceOut_, isFixed);
2113          }
2114          solution_[sequenceOut_] = valueOut_;
2115     } else {
2116          //if (objective_->type()<2)
2117          //assert (fabs(theta_)>1.0e-13);
2118          // flip from bound to bound
2119          // As Nonlinear costs may have moved bounds (to more feasible)
2120          // Redo using value
2121          if (fabs(valueIn_ - lower_[sequenceIn_]) < fabs(valueIn_ - upper_[sequenceIn_])) {
2122               // as if from upper bound
2123               setStatus(sequenceIn_, atLowerBound);
2124          } else {
2125               // as if from lower bound
2126               setStatus(sequenceIn_, atUpperBound);
2127          }
2128     }
2129
2130     // Update hidden stuff e.g. effective RHS and gub
2131     int invertNow=matrix_->updatePivot(this, oldIn, oldOut);
2132     objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
2133     if (handler_->logLevel() > 7) {
2134          //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
2135          handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
2136                    << numberIterations_ << objectiveValue()
2137                    << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
2138                    << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
2139          handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
2140          handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
2141          handler_->message() << CoinMessageEol;
2142     }
2143#if 0
2144     if (numberIterations_ > 10000)
2145          printf(" it %d %g %c%d %c%d\n"
2146                 , numberIterations_, objectiveValue()
2147                 , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
2148                 , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
2149#endif
2150     if (trustedUserPointer_ && trustedUserPointer_->typeStruct == 1) {
2151          if (algorithm_ > 0 && integerType_ && !nonLinearCost_->numberInfeasibilities()) {
2152               if (fabs(theta_) > 1.0e-6 || !numberIterations_) {
2153                    // For saving solutions
2154                    typedef struct {
2155                         int numberSolutions;
2156                         int maximumSolutions;
2157                         int numberColumns;
2158                         double ** solution;
2159                         int * numberUnsatisfied;
2160                    } clpSolution;
2161                    clpSolution * solution = reinterpret_cast<clpSolution *> (trustedUserPointer_->data);
2162                    if (solution->numberSolutions == solution->maximumSolutions) {
2163                         int n =  solution->maximumSolutions;
2164                         int n2 = (n * 3) / 2 + 10;
2165                         solution->maximumSolutions = n2;
2166                         double ** temp = new double * [n2];
2167                         for (int i = 0; i < n; i++)
2168                              temp[i] = solution->solution[i];
2169                         delete [] solution->solution;
2170                         solution->solution = temp;
2171                         int * tempN = new int [n2];
2172                         for (int i = 0; i < n; i++)
2173                              tempN[i] = solution->numberUnsatisfied[i];
2174                         delete [] solution->numberUnsatisfied;
2175                         solution->numberUnsatisfied = tempN;
2176                    }
2177                    assert (numberColumns_ == solution->numberColumns);
2178                    double * sol = new double [numberColumns_];
2179                    solution->solution[solution->numberSolutions] = sol;
2180                    int numberFixed = 0;
2181                    int numberUnsat = 0;
2182                    int numberSat = 0;
2183                    double sumUnsat = 0.0;
2184                    double tolerance = 10.0 * primalTolerance_;
2185                    double mostAway = 0.0;
2186                    for (int i = 0; i < numberColumns_; i++) {
2187                         // Save anyway
2188                         sol[i] = columnScale_ ? solution_[i] * columnScale_[i] : solution_[i];
2189                         // rest is optional
2190                         if (upper_[i] > lower_[i]) {
2191                              double value = solution_[i];
2192                              if (value > lower_[i] + tolerance &&
2193                                        value < upper_[i] - tolerance && integerType_[i]) {
2194                                   // may have to modify value if scaled
2195                                   if (columnScale_)
2196                                        value *= columnScale_[i];
2197                                   double closest = floor(value + 0.5);
2198                                   // problem may be perturbed so relax test
2199                                   if (fabs(value - closest) > 1.0e-4) {
2200                                        numberUnsat++;
2201                                        sumUnsat += fabs(value - closest);
2202                                        if (mostAway < fabs(value - closest)) {
2203                                             mostAway = fabs(value - closest);
2204                                        }
2205                                   } else {
2206                                        numberSat++;
2207                                   }
2208                              } else {
2209                                   numberSat++;
2210                              }
2211                         } else {
2212                              numberFixed++;
2213                         }
2214                    }
2215                    solution->numberUnsatisfied[solution->numberSolutions++] = numberUnsat;
2216                    COIN_DETAIL_PRINT(printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n",
2217                                             numberIterations_, numberUnsat, sumUnsat, mostAway, numberFixed, numberSat));
2218               }
2219          }
2220     }
2221     if (hitMaximumIterations())
2222          return 2;
2223#if 1
2224     //if (numberIterations_>14000)
2225     //handler_->setLogLevel(63);
2226     //if (numberIterations_>24000)
2227     //exit(77);
2228     // check for small cycles
2229     int in = sequenceIn_;
2230     int out = sequenceOut_;
2231     matrix_->correctSequence(this, in, out);
2232     int cycle = progress_.cycle(in, out,
2233                                 directionIn_, directionOut_);
2234     if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
2235          //if (cycle>0) {
2236          if (handler_->logLevel() >= 63)
2237               printf("Cycle of %d\n", cycle);
2238          // reset
2239          progress_.startCheck();
2240          double random = randomNumberGenerator_.randomDouble();
2241          int extra = static_cast<int> (9.999 * random);
2242          int off[] = {1, 1, 1, 1, 2, 2, 2, 3, 3, 4};
2243          if (factorization_->pivots() > cycle) {
2244               forceFactorization_ = CoinMax(1, cycle - off[extra]);
2245          } else {
2246            /* need to reject something
2247               should be better if don't reject incoming
2248               as it is in basis */
2249               int iSequence;
2250               //if (algorithm_ > 0)
2251               //   iSequence = sequenceIn_;
2252               //else
2253                    iSequence = sequenceOut_;
2254               char x = isColumn(iSequence) ? 'C' : 'R';
2255               if (handler_->logLevel() >= 63)
2256                    handler_->message(CLP_SIMPLEX_FLAG, messages_)
2257                              << x << sequenceWithin(iSequence)
2258                              << CoinMessageEol;
2259               setFlagged(iSequence);
2260               //printf("flagging %d\n",iSequence);
2261          }
2262          return 1;
2263     }
2264#endif
2265     // only time to re-factorize if one before real time
2266     // this is so user won't be surprised that maximumPivots has exact meaning
2267     int numberPivots = factorization_->pivots();
2268     int maximumPivots = factorization_->maximumPivots();
2269     int numberDense = factorization_->numberDense();
2270     bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 >
2271                        2 * maximumIterations());
2272     if (numberPivots == maximumPivots ||
2273               maximumPivots < 2) {
2274          // If dense then increase
2275          if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots
2276              && false) {
2277               factorization_->maximumPivots(numberDense);
2278               dualRowPivot_->maximumPivotsChanged();
2279               primalColumnPivot_->maximumPivotsChanged();
2280               // and redo arrays
2281               for (int iRow = 0; iRow < 4; iRow++) {
2282                    int length = rowArray_[iRow]->capacity() + numberDense - maximumPivots;
2283                    rowArray_[iRow]->reserve(length);
2284               }
2285          }
2286#if CLP_FACTORIZATION_NEW_TIMING>1
2287          factorization_->statsRefactor('M');
2288#endif
2289          return 1;
2290     } else if ((factorization_->timeToRefactorize() && !dontInvert)
2291                ||invertNow) {
2292          //printf("ret after %d pivots\n",factorization_->pivots());
2293#if CLP_FACTORIZATION_NEW_TIMING>1
2294          factorization_->statsRefactor('T');
2295#endif
2296          return 1;
2297     } else if (forceFactorization_ > 0 &&
2298                factorization_->pivots() == forceFactorization_) {
2299          // relax
2300          forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
2301          if (forceFactorization_ > factorization_->maximumPivots())
2302               forceFactorization_ = -1; //off
2303#if CLP_FACTORIZATION_NEW_TIMING>1
2304          factorization_->statsRefactor('F');
2305#endif
2306          return 1;
2307     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
2308          double random = randomNumberGenerator_.randomDouble();
2309          while (random<0.45)
2310            random *= 2.0;
2311          int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
2312          if (factorization_->pivots() >= random * maxNumber) {
2313               return 1;
2314          } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) &&
2315                     numberIterations_ < 1001000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
2316               return 1;
2317          } else {
2318               // carry on iterating
2319               return 0;
2320          }
2321     } else {
2322          // carry on iterating
2323          return 0;
2324     }
2325}
2326// Copy constructor.
2327ClpSimplex::ClpSimplex(const ClpSimplex &rhs, int scalingMode) :
2328     ClpModel(rhs, scalingMode),
2329     bestPossibleImprovement_(0.0),
2330     zeroTolerance_(1.0e-13),
2331     columnPrimalSequence_(-2),
2332     rowPrimalSequence_(-2),
2333     bestObjectiveValue_(rhs.bestObjectiveValue_),
2334     moreSpecialOptions_(2),
2335     baseIteration_(0),
2336     vectorMode_(0),
2337     primalToleranceToGetOptimal_(-1.0),
2338     largeValue_(1.0e15),
2339     largestPrimalError_(0.0),
2340     largestDualError_(0.0),
2341     alphaAccuracy_(-1.0),
2342     dualBound_(1.0e10),
2343     alpha_(0.0),
2344     theta_(0.0),
2345     lowerIn_(0.0),
2346     valueIn_(0.0),
2347     upperIn_(-COIN_DBL_MAX),
2348     dualIn_(0.0),
2349     lowerOut_(-1),
2350     valueOut_(-1),
2351     upperOut_(-1),
2352     dualOut_(-1),
2353     dualTolerance_(1.0e-7),
2354     primalTolerance_(1.0e-7),
2355     sumDualInfeasibilities_(0.0),
2356     sumPrimalInfeasibilities_(0.0),
2357     infeasibilityCost_(1.0e10),
2358     sumOfRelaxedDualInfeasibilities_(0.0),
2359     sumOfRelaxedPrimalInfeasibilities_(0.0),
2360     acceptablePivot_(1.0e-8),
2361     lower_(NULL),
2362     rowLowerWork_(NULL),
2363     columnLowerWork_(NULL),
2364     upper_(NULL),
2365     rowUpperWork_(NULL),
2366     columnUpperWork_(NULL),
2367     cost_(NULL),
2368     rowObjectiveWork_(NULL),
2369     objectiveWork_(NULL),
2370     sequenceIn_(-1),
2371     directionIn_(-1),
2372     sequenceOut_(-1),
2373     directionOut_(-1),
2374     pivotRow_(-1),
2375     lastGoodIteration_(-100),
2376     dj_(NULL),
2377     rowReducedCost_(NULL),
2378     reducedCostWork_(NULL),
2379     solution_(NULL),
2380     rowActivityWork_(NULL),
2381     columnActivityWork_(NULL),
2382     numberDualInfeasibilities_(0),
2383     numberDualInfeasibilitiesWithoutFree_(0),
2384     numberPrimalInfeasibilities_(100),
2385     numberRefinements_(0),
2386     pivotVariable_(NULL),
2387     factorization_(NULL),
2388     savedSolution_(NULL),
2389     numberTimesOptimal_(0),
2390     disasterArea_(NULL),
2391     changeMade_(1),
2392     algorithm_(0),
2393     forceFactorization_(-1),
2394     perturbation_(100),
2395     nonLinearCost_(NULL),
2396     lastBadIteration_(-999999),
2397     lastFlaggedIteration_(-999999),
2398     numberFake_(0),
2399     numberChanged_(0),
2400     progressFlag_(0),
2401     firstFree_(-1),
2402     numberExtraRows_(0),
2403     maximumBasic_(0),
2404     dontFactorizePivots_(0),
2405     incomingInfeasibility_(1.0),
2406     allowedInfeasibility_(10.0),
2407     automaticScale_(0),
2408     maximumPerturbationSize_(0),
2409     perturbationArray_(NULL),
2410     baseModel_(NULL)
2411#ifdef ABC_INHERIT
2412     ,abcSimplex_(NULL),
2413     abcState_(0)
2414#endif
2415{
2416     int i;
2417     for (i = 0; i < 6; i++) {
2418          rowArray_[i] = NULL;
2419          columnArray_[i] = NULL;
2420     }
2421     for (i = 0; i < 4; i++) {
2422          spareIntArray_[i] = 0;
2423          spareDoubleArray_[i] = 0.0;
2424     }
2425     saveStatus_ = NULL;
2426     factorization_ = NULL;
2427     dualRowPivot_ = NULL;
2428     primalColumnPivot_ = NULL;
2429     gutsOfDelete(0);
2430     delete nonLinearCost_;
2431     nonLinearCost_ = NULL;
2432     gutsOfCopy(rhs);
2433     solveType_ = 1; // say simplex based life form
2434}
2435// Copy constructor from model
2436ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
2437     ClpModel(rhs, scalingMode),
2438     bestPossibleImprovement_(0.0),
2439     zeroTolerance_(1.0e-13),
2440     columnPrimalSequence_(-2),
2441     rowPrimalSequence_(-2),
2442     bestObjectiveValue_(-COIN_DBL_MAX),
2443     moreSpecialOptions_(2),
2444     baseIteration_(0),
2445     vectorMode_(0),
2446     primalToleranceToGetOptimal_(-1.0),
2447     largeValue_(1.0e15),
2448     largestPrimalError_(0.0),
2449     largestDualError_(0.0),
2450     alphaAccuracy_(-1.0),
2451     dualBound_(1.0e10),
2452     alpha_(0.0),
2453     theta_(0.0),
2454     lowerIn_(0.0),
2455     valueIn_(0.0),
2456     upperIn_(-COIN_DBL_MAX),
2457     dualIn_(0.0),
2458     lowerOut_(-1),
2459     valueOut_(-1),
2460     upperOut_(-1),
2461     dualOut_(-1),
2462     dualTolerance_(1.0e-7),
2463     primalTolerance_(1.0e-7),
2464     sumDualInfeasibilities_(0.0),
2465     sumPrimalInfeasibilities_(0.0),
2466     infeasibilityCost_(1.0e10),
2467     sumOfRelaxedDualInfeasibilities_(0.0),
2468     sumOfRelaxedPrimalInfeasibilities_(0.0),
2469     acceptablePivot_(1.0e-8),
2470     lower_(NULL),
2471     rowLowerWork_(NULL),
2472     columnLowerWork_(NULL),
2473     upper_(NULL),
2474     rowUpperWork_(NULL),
2475     columnUpperWork_(NULL),
2476     cost_(NULL),
2477     rowObjectiveWork_(NULL),
2478     objectiveWork_(NULL),
2479     sequenceIn_(-1),
2480     directionIn_(-1),
2481     sequenceOut_(-1),
2482     directionOut_(-1),
2483     pivotRow_(-1),
2484     lastGoodIteration_(-100),
2485     dj_(NULL),
2486     rowReducedCost_(NULL),
2487     reducedCostWork_(NULL),
2488     solution_(NULL),
2489     rowActivityWork_(NULL),
2490     columnActivityWork_(NULL),
2491     numberDualInfeasibilities_(0),
2492     numberDualInfeasibilitiesWithoutFree_(0),
2493     numberPrimalInfeasibilities_(100),
2494     numberRefinements_(0),
2495     pivotVariable_(NULL),
2496     factorization_(NULL),
2497     savedSolution_(NULL),
2498     numberTimesOptimal_(0),
2499     disasterArea_(NULL),
2500     changeMade_(1),
2501     algorithm_(0),
2502     forceFactorization_(-1),
2503     perturbation_(100),
2504     nonLinearCost_(NULL),
2505     lastBadIteration_(-999999),
2506     lastFlaggedIteration_(-999999),
2507     numberFake_(0),
2508     numberChanged_(0),
2509     progressFlag_(0),
2510     firstFree_(-1),
2511     numberExtraRows_(0),
2512     maximumBasic_(0),
2513     dontFactorizePivots_(0),
2514     incomingInfeasibility_(1.0),
2515     allowedInfeasibility_(10.0),
2516     automaticScale_(0),
2517     maximumPerturbationSize_(0),
2518     perturbationArray_(NULL),
2519     baseModel_(NULL)
2520#ifdef ABC_INHERIT
2521     ,abcSimplex_(NULL),
2522     abcState_(0)
2523#endif
2524{
2525     int i;
2526     for (i = 0; i < 6; i++) {
2527          rowArray_[i] = NULL;
2528          columnArray_[i] = NULL;
2529     }
2530     for (i = 0; i < 4; i++) {
2531          spareIntArray_[i] = 0;
2532          spareDoubleArray_[i] = 0.0;
2533     }
2534     saveStatus_ = NULL;
2535     // get an empty factorization so we can set tolerances etc
2536     getEmptyFactorization();
2537     // say Steepest pricing
2538     dualRowPivot_ = new ClpDualRowSteepest();
2539     // say Steepest pricing
2540     primalColumnPivot_ = new ClpPrimalColumnSteepest();
2541     solveType_ = 1; // say simplex based life form
2542
2543}
2544// Assignment operator. This copies the data
2545ClpSimplex &
2546ClpSimplex::operator=(const ClpSimplex & rhs)
2547{
2548     if (this != &rhs) {
2549          gutsOfDelete(0);
2550          delete nonLinearCost_;
2551          nonLinearCost_ = NULL;
2552          ClpModel::operator=(rhs);
2553          gutsOfCopy(rhs);
2554     }
2555     return *this;
2556}
2557void
2558ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
2559{
2560     assert (numberRows_ == rhs.numberRows_);
2561     assert (numberColumns_ == rhs.numberColumns_);
2562     numberExtraRows_ = rhs.numberExtraRows_;
2563     maximumBasic_ = rhs.maximumBasic_;
2564     dontFactorizePivots_ = rhs.dontFactorizePivots_;
2565     int numberRows2 = numberRows_ + numberExtraRows_;
2566     moreSpecialOptions_ = rhs.moreSpecialOptions_;
2567     if ((whatsChanged_ & 1) != 0) {
2568          int numberTotal = numberColumns_ + numberRows2;
2569          if ((specialOptions_ & 65536) != 0 && maximumRows_ >= 0) {
2570               assert (maximumInternalRows_ >= numberRows2);
2571               assert (maximumInternalColumns_ >= numberColumns_);
2572               numberTotal = 2 * (maximumInternalColumns_ + maximumInternalRows_);
2573          }
2574          lower_ = ClpCopyOfArray(rhs.lower_, numberTotal);
2575          rowLowerWork_ = lower_ + numberColumns_;
2576          columnLowerWork_ = lower_;
2577          upper_ = ClpCopyOfArray(rhs.upper_, numberTotal);
2578          rowUpperWork_ = upper_ + numberColumns_;
2579          columnUpperWork_ = upper_;
2580          cost_ = ClpCopyOfArray(rhs.cost_, numberTotal);
2581          objectiveWork_ = cost_;
2582          rowObjectiveWork_ = cost_ + numberColumns_;
2583          dj_ = ClpCopyOfArray(rhs.dj_, numberTotal);
2584          if (dj_) {
2585               reducedCostWork_ = dj_;
2586               rowReducedCost_ = dj_ + numberColumns_;
2587          }
2588          solution_ = ClpCopyOfArray(rhs.solution_, numberTotal);
2589          if (solution_) {
2590               columnActivityWork_ = solution_;
2591               rowActivityWork_ = solution_ + numberColumns_;
2592          }
2593          if (rhs.pivotVariable_) {
2594               pivotVariable_ = new int[numberRows2];
2595               CoinMemcpyN ( rhs.pivotVariable_, numberRows2 , pivotVariable_);
2596          } else {
2597               pivotVariable_ = NULL;
2598          }
2599          savedSolution_ = ClpCopyOfArray(rhs.savedSolution_, numberTotal);
2600          int i;
2601          for (i = 0; i < 6; i++) {
2602               rowArray_[i] = NULL;
2603               if (rhs.rowArray_[i])
2604                    rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
2605               columnArray_[i] = NULL;
2606               if (rhs.columnArray_[i])
2607                    columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
2608          }
2609          if (rhs.saveStatus_) {
2610               saveStatus_ = ClpCopyOfArray( rhs.saveStatus_, numberTotal);
2611          }
2612     } else {
2613          lower_ = NULL;
2614          rowLowerWork_ = NULL;
2615          columnLowerWork_ = NULL;
2616          upper_ = NULL;
2617          rowUpperWork_ = NULL;
2618          columnUpperWork_ = NULL;
2619          cost_ = NULL;
2620          objectiveWork_ = NULL;
2621          rowObjectiveWork_ = NULL;
2622          dj_ = NULL;
2623          reducedCostWork_ = NULL;
2624          rowReducedCost_ = NULL;
2625          solution_ = NULL;
2626          columnActivityWork_ = NULL;
2627          rowActivityWork_ = NULL;
2628          pivotVariable_ = NULL;
2629          savedSolution_ = NULL;
2630          int i;
2631          for (i = 0; i < 6; i++) {
2632               rowArray_[i] = NULL;
2633               columnArray_[i] = NULL;
2634          }
2635          saveStatus_ = NULL;
2636     }
2637     if (rhs.factorization_) {
2638          setFactorization(*rhs.factorization_);
2639     } else {
2640          delete factorization_;
2641          factorization_ = NULL;
2642     }
2643     bestPossibleImprovement_ = rhs.bestPossibleImprovement_;
2644     columnPrimalSequence_ = rhs.columnPrimalSequence_;
2645     zeroTolerance_ = rhs.zeroTolerance_;
2646     rowPrimalSequence_ = rhs.rowPrimalSequence_;
2647     bestObjectiveValue_ = rhs.bestObjectiveValue_;
2648     baseIteration_ = rhs.baseIteration_;
2649     vectorMode_ = rhs.vectorMode_;
2650     primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
2651     largeValue_ = rhs.largeValue_;
2652     largestPrimalError_ = rhs.largestPrimalError_;
2653     largestDualError_ = rhs.largestDualError_;
2654     alphaAccuracy_ = rhs.alphaAccuracy_;
2655     dualBound_ = rhs.dualBound_;
2656     alpha_ = rhs.alpha_;
2657     theta_ = rhs.theta_;
2658     lowerIn_ = rhs.lowerIn_;
2659     valueIn_ = rhs.valueIn_;
2660     upperIn_ = rhs.upperIn_;
2661     dualIn_ = rhs.dualIn_;
2662     sequenceIn_ = rhs.sequenceIn_;
2663     directionIn_ = rhs.directionIn_;
2664     lowerOut_ = rhs.lowerOut_;
2665     valueOut_ = rhs.valueOut_;
2666     upperOut_ = rhs.upperOut_;
2667     dualOut_ = rhs.dualOut_;
2668     sequenceOut_ = rhs.sequenceOut_;
2669     directionOut_ = rhs.directionOut_;
2670     pivotRow_ = rhs.pivotRow_;
2671     lastGoodIteration_ = rhs.lastGoodIteration_;
2672     numberRefinements_ = rhs.numberRefinements_;
2673     dualTolerance_ = rhs.dualTolerance_;
2674     primalTolerance_ = rhs.primalTolerance_;
2675     sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
2676     numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
2677     numberDualInfeasibilitiesWithoutFree_ =
2678          rhs.numberDualInfeasibilitiesWithoutFree_;
2679     sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
2680     numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
2681     dualRowPivot_ = rhs.dualRowPivot_->clone(true);
2682     dualRowPivot_->setModel(this);
2683     primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
2684     primalColumnPivot_->setModel(this);
2685     numberTimesOptimal_ = rhs.numberTimesOptimal_;
2686     disasterArea_ = NULL;
2687     changeMade_ = rhs.changeMade_;
2688     algorithm_ = rhs.algorithm_;
2689     forceFactorization_ = rhs.forceFactorization_;
2690     perturbation_ = rhs.perturbation_;
2691     infeasibilityCost_ = rhs.infeasibilityCost_;
2692     lastBadIteration_ = rhs.lastBadIteration_;
2693     lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
2694     numberFake_ = rhs.numberFake_;
2695     numberChanged_ = rhs.numberChanged_;
2696     progressFlag_ = rhs.progressFlag_;
2697     firstFree_ = rhs.firstFree_;
2698     incomingInfeasibility_ = rhs.incomingInfeasibility_;
2699     allowedInfeasibility_ = rhs.allowedInfeasibility_;
2700     automaticScale_ = rhs.automaticScale_;
2701#ifdef ABC_INHERIT
2702     abcSimplex_ = NULL;
2703     abcState_ = rhs.abcState_;
2704#endif
2705     maximumPerturbationSize_ = rhs.maximumPerturbationSize_;
2706     if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) {
2707          perturbationArray_ = CoinCopyOfArray(rhs.perturbationArray_,
2708                                               maximumPerturbationSize_);
2709     } else {
2710          maximumPerturbationSize_ = 0;
2711          perturbationArray_ = NULL;
2712     }
2713     if (rhs.baseModel_) {
2714          baseModel_ = new ClpSimplex(*rhs.baseModel_);
2715     } else {
2716          baseModel_ = NULL;
2717     }
2718     progress_ = rhs.progress_;
2719     for (int i = 0; i < 4; i++) {
2720          spareIntArray_[i] = rhs.spareIntArray_[i];
2721          spareDoubleArray_[i] = rhs.spareDoubleArray_[i];
2722     }
2723     sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
2724     sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
2725     acceptablePivot_ = rhs.acceptablePivot_;
2726     if (rhs.nonLinearCost_ != NULL)
2727          nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
2728     else
2729          nonLinearCost_ = NULL;
2730     solveType_ = rhs.solveType_;
2731     eventHandler_->setSimplex(this);
2732}
2733// type == 0 do everything, most + pivot data, 2 factorization data as well
2734void
2735ClpSimplex::gutsOfDelete(int type)
2736{
2737     if (!type || (specialOptions_ & 65536) == 0) {
2738          maximumInternalColumns_ = -1;
2739          maximumInternalRows_ = -1;
2740          delete [] lower_;
2741          lower_ = NULL;
2742          rowLowerWork_ = NULL;
2743          columnLowerWork_ = NULL;
2744          delete [] upper_;
2745          upper_ = NULL;
2746          rowUpperWork_ = NULL;
2747          columnUpperWork_ = NULL;
2748          delete [] cost_;
2749          cost_ = NULL;
2750          objectiveWork_ = NULL;
2751          rowObjectiveWork_ = NULL;
2752          delete [] dj_;
2753          dj_ = NULL;
2754          reducedCostWork_ = NULL;
2755          rowReducedCost_ = NULL;
2756          delete [] solution_;
2757          solution_ = NULL;
2758          rowActivityWork_ = NULL;
2759          columnActivityWork_ = NULL;
2760          delete [] savedSolution_;
2761          savedSolution_ = NULL;
2762     }
2763     if ((specialOptions_ & 2) == 0) {
2764          delete nonLinearCost_;
2765          nonLinearCost_ = NULL;
2766     }
2767     int i;
2768     if ((specialOptions_ & 65536) == 0) {
2769          for (i = 0; i < 6; i++) {
2770               delete rowArray_[i];
2771               rowArray_[i] = NULL;
2772               delete columnArray_[i];
2773               columnArray_[i] = NULL;
2774          }
2775     }
2776     delete [] saveStatus_;
2777     saveStatus_ = NULL;
2778     if (type != 1) {
2779          delete rowCopy_;
2780          rowCopy_ = NULL;
2781     }
2782     if (!type) {
2783          // delete everything
2784          setEmptyFactorization();
2785          delete [] pivotVariable_;
2786          pivotVariable_ = NULL;
2787          delete dualRowPivot_;
2788          dualRowPivot_ = NULL;
2789          delete primalColumnPivot_;
2790          primalColumnPivot_ = NULL;
2791          delete baseModel_;
2792          baseModel_ = NULL;
2793          delete [] perturbationArray_;
2794          perturbationArray_ = NULL;
2795          maximumPerturbationSize_ = 0;
2796     } else {
2797          // delete any size information in methods
2798          if (type > 1) {
2799               //assert (factorization_);
2800               if (factorization_)
2801                    factorization_->clearArrays();
2802               delete [] pivotVariable_;
2803               pivotVariable_ = NULL;
2804          }
2805          dualRowPivot_->clearArrays();
2806          primalColumnPivot_->clearArrays();
2807     }
2808}
2809// This sets largest infeasibility and most infeasible
2810void
2811ClpSimplex::checkPrimalSolution(const double * rowActivities,
2812                                const double * columnActivities)
2813{
2814     double * solution;
2815     int iRow, iColumn;
2816
2817     objectiveValue_ = 0.0;
2818     // now look at primal solution
2819     solution = rowActivityWork_;
2820     sumPrimalInfeasibilities_ = 0.0;
2821     numberPrimalInfeasibilities_ = 0;
2822     double primalTolerance = primalTolerance_;
2823     double relaxedTolerance = primalTolerance_;
2824     // we can't really trust infeasibilities if there is primal error
2825     double error = CoinMin(1.0e-2, largestPrimalError_);
2826     // allow tolerance at least slightly bigger than standard
2827     relaxedTolerance = relaxedTolerance +  error;
2828     sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2829     for (iRow = 0; iRow < numberRows_; iRow++) {
2830          //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
2831          double infeasibility = 0.0;
2832          objectiveValue_ += solution[iRow] * rowObjectiveWork_[iRow];
2833          if (solution[iRow] > rowUpperWork_[iRow]) {
2834               infeasibility = solution[iRow] - rowUpperWork_[iRow];
2835          } else if (solution[iRow] < rowLowerWork_[iRow]) {
2836               infeasibility = rowLowerWork_[iRow] - solution[iRow];
2837          }
2838          if (infeasibility > primalTolerance) {
2839               sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2840               if (infeasibility > relaxedTolerance)
2841                    sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2842               numberPrimalInfeasibilities_ ++;
2843          }
2844          infeasibility = fabs(rowActivities[iRow] - solution[iRow]);
2845     }
2846     // Check any infeasibilities from dynamic rows
2847     matrix_->primalExpanded(this, 2);
2848     solution = columnActivityWork_;
2849     if (!matrix_->rhsOffset(this)) {
2850          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2851               //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2852               double infeasibility = 0.0;
2853               objectiveValue_ += objectiveWork_[iColumn] * solution[iColumn];
2854               if (solution[iColumn] > columnUpperWork_[iColumn]) {
2855                    infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2856               } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2857                    infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2858               }
2859               if (infeasibility > primalTolerance) {
2860                    sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2861                    if (infeasibility > relaxedTolerance)
2862                         sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2863                    numberPrimalInfeasibilities_ ++;
2864               }
2865               infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2866          }
2867     } else {
2868          // as we are using effective rhs we only check basics
2869          // But we do need to get objective
2870          objectiveValue_ += innerProduct(objectiveWork_, numberColumns_, solution);
2871          for (int j = 0; j < numberRows_; j++) {
2872               int iColumn = pivotVariable_[j];
2873               //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2874               double infeasibility = 0.0;
2875               if (solution[iColumn] > columnUpperWork_[iColumn]) {
2876                    infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2877               } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2878                    infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2879               }
2880               if (infeasibility > primalTolerance) {
2881                    sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2882                    if (infeasibility > relaxedTolerance)
2883                         sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2884                    numberPrimalInfeasibilities_ ++;
2885               }
2886               infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2887          }
2888     }
2889     objectiveValue_ += objective_->nonlinearOffset();
2890     objectiveValue_ /= (objectiveScale_ * rhsScale_);
2891}
2892void
2893ClpSimplex::checkDualSolution()
2894{
2895
2896     int iRow, iColumn;
2897     sumDualInfeasibilities_ = 0.0;
2898     numberDualInfeasibilities_ = 0;
2899     numberDualInfeasibilitiesWithoutFree_ = 0;
2900     if (matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) {
2901          // pretend we found dual infeasibilities
2902          sumOfRelaxedDualInfeasibilities_ = 1.0;
2903          sumDualInfeasibilities_ = 1.0;
2904          numberDualInfeasibilities_ = 1;
2905          return;
2906     }
2907     int firstFreePrimal = -1;
2908     int firstFreeDual = -1;
2909     int numberSuperBasicWithDj = 0;
2910     bestPossibleImprovement_ = 0.0;
2911     // we can't really trust infeasibilities if there is dual error
2912     double error = CoinMin(1.0e-2, largestDualError_);
2913     // allow tolerance at least slightly bigger than standard
2914     double relaxedTolerance = dualTolerance_ +  error;
2915     // allow bigger tolerance for possible improvement
2916     double possTolerance = 5.0 * relaxedTolerance;
2917     sumOfRelaxedDualInfeasibilities_ = 0.0;
2918
2919     // Check any djs from dynamic rows
2920     matrix_->dualExpanded(this, NULL, NULL, 3);
2921     numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_;
2922     objectiveValue_ = 0.0;
2923     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2924          objectiveValue_ += objectiveWork_[iColumn] * columnActivityWork_[iColumn];
2925          if (getColumnStatus(iColumn) != basic && !flagged(iColumn)) {
2926               // not basic
2927               double distanceUp = columnUpperWork_[iColumn] -
2928                                   columnActivityWork_[iColumn];
2929               double distanceDown = columnActivityWork_[iColumn] -
2930                                     columnLowerWork_[iColumn];
2931               if (distanceUp > primalTolerance_) {
2932                    double value = reducedCostWork_[iColumn];
2933                    // Check if "free"
2934                    if (distanceDown > primalTolerance_) {
2935                         if (fabs(value) > 1.0e2 * relaxedTolerance) {
2936                              numberSuperBasicWithDj++;
2937                              if (firstFreeDual < 0)
2938                                   firstFreeDual = iColumn;
2939                         }
2940                         if (firstFreePrimal < 0)
2941                              firstFreePrimal = iColumn;
2942                    }
2943                    // should not be negative
2944                    if (value < 0.0) {
2945                         value = - value;
2946                         if (value > dualTolerance_) {
2947                              if (getColumnStatus(iColumn) != isFree) {
2948                                   numberDualInfeasibilitiesWithoutFree_ ++;
2949                                   sumDualInfeasibilities_ += value - dualTolerance_;
2950                                   if (value > possTolerance)
2951                                        bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2952                                   if (value > relaxedTolerance)
2953                                        sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2954                                   numberDualInfeasibilities_ ++;
2955                              } else {
2956                                   // free so relax a lot
2957                                   value *= 0.01;
2958                                   if (value > dualTolerance_) {
2959                                        sumDualInfeasibilities_ += value - dualTolerance_;
2960                                        if (value > possTolerance)
2961                                             bestPossibleImprovement_ = 1.0e100;
2962                                        if (value > relaxedTolerance)
2963                                             sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2964                                        numberDualInfeasibilities_ ++;
2965                                   }
2966                              }
2967                         }
2968                    }
2969               }
2970               if (distanceDown > primalTolerance_) {
2971                    double value = reducedCostWork_[iColumn];
2972                    // should not be positive
2973                    if (value > 0.0) {
2974                         if (value > dualTolerance_) {
2975                              sumDualInfeasibilities_ += value - dualTolerance_;
2976                              if (value > possTolerance)
2977                                   bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2978                              if (value > relaxedTolerance)
2979                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2980                              numberDualInfeasibilities_ ++;
2981                              if (getColumnStatus(iColumn) != isFree)
2982                                   numberDualInfeasibilitiesWithoutFree_ ++;
2983                              // maybe we can make feasible by increasing tolerance
2984                         }
2985                    }
2986               }
2987          }
2988     }
2989     for (iRow = 0; iRow < numberRows_; iRow++) {
2990          objectiveValue_ += rowActivityWork_[iRow] * rowObjectiveWork_[iRow];
2991          if (getRowStatus(iRow) != basic && !flagged(iRow + numberColumns_)) {
2992               // not basic
2993               double distanceUp = rowUpperWork_[iRow] - rowActivityWork_[iRow];
2994               double distanceDown = rowActivityWork_[iRow] - rowLowerWork_[iRow];
2995               if (distanceUp > primalTolerance_) {
2996                    double value = rowReducedCost_[iRow];
2997                    // Check if "free"
2998                    if (distanceDown > primalTolerance_) {
2999                         if (fabs(value) > 1.0e2 * relaxedTolerance) {
3000                              numberSuperBasicWithDj++;
3001                              if (firstFreeDual < 0)
3002                                   firstFreeDual = iRow + numberColumns_;
3003                         }
3004                         if (firstFreePrimal < 0)
3005                              firstFreePrimal = iRow + numberColumns_;
3006                    }
3007                    // should not be negative
3008                    if (value < 0.0) {
3009                         value = - value;
3010                         if (value > dualTolerance_) {
3011                              sumDualInfeasibilities_ += value - dualTolerance_;
3012                              if (value > possTolerance)
3013                                   bestPossibleImprovement_ += value * CoinMin(distanceUp, 1.0e10);
3014                              if (value > relaxedTolerance)
3015                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
3016                              numberDualInfeasibilities_ ++;
3017                              if (getRowStatus(iRow) != isFree)
3018                                   numberDualInfeasibilitiesWithoutFree_ ++;
3019                         }
3020                    }
3021               }
3022               if (distanceDown > primalTolerance_) {
3023                    double value = rowReducedCost_[iRow];
3024                    // should not be positive
3025                    if (value > 0.0) {
3026                         if (value > dualTolerance_) {
3027                              sumDualInfeasibilities_ += value - dualTolerance_;
3028                              if (value > possTolerance)
3029                                   bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
3030                              if (value > relaxedTolerance)
3031                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
3032                              numberDualInfeasibilities_ ++;
3033                              if (getRowStatus(iRow) != isFree)
3034                                   numberDualInfeasibilitiesWithoutFree_ ++;
3035                              // maybe we can make feasible by increasing tolerance
3036                         }
3037                    }
3038               }
3039          }
3040     }
3041     if (algorithm_ < 0 && firstFreeDual >= 0) {
3042          // dual
3043          firstFree_ = firstFreeDual;
3044     } else if (numberSuperBasicWithDj ||
3045                (progress_.lastIterationNumber(0) <= 0)) {
3046          firstFree_ = firstFreePrimal;
3047     }
3048     objectiveValue_ += objective_->nonlinearOffset();
3049     objectiveValue_ /= (objectiveScale_ * rhsScale_);
3050}
3051/* This sets sum and number of infeasibilities (Dual and Primal) */
3052void
3053ClpSimplex::checkBothSolutions()
3054{
3055     if ((matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) ||
3056               matrix_->rhsOffset(this)) {
3057          // Say may be free or superbasic
3058          moreSpecialOptions_ &= ~8;
3059          // old way
3060          checkPrimalSolution(rowActivityWork_, columnActivityWork_);
3061          checkDualSolution();
3062          return;
3063     }
3064     int iSequence;
3065     assert (dualTolerance_ > 0.0 && dualTolerance_ < 1.0e10);
3066     assert (primalTolerance_ > 0.0 && primalTolerance_ < 1.0e10);
3067     objectiveValue_ = 0.0;
3068     sumPrimalInfeasibilities_ = 0.0;
3069     numberPrimalInfeasibilities_ = 0;
3070     double primalTolerance = primalTolerance_;
3071     double relaxedToleranceP = primalTolerance_;
3072     // we can't really trust infeasibilities if there is primal error
3073     double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_,0.0*primalTolerance_));
3074     // allow tolerance at least slightly bigger than standard
3075     relaxedToleranceP = relaxedToleranceP +  error;
3076     sumOfRelaxedPrimalInfeasibilities_ = 0.0;
3077     sumDualInfeasibilities_ = 0.0;
3078     numberDualInfeasibilities_ = 0;
3079     double dualTolerance = dualTolerance_;
3080     double relaxedToleranceD = dualTolerance;
3081     // we can't really trust infeasibilities if there is dual error
3082     error = CoinMin(1.0e-2, CoinMax(largestDualError_,5.0*dualTolerance_));
3083     // allow tolerance at least slightly bigger than standard
3084     relaxedToleranceD = relaxedToleranceD +  error;
3085     // allow bigger tolerance for possible improvement
3086     double possTolerance = 5.0 * relaxedToleranceD;
3087     sumOfRelaxedDualInfeasibilities_ = 0.0;
3088     bestPossibleImprovement_ = 0.0;
3089
3090     // Check any infeasibilities from dynamic rows
3091     matrix_->primalExpanded(this, 2);
3092     // Check any djs from dynamic rows
3093     matrix_->dualExpanded(this, NULL, NULL, 3);
3094     int numberDualInfeasibilitiesFree = 0;
3095     int firstFreePrimal = -1;
3096     int firstFreeDual = -1;
3097     int numberSuperBasicWithDj = 0;
3098
3099     int numberTotal = numberRows_ + numberColumns_;
3100     // Say no free or superbasic
3101     moreSpecialOptions_ |= 8;
3102     //#define PRINT_INFEAS
3103#ifdef PRINT_INFEAS
3104     int seqInf[10];
3105#endif
3106     for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3107          double value = solution_[iSequence];
3108#ifdef COIN_DEBUG
3109          if (fabs(value) > 1.0e20)
3110               printf("%d values %g %g %g - status %d\n", iSequence, lower_[iSequence],
3111                      solution_[iSequence], upper_[iSequence], status_[iSequence]);
3112#endif
3113          objectiveValue_ += value * cost_[iSequence];
3114          double distanceUp = upper_[iSequence] - value;
3115          double distanceDown = value - lower_[iSequence];
3116          if (distanceUp < -primalTolerance) {
3117               double infeasibility = -distanceUp;
3118               sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
3119               if (infeasibility > relaxedToleranceP)
3120                    sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
3121#ifdef PRINT_INFEAS
3122               if (numberPrimalInfeasibilities_<10) {
3123                 seqInf[numberPrimalInfeasibilities_]=iSequence;
3124               }
3125#endif
3126               numberPrimalInfeasibilities_ ++;
3127          } else if (distanceDown < -primalTolerance) {
3128               double infeasibility = -distanceDown;
3129               sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
3130               if (infeasibility > relaxedToleranceP)
3131                    sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
3132#ifdef PRINT_INFEAS
3133               if (numberPrimalInfeasibilities_<10) {
3134                 seqInf[numberPrimalInfeasibilities_]=iSequence;
3135               }
3136#endif
3137               numberPrimalInfeasibilities_ ++;
3138          } else {
3139               // feasible (so could be free)
3140               if (getStatus(iSequence) != basic && !flagged(iSequence)) {
3141                    // not basic
3142                    double djValue = dj_[iSequence];
3143                    if (distanceDown < primalTolerance) {
3144                         if (distanceUp > primalTolerance && djValue < -dualTolerance) {
3145                              sumDualInfeasibilities_ -= djValue + dualTolerance;
3146                              if (djValue < -possTolerance)
3147                                   bestPossibleImprovement_ -= distanceUp * djValue;
3148                              if (djValue < -relaxedToleranceD)
3149                                   sumOfRelaxedDualInfeasibilities_ -= djValue + relaxedToleranceD;
3150                              numberDualInfeasibilities_ ++;
3151                         }
3152                    } else if (distanceUp < primalTolerance) {
3153                         if (djValue > dualTolerance) {
3154                              sumDualInfeasibilities_ += djValue - dualTolerance;
3155                              if (djValue > possTolerance)
3156                                   bestPossibleImprovement_ += distanceDown * djValue;
3157                              if (djValue > relaxedToleranceD)
3158                                   sumOfRelaxedDualInfeasibilities_ += djValue - relaxedToleranceD;
3159                              numberDualInfeasibilities_ ++;
3160                         }
3161                    } else {
3162                         // may be free
3163                         // Say free or superbasic
3164                         moreSpecialOptions_ &= ~8;
3165                         djValue *= 0.01;
3166                         if (fabs(djValue) > dualTolerance) {
3167                              if (getStatus(iSequence) == isFree)
3168                                   numberDualInfeasibilitiesFree++;
3169                              sumDualInfeasibilities_ += fabs(djValue) - dualTolerance;
3170                              bestPossibleImprovement_ = 1.0e100;
3171                              numberDualInfeasibilities_ ++;
3172                              if (fabs(djValue) > relaxedToleranceD) {
3173                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedToleranceD;
3174                                   numberSuperBasicWithDj++;
3175                                   if (firstFreeDual < 0)
3176                                        firstFreeDual = iSequence;
3177                              }
3178                         }
3179                         if (firstFreePrimal < 0)
3180                              firstFreePrimal = iSequence;
3181                    }
3182               }
3183          }
3184     }
3185     objectiveValue_ += objective_->nonlinearOffset();
3186     objectiveValue_ /= (objectiveScale_ * rhsScale_);
3187     numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ -
3188                                             numberDualInfeasibilitiesFree;
3189#ifdef PRINT_INFEAS
3190     if (numberPrimalInfeasibilities_<=10) {
3191       printf("---------------start-----------\n");
3192       if (!rowScale_) {
3193         for (int i=0;i<numberPrimalInfeasibilities_;i++) {
3194           int iSeq = seqInf[i];
3195           double infeas;
3196           if (solution_[iSeq]<lower_[iSeq])
3197             infeas = lower_[iSeq]-solution_[iSeq];
3198           else
3199             infeas = solution_[iSeq]-upper_[iSeq];
3200           if (iSeq<numberColumns_) {
3201             printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g\n",
3202                    iSeq,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas);
3203           } else {
3204             printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g\n",
3205                    iSeq-numberColumns_,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas);
3206           }
3207         }
3208       } else {
3209         for (int i=0;i<numberPrimalInfeasibilities_;i++) {
3210           int iSeq = seqInf[i];
3211           double infeas;
3212           if (solution_[iSeq]<lower_[iSeq])
3213             infeas = lower_[iSeq]-solution_[iSeq];
3214           else
3215             infeas = solution_[iSeq]-upper_[iSeq];
3216           double unscaled = infeas;
3217           if (iSeq<numberColumns_) {
3218             unscaled *= columnScale_[iSeq];
3219             printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3220                    iSeq,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas,unscaled);
3221           } else {
3222             unscaled /= rowScale_[iSeq-numberColumns_];
3223             printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3224                    iSeq-numberColumns_,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas,unscaled);
3225           }
3226         }
3227       }
3228     }
3229#endif
3230     if (algorithm_ < 0 && firstFreeDual >= 0) {
3231          // dual
3232          firstFree_ = firstFreeDual;
3233     } else if (numberSuperBasicWithDj ||
3234                (progress_.lastIterationNumber(0) <= 0)) {
3235          firstFree_ = firstFreePrimal;
3236     }
3237}
3238/* Adds multiple of a column into an array */
3239void
3240ClpSimplex::add(double * array,
3241                int sequence, double multiplier) const
3242{
3243     if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3244          //slack
3245          array [sequence-numberColumns_] -= multiplier;
3246     } else {
3247          // column
3248          matrix_->add(this, array, sequence, multiplier);
3249     }
3250}
3251/*
3252  Unpacks one column of the matrix into indexed array
3253*/
3254void
3255ClpSimplex::unpack(CoinIndexedVector * rowArray) const
3256{
3257     rowArray->clear();
3258     if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3259          //slack
3260          rowArray->insert(sequenceIn_ - numberColumns_, -1.0);
3261     } else {
3262          // column
3263          matrix_->unpack(this, rowArray, sequenceIn_);
3264     }
3265}
3266void
3267ClpSimplex::unpack(CoinIndexedVector * rowArray, int sequence) const
3268{
3269     rowArray->clear();
3270     if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3271          //slack
3272          rowArray->insert(sequence - numberColumns_, -1.0);
3273     } else {
3274          // column
3275          matrix_->unpack(this, rowArray, sequence);
3276     }
3277}
3278/*
3279  Unpacks one column of the matrix into indexed array
3280*/
3281void
3282ClpSimplex::unpackPacked(CoinIndexedVector * rowArray)
3283{
3284     rowArray->clear();
3285     if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3286          //slack
3287          int * index = rowArray->getIndices();
3288          double * array = rowArray->denseVector();
3289          array[0] = -1.0;
3290          index[0] = sequenceIn_ - numberColumns_;
3291          rowArray->setNumElements(1);
3292          rowArray->setPackedMode(true);
3293     } else {
3294          // column
3295          matrix_->unpackPacked(this, rowArray, sequenceIn_);
3296     }
3297}
3298void
3299ClpSimplex::unpackPacked(CoinIndexedVector * rowArray, int sequence)
3300{
3301     rowArray->clear();
3302     if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3303          //slack
3304          int * index = rowArray->getIndices();
3305          double * array = rowArray->denseVector();
3306          array[0] = -1.0;
3307          index[0] = sequence - numberColumns_;
3308          rowArray->setNumElements(1);
3309          rowArray->setPackedMode(true);
3310     } else {
3311          // column
3312          matrix_->unpackPacked(this, rowArray, sequence);
3313     }
3314}
3315//static int x_gaps[4]={0,0,0,0};
3316//static int scale_times[]={0,0,0,0};
3317bool
3318ClpSimplex::createRim(int what, bool makeRowCopy, int startFinishOptions)
3319{
3320     bool goodMatrix = true;
3321     int saveLevel = handler_->logLevel();
3322     spareIntArray_[0] = 0;
3323     if (!matrix_->canGetRowCopy())
3324          makeRowCopy = false; // switch off row copy if can't produce
3325     // Arrays will be there and correct size unless what is 63
3326     bool newArrays = (what == 63);
3327     // We may be restarting with same size
3328     bool keepPivots = false;
3329     if (startFinishOptions == -1) {
3330          startFinishOptions = 0;
3331          keepPivots = true;
3332     }
3333     bool oldMatrix = ((startFinishOptions & 4) != 0 && (whatsChanged_ & 1) != 0);
3334     if (what == 63) {
3335          pivotRow_ = -1;
3336          if (!status_)
3337               createStatus();
3338          if (oldMatrix)
3339               newArrays = false;
3340          if (problemStatus_ == 10) {
3341               handler_->setLogLevel(0); // switch off messages
3342               if (rowArray_[0]) {
3343                    // stuff is still there
3344                    oldMatrix = true;
3345                    newArrays = false;
3346                    keepPivots = true;
3347                    for (int iRow = 0; iRow < 4; iRow++) {
3348                         rowArray_[iRow]->clear();
3349                    }
3350                    for (int iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
3351                         columnArray_[iColumn]->clear();
3352                    }
3353               }
3354          } else if (factorization_) {
3355               // match up factorization messages
3356               if (handler_->logLevel() < 3)
3357                    factorization_->messageLevel(0);
3358               else
3359                    factorization_->messageLevel(CoinMax(3, factorization_->messageLevel()));
3360               /* Faster to keep pivots rather than re-scan matrix.  Matrix may have changed
3361                  i.e. oldMatrix false but okay as long as same number rows and status array exists
3362               */
3363               if ((startFinishOptions & 2) != 0 && factorization_->numberRows() == numberRows_ && status_)
3364                    keepPivots = true;
3365          }
3366          numberExtraRows_ = matrix_->generalExpanded(this, 2, maximumBasic_);
3367          if (numberExtraRows_ && newArrays) {
3368               // make sure status array large enough
3369               assert (status_);
3370               int numberOld = numberRows_ + numberColumns_;
3371               int numberNew = numberRows_ + numberColumns_ + numberExtraRows_;
3372               unsigned char * newStatus = new unsigned char [numberNew];
3373               memset(newStatus + numberOld, 0, numberExtraRows_);
3374               CoinMemcpyN(status_, numberOld, newStatus);
3375               delete [] status_;
3376               status_ = newStatus;
3377          }
3378     }
3379     int numberRows2 = numberRows_ + numberExtraRows_;
3380     int numberTotal = numberRows2 + numberColumns_;
3381     if ((specialOptions_ & 65536) != 0) {
3382          assert (!numberExtraRows_);
3383          if (!cost_ || numberRows2 > maximumInternalRows_ ||
3384                    numberColumns_ > maximumInternalColumns_) {
3385               newArrays = true;
3386               keepPivots = false;
3387               COIN_DETAIL_PRINT(printf("createrim a %d rows, %d maximum rows %d maxinternal\n",
3388                                        numberRows_, maximumRows_, maximumInternalRows_));
3389               int oldMaximumRows = maximumInternalRows_;
3390               int oldMaximumColumns = maximumInternalColumns_;
3391               if (cost_) {
3392                    if (numberRows2 > maximumInternalRows_)
3393                         maximumInternalRows_ = numberRows2;
3394                    if (numberColumns_ > maximumInternalColumns_)
3395                         maximumInternalColumns_ = numberColumns_;
3396               } else {
3397                    maximumInternalRows_ = numberRows2;
3398                    maximumInternalColumns_ = numberColumns_;
3399               }
3400               //maximumRows_=CoinMax(maximumInternalRows_,maximumRows_);
3401               //maximumColumns_=CoinMax(maximumInternalColumns_,maximumColumns_);
3402               assert(maximumInternalRows_ == maximumRows_);
3403               assert(maximumInternalColumns_ == maximumColumns_);
3404               COIN_DETAIL_PRINT(printf("createrim b %d rows, %d maximum rows, %d maxinternal\n",
3405                                        numberRows_, maximumRows_, maximumInternalRows_));
3406               int numberTotal2 = (maximumInternalRows_ + maximumInternalColumns_) * 2;
3407               delete [] cost_;
3408               cost_ = new double[numberTotal2];
3409               delete [] lower_;
3410               delete [] upper_;
3411               lower_ = new double[numberTotal2];
3412               upper_ = new double[numberTotal2];
3413               delete [] dj_;
3414               dj_ = new double[numberTotal2];
3415               delete [] solution_;
3416               solution_ = new double[numberTotal2];
3417               // ***** should be non NULL but seems to be too much
3418               //printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_);
3419               if (savedRowScale_) {
3420                    assert (oldMaximumRows > 0);
3421                    double * temp;
3422                    temp = new double [4*maximumRows_];
3423                    CoinFillN(temp, 4 * maximumRows_, 1.0);
3424                    CoinMemcpyN(savedRowScale_, numberRows_, temp);
3425                    CoinMemcpyN(savedRowScale_ + oldMaximumRows, numberRows_, temp + maximumRows_);
3426                    CoinMemcpyN(savedRowScale_ + 2 * oldMaximumRows, numberRows_, temp + 2 * maximumRows_);
3427                    CoinMemcpyN(savedRowScale_ + 3 * oldMaximumRows, numberRows_, temp + 3 * maximumRows_);
3428                    delete [] savedRowScale_;
3429                    savedRowScale_ = temp;
3430                    temp = new double [4*maximumColumns_];
3431                    CoinFillN(temp, 4 * maximumColumns_, 1.0);
3432                    CoinMemcpyN(savedColumnScale_, numberColumns_, temp);
3433                    CoinMemcpyN(savedColumnScale_ + oldMaximumColumns, numberColumns_, temp + maximumColumns_);
3434                    CoinMemcpyN(savedColumnScale_ + 2 * oldMaximumColumns, numberColumns_, temp + 2 * maximumColumns_);
3435                    CoinMemcpyN(savedColumnScale_ + 3 * oldMaximumColumns, numberColumns_, temp + 3 * maximumColumns_);
3436                    delete [] savedColumnScale_;
3437                    savedColumnScale_ = temp;
3438               }
3439          }
3440     }
3441     int i;
3442     bool doSanityCheck = true;
3443     if (what == 63) {
3444          // We may want to switch stuff off for speed
3445          if ((specialOptions_ & 256) != 0)
3446               makeRowCopy = false; // no row copy
3447          if ((specialOptions_ & 128) != 0)
3448               doSanityCheck = false; // no sanity check
3449          //check matrix
3450          if (!matrix_)
3451               matrix_ = new ClpPackedMatrix();
3452          int checkType = (doSanityCheck) ? 15 : 14;
3453          if (oldMatrix)
3454               checkType = 14;
3455          bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
3456          if (inCbcOrOther)
3457               checkType -= 4; // don't check for duplicates
3458          if (!matrix_->allElementsInRange(this, smallElement_, 1.0e20, checkType)) {
3459               problemStatus_ = 4;
3460               secondaryStatus_ = 8;
3461               //goodMatrix= false;
3462               return false;
3463          }
3464          bool rowCopyIsScaled;
3465          if (makeRowCopy) {
3466               if(!oldMatrix || !rowCopy_) {
3467                    delete rowCopy_;
3468                    // may return NULL if can't give row copy
3469                    rowCopy_ = matrix_->reverseOrderedCopy();
3470                    rowCopyIsScaled = false;
3471               } else {
3472                    rowCopyIsScaled = true;
3473               }
3474          }
3475#if 0
3476          if (what == 63) {
3477               int k = rowScale_ ? 1 : 0;
3478               if (oldMatrix)
3479                    k += 2;
3480               scale_times[k]++;
3481               if ((scale_times[0] + scale_times[1] + scale_times[2] + scale_times[3]) % 1000 == 0)
3482                    printf("scale counts %d %d %d %d\n",
3483                           scale_times[0], scale_times[1], scale_times[2], scale_times[3]);
3484          }
3485#endif
3486          // do scaling if needed
3487          if (!oldMatrix && scalingFlag_ < 0) {
3488               if (scalingFlag_ < 0 && rowScale_) {
3489                    //if (handler_->logLevel()>0)
3490                    printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
3491                           scalingFlag_);
3492                    scalingFlag_ = 0;
3493               }
3494               delete [] rowScale_;
3495               delete [] columnScale_;
3496               rowScale_ = NULL;
3497               columnScale_ = NULL;
3498          }
3499          inverseRowScale_ = NULL;
3500          inverseColumnScale_ = NULL;
3501          if (scalingFlag_ > 0 &&(specialOptions_ & 65536) != 0&&
3502              rowScale_&&rowScale_==savedRowScale_) 
3503            rowScale_=NULL;
3504          if (scalingFlag_ > 0 && !rowScale_) {
3505               if ((specialOptions_ & 65536) != 0) {
3506                    assert (!rowScale_);
3507                    rowScale_ = savedRowScale_;
3508                    columnScale_ = savedColumnScale_;
3509                    // put back original
3510                    if (savedRowScale_) {
3511                         inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
3512                         inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
3513                         CoinMemcpyN(savedRowScale_ + 2 * maximumInternalRows_,
3514                                     numberRows2, savedRowScale_);
3515                         CoinMemcpyN(savedRowScale_ + 3 * maximumInternalRows_,
3516                                     numberRows2, inverseRowScale_);
3517                         CoinMemcpyN(savedColumnScale_ + 2 * maximumColumns_,
3518                                     numberColumns_, savedColumnScale_);
3519                         CoinMemcpyN(savedColumnScale_ + 3 * maximumColumns_,
3520                                     numberColumns_, inverseColumnScale_);
3521                    }
3522               }
3523               if (matrix_->scale(this,this))
3524                    scalingFlag_ = -scalingFlag_; // not scaled after all
3525               if (rowScale_ && automaticScale_) {
3526                    if (!savedRowScale_) {
3527                      inverseRowScale_ = rowScale_ + numberRows2;
3528                      inverseColumnScale_ = columnScale_ + numberColumns_;
3529                    }
3530                    // try automatic scaling
3531                    double smallestObj = 1.0e100;
3532                    double largestObj = 0.0;
3533                    double largestRhs = 0.0;
3534                    const double * obj = objective();
3535                    for (i = 0; i < numberColumns_; i++) {
3536                         double value = fabs(obj[i]);
3537                         value *= columnScale_[i];
3538                         if (value && columnLower_[i] != columnUpper_[i]) {
3539                              smallestObj = CoinMin(smallestObj, value);
3540                              largestObj = CoinMax(largestObj, value);
3541                         }
3542                         if (columnLower_[i] > 0.0 || columnUpper_[i] < 0.0) {
3543                              double scale = 1.0 * inverseColumnScale_[i];
3544                              //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3545                              if (columnLower_[i] > 0)
3546                                   largestRhs = CoinMax(largestRhs, columnLower_[i] * scale);
3547                              if (columnUpper_[i] < 0.0)
3548                                   largestRhs = CoinMax(largestRhs, -columnUpper_[i] * scale);
3549                         }
3550                    }
3551                    for (i = 0; i < numberRows_; i++) {
3552                         if (rowLower_[i] > 0.0 || rowUpper_[i] < 0.0) {
3553                              double scale = rowScale_[i];
3554                              //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3555                              if (rowLower_[i] > 0)
3556                                   largestRhs = CoinMax(largestRhs, rowLower_[i] * scale);
3557                              if (rowUpper_[i] < 0.0)
3558                                   largestRhs = CoinMax(largestRhs, -rowUpper_[i] * scale);
3559                         }
3560                    }
3561                    COIN_DETAIL_PRINT(printf("small obj %g, large %g - rhs %g\n", smallestObj, largestObj, largestRhs));
3562                    bool scalingDone = false;
3563                    // look at element range
3564                    double smallestNegative;
3565                    double largestNegative;
3566                    double smallestPositive;
3567                    double largestPositive;
3568                    matrix_->rangeOfElements(smallestNegative, largestNegative,
3569                                             smallestPositive, largestPositive);
3570                    smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
3571                    largestPositive = CoinMax(fabs(largestNegative), largestPositive);
3572                    if (largestObj) {
3573                         double ratio = largestObj / smallestObj;
3574                         double scale = 1.0;
3575                         if (ratio < 1.0e8) {
3576                              // reasonable
3577                              if (smallestObj < 1.0e-4) {
3578                                   // may as well scale up
3579                                   scalingDone = true;
3580                                   scale = 1.0e-3 / smallestObj;
3581                              } else if (largestObj < 1.0e6 || (algorithm_ > 0 && largestObj < 1.0e-4 * infeasibilityCost_)) {
3582                                   //done=true;
3583                              } else {
3584                                   scalingDone = true;
3585                                   if (algorithm_ < 0) {
3586                                        scale = 1.0e6 / largestObj;
3587                                   } else {
3588                                        scale = CoinMax(1.0e6, 1.0e-4 * infeasibilityCost_) / largestObj;
3589                                   }
3590                              }
3591                         } else if (ratio < 1.0e12) {
3592                              // not so good
3593                              if (smallestObj < 1.0e-7) {
3594                                   // may as well scale up
3595                                   scalingDone = true;
3596                                   scale = 1.0e-6 / smallestObj;
3597                              } else if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3598                                   //done=true;
3599                              } else {
3600                                   scalingDone = true;
3601                                   if (algorithm_ < 0) {
3602                                        scale = 1.0e7 / largestObj;
3603                                   } else {
3604                                        scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3605                                   }
3606                              }
3607                         } else {
3608                              // Really nasty problem
3609                              if (smallestObj < 1.0e-8) {
3610                                   // may as well scale up
3611                                   scalingDone = true;
3612                                   scale = 1.0e-7 / smallestObj;
3613                                   largestObj *= scale;
3614                              }
3615                              if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3616                                   //done=true;
3617                              } else {
3618                                   scalingDone = true;
3619                                   if (algorithm_ < 0) {
3620                                        scale = 1.0e7 / largestObj;
3621                                   } else {
3622                                        scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3623                                   }
3624                              }
3625                         }
3626                         objectiveScale_ = scale;
3627                    }
3628                    if (largestRhs > 1.0e12) {
3629                         scalingDone = true;
3630                         rhsScale_ = 1.0e9 / largestRhs;
3631                    } else if (largestPositive > 1.0e-14 * smallestPositive && largestRhs > 1.0e6) {
3632                         scalingDone = true;
3633                         rhsScale_ = 1.0e6 / largestRhs;
3634                    } else {
3635                         rhsScale_ = 1.0;
3636                    }
3637                    if (scalingDone) {
3638                         handler_->message(CLP_RIM_SCALE, messages_)
3639                                   << objectiveScale_ << rhsScale_
3640                                   << CoinMessageEol;
3641                    }
3642               }
3643          } else if (makeRowCopy && scalingFlag_ > 0 && !rowCopyIsScaled) {
3644               matrix_->scaleRowCopy(this);
3645          }
3646          if (rowScale_ && !savedRowScale_) {
3647               inverseRowScale_ = rowScale_ + numberRows2;
3648               inverseColumnScale_ = columnScale_ + numberColumns_;
3649          }
3650          // See if we can try for faster row copy
3651          if (makeRowCopy && !oldMatrix) {
3652               ClpPackedMatrix* clpMatrix =
3653                    dynamic_cast< ClpPackedMatrix*>(matrix_);
3654               if (clpMatrix && numberThreads_)
3655                    clpMatrix->specialRowCopy(this, rowCopy_);
3656               if (clpMatrix)
3657                    clpMatrix->specialColumnCopy(this);
3658          }
3659     }
3660     if (what == 63) {
3661#if 0
3662          {
3663               x_gaps[0]++;
3664               ClpPackedMatrix* clpMatrix =
3665               dynamic_cast< ClpPackedMatrix*>(matrix_);
3666               if (clpMatrix) {
3667                    if (!clpMatrix->getPackedMatrix()->hasGaps())
3668                         x_gaps[1]++;
3669                    if ((clpMatrix->flags() & 2) == 0)
3670                         x_gaps[3]++;
3671               } else {
3672                    x_gaps[2]++;
3673               }
3674               if ((x_gaps[0] % 1000) == 0)
3675                    printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n",
3676                           x_gaps[0], x_gaps[1], x_gaps[2], x_gaps[3]);
3677          }
3678#endif
3679          if (newArrays && (specialOptions_ & 65536) == 0) {
3680               delete [] cost_;
3681               cost_ = new double[2*numberTotal];
3682               delete [] lower_;
3683               delete [] upper_;
3684               lower_ = new double[numberTotal];
3685               upper_ = new double[numberTotal];
3686               delete [] dj_;
3687               dj_ = new double[numberTotal];
3688               delete [] solution_;
3689               solution_ = new double[numberTotal];
3690          }
3691          reducedCostWork_ = dj_;
3692          rowReducedCost_ = dj_ + numberColumns_;
3693          columnActivityWork_ = solution_;
3694          rowActivityWork_ = solution_ + numberColumns_;
3695          objectiveWork_ = cost_;
3696          rowObjectiveWork_ = cost_ + numberColumns_;
3697          rowLowerWork_ = lower_ + numberColumns_;
3698          columnLowerWork_ = lower_;
3699          rowUpperWork_ = upper_ + numberColumns_;
3700          columnUpperWork_ = upper_;
3701     }
3702     if ((what & 4) != 0) {
3703          double direction = optimizationDirection_ * objectiveScale_;
3704          const double * obj = objective();
3705          const double * rowScale = rowScale_;
3706          const double * columnScale = columnScale_;
3707          // and also scale by scale factors
3708          if (rowScale) {
3709               if (rowObjective_) {
3710                    for (i = 0; i < numberRows_; i++)
3711                         rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
3712               } else {
3713                    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3714               }
3715               // If scaled then do all columns later in one loop
3716               if (what != 63) {
3717                    for (i = 0; i < numberColumns_; i++) {
3718                         CoinAssert(fabs(obj[i]) < 1.0e25);
3719                         objectiveWork_[i] = obj[i] * direction * columnScale[i];
3720                    }
3721               }
3722          } else {
3723               if (rowObjective_) {
3724                    for (i = 0; i < numberRows_; i++)
3725                         rowObjectiveWork_[i] = rowObjective_[i] * direction;
3726               } else {
3727                    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3728               }
3729               for (i = 0; i < numberColumns_; i++) {
3730                    CoinAssert(fabs(obj[i]) < 1.0e25);
3731                    objectiveWork_[i] = obj[i] * direction;
3732               }
3733          }
3734     }
3735     if ((what & 1) != 0) {
3736          const double * rowScale = rowScale_;
3737          // clean up any mismatches on infinity
3738          // and fix any variables with tiny gaps
3739          double primalTolerance = dblParam_[ClpPrimalTolerance];
3740          if(rowScale) {
3741               // If scaled then do all columns later in one loop
3742               if (what != 63) {
3743                    const double * inverseScale = inverseColumnScale_;
3744                    for (i = 0; i < numberColumns_; i++) {
3745                         double multiplier = rhsScale_ * inverseScale[i];
3746                         double lowerValue = columnLower_[i];
3747                         double upperValue = columnUpper_[i];
3748                         if (lowerValue > -1.0e20) {
3749                              columnLowerWork_[i] = lowerValue * multiplier;
3750                              if (upperValue >= 1.0e20) {
3751                                   columnUpperWork_[i] = COIN_DBL_MAX;
3752                              } else {
3753                                   columnUpperWork_[i] = upperValue * multiplier;
3754                                   if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3755                                        if (columnLowerWork_[i] >= 0.0) {
3756                                             columnUpperWork_[i] = columnLowerWork_[i];
3757                                        } else if (columnUpperWork_[i] <= 0.0) {
3758                                             columnLowerWork_[i] = columnUpperWork_[i];
3759                                        } else {
3760                                             columnUpperWork_[i] = 0.0;
3761                                             columnLowerWork_[i] = 0.0;
3762                                        }
3763                                   }
3764                              }
3765                         } else if (upperValue < 1.0e20) {
3766                              columnLowerWork_[i] = -COIN_DBL_MAX;
3767                              columnUpperWork_[i] = upperValue * multiplier;
3768                         } else {
3769                              // free
3770                              columnLowerWork_[i] = -COIN_DBL_MAX;
3771                              columnUpperWork_[i] = COIN_DBL_MAX;
3772                         }
3773                    }
3774               }
3775               for (i = 0; i < numberRows_; i++) {
3776                    double multiplier = rhsScale_ * rowScale[i];
3777                    double lowerValue = rowLower_[i];
3778                    double upperValue = rowUpper_[i];
3779                    if (lowerValue > -1.0e20) {
3780                         rowLowerWork_[i] = lowerValue * multiplier;
3781                         if (upperValue >= 1.0e20) {
3782                              rowUpperWork_[i] = COIN_DBL_MAX;
3783                         } else {
3784                              rowUpperWork_[i] = upperValue * multiplier;
3785                              if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3786                                   if (rowLowerWork_[i] >= 0.0) {
3787                                        rowUpperWork_[i] = rowLowerWork_[i];
3788                                   } else if (rowUpperWork_[i] <= 0.0) {
3789                                        rowLowerWork_[i] = rowUpperWork_[i];
3790                                   } else {
3791                                        rowUpperWork_[i] = 0.0;
3792                                        rowLowerWork_[i] = 0.0;
3793                                   }
3794                              }
3795                         }
3796                    } else if (upperValue < 1.0e20) {
3797                         rowLowerWork_[i] = -COIN_DBL_MAX;
3798                         rowUpperWork_[i] = upperValue * multiplier;
3799                    } else {
3800                         // free
3801                         rowLowerWork_[i] = -COIN_DBL_MAX;
3802                         rowUpperWork_[i] = COIN_DBL_MAX;
3803                    }
3804               }
3805          } else if (rhsScale_ != 1.0) {
3806               for (i = 0; i < numberColumns_; i++) {
3807                    double lowerValue = columnLower_[i];
3808                    double upperValue = columnUpper_[i];
3809                    if (lowerValue > -1.0e20) {
3810                         columnLowerWork_[i] = lowerValue * rhsScale_;
3811                         if (upperValue >= 1.0e20) {
3812                              columnUpperWork_[i] = COIN_DBL_MAX;
3813                         } else {
3814                              columnUpperWork_[i] = upperValue * rhsScale_;
3815                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3816                                   if (columnLowerWork_[i] >= 0.0) {
3817                                        columnUpperWork_[i] = columnLowerWork_[i];
3818                                   } else if (columnUpperWork_[i] <= 0.0) {
3819                                        columnLowerWork_[i] = columnUpperWork_[i];
3820                                   } else {
3821                                        columnUpperWork_[i] = 0.0;
3822                                        columnLowerWork_[i] = 0.0;
3823                                   }
3824                              }
3825                         }
3826                    } else if (upperValue < 1.0e20) {
3827                         columnLowerWork_[i] = -COIN_DBL_MAX;
3828                         columnUpperWork_[i] = upperValue * rhsScale_;
3829                    } else {
3830                         // free
3831                         columnLowerWork_[i] = -COIN_DBL_MAX;
3832                         columnUpperWork_[i] = COIN_DBL_MAX;
3833                    }
3834               }
3835               for (i = 0; i < numberRows_; i++) {
3836                    double lowerValue = rowLower_[i];
3837                    double upperValue = rowUpper_[i];
3838                    if (lowerValue > -1.0e20) {
3839                         rowLowerWork_[i] = lowerValue * rhsScale_;
3840                         if (upperValue >= 1.0e20) {
3841                              rowUpperWork_[i] = COIN_DBL_MAX;
3842                         } else {
3843                              rowUpperWork_[i] = upperValue * rhsScale_;
3844                              if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3845                                   if (rowLowerWork_[i] >= 0.0) {
3846                                        rowUpperWork_[i] = rowLowerWork_[i];
3847                                   } else if (rowUpperWork_[i] <= 0.0) {
3848                                        rowLowerWork_[i] = rowUpperWork_[i];
3849                                   } else {
3850                                        rowUpperWork_[i] = 0.0;
3851                                        rowLowerWork_[i] = 0.0;
3852                                   }
3853                              }
3854                         }
3855                    } else if (upperValue < 1.0e20) {
3856                         rowLowerWork_[i] = -COIN_DBL_MAX;
3857                         rowUpperWork_[i] = upperValue * rhsScale_;
3858                    } else {
3859                         // free
3860                         rowLowerWork_[i] = -COIN_DBL_MAX;
3861                         rowUpperWork_[i] = COIN_DBL_MAX;
3862                    }
3863               }
3864          } else {
3865               for (i = 0; i < numberColumns_; i++) {
3866                    double lowerValue = columnLower_[i];
3867                    double upperValue = columnUpper_[i];
3868                    if (lowerValue > -1.0e20) {
3869                         columnLowerWork_[i] = lowerValue;
3870                         if (upperValue >= 1.0e20) {
3871                              columnUpperWork_[i] = COIN_DBL_MAX;
3872                         } else {
3873                              columnUpperWork_[i] = upperValue;
3874                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3875                                   if (columnLowerWork_[i] >= 0.0) {
3876                                        columnUpperWork_[i] = columnLowerWork_[i];
3877                                   } else if (columnUpperWork_[i] <= 0.0) {
3878                                        columnLowerWork_[i] = columnUpperWork_[i];
3879                                   } else {
3880                                        columnUpperWork_[i] = 0.0;
3881                                        columnLowerWork_[i] = 0.0;
3882                                   }
3883                              }
3884                         }
3885                    } else if (upperValue < 1.0e20) {
3886                         columnLowerWork_[i] = -COIN_DBL_MAX;
3887                         columnUpperWork_[i] = upperValue;
3888                    } else {
3889                         // free
3890                         columnLowerWork_[i] = -COIN_DBL_MAX;
3891                         columnUpperWork_[i] = COIN_DBL_MAX;
3892                    }
3893               }
3894               for (i = 0; i < numberRows_; i++) {
3895                    double lowerValue = rowLower_[i];
3896                    double upperValue = rowUpper_[i];
3897                    if (lowerValue > -1.0e20) {
3898                         rowLowerWork_[i] = lowerValue;
3899                         if (upperValue >= 1.0e20) {
3900                              rowUpperWork_[i] = COIN_DBL_MAX;
3901                         } else {
3902                              rowUpperWork_[i] = upperValue;
3903                              if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3904                                   if (rowLowerWork_[i] >= 0.0) {
3905                                        rowUpperWork_[i] = rowLowerWork_[i];
3906                                   } else if (rowUpperWork_[i] <= 0.0) {
3907                                        rowLowerWork_[i] = rowUpperWork_[i];
3908                                   } else {
3909                                        rowUpperWork_[i] = 0.0;
3910                                        rowLowerWork_[i] = 0.0;
3911                                   }
3912                              }
3913                         }
3914                    } else if (upperValue < 1.0e20) {
3915                         rowLowerWork_[i] = -COIN_DBL_MAX;
3916                         rowUpperWork_[i] = upperValue;
3917                    } else {
3918                         // free
3919                         rowLowerWork_[i] = -COIN_DBL_MAX;
3920                         rowUpperWork_[i] = COIN_DBL_MAX;
3921                    }
3922               }
3923          }
3924     }
3925     if (what == 63) {
3926          // move information to work arrays
3927          double direction = optimizationDirection_;
3928          // direction is actually scale out not scale in
3929          if (direction)
3930               direction = 1.0 / direction;
3931          if (direction != 1.0) {
3932               // reverse all dual signs
3933               for (i = 0; i < numberColumns_; i++)
3934                    reducedCost_[i] *= direction;
3935               for (i = 0; i < numberRows_; i++)
3936                    dual_[i] *= direction;
3937          }
3938          for (i = 0; i < numberRows_ + numberColumns_; i++) {
3939               setFakeBound(i, noFake);
3940          }
3941          if (rowScale_) {
3942               const double * obj = objective();
3943               double direction = optimizationDirection_ * objectiveScale_;
3944               // clean up any mismatches on infinity
3945               // and fix any variables with tiny gaps
3946               double primalTolerance = dblParam_[ClpPrimalTolerance];
3947               // on entry
3948               const double * inverseScale = inverseColumnScale_;
3949               for (i = 0; i < numberColumns_; i++) {
3950                    CoinAssert(fabs(obj[i]) < 1.0e25);
3951                    double scaleFactor = columnScale_[i];
3952                    double multiplier = rhsScale_ * inverseScale[i];
3953                    scaleFactor *= direction;
3954                    objectiveWork_[i] = obj[i] * scaleFactor;
3955                    reducedCostWork_[i] = reducedCost_[i] * scaleFactor;
3956                    double lowerValue = columnLower_[i];
3957                    double upperValue = columnUpper_[i];
3958                    if (lowerValue > -1.0e20) {
3959                         columnLowerWork_[i] = lowerValue * multiplier;
3960                         if (upperValue >= 1.0e20) {
3961                              columnUpperWork_[i] = COIN_DBL_MAX;
3962                         } else {
3963                              columnUpperWork_[i] = upperValue * multiplier;
3964                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3965                                   if (columnLowerWork_[i] >= 0.0) {
3966                                        columnUpperWork_[i] = columnLowerWork_[i];
3967                                   } else if (columnUpperWork_[i] <= 0.0) {
3968                                        columnLowerWork_[i] = columnUpperWork_[i];
3969                                   } else {
3970                                        columnUpperWork_[i] = 0.0;
3971                                        columnLowerWork_[i] = 0.0;
3972                                   }
3973                              }
3974                         }
3975                    } else if (upperValue < 1.0e20) {
3976                         columnLowerWork_[i] = -COIN_DBL_MAX;
3977                         columnUpperWork_[i] = upperValue * multiplier;
3978                    } else {
3979                         // free
3980                         columnLowerWork_[i] = -COIN_DBL_MAX;
3981                         columnUpperWork_[i] = COIN_DBL_MAX;
3982                    }
3983                    double value = columnActivity_[i] * multiplier;
3984                    if (fabs(value) > 1.0e20) {
3985                         //printf("bad value of %g for column %d\n",value,i);
3986                         setColumnStatus(i, superBasic);
3987                         if (columnUpperWork_[i] < 0.0) {
3988                              value = columnUpperWork_[i];
3989                         } else if (columnLowerWork_[i] > 0.0) {
3990                              value = columnLowerWork_[i];
3991                         } else {
3992                              value = 0.0;
3993                         }
3994                    }
3995                    columnActivityWork_[i] = value;
3996               }
3997               inverseScale = inverseRowScale_;
3998               for (i = 0; i < numberRows_; i++) {
3999                    dual_[i] *= inverseScale[i];
4000                    dual_[i] *= objectiveScale_;
4001                    rowReducedCost_[i] = dual_[i];
4002                    double multiplier = rhsScale_ * rowScale_[i];
4003                    double value = rowActivity_[i] * multiplier;
4004                    if (fabs(value) > 1.0e20) {
4005                         //printf("bad value of %g for row %d\n",value,i);
4006                         setRowStatus(i, superBasic);
4007                         if (rowUpperWork_[i] < 0.0) {
4008                              value = rowUpperWork_[i];
4009                         } else if (rowLowerWork_[i] > 0.0) {
4010                              value = rowLowerWork_[i];
4011                         } else {
4012                              value = 0.0;
4013                         }
4014                    }
4015                    rowActivityWork_[i] = value;
4016               }
4017          } else if (objectiveScale_ != 1.0 || rhsScale_ != 1.0) {
4018               // on entry
4019               for (i = 0; i < numberColumns_; i++) {
4020                    double value = columnActivity_[i];
4021                    value *= rhsScale_;
4022                    if (fabs(value) > 1.0e20) {
4023                         //printf("bad value of %g for column %d\n",value,i);
4024                         setColumnStatus(i, superBasic);
4025                         if (columnUpperWork_[i] < 0.0) {
4026                              value = columnUpperWork_[i];
4027                         } else if (columnLowerWork_[i] > 0.0) {
4028                              value = columnLowerWork_[i];
4029                         } else {
4030                              value = 0.0;
4031                         }
4032                    }
4033                    columnActivityWork_[i] = value;
4034                    reducedCostWork_[i] = reducedCost_[i] * objectiveScale_;
4035               }
4036               for (i = 0; i < numberRows_; i++) {
4037                    double value = rowActivity_[i];
4038                    value *= rhsScale_;
4039                    if (fabs(value) > 1.0e20) {
4040                         //printf("bad value of %g for row %d\n",value,i);
4041                         setRowStatus(i, superBasic);
4042                         if (rowUpperWork_[i] < 0.0) {
4043                              value = rowUpperWork_[i];
4044                         } else if (rowLowerWork_[i] > 0.0) {
4045                              value = rowLowerWork_[i];
4046                         } else {
4047                              value = 0.0;
4048                         }
4049                    }
4050                    rowActivityWork_[i] = value;
4051                    dual_[i] *= objectiveScale_;
4052                    rowReducedCost_[i] = dual_[i];
4053               }
4054          } else {
4055               // on entry
4056               for (i = 0; i < numberColumns_; i++) {
4057                    double value = columnActivity_[i];
4058                    if (fabs(value) > 1.0e20) {
4059                         //printf("bad value of %g for column %d\n",value,i);
4060                         setColumnStatus(i, superBasic);
4061                         if (columnUpperWork_[i] < 0.0) {
4062                              value = columnUpperWork_[i];
4063                         } else if (columnLowerWork_[i] > 0.0) {
4064                              value = columnLowerWork_[i];
4065                         } else {
4066                              value = 0.0;
4067                         }
4068                    }
4069                    columnActivityWork_[i] = value;
4070                    reducedCostWork_[i] = reducedCost_[i];
4071               }
4072               for (i = 0; i < numberRows_; i++) {
4073                    double value = rowActivity_[i];
4074                    if (fabs(value) > 1.0e20) {
4075                         //printf("bad value of %g for row %d\n",value,i);
4076                         setRowStatus(i, superBasic);
4077                         if (rowUpperWork_[i] < 0.0) {
4078                              value = rowUpperWork_[i];
4079                         } else if (rowLowerWork_[i] > 0.0) {
4080                              value = rowLowerWork_[i];
4081                         } else {
4082                              value = 0.0;
4083                         }
4084                    }
4085                    rowActivityWork_[i] = value;
4086                    rowReducedCost_[i] = dual_[i];
4087               }
4088          }
4089     }
4090
4091     if (what == 63 && doSanityCheck) {
4092          // check rim of problem okay
4093          if (!sanityCheck())
4094               goodMatrix = false;
4095     }
4096     // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
4097     // maybe we need to move scales to SimplexModel for factorization?
4098     if ((what == 63 && !pivotVariable_) || (newArrays && !keepPivots)) {
4099          delete [] pivotVariable_;
4100          pivotVariable_ = new int[numberRows2];
4101          for (int i = 0; i < numberRows2; i++)
4102               pivotVariable_[i] = -1;
4103     } else if (what == 63 && !keepPivots) {
4104          // just reset
4105          for (int i = 0; i < numberRows2; i++)
4106               pivotVariable_[i] = -1;
4107     } else if (what == 63) {
4108          // check pivots
4109          for (int i = 0; i < numberRows2; i++) {
4110               int iSequence = pivotVariable_[i];
4111               if (iSequence<0||(iSequence < numberRows_ + numberColumns_ &&
4112                                 getStatus(iSequence) != basic)) {
4113                    keepPivots = false;
4114                    break;
4115               }
4116          }
4117          if (!keepPivots) {
4118               // reset
4119               for (int i = 0; i < numberRows2; i++)
4120                    pivotVariable_[i] = -1;
4121          } else {
4122               // clean
4123               for (int i = 0; i < numberColumns_ + numberRows_; i++) {
4124                    Status status = getStatus(i);
4125                    if (status != basic) {
4126                         if (upper_[i] == lower_[i]) {
4127                              setStatus(i, isFixed);
4128                              solution_[i] = lower_[i];
4129                         } else if (status == atLowerBound) {
4130                              if (lower_[i] > -1.0e20) {
4131                                   solution_[i] = lower_[i];
4132                              } else {
4133                                   //printf("seq %d at lower of %g\n",i,lower_[i]);
4134                                   if (upper_[i] < 1.0e20) {
4135                                        solution_[i] = upper_[i];
4136                                        setStatus(i, atUpperBound);
4137                                   } else {
4138                                        setStatus(i, isFree);
4139                                   }
4140                              }
4141                         } else if (status == atUpperBound) {
4142                              if (upper_[i] < 1.0e20) {
4143                                   solution_[i] = upper_[i];
4144                              } else {
4145                                   //printf("seq %d at upper of %g\n",i,upper_[i]);
4146                                   if (lower_[i] > -1.0e20) {
4147                                        solution_[i] = lower_[i];
4148                                        setStatus(i, atLowerBound);
4149                                   } else {
4150                                        setStatus(i, isFree);
4151                                   }
4152                              }
4153                         } else if (status == isFixed && upper_[i] > lower_[i]) {
4154                              // was fixed - not now
4155                              if (solution_[i] <= lower_[i]+primalTolerance_) {
4156                                   setStatus(i, atLowerBound);
4157                              } else if (solution_[i] >= upper_[i]-primalTolerance_) {
4158                                   setStatus(i, atUpperBound);
4159                              } else {
4160                                   setStatus(i, superBasic);
4161                              }
4162                         }
4163                    }
4164               }
4165          }
4166     }
4167
4168     if (what == 63) {
4169          if (newArrays) {
4170               // get some arrays
4171               int iRow, iColumn;
4172               // these are "indexed" arrays so we always know where nonzeros are
4173               /**********************************************************
4174               rowArray_[3] is long enough for rows+columns (2 also maybe)
4175               rowArray_[1] is long enough for max(rows,columns)
4176               *********************************************************/
4177               for (iRow = 0; iRow < 4; iRow++) {
4178                    int length = numberRows2 + factorization_->maximumPivots();
4179                    if (iRow > SHORT_REGION || objective_->type() > 1)
4180                         length += numberColumns_;
4181                    else if (iRow == 1)
4182                         length = CoinMax(length, numberColumns_);
4183                    if ((specialOptions_ & 65536) == 0 || !rowArray_[iRow]) {
4184                         delete rowArray_[iRow];
4185                         rowArray_[iRow] = new CoinIndexedVector();
4186                    }
4187                    rowArray_[iRow]->reserve(length);
4188               }
4189
4190               for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
4191                    if ((specialOptions_ & 65536) == 0 || !columnArray_[iColumn]) {
4192                         delete columnArray_[iColumn];
4193                         columnArray_[iColumn] = new CoinIndexedVector();
4194                    }
4195                    columnArray_[iColumn]->reserve(numberColumns_+numberRows2);
4196               }
4197          } else {
4198               int iRow, iColumn;
4199               for (iRow = 0; iRow < 4; iRow++) {
4200                    int length = numberRows2 + factorization_->maximumPivots();
4201                    if (iRow > SHORT_REGION || objective_->type() > 1)
4202                         length += numberColumns_;
4203                    if(rowArray_[iRow]->capacity() >= length) {
4204                         rowArray_[iRow]->clear();
4205                    } else {
4206                         // model size or maxinv changed
4207                         rowArray_[iRow]->reserve(length);
4208                    }
4209#ifndef NDEBUG
4210                    rowArray_[iRow]->checkClear();
4211#endif
4212               }
4213
4214               for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) {
4215                    int length = numberColumns_;
4216                    if (iColumn)
4217                         length = CoinMax(numberRows2, numberColumns_);
4218                    if(columnArray_[iColumn]->capacity() >= length) {
4219                         columnArray_[iColumn]->clear();
4220                    } else {
4221                         // model size or maxinv changed
4222                         columnArray_[iColumn]->reserve(length);
4223                    }
4224#ifndef NDEBUG
4225                    columnArray_[iColumn]->checkClear();
4226#endif
4227               }
4228          }
4229     }
4230     if (problemStatus_ == 10) {
4231          problemStatus_ = -1;
4232          handler_->setLogLevel(saveLevel); // switch back messages
4233     }
4234     if ((what & 5) != 0)
4235          matrix_->generalExpanded(this, 9, what); // update costs and bounds if necessary
4236     if (goodMatrix && (specialOptions_ & 65536) != 0) {
4237          int save = maximumColumns_ + maximumRows_;
4238          CoinMemcpyN(cost_, numberTotal, cost_ + save);
4239          CoinMemcpyN(lower_, numberTotal, lower_ + save);
4240          CoinMemcpyN(upper_, numberTotal, upper_ + save);
4241          CoinMemcpyN(dj_, numberTotal, dj_ + save);
4242          CoinMemcpyN(solution_, numberTotal, solution_ + save);
4243          if (rowScale_ && !savedRowScale_) {
4244               double * temp;
4245               temp = new double [4*maximumRows_];
4246               CoinFillN(temp, 4 * maximumRows_, 1.0);
4247               CoinMemcpyN(rowScale_, numberRows2, temp);
4248               CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + maximumRows_);
4249               CoinMemcpyN(rowScale_, numberRows2, temp + 2 * maximumRows_);
4250               CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + 3 * maximumRows_);
4251               delete [] rowScale_;
4252               savedRowScale_ = temp;
4253               rowScale_ = savedRowScale_;
4254               inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
4255               temp = new double [4*maximumColumns_];
4256               CoinFillN(temp, 4 * maximumColumns_, 1.0);
4257               CoinMemcpyN(columnScale_, numberColumns_, temp);
4258               CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + maximumColumns_);
4259               CoinMemcpyN(columnScale_, numberColumns_, temp + 2 * maximumColumns_);
4260               CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + 3 * maximumColumns_);
4261               delete [] columnScale_;
4262               savedColumnScale_ = temp;
4263               columnScale_ = savedColumnScale_;
4264               inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
4265          }
4266     }
4267#ifdef CLP_USER_DRIVEN
4268     eventHandler_->event(ClpEventHandler::endOfCreateRim);
4269#endif
4270     return goodMatrix;
4271}
4272// Does rows and columns
4273void
4274ClpSimplex::createRim1(bool initial)
4275{
4276     int i;
4277     int numberRows2 = numberRows_ + numberExtraRows_;
4278     int numberTotal = numberRows2 + numberColumns_;
4279     if ((specialOptions_ & 65536) != 0 && true) {
4280          assert (!initial);
4281          int save = maximumColumns_ + maximumRows_;
4282          CoinMemcpyN(lower_ + save, numberTotal, lower_);
4283          CoinMemcpyN(upper_ + save, numberTotal, upper_);
4284          return;
4285     }
4286     const double * rowScale = rowScale_;
4287     // clean up any mismatches on infinity
4288     // and fix any variables with tiny gaps
4289     double primalTolerance = dblParam_[ClpPrimalTolerance];
4290     if(rowScale) {
4291          // If scaled then do all columns later in one loop
4292          if (!initial) {
4293               const double * inverseScale = inverseColumnScale_;
4294               for (i = 0; i < numberColumns_; i++) {
4295                    double multiplier = rhsScale_ * inverseScale[i];
4296                    double lowerValue = columnLower_[i];
4297                    double upperValue = columnUpper_[i];
4298                    if (lowerValue > -1.0e20) {
4299                         columnLowerWork_[i] = lowerValue * multiplier;
4300                         if (upperValue >= 1.0e20) {
4301                              columnUpperWork_[i] = COIN_DBL_MAX;
4302                         } else {
4303                              columnUpperWork_[i] = upperValue * multiplier;
4304                              if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4305                                   if (columnLowerWork_[i] >= 0.0) {
4306                                        columnUpperWork_[i] = columnLowerWork_[i];
4307                                   } else if (columnUpperWork_[i] <= 0.0) {
4308                                        columnLowerWork_[i] = columnUpperWork_[i];
4309                                   } else {
4310                                        columnUpperWork_[i] = 0.0;
4311                                        columnLowerWork_[i] = 0.0;
4312                                   }
4313                              }
4314                         }
4315                    } else if (upperValue < 1.0e20) {
4316                         columnLowerWork_[i] = -COIN_DBL_MAX;
4317                         columnUpperWork_[i] = upperValue * multiplier;
4318                    } else {
4319                         // free
4320                         columnLowerWork_[i] = -COIN_DBL_MAX;
4321                         columnUpperWork_[i] = COIN_DBL_MAX;
4322                    }
4323               }
4324          }
4325          for (i = 0; i < numberRows_; i++) {
4326               double multiplier = rhsScale_ * rowScale[i];
4327               double lowerValue = rowLower_[i];
4328               double upperValue = rowUpper_[i];
4329               if (lowerValue > -1.0e20) {
4330                    rowLowerWork_[i] = lowerValue * multiplier;
4331                    if (upperValue >= 1.0e20) {
4332                         rowUpperWork_[i] = COIN_DBL_MAX;
4333                    } else {
4334                         rowUpperWork_[i] = upperValue * multiplier;
4335                         if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4336                              if (rowLowerWork_[i] >= 0.0) {
4337                                   rowUpperWork_[i] = rowLowerWork_[i];
4338                              } else if (rowUpperWork_[i] <= 0.0) {
4339                                   rowLowerWork_[i] = rowUpperWork_[i];
4340                              } else {
4341                                   rowUpperWork_[i] = 0.0;
4342                                   rowLowerWork_[i] = 0.0;
4343                              }
4344                         }
4345                    }
4346               } else if (upperValue < 1.0e20) {
4347                    rowLowerWork_[i] = -COIN_DBL_MAX;
4348                    rowUpperWork_[i] = upperValue * multiplier;
4349               } else {
4350                    // free
4351                    rowLowerWork_[i] = -COIN_DBL_MAX;
4352                    rowUpperWork_[i] = COIN_DBL_MAX;
4353               }
4354          }
4355     } else if (rhsScale_ != 1.0) {
4356          for (i = 0; i < numberColumns_; i++) {
4357               double lowerValue = columnLower_[i];
4358               double upperValue = columnUpper_[i];
4359               if (lowerValue > -1.0e20) {
4360                    columnLowerWork_[i] = lowerValue * rhsScale_;
4361                    if (upperValue >= 1.0e20) {
4362                         columnUpperWork_[i] = COIN_DBL_MAX;
4363                    } else {
4364                         columnUpperWork_[i] = upperValue * rhsScale_;
4365                         if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4366                              if (columnLowerWork_[i] >= 0.0) {
4367                                   columnUpperWork_[i] = columnLowerWork_[i];
4368                              } else if (columnUpperWork_[i] <= 0.0) {
4369                                   columnLowerWork_[i] = columnUpperWork_[i];
4370                              } else {
4371                                   columnUpperWork_[i] = 0.0;
4372                                   columnLowerWork_[i] = 0.0;
4373                              }
4374                         }
4375                    }
4376               } else if (upperValue < 1.0e20) {
4377                    columnLowerWork_[i] = -COIN_DBL_MAX;
4378                    columnUpperWork_[i] = upperValue * rhsScale_;
4379               } else {
4380                    // free
4381                    columnLowerWork_[i] = -COIN_DBL_MAX;
4382                    columnUpperWork_[i] = COIN_DBL_MAX;
4383               }
4384          }
4385          for (i = 0; i < numberRows_; i++) {
4386               double lowerValue = rowLower_[i];
4387               double upperValue = rowUpper_[i];
4388               if (lowerValue > -1.0e20) {
4389                    rowLowerWork_[i] = lowerValue * rhsScale_;
4390                    if (upperValue >= 1.0e20) {
4391                         rowUpperWork_[i] = COIN_DBL_MAX;
4392                    } else {
4393                         rowUpperWork_[i] = upperValue * rhsScale_;
4394                         if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4395                              if (rowLowerWork_[i] >= 0.0) {
4396                                   rowUpperWork_[i] = rowLowerWork_[i];
4397                              } else if (rowUpperWork_[i] <= 0.0) {
4398                                   rowLowerWork_[i] = rowUpperWork_[i];
4399                              } else {
4400                                   rowUpperWork_[i] = 0.0;
4401                                   rowLowerWork_[i] = 0.0;
4402                              }
4403                         }
4404                    }
4405               } else if (upperValue < 1.0e20) {
4406                    rowLowerWork_[i] = -COIN_DBL_MAX;
4407                    rowUpperWork_[i] = upperValue * rhsScale_;
4408               } else {
4409                    // free
4410                    rowLowerWork_[i] = -COIN_DBL_MAX;
4411                    rowUpperWork_[i] = COIN_DBL_MAX;
4412               }
4413          }
4414     } else {
4415          for (i = 0; i < numberColumns_; i++) {
4416               double lowerValue = columnLower_[i];
4417               double upperValue = columnUpper_[i];
4418               if (lowerValue > -1.0e20) {
4419                    columnLowerWork_[i] = lowerValue;
4420                    if (upperValue >= 1.0e20) {
4421                         columnUpperWork_[i] = COIN_DBL_MAX;
4422                    } else {
4423                         columnUpperWork_[i] = upperValue;
4424                         if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4425                              if (columnLowerWork_[i] >= 0.0) {
4426                                   columnUpperWork_[i] = columnLowerWork_[i];
4427                              } else if (columnUpperWork_[i] <= 0.0) {
4428                                   columnLowerWork_[i] = columnUpperWork_[i];
4429                              } else {
4430                                   columnUpperWork_[i] = 0.0;
4431                                   columnLowerWork_[i] = 0.0;
4432                              }
4433                         }
4434                    }
4435               } else if (upperValue < 1.0e20) {
4436                    columnLowerWork_[i] = -COIN_DBL_MAX;
4437                    columnUpperWork_[i] = upperValue;
4438               } else {
4439                    // free
4440                    columnLowerWork_[i] = -COIN_DBL_MAX;
4441                    columnUpperWork_[i] = COIN_DBL_MAX;
4442               }
4443          }
4444          for (i = 0; i < numberRows_; i++) {
4445               double lowerValue = rowLower_[i];
4446               double upperValue = rowUpper_[i];
4447               if (lowerValue > -1.0e20) {
4448                    rowLowerWork_[i] = lowerValue;
4449                    if (upperValue >= 1.0e20) {
4450                         rowUpperWork_[i] = COIN_DBL_MAX;
4451                    } else {
4452                         rowUpperWork_[i] = upperValue;
4453                         if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4454                              if (rowLowerWork_[i] >= 0.0) {
4455                                   rowUpperWork_[i] = rowLowerWork_[i];
4456                              } else if (rowUpperWork_[i] <= 0.0) {
4457                                   rowLowerWork_[i] = rowUpperWork_[i];
4458                              } else {
4459                                   rowUpperWork_[i] = 0.0;
4460                                   rowLowerWork_[i] = 0.0;
4461                              }
4462                         }
4463                    }
4464               } else if (upperValue < 1.0e20) {
4465                    rowLowerWork_[i] = -COIN_DBL_MAX;
4466                    rowUpperWork_[i] = upperValue;
4467               } else {
4468                    // free
4469                    rowLowerWork_[i] = -COIN_DBL_MAX;
4470                    rowUpperWork_[i] = COIN_DBL_MAX;
4471               }
4472          }
4473     }
4474#ifndef NDEBUG
4475     if ((specialOptions_ & 65536) != 0 && false) {
4476          assert (!initial);
4477          int save = maximumColumns_ + maximumRows_;
4478          for (int i = 0; i < numberTotal; i++) {
4479               assert (fabs(lower_[i] - lower_[i+save]) < 1.0e-5);
4480               assert (fabs(upper_[i] - upper_[i+save]) < 1.0e-5);
4481          }
4482     }
4483#endif
4484}
4485// Does objective
4486void
4487ClpSimplex::createRim4(bool initial)
4488{
4489     int i;
4490     int numberRows2 = numberRows_ + numberExtraRows_;
4491     int numberTotal = numberRows2 + numberColumns_;
4492     if ((specialOptions_ & 65536) != 0 && true) {
4493          assert (!initial);
4494          int save = maximumColumns_ + maximumRows_;
4495          CoinMemcpyN(cost_ + save, numberTotal, cost_);
4496          return;
4497     }
4498     double direction = optimizationDirection_ * objectiveScale_;
4499     const double * obj = objective();
4500     const double * rowScale = rowScale_;
4501     const double * columnScale = columnScale_;
4502     // and also scale by scale factors
4503     if (rowScale) {
4504          if (rowObjective_) {
4505               for (i = 0; i < numberRows_; i++)
4506                    rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
4507          } else {
4508               memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4509          }
4510          // If scaled then do all columns later in one loop
4511          if (!initial) {
4512               for (i = 0; i < numberColumns_; i++) {
4513                    CoinAssert(fabs(obj[i]) < 1.0e25);
4514                    objectiveWork_[i] = obj[i] * direction * columnScale[i];
4515               }
4516          }
4517     } else {
4518          if (rowObjective_) {
4519               for (i = 0; i < numberRows_; i++)
4520                    rowObjectiveWork_[i] = rowObjective_[i] * direction;
4521          } else {
4522               memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4523          }
4524          for (i = 0; i < numberColumns_; i++) {
4525               CoinAssert(fabs(obj[i]) < 1.0e25);
4526               objectiveWork_[i] = obj[i] * direction;
4527          }
4528     }
4529}
4530// Does rows and columns and objective
4531void
4532ClpSimplex::createRim5(bool initial)
4533{
4534     createRim4(initial);
4535     createRim1(initial);
4536}
4537void
4538ClpSimplex::deleteRim(int getRidOfFactorizationData)
4539{
4540#ifdef CLP_USER_DRIVEN
4541     eventHandler_->event(ClpEventHandler::beforeDeleteRim);
4542#endif
4543     // Just possible empty problem
4544     int numberRows = numberRows_;
4545     int numberColumns = numberColumns_;
4546     if (!numberRows || !numberColumns) {
4547          numberRows = 0;
4548          if (objective_->type() < 2)
4549               numberColumns = 0;
4550     }
4551     int i;
4552     if (problemStatus_ != 1 && problemStatus_ != 2) {
4553          delete [] ray_;
4554          ray_ = NULL;
4555     }
4556     // set upperOut_ to furthest away from bound so can use in dual for dualBound_
4557     upperOut_ = 1.0;
4558#if 0
4559     {
4560          int nBad = 0;
4561          for (i = 0; i < numberColumns; i++) {
4562               if (lower_[i] == upper_[i] && getColumnStatus(i) == basic)
4563                    nBad++;
4564          }
4565          if (nBad)
4566               printf("yy %d basic fixed\n", nBad);
4567     }
4568#endif
4569     if ((moreSpecialOptions_&4194304)!=0) {
4570       // preset tolerances were changed
4571       moreSpecialOptions_ &= ~4194304;
4572       primalTolerance_=1.0e-7;
4573       dblParam_[ClpPrimalTolerance]=primalTolerance_;
4574       dualTolerance_=1.0e-7;
4575       dblParam_[ClpDualTolerance]=dualTolerance_;
4576     }
4577     // ray may be null if in branch and bound
4578     if (rowScale_ && solution_) {
4579          // Collect infeasibilities
4580          int numberPrimalScaled = 0;
4581          int numberPrimalUnscaled = 0;
4582          int numberDualScaled = 0;
4583          int numberDualUnscaled = 0;
4584          double scaleC = 1.0 / objectiveScale_;
4585          double scaleR = 1.0 / rhsScale_;
4586          const double * inverseScale = inverseColumnScale_;
4587          for (i = 0; i < numberColumns; i++) {
4588               double scaleFactor = columnScale_[i];
4589               double valueScaled = columnActivityWork_[i];
4590               double lowerScaled = columnLowerWork_[i];
4591               double upperScaled = columnUpperWork_[i];
4592               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4593                    if (valueScaled < lowerScaled - primalTolerance_ ||
4594                              valueScaled > upperScaled + primalTolerance_)
4595                         numberPrimalScaled++;
4596                    else
4597                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4598               }
4599               columnActivity_[i] = valueScaled * scaleFactor * scaleR;
4600               double value = columnActivity_[i];
4601               if (value < columnLower_[i] - primalTolerance_)
4602                    numberPrimalUnscaled++;
4603               else if (value > columnUpper_[i] + primalTolerance_)
4604                    numberPrimalUnscaled++;
4605               double valueScaledDual = reducedCostWork_[i];
4606               if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4607                    numberDualScaled++;
4608               if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4609                    numberDualScaled++;
4610               reducedCost_[i] = (valueScaledDual * scaleC) * inverseScale[i];
4611               double valueDual = reducedCost_[i];
4612               if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4613                    numberDualUnscaled++;
4614               if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4615                    numberDualUnscaled++;
4616          }
4617          inverseScale = inverseRowScale_;
4618          for (i = 0; i < numberRows; i++) {
4619               double scaleFactor = rowScale_[i];
4620               double valueScaled = rowActivityWork_[i];
4621               double lowerScaled = rowLowerWork_[i];
4622               double upperScaled = rowUpperWork_[i];
4623               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4624                    if (valueScaled < lowerScaled - primalTolerance_ ||
4625                              valueScaled > upperScaled + primalTolerance_)
4626                         numberPrimalScaled++;
4627                    else
4628                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4629               }
4630               rowActivity_[i] = (valueScaled * scaleR) * inverseScale[i];
4631               double value = rowActivity_[i];
4632               if (value < rowLower_[i] - primalTolerance_)
4633                    numberPrimalUnscaled++;
4634               else if (value > rowUpper_[i] + primalTolerance_)
4635                    numberPrimalUnscaled++;
4636               double valueScaledDual = dual_[i] + rowObjectiveWork_[i];;
4637               if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4638                    numberDualScaled++;
4639               if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4640                    numberDualScaled++;
4641               dual_[i] *= scaleFactor * scaleC;
4642               double valueDual = dual_[i];
4643               if (rowObjective_)
4644                    valueDual += rowObjective_[i];
4645               if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4646                    numberDualUnscaled++;
4647               if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4648                    numberDualUnscaled++;
4649          }
4650          if (!problemStatus_ && !secondaryStatus_) {
4651               // See if we need to set secondary status
4652               if (numberPrimalUnscaled) {
4653                    if (numberDualUnscaled)
4654                         secondaryStatus_ = 4;
4655                    else
4656                         secondaryStatus_ = 2;
4657               } else {
4658                    if (numberDualUnscaled)
4659                         secondaryStatus_ = 3;
4660               }
4661          }
4662          if (problemStatus_ == 2 && ray_) {
4663               for (i = 0; i < numberColumns; i++) {
4664                    ray_[i] *= columnScale_[i];
4665               }
4666          } else if (problemStatus_ == 1 && ray_) {
4667               for (i = 0; i < numberRows; i++) {
4668                    ray_[i] *= rowScale_[i];
4669               }
4670          }
4671     } else if (rhsScale_ != 1.0 || objectiveScale_ != 1.0) {
4672          // Collect infeasibilities
4673          int numberPrimalScaled = 0;
4674          int numberPrimalUnscaled = 0;
4675          int numberDualScaled = 0;
4676          int numberDualUnscaled = 0;
4677          double scaleC = 1.0 / objectiveScale_;
4678          double scaleR = 1.0 / rhsScale_;
4679          for (i = 0; i < numberColumns; i++) {
4680               double valueScaled = columnActivityWork_[i];
4681               double lowerScaled = columnLowerWork_[i];
4682               double upperScaled = columnUpperWork_[i];
4683               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4684                    if (valueScaled < lowerScaled - primalTolerance_ ||
4685                              valueScaled > upperScaled + primalTolerance_)
4686                         numberPrimalScaled++;
4687                    else
4688                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4689               }
4690               columnActivity_[i] = valueScaled * scaleR;
4691               double value = columnActivity_[i];
4692               if (value < columnLower_[i] - primalTolerance_)
4693                    numberPrimalUnscaled++;
4694               else if (value > columnUpper_[i] + primalTolerance_)
4695                    numberPrimalUnscaled++;
4696               double valueScaledDual = reducedCostWork_[i];
4697               if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4698                    numberDualScaled++;
4699               if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4700                    numberDualScaled++;
4701               reducedCost_[i] = valueScaledDual * scaleC;
4702               double valueDual = reducedCost_[i];
4703               if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4704                    numberDualUnscaled++;
4705               if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4706                    numberDualUnscaled++;
4707          }
4708          for (i = 0; i < numberRows; i++) {
4709               double valueScaled = rowActivityWork_[i];
4710               double lowerScaled = rowLowerWork_[i];
4711               double upperScaled = rowUpperWork_[i];
4712               if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
4713                    if (valueScaled < lowerScaled - primalTolerance_ ||
4714                              valueScaled > upperScaled + primalTolerance_)
4715                         numberPrimalScaled++;
4716                    else
4717                         upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4718               }
4719               rowActivity_[i] = valueScaled * scaleR;
4720               double value = rowActivity_[i];
4721               if (value < rowLower_[i] - primalTolerance_)
4722                    numberPrimalUnscaled++;
4723               else if (value > rowUpper_[i] + primalTolerance_)
4724                    numberPrimalUnscaled++;
4725               double valueScaledDual = dual_[i] + rowObjectiveWork_[i];;
4726               if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4727                    numberDualScaled++;
4728               if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4729                    numberDualScaled++;
4730               dual_[i] *= scaleC;
4731               double valueDual = dual_[i];
4732               if (rowObjective_)
4733                    valueDual += rowObjective_[i];
4734               if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4735                    numberDualUnscaled++;
4736               if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4737                    numberDualUnscaled++;
4738          }
4739          if (!problemStatus_ && !secondaryStatus_) {
4740               // See if we need to set secondary status
4741               if (numberPrimalUnscaled) {
4742                    if (numberDualUnscaled)
4743                         secondaryStatus_ = 4;
4744                    else
4745                         secondaryStatus_ = 2;
4746               } else {
4747                    if (numberDualUnscaled)
4748                         secondaryStatus_ = 3;
4749               }
4750          }
4751     } else {
4752          if (columnActivityWork_) {
4753               for (i = 0; i < numberColumns; i++) {
4754                    double value = columnActivityWork_[i];
4755                    double lower = columnLowerWork_[i];
4756                    double upper = columnUpperWork_[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                    columnActivity_[i] = columnActivityWork_[i];
4762                    reducedCost_[i] = reducedCostWork_[i];
4763               }
4764               for (i = 0; i < numberRows; i++) {
4765                    double value = rowActivityWork_[i];
4766                    double lower = rowLowerWork_[i];
4767                    double upper = rowUpperWork_[i];
4768                    if (lower > -1.0e20 || upper < 1.0e20) {
4769                         if (value > lower && value < upper)
4770                              upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4771                    }
4772                    rowActivity_[i] = rowActivityWork_[i];
4773               }
4774          }
4775     }
4776     // switch off scalefactor if auto
4777     if (automaticScale_) {
4778          rhsScale_ = 1.0;
4779          objectiveScale_ = 1.0;
4780     }
4781     if (optimizationDirection_ != 1.0) {
4782          // and modify all dual signs
4783          for (i = 0; i < numberColumns; i++)
4784               reducedCost_[i] *= optimizationDirection_;
4785          for (i = 0; i < numberRows; i++)
4786               dual_[i] *= optimizationDirection_;
4787     }
4788     // scaling may have been turned off
4789     scalingFlag_ = abs(scalingFlag_);
4790     if(getRidOfFactorizationData > 0) {
4791          gutsOfDelete(getRidOfFactorizationData + 1);
4792     } else {
4793          // at least get rid of nonLinearCost_
4794          delete nonLinearCost_;
4795          nonLinearCost_ = NULL;
4796     }
4797     if (!rowObjective_ && problemStatus_ == 0 && objective_->type() == 1 &&
4798               numberRows && numberColumns) {
4799          // Redo objective value
4800          double objectiveValue = 0.0;
4801          const double * cost = objective();
4802          for (int i = 0; i < numberColumns; i++) {
4803               double value = columnActivity_[i];
4804               objectiveValue += value * cost[i];
4805          }
4806          //if (fabs(objectiveValue_ -objectiveValue*optimizationDirection())>1.0e-5)
4807          //printf("old obj %g new %g\n",objectiveValue_, objectiveValue*optimizationDirection());
4808          objectiveValue_ = objectiveValue * optimizationDirection();
4809     }
4810     // get rid of data
4811     matrix_->generalExpanded(this, 13, scalingFlag_);
4812}
4813void
4814ClpSimplex::setDualBound(double value)
4815{
4816     if (value > 0.0)
4817          dualBound_ = value;
4818}
4819void
4820ClpSimplex::setInfeasibilityCost(double value)
4821{
4822     if (value > 0.0)
4823          infeasibilityCost_ = value;
4824}
4825void ClpSimplex::setNumberRefinements( int value)
4826{
4827     if (value >= 0 && value < 10)
4828          numberRefinements_ = value;
4829}
4830// Sets row pivot choice algorithm in dual
4831void
4832ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
4833{
4834     delete dualRowPivot_;
4835     dualRowPivot_ = choice.clone(true);
4836     dualRowPivot_->setModel(this);
4837}
4838// Sets column pivot choice algorithm in primal
4839void
4840ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
4841{
4842     delete primalColumnPivot_;
4843     primalColumnPivot_ = choice.clone(true);
4844     primalColumnPivot_->setModel(this);
4845}
4846void
4847ClpSimplex::setFactorization( ClpFactorization & factorization)
4848{
4849     if (factorization_)
4850          factorization_->setFactorization(factorization);
4851     else
4852          factorization_ = new ClpFactorization(factorization,
4853                                                numberRows_);
4854}
4855
4856// Swaps factorization
4857ClpFactorization *
4858ClpSimplex::swapFactorization( ClpFactorization * factorization)
4859{
4860     ClpFactorization * swap = factorization_;
4861     factorization_ = factorization;
4862     return swap;
4863}
4864// Copies in factorization to existing one
4865void
4866ClpSimplex::copyFactorization( ClpFactorization & factorization)
4867{
4868     *factorization_ = factorization;
4869}
4870/* Perturbation:
4871   -50 to +50 - perturb by this power of ten (-6 sounds good)
4872   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
4873   101 - we are perturbed
4874   102 - don't try perturbing again
4875   default is 100
4876*/
4877void
4878ClpSimplex::setPerturbation(int value)
4879{
4880     if(value <= 102 && value >= -1000) {
4881          perturbation_ = value;
4882     }
4883}
4884// Sparsity on or off
4885bool
4886ClpSimplex::sparseFactorization() const
4887{
4888     return factorization_->sparseThreshold() != 0;
4889}
4890void
4891ClpSimplex::setSparseFactorization(bool value)
4892{
4893     if (value) {
4894          if (!factorization_->sparseThreshold())
4895               factorization_->goSparse();
4896     } else {
4897          factorization_->sparseThreshold(0);
4898     }
4899}
4900void checkCorrect(ClpSimplex * /*model*/, int iRow,
4901                  const double * element, const int * rowStart, const int * rowLength,
4902                  const int * column,
4903                  const double * columnLower_, const double * columnUpper_,
4904                  int /*infiniteUpperC*/,
4905                  int /*infiniteLowerC*/,
4906                  double &maximumUpC,
4907                  double &maximumDownC)
4908{
4909     int infiniteUpper = 0;
4910     int infiniteLower = 0;
4911     double maximumUp = 0.0;
4912     double maximumDown = 0.0;
4913     CoinBigIndex rStart = rowStart[iRow];
4914     CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4915     CoinBigIndex j;
4916     double large = 1.0e15;
4917     int iColumn;
4918     // Compute possible lower and upper ranges
4919
4920     for (j = rStart; j < rEnd; ++j) {
4921          double value = element[j];
4922          iColumn = column[j];
4923          if (value > 0.0) {
4924               if (columnUpper_[iColumn] >= large) {
4925                    ++infiniteUpper;
4926               } else {
4927                    maximumUp += columnUpper_[iColumn] * value;
4928               }
4929               if (columnLower_[iColumn] <= -large) {
4930                    ++infiniteLower;
4931               } else {
4932                    maximumDown += columnLower_[iColumn] * value;
4933               }
4934          } else if (value < 0.0) {
4935               if (columnUpper_[iColumn] >= large) {
4936                    ++infiniteLower;
4937               } else {
4938                    maximumDown += columnUpper_[iColumn] * value;
4939               }
4940               if (columnLower_[iColumn] <= -large) {
4941                    ++infiniteUpper;
4942               } else {
4943                    maximumUp += columnLower_[iColumn] * value;
4944               }
4945          }
4946     }
4947     //assert (infiniteLowerC==infiniteLower);
4948     //assert (infiniteUpperC==infiniteUpper);
4949     if (fabs(maximumUp - maximumUpC) > 1.0e-12 * CoinMax(fabs(maximumUp), fabs(maximumUpC)))
4950          COIN_DETAIL_PRINT(printf("row %d comp up %g, true up %g\n", iRow,
4951                                   maximumUpC, maximumUp));
4952     if (fabs(maximumDown - maximumDownC) > 1.0e-12 * CoinMax(fabs(maximumDown), fabs(maximumDownC)))
4953          COIN_DETAIL_PRINT(printf("row %d comp down %g, true down %g\n", iRow,
4954                 maximumDownC, maximumDown));
4955     maximumUpC = maximumUp;
4956     maximumDownC = maximumDown;
4957}
4958
4959/* Tightens primal bounds to make dual faster.  Unless
4960   fixed, bounds are slightly looser than they could be.
4961   This is to make dual go faster and is probably not needed
4962   with a presolve.  Returns non-zero if problem infeasible
4963
4964   Fudge for branch and bound - put bounds on columns of factor *
4965   largest value (at continuous) - should improve stability
4966   in branch and bound on infeasible branches (0.0 is off)
4967*/
4968int
4969ClpSimplex::tightenPrimalBounds(double factor, int doTight, bool tightIntegers)
4970{
4971
4972     // Get a row copy in standard format
4973     CoinPackedMatrix copy;
4974     copy.setExtraGap(0.0);
4975     copy.setExtraMajor(0.0);
4976     copy.reverseOrderedCopyOf(*matrix());
4977     // Matrix may have been created so get rid of it
4978     matrix_->releasePackedMatrix();
4979     // get matrix data pointers
4980     const int * column = copy.getIndices();
4981     const CoinBigIndex * rowStart = copy.getVectorStarts();
4982     const int * rowLength = copy.getVectorLengths();
4983     const double * element = copy.getElements();
4984     int numberChanged = 1, iPass = 0;
4985     double large = largeValue(); // treat bounds > this as infinite
4986#ifndef NDEBUG
4987     double large2 = 1.0e10 * large;
4988#endif
4989     int numberInfeasible = 0;
4990     int totalTightened = 0;
4991
4992     double tolerance = primalTolerance();
4993
4994
4995     // Save column bounds
4996     double * saveLower = new double [numberColumns_];
4997     CoinMemcpyN(columnLower_, numberColumns_, saveLower);
4998     double * saveUpper = new double [numberColumns_];
4999     CoinMemcpyN(columnUpper_, numberColumns_, saveUpper);
5000
5001     int iRow, iColumn;
5002     // If wanted compute a reasonable dualBound_
5003     if (factor == COIN_DBL_MAX) {
5004          factor = 0.0;
5005          if (dualBound_ == 1.0e10) {
5006               // get largest scaled away from bound
5007               double largest = 1.0e-12;
5008               double largestScaled = 1.0e-12;
5009               int iRow;
5010               for (iRow = 0; iRow < numberRows_; iRow++) {
5011                    double value = rowActivity_[iRow];
5012                    double above = value - rowLower_[iRow];
5013                    double below = rowUpper_[iRow] - value;
5014                    if (above < 1.0e12) {
5015                         largest = CoinMax(largest, above);
5016                    }
5017                    if (below < 1.0e12) {
5018                         largest = CoinMax(largest, below);
5019                    }
5020                    if (rowScale_) {
5021                         double multiplier = rowScale_[iRow];
5022                         above *= multiplier;
5023                         below *= multiplier;
5024                    }
5025                    if (above < 1.0e12) {
5026                         largestScaled = CoinMax(largestScaled, above);
5027                    }
5028                    if (below < 1.0e12) {
5029                         largestScaled = CoinMax(largestScaled, below);
5030                    }
5031               }
5032
5033               int iColumn;
5034               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5035                    double value = columnActivity_[iColumn];
5036                    double above = value - columnLower_[iColumn];
5037                    double below = columnUpper_[iColumn] - value;
5038                    if (above < 1.0e12) {
5039                         largest = CoinMax(largest, above);
5040                    }
5041                    if (below < 1.0e12) {
5042                         largest = CoinMax(largest, below);
5043                    }
5044                    if (columnScale_) {
5045                         double multiplier = 1.0 / columnScale_[iColumn];
5046                         above *= multiplier;
5047                         below *= multiplier;
5048                    }
5049                    if (above < 1.0e12) {
5050                         largestScaled = CoinMax(largestScaled, above);
5051                    }
5052                    if (below < 1.0e12) {
5053                         largestScaled = CoinMax(largestScaled, below);
5054                    }
5055               }
5056               std::cout << "Largest (scaled) away from bound " << largestScaled
5057                         << " unscaled " << largest << std::endl;
5058               dualBound_ = CoinMax(1.0001e7, CoinMin(100.0 * largest, 1.00001e10));
5059          }
5060     }
5061
5062     // If wanted - tighten column bounds using solution
5063     if (factor) {
5064          double largest = 0.0;
5065          if (factor > 0.0) {
5066               assert (factor > 1.0);
5067               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5068                    if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
5069                         largest = CoinMax(largest, fabs(columnActivity_[iColumn]));
5070                    }
5071               }
5072               largest *= factor;
5073          } else {
5074               // absolute
5075               largest = - factor;
5076          }
5077          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5078               if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
5079                    columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn], largest);
5080                    columnLower_[iColumn] = CoinMax(columnLower_[iColumn], -largest);
5081               }
5082          }
5083     }
5084#define MAXPASS 10
5085
5086     // Loop round seeing if we can tighten bounds
5087     // Would be faster to have a stack of possible rows
5088     // and we put altered rows back on stack
5089     int numberCheck = -1;
5090     while(numberChanged > numberCheck) {
5091
5092          numberChanged = 0; // Bounds tightened this pass
5093
5094          if (iPass == MAXPASS) break;
5095          iPass++;
5096
5097          for (iRow = 0; iRow < numberRows_; iRow++) {
5098
5099               if (rowLower_[iRow] > -large || rowUpper_[iRow] < large) {
5100
5101                    // possible row
5102                    int infiniteUpper = 0;
5103                    int infiniteLower = 0;
5104                    double maximumUp = 0.0;
5105                    double maximumDown = 0.0;
5106                    double newBound;
5107                    CoinBigIndex rStart = rowStart[iRow];
5108                    CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
5109                    CoinBigIndex j;
5110                    // Compute possible lower and upper ranges
5111
5112                    for (j = rStart; j < rEnd; ++j) {
5113                         double value = element[j];
5114                         iColumn = column[j];
5115                         if (value > 0.0) {
5116                              if (columnUpper_[iColumn] >= large) {
5117                                   ++infiniteUpper;
5118                              } else {
5119                                   maximumUp += columnUpper_[iColumn] * value;
5120                              }
5121                              if (columnLower_[iColumn] <= -large) {
5122                                   ++infiniteLower;
5123                              } else {
5124                                   maximumDown += columnLower_[iColumn] * value;
5125                              }
5126                         } else if (value < 0.0) {
5127                              if (columnUpper_[iColumn] >= large) {
5128                                   ++infiniteLower;
5129                              } else {
5130                                   maximumDown += columnUpper_[iColumn] * value;
5131                              }
5132                              if (columnLower_[iColumn] <= -large) {
5133                                   ++infiniteUpper;
5134                              } else {
5135                                   maximumUp += columnLower_[iColumn] * value;
5136                              }
5137                         }
5138                    }
5139                    // Build in a margin of error
5140                    maximumUp += 1.0e-8 * fabs(maximumUp);
5141                    maximumDown -= 1.0e-8 * fabs(maximumDown);
5142                    double maxUp = maximumUp + infiniteUpper * 1.0e31;
5143                    double maxDown = maximumDown - infiniteLower * 1.0e31;
5144                    if (maxUp <= rowUpper_[iRow] + tolerance &&
5145                              maxDown >= rowLower_[iRow] - tolerance) {
5146
5147                         // Row is redundant - make totally free
5148                         // NO - can't do this for postsolve
5149                         // rowLower_[iRow]=-COIN_DBL_MAX;
5150                         // rowUpper_[iRow]=COIN_DBL_MAX;
5151                         //printf("Redundant row in presolveX %d\n",iRow);
5152
5153                    } else {
5154                         if (maxUp < rowLower_[iRow] - 100.0 * tolerance ||
5155                                   maxDown > rowUpper_[iRow] + 100.0 * tolerance) {
5156                              // problem is infeasible - exit at once
5157                              numberInfeasible++;
5158                              break;
5159                         }
5160                         double lower = rowLower_[iRow];
5161                         double upper = rowUpper_[iRow];
5162                         for (j = rStart; j < rEnd; ++j) {
5163                              double value = element[