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

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

changes for advanced use of Clp

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