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

Last change on this file since 1723 was 1723, checked in by forrest, 8 years ago

out some printf statements

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