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

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

patches to try and make parametrics faster

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