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

Last change on this file since 1696 was 1696, checked in by forrest, 9 years ago

for integer problems with large values

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