source: stable/1.14/Clp/src/ClpSimplex.cpp @ 1866

Last change on this file since 1866 was 1866, checked in by stefan, 8 years ago

recompute reduced costs if factorization is reused but objective has changed

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