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

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

various fixes plus slightly weighted pricing plus lagomory switches

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