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

Last change on this file since 1691 was 1691, checked in by stefan, 9 years ago

check for AMD and CHOLMOD and GLPK's AMD in configure; compile ClpCholeskyUfl? only if AMD, CHOLMOD, or GLPK available

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