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

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

adjust to changes in CoinUtils? header files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 495.4 KB
<
RevLine 
[1370]1/* $Id: ClpSimplex.cpp 1722 2011-04-17 09:58:37Z stefan $ */
[2]2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1665]4// This code is licensed under the terms of the Eclipse Public License (EPL).
[2]5
[655]6//#undef NDEBUG
[2]7
[772]8#include "ClpConfig.h"
[760]9
[63]10#include "CoinPragma.hpp"
[2]11#include <math.h>
12
[651]13#if SLIM_CLP==2
14#define SLIM_NOIO
15#endif
[2]16#include "CoinHelperFunctions.hpp"
[1722]17#include "CoinFloatEqual.hpp"
[2]18#include "ClpSimplex.hpp"
19#include "ClpFactorization.hpp"
[50]20#include "ClpPackedMatrix.hpp"
21#include "CoinIndexedVector.hpp"
22#include "ClpDualRowDantzig.hpp"
[2]23#include "ClpDualRowSteepest.hpp"
24#include "ClpPrimalColumnDantzig.hpp"
25#include "ClpPrimalColumnSteepest.hpp"
26#include "ClpNonLinearCost.hpp"
27#include "ClpMessage.hpp"
[344]28#include "ClpEventHandler.hpp"
[119]29#include "ClpLinearObjective.hpp"
[286]30#include "ClpHelperFunctions.hpp"
[546]31#include "CoinModel.hpp"
[1034]32#include "CoinLpIO.hpp"
[2]33#include <cfloat>
34
35#include <string>
36#include <stdio.h>
37#include <iostream>
38//#############################################################################
39
[1034]40ClpSimplex::ClpSimplex (bool emptyMessages) :
[2]41
[1525]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)
[2]124{
[1525]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
[2]145}
146
[225]147// Subproblem constructor
148ClpSimplex::ClpSimplex ( const ClpModel * rhs,
[1525]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)
[225]235{
[1525]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     }
[225]294}
[1154]295// Subproblem constructor
296ClpSimplex::ClpSimplex ( const ClpSimplex * rhs,
[1525]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)
[1154]383{
[1525]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     }
[1154]457}
[618]458// Puts solution back small model
[1525]459void
460ClpSimplex::getbackSolution(const ClpSimplex & smallModel, const int * whichRow, const int * whichColumn)
[618]461{
[1525]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_);
[1286]488#if 0
[1525]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     }
[1286]503#endif
[1525]504     matrix()->times(columnActivity_, rowActivity_) ;
[618]505}
[2]506
507//-----------------------------------------------------------------------------
508
509ClpSimplex::~ClpSimplex ()
510{
[1525]511     setPersistenceFlag(0);
512     gutsOfDelete(0);
513     delete nonLinearCost_;
[2]514}
515//#############################################################################
[1525]516void ClpSimplex::setLargeValue( double value)
[2]517{
[1525]518     if (value > 0.0 && value < COIN_DBL_MAX)
519          largeValue_ = value;
[2]520}
521int
[225]522ClpSimplex::gutsOfSolution ( double * givenDuals,
[1525]523                             const double * givenPrimals,
524                             bool valuesPass)
[2]525{
526
527
[1525]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();
[50]573#ifdef CLP_DEBUG
[1525]574          std::cout << "Largest given infeasibility " << oldValue
575                    << " now " << nonLinearCost_->largestInfeasibility() << std::endl;
[50]576#endif
[1525]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];
[2]601
[1525]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);
[1439]615#if 0
[1525]616                    allSlackBasis(true);
617                    CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
[1439]618#else
[1525]619                    // allow
620                    numberOut = 0;
[1439]621#endif
[1525]622               }
623               CoinSort_2(save, save + numberOut, sort);
624               numberOut = CoinMin(1000, numberOut);
625               for (iRow = 0; iRow < numberOut; iRow++) {
626                    int jRow = sort[iRow];
627                    int iColumn = pivotVariable_[jRow];
628                    setColumnStatus(iColumn, superBasic);
629                    setRowStatus(jRow, basic);
630                    pivotVariable_[jRow] = jRow + numberColumns_;
631                    if (fabs(solution_[iColumn]) > 1.0e10) {
632                         if (upper_[iColumn] < 0.0) {
633                              solution_[iColumn] = upper_[iColumn];
634                         } else if (lower_[iColumn] > 0.0) {
635                              solution_[iColumn] = lower_[iColumn];
636                         } else {
637                              solution_[iColumn] = 0.0;
638                         }
639                    }
640               }
641               delete [] sort;
642          }
643          delete [] save;
644          if (numberOut)
645               return numberOut;
646     }
647     if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
648          //printf("trying feas pump\n");
649          const char * integerType = integerInformation();
650          assert (integerType);
651          assert (perturbationArray_);
652          CoinZeroN(cost_, numberRows_ + numberColumns_);
653          for (int i = 0; i < numberRows_ - numberRows_; i++) {
654               int iSequence = pivotVariable_[i];
655               if (iSequence < numberColumns_ && integerType[iSequence]) {
656                    double lower = lower_[iSequence];
657                    double upper = upper_[iSequence];
658                    double value = solution_[iSequence];
659                    if (value >= lower - primalTolerance_ &&
660                              value <= upper + primalTolerance_) {
661                         double sign;
662                         if (value - lower < upper - value)
663                              sign = 1.0;
664                         else
665                              sign = -1.0;
666                         cost_[iSequence] = sign * perturbationArray_[iSequence];
667                    }
668               }
669          }
670     }
671     computeDuals(givenDuals);
672     if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
673          const char * integerType = integerInformation();
674          // Need to do columns and rows to stay dual feasible
675          for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
676               if (integerType[iSequence] && getStatus(iSequence) != basic) {
677                    double djValue = dj_[iSequence];
678                    double change = 0.0;
679                    if (getStatus(iSequence) == atLowerBound)
680                         change = CoinMax(-djValue, 10.0 * perturbationArray_[iSequence]);
681                    else if (getStatus(iSequence) == atUpperBound)
682                         change = CoinMin(-djValue, -10.0 * perturbationArray_[iSequence]);
683                    cost_[iSequence] = change;
684                    dj_[iSequence] += change;
685               }
686          }
687     }
[2]688
[1525]689     // now check solutions
690     //checkPrimalSolution( rowActivityWork_, columnActivityWork_);
691     //checkDualSolution();
692     checkBothSolutions();
693     objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_);
694     if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 ||
695                                      largestDualError_ > 1.0e-2))
696          handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
697                    << largestPrimalError_
698                    << largestDualError_
699                    << CoinMessageEol;
700     if (largestPrimalError_ > 1.0e-1 && numberRows_ > 100 && numberIterations_) {
701          // Change factorization tolerance
702          if (factorization_->zeroTolerance() > 1.0e-18)
703               factorization_->zeroTolerance(1.0e-18);
704     }
705     // Switch off false values pass indicator
706     if (!valuesPass && algorithm_ > 0)
707          firstFree_ = -1;
708     return 0;
[2]709}
710void
711ClpSimplex::computePrimals ( const double * rowActivities,
[1525]712                             const double * columnActivities)
[2]713{
714
[1525]715     //work space
716     CoinIndexedVector  * workSpace = rowArray_[0];
[2]717
[1525]718     CoinIndexedVector * arrayVector = rowArray_[1];
719     arrayVector->clear();
720     CoinIndexedVector * previousVector = rowArray_[2];
721     previousVector->clear();
722     // accumulate non basic stuff
[2]723
[1525]724     int iRow;
725     // order is this way for scaling
726     if (columnActivities != columnActivityWork_)
727          ClpDisjointCopyN(columnActivities, numberColumns_, columnActivityWork_);
728     if (rowActivities != rowActivityWork_)
729          ClpDisjointCopyN(rowActivities, numberRows_, rowActivityWork_);
730     double * array = arrayVector->denseVector();
731     int * index = arrayVector->getIndices();
732     int number = 0;
733     const double * rhsOffset = matrix_->rhsOffset(this, false, true);
734     if (!rhsOffset) {
735          // Use whole matrix every time to make it easier for ClpMatrixBase
736          // So zero out basic
737          for (iRow = 0; iRow < numberRows_; iRow++) {
738               int iPivot = pivotVariable_[iRow];
739               assert (iPivot >= 0);
740               solution_[iPivot] = 0.0;
[1360]741#ifdef CLP_INVESTIGATE
[1525]742               assert (getStatus(iPivot) == basic);
[1360]743#endif
[1525]744          }
745          // Extended solution before "update"
746          matrix_->primalExpanded(this, 0);
747          times(-1.0, columnActivityWork_, array);
748          for (iRow = 0; iRow < numberRows_; iRow++) {
749               double value = array[iRow] + rowActivityWork_[iRow];
750               if (value) {
751                    array[iRow] = value;
752                    index[number++] = iRow;
753               } else {
754                    array[iRow] = 0.0;
755               }
756          }
757     } else {
758          // we have an effective rhs lying around
759          // zero out basic (really just for slacks)
760          for (iRow = 0; iRow < numberRows_; iRow++) {
761               int iPivot = pivotVariable_[iRow];
762               solution_[iPivot] = 0.0;
763          }
764          for (iRow = 0; iRow < numberRows_; iRow++) {
765               double value = rhsOffset[iRow] + rowActivityWork_[iRow];
766               if (value) {
767                    array[iRow] = value;
768                    index[number++] = iRow;
769               } else {
770                    array[iRow] = 0.0;
771               }
772          }
773     }
774     arrayVector->setNumElements(number);
[336]775#ifdef CLP_DEBUG
[1525]776     if (numberIterations_ == -3840) {
777          int i;
778          for (i = 0; i < numberRows_ + numberColumns_; i++)
779               printf("%d status %d\n", i, status_[i]);
780          printf("xxxxx1\n");
781          for (i = 0; i < numberRows_; i++)
782               if (array[i])
783                    printf("%d rhs %g\n", i, array[i]);
784          printf("xxxxx2\n");
785          for (i = 0; i < numberRows_ + numberColumns_; i++)
786               if (getStatus(i) != basic)
787                    printf("%d non basic %g %g %g\n", i, lower_[i], solution_[i], upper_[i]);
788          printf("xxxxx3\n");
789     }
[336]790#endif
[1525]791     // Ftran adjusted RHS and iterate to improve accuracy
792     double lastError = COIN_DBL_MAX;
793     int iRefine;
794     CoinIndexedVector * thisVector = arrayVector;
795     CoinIndexedVector * lastVector = previousVector;
796     if (number)
797          factorization_->updateColumn(workSpace, thisVector);
798     double * work = workSpace->denseVector();
[336]799#ifdef CLP_DEBUG
[1525]800     if (numberIterations_ == -3840) {
801          int i;
802          for (i = 0; i < numberRows_; i++)
803               if (array[i])
804                    printf("%d after rhs %g\n", i, array[i]);
805          printf("xxxxx4\n");
806     }
[336]807#endif
[1525]808     bool goodSolution = true;
809     for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
[2]810
[1525]811          int numberIn = thisVector->getNumElements();
812          int * indexIn = thisVector->getIndices();
813          double * arrayIn = thisVector->denseVector();
814          // put solution in correct place
815          if (!rhsOffset) {
816               int j;
817               for (j = 0; j < numberIn; j++) {
818                    iRow = indexIn[j];
819                    int iPivot = pivotVariable_[iRow];
820                    solution_[iPivot] = arrayIn[iRow];
821                    //assert (fabs(solution_[iPivot])<1.0e100);
822               }
823          } else {
824               for (iRow = 0; iRow < numberRows_; iRow++) {
825                    int iPivot = pivotVariable_[iRow];
826                    solution_[iPivot] = arrayIn[iRow];
827                    //assert (fabs(solution_[iPivot])<1.0e100);
828               }
829          }
830          // Extended solution after "update"
831          matrix_->primalExpanded(this, 1);
832          // check Ax == b  (for all)
833          // signal column generated matrix to just do basic (and gub)
834          unsigned int saveOptions = specialOptions();
835          setSpecialOptions(16);
836          times(-1.0, columnActivityWork_, work);
837          setSpecialOptions(saveOptions);
838          largestPrimalError_ = 0.0;
839          double multiplier = 131072.0;
840          for (iRow = 0; iRow < numberRows_; iRow++) {
841               double value = work[iRow] + rowActivityWork_[iRow];
842               work[iRow] = value * multiplier;
843               if (fabs(value) > largestPrimalError_) {
844                    largestPrimalError_ = fabs(value);
845               }
846          }
847          if (largestPrimalError_ >= lastError) {
848               // restore
849               CoinIndexedVector * temp = thisVector;
850               thisVector = lastVector;
851               lastVector = temp;
852               goodSolution = false;
853               break;
854          }
855          if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
856               // try and make better
857               // save this
858               CoinIndexedVector * temp = thisVector;
859               thisVector = lastVector;
860               lastVector = temp;
861               int * indexOut = thisVector->getIndices();
862               int number = 0;
863               array = thisVector->denseVector();
864               thisVector->clear();
865               for (iRow = 0; iRow < numberRows_; iRow++) {
866                    double value = work[iRow];
867                    if (value) {
868                         array[iRow] = value;
869                         indexOut[number++] = iRow;
870                         work[iRow] = 0.0;
871                    }
872               }
873               thisVector->setNumElements(number);
874               lastError = largestPrimalError_;
875               factorization_->updateColumn(workSpace, thisVector);
876               multiplier = 1.0 / multiplier;
877               double * previous = lastVector->denseVector();
878               number = 0;
879               for (iRow = 0; iRow < numberRows_; iRow++) {
880                    double value = previous[iRow] + multiplier * array[iRow];
881                    if (value) {
882                         array[iRow] = value;
883                         indexOut[number++] = iRow;
884                    } else {
885                         array[iRow] = 0.0;
886                    }
887               }
888               thisVector->setNumElements(number);
889          } else {
890               break;
891          }
892     }
[2]893
[1525]894     // solution as accurate as we are going to get
895     ClpFillN(work, numberRows_, 0.0);
896     if (!goodSolution) {
897          array = thisVector->denseVector();
898          // put solution in correct place
899          for (iRow = 0; iRow < numberRows_; iRow++) {
900               int iPivot = pivotVariable_[iRow];
901               solution_[iPivot] = array[iRow];
902               //assert (fabs(solution_[iPivot])<1.0e100);
903          }
904     }
905     arrayVector->clear();
906     previousVector->clear();
[336]907#ifdef CLP_DEBUG
[1525]908     if (numberIterations_ == -3840) {
909          exit(77);
910     }
[336]911#endif
[2]912}
913// now dual side
914void
[137]915ClpSimplex::computeDuals(double * givenDjs)
[2]916{
[651]917#ifndef SLIM_CLP
[1525]918     if (objective_->type() == 1 || !objective_->activated()) {
[651]919#endif
[1525]920          // Linear
921          //work space
922          CoinIndexedVector  * workSpace = rowArray_[0];
923
924          CoinIndexedVector * arrayVector = rowArray_[1];
925          arrayVector->clear();
926          CoinIndexedVector * previousVector = rowArray_[2];
927          previousVector->clear();
928          int iRow;
[2]929#ifdef CLP_DEBUG
[1525]930          workSpace->checkClear();
[2]931#endif
[1525]932          double * array = arrayVector->denseVector();
933          int * index = arrayVector->getIndices();
934          int number = 0;
935          if (!givenDjs) {
936               for (iRow = 0; iRow < numberRows_; iRow++) {
937                    int iPivot = pivotVariable_[iRow];
938                    double value = cost_[iPivot];
939                    if (value) {
940                         array[iRow] = value;
941                         index[number++] = iRow;
942                    }
943               }
944          } else {
945               // dual values pass - djs may not be zero
946               for (iRow = 0; iRow < numberRows_; iRow++) {
947                    int iPivot = pivotVariable_[iRow];
948                    // make sure zero if done
949                    if (!pivoted(iPivot))
950                         givenDjs[iPivot] = 0.0;
951                    double value = cost_[iPivot] - givenDjs[iPivot];
952                    if (value) {
953                         array[iRow] = value;
954                         index[number++] = iRow;
955                    }
956               }
957          }
958          arrayVector->setNumElements(number);
959          // Extended duals before "updateTranspose"
960          matrix_->dualExpanded(this, arrayVector, givenDjs, 0);
961
962          // Btran basic costs and get as accurate as possible
963          double lastError = COIN_DBL_MAX;
964          int iRefine;
965          double * work = workSpace->denseVector();
966          CoinIndexedVector * thisVector = arrayVector;
967          CoinIndexedVector * lastVector = previousVector;
968          factorization_->updateColumnTranspose(workSpace, thisVector);
969
970          for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
971               // check basic reduced costs zero
972               largestDualError_ = 0.0;
973               if (!numberExtraRows_) {
974                    // Just basic
975                    int * index2 = workSpace->getIndices();
976                    // use reduced costs for slacks as work array
977                    double * work2 = reducedCostWork_ + numberColumns_;
978                    int numberStructurals = 0;
979                    for (iRow = 0; iRow < numberRows_; iRow++) {
980                         int iPivot = pivotVariable_[iRow];
981                         if (iPivot < numberColumns_)
982                              index2[numberStructurals++] = iPivot;
983                    }
984                    matrix_->listTransposeTimes(this, array, index2, numberStructurals, work2);
985                    numberStructurals = 0;
986                    if (!givenDjs) {
987                         for (iRow = 0; iRow < numberRows_; iRow++) {
988                              int iPivot = pivotVariable_[iRow];
989                              double value;
990                              if (iPivot >= numberColumns_) {
991                                   // slack
992                                   value = rowObjectiveWork_[iPivot-numberColumns_]
993                                           + array[iPivot-numberColumns_];
994                              } else {
995                                   // column
996                                   value = objectiveWork_[iPivot] - work2[numberStructurals++];
997                              }
998                              work[iRow] = value;
999                              if (fabs(value) > largestDualError_) {
1000                                   largestDualError_ = fabs(value);
1001                              }
1002                         }
1003                    } else {
1004                         for (iRow = 0; iRow < numberRows_; iRow++) {
1005                              int iPivot = pivotVariable_[iRow];
1006                              if (iPivot >= numberColumns_) {
1007                                   // slack
1008                                   work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
1009                                                + array[iPivot-numberColumns_] - givenDjs[iPivot];
1010                              } else {
1011                                   // column
1012                                   work[iRow] = objectiveWork_[iPivot] - work2[numberStructurals++]
1013                                                - givenDjs[iPivot];
1014                              }
1015                              if (fabs(work[iRow]) > largestDualError_) {
1016                                   largestDualError_ = fabs(work[iRow]);
1017                                   //assert (largestDualError_<1.0e-7);
1018                                   //if (largestDualError_>1.0e-7)
1019                                   //printf("large dual error %g\n",largestDualError_);
1020                              }
1021                         }
1022                    }
1023               } else {
1024                    // extra rows - be more careful
[451]1025#if 1
[1525]1026                    // would be faster to do just for basic but this reduces code
1027                    ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
1028                    transposeTimes(-1.0, array, reducedCostWork_);
[451]1029#else
[1525]1030                    // Just basic
1031                    int * index2 = workSpace->getIndices();
1032                    int numberStructurals = 0;
1033                    for (iRow = 0; iRow < numberRows_; iRow++) {
1034                         int iPivot = pivotVariable_[iRow];
1035                         if (iPivot < numberColumns_)
1036                              index2[numberStructurals++] = iPivot;
1037                    }
1038                    matrix_->listTransposeTimes(this, array, index2, numberStructurals, work);
1039                    for (iRow = 0; iRow < numberStructurals; iRow++) {
1040                         int iPivot = index2[iRow];
1041                         reducedCostWork_[iPivot] = objectiveWork_[iPivot] - work[iRow];
1042                    }
[451]1043#endif
[1525]1044                    // update by duals on sets
1045                    matrix_->dualExpanded(this, NULL, NULL, 1);
1046                    if (!givenDjs) {
1047                         for (iRow = 0; iRow < numberRows_; iRow++) {
1048                              int iPivot = pivotVariable_[iRow];
1049                              double value;
1050                              if (iPivot >= numberColumns_) {
1051                                   // slack
1052                                   value = rowObjectiveWork_[iPivot-numberColumns_]
1053                                           + array[iPivot-numberColumns_];
1054                              } else {
1055                                   // column
1056                                   value = reducedCostWork_[iPivot];
1057                              }
1058                              work[iRow] = value;
1059                              if (fabs(value) > largestDualError_) {
1060                                   largestDualError_ = fabs(value);
1061                              }
1062                         }
1063                    } else {
1064                         for (iRow = 0; iRow < numberRows_; iRow++) {
1065                              int iPivot = pivotVariable_[iRow];
1066                              if (iPivot >= numberColumns_) {
1067                                   // slack
1068                                   work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
1069                                                + array[iPivot-numberColumns_] - givenDjs[iPivot];
1070                              } else {
1071                                   // column
1072                                   work[iRow] = reducedCostWork_[iPivot] - givenDjs[iPivot];
1073                              }
1074                              if (fabs(work[iRow]) > largestDualError_) {
1075                                   largestDualError_ = fabs(work[iRow]);
1076                                   //assert (largestDualError_<1.0e-7);
1077                                   //if (largestDualError_>1.0e-7)
1078                                   //printf("large dual error %g\n",largestDualError_);
1079                              }
1080                         }
1081                    }
1082               }
1083               if (largestDualError_ >= lastError) {
1084                    // restore
1085                    CoinIndexedVector * temp = thisVector;
1086                    thisVector = lastVector;
1087                    lastVector = temp;
1088                    break;
1089               }
1090               if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
1091                         && !givenDjs) {
1092                    // try and make better
1093                    // save this
1094                    CoinIndexedVector * temp = thisVector;
1095                    thisVector = lastVector;
1096                    lastVector = temp;
1097                    int * indexOut = thisVector->getIndices();
1098                    int number = 0;
1099                    array = thisVector->denseVector();
1100                    thisVector->clear();
1101                    double multiplier = 131072.0;
1102                    for (iRow = 0; iRow < numberRows_; iRow++) {
1103                         double value = multiplier * work[iRow];
1104                         if (value) {
1105                              array[iRow] = value;
1106                              indexOut[number++] = iRow;
1107                              work[iRow] = 0.0;
1108                         }
1109                         work[iRow] = 0.0;
1110                    }
1111                    thisVector->setNumElements(number);
1112                    lastError = largestDualError_;
1113                    factorization_->updateColumnTranspose(workSpace, thisVector);
1114                    multiplier = 1.0 / multiplier;
1115                    double * previous = lastVector->denseVector();
1116                    number = 0;
1117                    for (iRow = 0; iRow < numberRows_; iRow++) {
1118                         double value = previous[iRow] + multiplier * array[iRow];
1119                         if (value) {
1120                              array[iRow] = value;
1121                              indexOut[number++] = iRow;
1122                         } else {
1123                              array[iRow] = 0.0;
1124                         }
1125                    }
1126                    thisVector->setNumElements(number);
1127               } else {
1128                    break;
1129               }
1130          }
1131          // now look at dual solution
1132          array = thisVector->denseVector();
1133          for (iRow = 0; iRow < numberRows_; iRow++) {
1134               // slack
1135               double value = array[iRow];
1136               dual_[iRow] = value;
1137               value += rowObjectiveWork_[iRow];
1138               rowReducedCost_[iRow] = value;
1139          }
1140          // can use work if problem scaled (for better cache)
1141          ClpPackedMatrix* clpMatrix =
1142               dynamic_cast< ClpPackedMatrix*>(matrix_);
1143          double * saveRowScale = rowScale_;
1144          //double * saveColumnScale = columnScale_;
1145          if (scaledMatrix_) {
1146               rowScale_ = NULL;
1147               clpMatrix = scaledMatrix_;
1148          }
1149          if (clpMatrix && (clpMatrix->flags() & 2) == 0) {
1150               CoinIndexedVector * cVector = columnArray_[0];
1151               int * whichColumn = cVector->getIndices();
1152               assert (!cVector->getNumElements());
1153               int n = 0;
1154               for (int i = 0; i < numberColumns_; i++) {
1155                    if (getColumnStatus(i) != basic) {
1156                         whichColumn[n++] = i;
1157                         reducedCostWork_[i] = objectiveWork_[i];
1158                    } else {
1159                         reducedCostWork_[i] = 0.0;
1160                    }
1161               }
1162               if (numberRows_ > 4000)
1163                    clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
1164                                                    rowScale_, columnScale_, work);
1165               else
1166                    clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
1167                                                    rowScale_, columnScale_, NULL);
1168          } else {
1169               ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
1170               if (numberRows_ > 4000)
1171                    matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
1172                                            rowScale_, columnScale_, work);
1173               else
1174                    matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
1175                                            rowScale_, columnScale_, NULL);
1176          }
1177          rowScale_ = saveRowScale;
1178          //columnScale_ = saveColumnScale;
1179          ClpFillN(work, numberRows_, 0.0);
1180          // Extended duals and check dual infeasibility
1181          if (!matrix_->skipDualCheck() || algorithm_ < 0 || problemStatus_ != -2)
1182               matrix_->dualExpanded(this, NULL, NULL, 2);
1183          // If necessary - override results
1184          if (givenDjs) {
1185               // restore accurate duals
1186               CoinMemcpyN(dj_, (numberRows_ + numberColumns_), givenDjs);
1187          }
1188          arrayVector->clear();
1189          previousVector->clear();
[651]1190#ifndef SLIM_CLP
[1525]1191     } else {
1192          // Nonlinear
1193          objective_->reducedGradient(this, dj_, false);
1194          // get dual_ by moving from reduced costs for slacks
1195          CoinMemcpyN(dj_ + numberColumns_, numberRows_, dual_);
1196     }
[651]1197#endif
[2]1198}
[1525]1199/* Given an existing factorization computes and checks
[2]1200   primal and dual solutions.  Uses input arrays for variables at
1201   bounds.  Returns feasibility states */
[1402]1202int ClpSimplex::getSolution ( const double * /*rowActivities*/,
[1525]1203                              const double * /*columnActivities*/)
[2]1204{
[1525]1205     if (!factorization_->status()) {
1206          // put in standard form
1207          createRim(7 + 8 + 16 + 32, false, -1);
1208          if (pivotVariable_[0] < 0)
1209               internalFactorize(0);
1210          // do work
1211          gutsOfSolution ( NULL, NULL);
1212          // release extra memory
1213          deleteRim(0);
1214     }
1215     return factorization_->status();
[2]1216}
[1525]1217/* Given an existing factorization computes and checks
[2]1218   primal and dual solutions.  Uses current problem arrays for
1219   bounds.  Returns feasibility states */
1220int ClpSimplex::getSolution ( )
1221{
[1525]1222     double * rowActivities = new double[numberRows_];
1223     double * columnActivities = new double[numberColumns_];
1224     ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
1225     ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
1226     int status = getSolution( rowActivities, columnActivities);
1227     delete [] rowActivities;
1228     delete [] columnActivities;
1229     return status;
[2]1230}
1231// Factorizes using current basis.  This is for external use
1232// Return codes are as from ClpFactorization
1233int ClpSimplex::factorize ()
1234{
[1525]1235     // put in standard form
1236     createRim(7 + 8 + 16 + 32, false);
1237     // do work
1238     int status = internalFactorize(-1);
1239     // release extra memory
1240     deleteRim(0);
[2]1241
[1525]1242     return status;
[2]1243}
[225]1244// Clean up status
[1525]1245void
[225]1246ClpSimplex::cleanStatus()
1247{
[1525]1248     int iRow, iColumn;
1249     int numberBasic = 0;
1250     // make row activities correct
1251     memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
1252     times(1.0, columnActivityWork_, rowActivityWork_);
1253     if (!status_)
1254          createStatus();
1255     for (iRow = 0; iRow < numberRows_; iRow++) {
1256          if (getRowStatus(iRow) == basic)
1257               numberBasic++;
1258          else {
1259               setRowStatus(iRow, superBasic);
1260               // but put to bound if close
1261               if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
1262                         <= primalTolerance_) {
1263                    rowActivityWork_[iRow] = rowLowerWork_[iRow];
1264                    setRowStatus(iRow, atLowerBound);
1265               } else if (fabs(rowActivityWork_[iRow] - rowUpperWork_[iRow])
1266                          <= primalTolerance_) {
1267                    rowActivityWork_[iRow] = rowUpperWork_[iRow];
1268                    setRowStatus(iRow, atUpperBound);
1269               }
1270          }
1271     }
1272     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1273          if (getColumnStatus(iColumn) == basic) {
1274               if (numberBasic == numberRows_) {
1275                    // take out of basis
1276                    setColumnStatus(iColumn, superBasic);
1277                    // but put to bound if close
1278                    if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1279                              <= primalTolerance_) {
1280                         columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1281                         setColumnStatus(iColumn, atLowerBound);
1282                    } else if (fabs(columnActivityWork_[iColumn]
1283                                    - columnUpperWork_[iColumn])
1284                               <= primalTolerance_) {
1285                         columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1286                         setColumnStatus(iColumn, atUpperBound);
1287                    }
1288               } else
1289                    numberBasic++;
1290          } else {
1291               setColumnStatus(iColumn, superBasic);
1292               // but put to bound if close
1293               if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1294                         <= primalTolerance_) {
1295                    columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1296                    setColumnStatus(iColumn, atLowerBound);
1297               } else if (fabs(columnActivityWork_[iColumn]
1298                               - columnUpperWork_[iColumn])
1299                          <= primalTolerance_) {
1300                    columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1301                    setColumnStatus(iColumn, atUpperBound);
1302               }
1303          }
1304     }
[225]1305}
[2]1306
[1525]1307/* Factorizes using current basis.
1308   solveType - 1 iterating, 0 initial, -1 external
[509]1309   - 2 then iterating but can throw out of basis
1310   If 10 added then in primal values pass
1311   Return codes are as from ClpFactorization unless initial factorization
1312   when total number of singularities is returned.
1313   Special case is numberRows_+1 -> all slack basis.
[2]1314*/
1315int ClpSimplex::internalFactorize ( int solveType)
1316{
[1525]1317     int iRow, iColumn;
1318     int totalSlacks = numberRows_;
1319     if (!status_)
1320          createStatus();
[2]1321
[1525]1322     bool valuesPass = false;
1323     if (solveType >= 10) {
1324          valuesPass = true;
1325          solveType -= 10;
1326     }
[92]1327#ifdef CLP_DEBUG
[1525]1328     if (solveType > 0) {
1329          int numberFreeIn = 0, numberFreeOut = 0;
1330          double biggestDj = 0.0;
1331          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1332               switch(getColumnStatus(iColumn)) {
[225]1333
[1525]1334               case basic:
1335                    if (columnLower_[iColumn] < -largeValue_
1336                              && columnUpper_[iColumn] > largeValue_)
1337                         numberFreeIn++;
1338                    break;
1339               default:
1340                    if (columnLower_[iColumn] < -largeValue_
1341                              && columnUpper_[iColumn] > largeValue_) {
1342                         numberFreeOut++;
1343                         biggestDj = CoinMax(fabs(dj_[iColumn]), biggestDj);
1344                    }
1345                    break;
1346               }
1347          }
1348          if (numberFreeIn + numberFreeOut)
1349               printf("%d in basis, %d out - largest dj %g\n",
1350                      numberFreeIn, numberFreeOut, biggestDj);
1351     }
[92]1352#endif
[1525]1353     if (solveType <= 0) {
1354          // Make sure everything is clean
1355          for (iRow = 0; iRow < numberRows_; iRow++) {
1356               if(getRowStatus(iRow) == isFixed) {
1357                    // double check fixed
1358                    if (rowUpperWork_[iRow] > rowLowerWork_[iRow])
1359                         setRowStatus(iRow, atLowerBound);
1360               } else if (getRowStatus(iRow) == isFree) {
1361                    // may not be free after all
1362                    if (rowLowerWork_[iRow] > -largeValue_ || rowUpperWork_[iRow] < largeValue_)
1363                         setRowStatus(iRow, superBasic);
1364               }
1365          }
1366          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1367               if(getColumnStatus(iColumn) == isFixed) {
1368                    // double check fixed
1369                    if (columnUpperWork_[iColumn] > columnLowerWork_[iColumn])
1370                         setColumnStatus(iColumn, atLowerBound);
1371               } else if (getColumnStatus(iColumn) == isFree) {
1372                    // may not be free after all
1373                    if (columnLowerWork_[iColumn] > -largeValue_ || columnUpperWork_[iColumn] < largeValue_)
1374                         setColumnStatus(iColumn, superBasic);
1375               }
1376          }
1377          if (!valuesPass) {
1378               // not values pass so set to bounds
1379               bool allSlack = true;
1380               if (status_) {
1381                    for (iRow = 0; iRow < numberRows_; iRow++) {
1382                         if (getRowStatus(iRow) != basic) {
1383                              allSlack = false;
1384                              break;
1385                         }
1386                    }
1387               }
1388               if (!allSlack) {
1389                    //#define CLP_INVESTIGATE2
[1449]1390#ifdef CLP_INVESTIGATE3
[1525]1391                    int numberTotal = numberRows_ + numberColumns_;
1392                    double * saveSol = valuesPass ?
1393                                       CoinCopyOfArray(solution_, numberTotal) : NULL;
[1429]1394#endif
[1525]1395                    // set values from warm start (if sensible)
1396                    int numberBasic = 0;
1397                    for (iRow = 0; iRow < numberRows_; iRow++) {
1398                         switch(getRowStatus(iRow)) {
[92]1399
[1525]1400                         case basic:
1401                              numberBasic++;
1402                              break;
1403                         case atUpperBound:
1404                              rowActivityWork_[iRow] = rowUpperWork_[iRow];
1405                              if (rowActivityWork_[iRow] > largeValue_) {
1406                                   if (rowLowerWork_[iRow] > -largeValue_) {
1407                                        rowActivityWork_[iRow] = rowLowerWork_[iRow];
1408                                        setRowStatus(iRow, atLowerBound);
1409                                   } else {
1410                                        // say free
1411                                        setRowStatus(iRow, isFree);
1412                                        rowActivityWork_[iRow] = 0.0;
1413                                   }
1414                              }
1415                              break;
1416                         case ClpSimplex::isFixed:
1417                         case atLowerBound:
1418                              rowActivityWork_[iRow] = rowLowerWork_[iRow];
1419                              if (rowActivityWork_[iRow] < -largeValue_) {
1420                                   if (rowUpperWork_[iRow] < largeValue_) {
1421                                        rowActivityWork_[iRow] = rowUpperWork_[iRow];
1422                                        setRowStatus(iRow, atUpperBound);
1423                                   } else {
1424                                        // say free
1425                                        setRowStatus(iRow, isFree);
1426                                        rowActivityWork_[iRow] = 0.0;
1427                                   }
1428                              }
1429                              break;
1430                         case isFree:
1431                              break;
1432                              // not really free - fall through to superbasic
1433                         case superBasic:
1434                              if (rowUpperWork_[iRow] > largeValue_) {
1435                                   if (rowLowerWork_[iRow] > -largeValue_) {
1436                                        rowActivityWork_[iRow] = rowLowerWork_[iRow];
1437                                        setRowStatus(iRow, atLowerBound);
1438                                   } else {
1439                                        // say free
1440                                        setRowStatus(iRow, isFree);
1441                                        rowActivityWork_[iRow] = 0.0;
1442                                   }
1443                              } else {
1444                                   if (rowLowerWork_[iRow] > -largeValue_) {
1445                                        // set to nearest
1446                                        if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
1447                                                  < fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])) {
1448                                             rowActivityWork_[iRow] = rowLowerWork_[iRow];
1449                                             setRowStatus(iRow, atLowerBound);
1450                                        } else {
1451                                             rowActivityWork_[iRow] = rowUpperWork_[iRow];
1452                                             setRowStatus(iRow, atUpperBound);
1453                                        }
1454                                   } else {
1455                                        rowActivityWork_[iRow] = rowUpperWork_[iRow];
1456                                        setRowStatus(iRow, atUpperBound);
1457                                   }
1458                              }
1459                              break;
1460                         }
1461                    }
1462                    totalSlacks = numberBasic;
1463
1464                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1465                         switch(getColumnStatus(iColumn)) {
1466
1467                         case basic:
1468                              if (numberBasic == maximumBasic_) {
1469                                   // take out of basis
1470                                   if (columnLowerWork_[iColumn] > -largeValue_) {
1471                                        if (columnActivityWork_[iColumn] - columnLowerWork_[iColumn] <
1472                                                  columnUpperWork_[iColumn] - columnActivityWork_[iColumn]) {
1473                                             columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1474                                             setColumnStatus(iColumn, atLowerBound);
1475                                        } else {
1476                                             columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1477                                             setColumnStatus(iColumn, atUpperBound);
1478                                        }
1479                                   } else if (columnUpperWork_[iColumn] < largeValue_) {
1480                                        columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1481                                        setColumnStatus(iColumn, atUpperBound);
1482                                   } else {
1483                                        columnActivityWork_[iColumn] = 0.0;
1484                                        setColumnStatus(iColumn, isFree);
1485                                   }
1486                              } else {
1487                                   numberBasic++;
1488                              }
1489                              break;
1490                         case atUpperBound:
1491                              columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1492                              if (columnActivityWork_[iColumn] > largeValue_) {
1493                                   if (columnLowerWork_[iColumn] < -largeValue_) {
1494                                        columnActivityWork_[iColumn] = 0.0;
1495                                        setColumnStatus(iColumn, isFree);
1496                                   } else {
1497                                        columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1498                                        setColumnStatus(iColumn, atLowerBound);
1499                                   }
1500                              }
1501                              break;
1502                         case isFixed:
1503                         case atLowerBound:
1504                              columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1505                              if (columnActivityWork_[iColumn] < -largeValue_) {
1506                                   if (columnUpperWork_[iColumn] > largeValue_) {
1507                                        columnActivityWork_[iColumn] = 0.0;
1508                                        setColumnStatus(iColumn, isFree);
1509                                   } else {
1510                                        columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1511                                        setColumnStatus(iColumn, atUpperBound);
1512                                   }
1513                              }
1514                              break;
1515                         case isFree:
1516                              break;
1517                              // not really free - fall through to superbasic
1518                         case superBasic:
1519                              if (columnUpperWork_[iColumn] > largeValue_) {
1520                                   if (columnLowerWork_[iColumn] > -largeValue_) {
1521                                        columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1522                                        setColumnStatus(iColumn, atLowerBound);
1523                                   } else {
1524                                        // say free
1525                                        setColumnStatus(iColumn, isFree);
1526                                        columnActivityWork_[iColumn] = 0.0;
1527                                   }
1528                              } else {
1529                                   if (columnLowerWork_[iColumn] > -largeValue_) {
1530                                        // set to nearest
1531                                        if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1532                                                  < fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])) {
1533                                             columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1534                                             setColumnStatus(iColumn, atLowerBound);
1535                                        } else {
1536                                             columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1537                                             setColumnStatus(iColumn, atUpperBound);
1538                                        }
1539                                   } else {
1540                                        columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1541                                        setColumnStatus(iColumn, atUpperBound);
1542                                   }
1543                              }
1544                              break;
1545                         }
1546                    }
[1449]1547#ifdef CLP_INVESTIGATE3
[1525]1548                    if (saveSol) {
1549                         int numberChanged = 0;
1550                         double largestChanged = 0.0;
1551                         for (int i = 0; i < numberTotal; i++) {
1552                              double difference = fabs(solution_[i] - saveSol[i]);
1553                              if (difference > 1.0e-7) {
1554                                   numberChanged++;
1555                                   if (difference > largestChanged)
1556                                        largestChanged = difference;
1557                              }
1558                         }
1559                         if (numberChanged)
1560                              printf("%d changed, largest %g\n", numberChanged, largestChanged);
1561                         delete [] saveSol;
1562                    }
[1429]1563#endif
[1428]1564#if 0
[1525]1565                    if (numberBasic < numberRows_) {
1566                         // add some slacks in case odd warmstart
[1427]1567#ifdef CLP_INVESTIGATE
[1525]1568                         printf("BAD %d basic, %d rows %d slacks\n",
1569                                numberBasic, numberRows_, totalSlacks);
[1427]1570#endif
[1525]1571                         int iRow = numberRows_ - 1;
1572                         while (numberBasic < numberRows_) {
1573                              if (getRowStatus(iRow) != basic) {
1574                                   setRowStatus(iRow, basic);
1575                                   numberBasic++;
1576                                   totalSlacks++;
1577                                   iRow--;
1578                              } else {
1579                                   break;
1580                              }
1581                         }
1582                    }
[1428]1583#endif
[1525]1584               } else {
1585                    // all slack basis
1586                    int numberBasic = 0;
1587                    if (!status_) {
1588                         createStatus();
1589                    }
1590                    for (iRow = 0; iRow < numberRows_; iRow++) {
1591                         double lower = rowLowerWork_[iRow];
1592                         double upper = rowUpperWork_[iRow];
1593                         if (lower > -largeValue_ || upper < largeValue_) {
1594                              if (fabs(lower) <= fabs(upper)) {
1595                                   rowActivityWork_[iRow] = lower;
1596                              } else {
1597                                   rowActivityWork_[iRow] = upper;
1598                              }
1599                         } else {
1600                              rowActivityWork_[iRow] = 0.0;
1601                         }
1602                         setRowStatus(iRow, basic);
1603                         numberBasic++;
1604                    }
1605                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1606                         double lower = columnLowerWork_[iColumn];
1607                         double upper = columnUpperWork_[iColumn];
1608                         double big_bound = largeValue_;
1609                         if (lower > -big_bound || upper < big_bound) {
1610                              if ((getColumnStatus(iColumn) == atLowerBound &&
1611                                        columnActivityWork_[iColumn] == lower) ||
1612                                        (getColumnStatus(iColumn) == atUpperBound &&
1613                                         columnActivityWork_[iColumn] == upper)) {
1614                                   // status looks plausible
1615                              } else {
1616                                   // set to sensible
1617                                   if (fabs(lower) <= fabs(upper)) {
1618                                        setColumnStatus(iColumn, atLowerBound);
1619                                        columnActivityWork_[iColumn] = lower;
1620                                   } else {
1621                                        setColumnStatus(iColumn, atUpperBound);
1622                                        columnActivityWork_[iColumn] = upper;
1623                                   }
1624                              }
1625                         } else {
1626                              setColumnStatus(iColumn, isFree);
1627                              columnActivityWork_[iColumn] = 0.0;
1628                         }
1629                    }
1630               }
1631          } else {
1632               // values pass has less coding
1633               // make row activities correct and clean basis a bit
1634               cleanStatus();
1635               if (status_) {
1636                    int numberBasic = 0;
1637                    for (iRow = 0; iRow < numberRows_; iRow++) {
1638                         if (getRowStatus(iRow) == basic)
1639                              numberBasic++;
1640                    }
1641                    totalSlacks = numberBasic;
[1321]1642#if 0
[1525]1643                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1644                         if (getColumnStatus(iColumn) == basic)
1645                              numberBasic++;
1646                    }
[1321]1647#endif
[1525]1648               } else {
1649                    // all slack basis
1650                    int numberBasic = 0;
1651                    if (!status_) {
1652                         createStatus();
1653                    }
1654                    for (iRow = 0; iRow < numberRows_; iRow++) {
1655                         setRowStatus(iRow, basic);
1656                         numberBasic++;
1657                    }
1658                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1659                         setColumnStatus(iColumn, superBasic);
1660                         // but put to bound if close
1661                         if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1662                                   <= primalTolerance_) {
1663                              columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1664                              setColumnStatus(iColumn, atLowerBound);
1665                         } else if (fabs(columnActivityWork_[iColumn]
1666                                         - columnUpperWork_[iColumn])
1667                                    <= primalTolerance_) {
1668                              columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1669                              setColumnStatus(iColumn, atUpperBound);
1670                         }
1671                    }
1672               }
1673          }
1674          numberRefinements_ = 1;
1675          // set fixed if they are
1676          for (iRow = 0; iRow < numberRows_; iRow++) {
1677               if (getRowStatus(iRow) != basic ) {
1678                    if (rowLowerWork_[iRow] == rowUpperWork_[iRow]) {
1679                         rowActivityWork_[iRow] = rowLowerWork_[iRow];
1680                         setRowStatus(iRow, isFixed);
1681                    }
1682               }
1683          }
1684          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1685               if (getColumnStatus(iColumn) != basic ) {
1686                    if (columnLowerWork_[iColumn] == columnUpperWork_[iColumn]) {
1687                         columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1688                         setColumnStatus(iColumn, isFixed);
1689                    }
1690               }
1691          }
1692     }
1693     //for (iRow=0;iRow<numberRows_+numberColumns_;iRow++) {
1694     //if (fabs(solution_[iRow])>1.0e10) {
1695     //  printf("large %g at %d - status %d\n",
1696     //         solution_[iRow],iRow,status_[iRow]);
1697     //}
1698     //}
[1613]1699#    ifndef _MSC_VER
1700         // The local static var k is a problem when trying to build a DLL. Since this is
1701         // just for debugging (likely done on *nix), just hide it from Windows
1702         // -- lh, 101016 --
[1525]1703     if (0)  {
1704          static int k = 0;
1705          printf("start basis\n");
1706          int i;
1707          for (i = 0; i < numberRows_; i++)
1708               printf ("xx %d %d\n", i, pivotVariable_[i]);
1709          for (i = 0; i < numberRows_ + numberColumns_; i++)
1710               if (getColumnStatus(i) == basic)
1711                    printf ("yy %d basic\n", i);
1712          if (k > 20)
1713               exit(0);
1714          k++;
1715     }
[1613]1716#    endif
[1525]1717     int status = factorization_->factorize(this, solveType, valuesPass);
1718     if (status) {
1719          handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
1720                    << status
1721                    << CoinMessageEol;
1722          return -1;
1723     } else if (!solveType) {
1724          // Initial basis - return number of singularities
1725          int numberSlacks = 0;
1726          for (iRow = 0; iRow < numberRows_; iRow++) {
1727               if (getRowStatus(iRow) == basic)
1728                    numberSlacks++;
1729          }
1730          status = CoinMax(numberSlacks - totalSlacks, 0);
1731          // special case if all slack
1732          if (numberSlacks == numberRows_) {
1733               status = numberRows_ + 1;
1734          }
1735     }
[2]1736
[1525]1737     // sparse methods
1738     //if (factorization_->sparseThreshold()) {
1739     // get default value
1740     factorization_->sparseThreshold(0);
[1660]1741     if (!(moreSpecialOptions_&1024))
1742       factorization_->goSparse();
[1525]1743     //}
1744
1745     return status;
[2]1746}
[1525]1747/*
[2]1748   This does basis housekeeping and does values for in/out variables.
1749   Can also decide to re-factorize
1750*/
[1525]1751int
[2]1752ClpSimplex::housekeeping(double objectiveChange)
1753{
[1525]1754     // save value of incoming and outgoing
1755     double oldIn = solution_[sequenceIn_];
1756     double oldOut = solution_[sequenceOut_];
1757     numberIterations_++;
1758     changeMade_++; // something has happened
1759     // incoming variable
1760     if (handler_->logLevel() > 7) {
1761          //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
1762          handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
1763                    << directionOut_
1764                    << directionIn_ << theta_
1765                    << dualOut_ << dualIn_ << alpha_
1766                    << CoinMessageEol;
1767          if (getStatus(sequenceIn_) == isFree) {
1768               handler_->message(CLP_SIMPLEX_FREEIN, messages_)
1769                         << sequenceIn_
1770                         << CoinMessageEol;
1771          }
1772     }
[1376]1773#if 0
[1525]1774     printf("h1 %d %d %g %g %g %g",
1775            directionOut_
1776            , directionIn_, theta_
1777            , dualOut_, dualIn_, alpha_);
[1376]1778#endif
[1525]1779     // change of incoming
1780     char rowcol[] = {'R', 'C'};
1781     if (pivotRow_ >= 0)
1782          pivotVariable_[pivotRow_] = sequenceIn();
1783     if (upper_[sequenceIn_] > 1.0e20 && lower_[sequenceIn_] < -1.0e20)
1784          progressFlag_ |= 2; // making real progress
1785     solution_[sequenceIn_] = valueIn_;
1786     if (upper_[sequenceOut_] - lower_[sequenceOut_] < 1.0e-12)
1787          progressFlag_ |= 1; // making real progress
1788     if (sequenceIn_ != sequenceOut_) {
1789          if (alphaAccuracy_ > 0.0) {
1790               double value = fabs(alpha_);
1791               if (value > 1.0)
1792                    alphaAccuracy_ *= value;
1793               else
1794                    alphaAccuracy_ /= value;
1795          }
1796          //assert( getStatus(sequenceOut_)== basic);
1797          setStatus(sequenceIn_, basic);
1798          if (upper_[sequenceOut_] - lower_[sequenceOut_] > 0) {
1799               // As Nonlinear costs may have moved bounds (to more feasible)
1800               // Redo using value
1801               if (fabs(valueOut_ - lower_[sequenceOut_]) < fabs(valueOut_ - upper_[sequenceOut_])) {
1802                    // going to lower
1803                    setStatus(sequenceOut_, atLowerBound);
1804                    oldOut = lower_[sequenceOut_];
1805               } else {
1806                    // going to upper
1807                    setStatus(sequenceOut_, atUpperBound);
1808                    oldOut = upper_[sequenceOut_];
1809               }
1810          } else {
1811               // fixed
1812               setStatus(sequenceOut_, isFixed);
1813          }
1814          solution_[sequenceOut_] = valueOut_;
1815     } else {
1816          //if (objective_->type()<2)
1817          //assert (fabs(theta_)>1.0e-13);
1818          // flip from bound to bound
1819          // As Nonlinear costs may have moved bounds (to more feasible)
1820          // Redo using value
1821          if (fabs(valueIn_ - lower_[sequenceIn_]) < fabs(valueIn_ - upper_[sequenceIn_])) {
1822               // as if from upper bound
1823               setStatus(sequenceIn_, atLowerBound);
1824          } else {
1825               // as if from lower bound
1826               setStatus(sequenceIn_, atUpperBound);
1827          }
1828     }
1829
1830     // Update hidden stuff e.g. effective RHS and gub
[1678]1831     int invertNow=matrix_->updatePivot(this, oldIn, oldOut);
[1525]1832     objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
1833     if (handler_->logLevel() > 7) {
1834          //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
1835          handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
1836                    << numberIterations_ << objectiveValue()
1837                    << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
1838                    << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
1839          handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
1840          handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
1841          handler_->message() << CoinMessageEol;
1842     }
[1376]1843#if 0
[1525]1844     if (numberIterations_ > 10000)
1845          printf(" it %d %g %c%d %c%d\n"
1846                 , numberIterations_, objectiveValue()
1847                 , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
1848                 , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
[1376]1849#endif
[1525]1850     if (trustedUserPointer_ && trustedUserPointer_->typeStruct == 1) {
1851          if (algorithm_ > 0 && integerType_ && !nonLinearCost_->numberInfeasibilities()) {
1852               if (fabs(theta_) > 1.0e-6 || !numberIterations_) {
1853                    // For saving solutions
1854                    typedef struct {
1855                         int numberSolutions;
1856                         int maximumSolutions;
1857                         int numberColumns;
1858                         double ** solution;
1859                         int * numberUnsatisfied;
1860                    } clpSolution;
1861                    clpSolution * solution = reinterpret_cast<clpSolution *> (trustedUserPointer_->data);
1862                    if (solution->numberSolutions == solution->maximumSolutions) {
1863                         int n =  solution->maximumSolutions;
1864                         int n2 = (n * 3) / 2 + 10;
1865                         solution->maximumSolutions = n2;
1866                         double ** temp = new double * [n2];
1867                         for (int i = 0; i < n; i++)
1868                              temp[i] = solution->solution[i];
1869                         delete [] solution->solution;
1870                         solution->solution = temp;
1871                         int * tempN = new int [n2];
1872                         for (int i = 0; i < n; i++)
1873                              tempN[i] = solution->numberUnsatisfied[i];
1874                         delete [] solution->numberUnsatisfied;
1875                         solution->numberUnsatisfied = tempN;
1876                    }
1877                    assert (numberColumns_ == solution->numberColumns);
1878                    double * sol = new double [numberColumns_];
1879                    solution->solution[solution->numberSolutions] = sol;
1880                    int numberFixed = 0;
1881                    int numberUnsat = 0;
1882                    int numberSat = 0;
1883                    double sumUnsat = 0.0;
1884                    double tolerance = 10.0 * primalTolerance_;
1885                    double mostAway = 0.0;
1886                    int iAway = -1;
1887                    for (int i = 0; i < numberColumns_; i++) {
1888                         // Save anyway
1889                         sol[i] = columnScale_ ? solution_[i] * columnScale_[i] : solution_[i];
1890                         // rest is optional
1891                         if (upper_[i] > lower_[i]) {
1892                              double value = solution_[i];
1893                              if (value > lower_[i] + tolerance &&
1894                                        value < upper_[i] - tolerance && integerType_[i]) {
1895                                   // may have to modify value if scaled
1896                                   if (columnScale_)
1897                                        value *= columnScale_[i];
1898                                   double closest = floor(value + 0.5);
1899                                   // problem may be perturbed so relax test
1900                                   if (fabs(value - closest) > 1.0e-4) {
1901                                        numberUnsat++;
1902                                        sumUnsat += fabs(value - closest);
1903                                        if (mostAway < fabs(value - closest)) {
1904                                             mostAway = fabs(value - closest);
1905                                             iAway = i;
1906                                        }
1907                                   } else {
1908                                        numberSat++;
1909                                   }
1910                              } else {
1911                                   numberSat++;
1912                              }
1913                         } else {
1914                              numberFixed++;
1915                         }
1916                    }
1917                    solution->numberUnsatisfied[solution->numberSolutions++] = numberUnsat;
1918                    printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n",
1919                           numberIterations_, numberUnsat, sumUnsat, mostAway, numberFixed, numberSat);
1920               }
1921          }
1922     }
1923     if (hitMaximumIterations())
1924          return 2;
[266]1925#if 1
[1525]1926     //if (numberIterations_>14000)
1927     //handler_->setLogLevel(63);
1928     //if (numberIterations_>24000)
1929     //exit(77);
1930     // check for small cycles
1931     int in = sequenceIn_;
1932     int out = sequenceOut_;
1933     matrix_->correctSequence(this, in, out);
1934     int cycle = progress_.cycle(in, out,
1935                                 directionIn_, directionOut_);
[1676]1936     if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
[1525]1937          //if (cycle>0) {
1938          if (handler_->logLevel() >= 63)
1939               printf("Cycle of %d\n", cycle);
1940          // reset
1941          progress_.startCheck();
1942          double random = randomNumberGenerator_.randomDouble();
1943          int extra = static_cast<int> (9.999 * random);
1944          int off[] = {1, 1, 1, 1, 2, 2, 2, 3, 3, 4};
1945          if (factorization_->pivots() > cycle) {
1946               forceFactorization_ = CoinMax(1, cycle - off[extra]);
1947          } else {
[1676]1948            /* need to reject something
1949               should be better if don't reject incoming
1950               as it is in basis */
[1525]1951               int iSequence;
[1676]1952               //if (algorithm_ > 0)
1953               //   iSequence = sequenceIn_;
1954               //else
[1525]1955                    iSequence = sequenceOut_;
1956               char x = isColumn(iSequence) ? 'C' : 'R';
1957               if (handler_->logLevel() >= 63)
1958                    handler_->message(CLP_SIMPLEX_FLAG, messages_)
1959                              << x << sequenceWithin(iSequence)
1960                              << CoinMessageEol;
1961               setFlagged(iSequence);
1962               //printf("flagging %d\n",iSequence);
1963          }
1964          return 1;
1965     }
[225]1966#endif
[1525]1967     // only time to re-factorize if one before real time
1968     // this is so user won't be surprised that maximumPivots has exact meaning
1969     int numberPivots = factorization_->pivots();
1970     int maximumPivots = factorization_->maximumPivots();
1971     int numberDense = factorization_->numberDense();
1972     bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 >
1973                        2 * maximumIterations());
1974     if (numberPivots == maximumPivots ||
1975               maximumPivots < 2) {
1976          // If dense then increase
1977          if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots) {
1978               factorization_->maximumPivots(numberDense);
1979               dualRowPivot_->maximumPivotsChanged();
1980               primalColumnPivot_->maximumPivotsChanged();
1981               // and redo arrays
1982               for (int iRow = 0; iRow < 4; iRow++) {
1983                    int length = rowArray_[iRow]->capacity() + numberDense - maximumPivots;
1984                    rowArray_[iRow]->reserve(length);
1985               }
1986          }
1987          return 1;
[1678]1988     } else if ((factorization_->timeToRefactorize() && !dontInvert)
1989                ||invertNow) {
[1525]1990          //printf("ret after %d pivots\n",factorization_->pivots());
1991          return 1;
1992     } else if (forceFactorization_ > 0 &&
1993                factorization_->pivots() == forceFactorization_) {
1994          // relax
1995          forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
1996          if (forceFactorization_ > factorization_->maximumPivots())
1997               forceFactorization_ = -1; //off
1998          return 1;
[1676]1999     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
[1525]2000          double random = randomNumberGenerator_.randomDouble();
2001          int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
2002          if (factorization_->pivots() >= random * maxNumber) {
2003               return 1;
2004          } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) &&
2005                     numberIterations_ < 1001000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
2006               return 1;
2007          } else {
2008               // carry on iterating
2009               return 0;
2010          }
2011     } else {
2012          // carry on iterating
2013          return 0;
2014     }
[2]2015}
[1525]2016// Copy constructor.
2017ClpSimplex::ClpSimplex(const ClpSimplex &rhs, int scalingMode) :
2018     ClpModel(rhs, scalingMode),
2019     bestPossibleImprovement_(0.0),
2020     zeroTolerance_(1.0e-13),
2021     columnPrimalSequence_(-2),
2022     rowPrimalSequence_(-2),
2023     bestObjectiveValue_(rhs.bestObjectiveValue_),
2024     moreSpecialOptions_(2),
2025     baseIteration_(0),
2026     primalToleranceToGetOptimal_(-1.0),
2027     largeValue_(1.0e15),
2028     largestPrimalError_(0.0),
2029     largestDualError_(0.0),
2030     alphaAccuracy_(-1.0),
2031     dualBound_(1.0e10),
2032     alpha_(0.0),
2033     theta_(0.0),
2034     lowerIn_(0.0),
2035     valueIn_(0.0),
2036     upperIn_(-COIN_DBL_MAX),
2037     dualIn_(0.0),
2038     lowerOut_(-1),
2039     valueOut_(-1),
2040     upperOut_(-1),
2041     dualOut_(-1),
2042     dualTolerance_(1.0e-7),
2043     primalTolerance_(1.0e-7),
2044     sumDualInfeasibilities_(0.0),
2045     sumPrimalInfeasibilities_(0.0),
2046     infeasibilityCost_(1.0e10),
2047     sumOfRelaxedDualInfeasibilities_(0.0),
2048     sumOfRelaxedPrimalInfeasibilities_(0.0),
2049     acceptablePivot_(1.0e-8),
2050     lower_(NULL),
2051     rowLowerWork_(NULL),
2052     columnLowerWork_(NULL),
2053     upper_(NULL),
2054     rowUpperWork_(NULL),
2055     columnUpperWork_(NULL),
2056     cost_(NULL),
2057     rowObjectiveWork_(NULL),
2058     objectiveWork_(NULL),
2059     sequenceIn_(-1),
2060     directionIn_(-1),
2061     sequenceOut_(-1),
2062     directionOut_(-1),
2063     pivotRow_(-1),
2064     lastGoodIteration_(-100),
2065     dj_(NULL),
2066     rowReducedCost_(NULL),
2067     reducedCostWork_(NULL),
2068     solution_(NULL),
2069     rowActivityWork_(NULL),
2070     columnActivityWork_(NULL),
2071     numberDualInfeasibilities_(0),
2072     numberDualInfeasibilitiesWithoutFree_(0),
2073     numberPrimalInfeasibilities_(100),
2074     numberRefinements_(0),
2075     pivotVariable_(NULL),
2076     factorization_(NULL),
2077     savedSolution_(NULL),
2078     numberTimesOptimal_(0),
2079     disasterArea_(NULL),
2080     changeMade_(1),
2081     algorithm_(0),
2082     forceFactorization_(-1),
2083     perturbation_(100),
2084     nonLinearCost_(NULL),
2085     lastBadIteration_(-999999),
2086     lastFlaggedIteration_(-999999),
2087     numberFake_(0),
2088     numberChanged_(0),
2089     progressFlag_(0),
2090     firstFree_(-1),
2091     numberExtraRows_(0),
2092     maximumBasic_(0),
2093     dontFactorizePivots_(0),
2094     incomingInfeasibility_(1.0),
2095     allowedInfeasibility_(10.0),
2096     automaticScale_(0),
2097     maximumPerturbationSize_(0),
2098     perturbationArray_(NULL),
2099     baseModel_(NULL)
2100{
2101     int i;
2102     for (i = 0; i < 6; i++) {
2103          rowArray_[i] = NULL;
2104          columnArray_[i] = NULL;
2105     }
2106     for (i = 0; i < 4; i++) {
2107          spareIntArray_[i] = 0;
2108          spareDoubleArray_[i] = 0.0;
2109     }
2110     saveStatus_ = NULL;
2111     factorization_ = NULL;
2112     dualRowPivot_ = NULL;
2113     primalColumnPivot_ = NULL;
2114     gutsOfDelete(0);
2115     delete nonLinearCost_;
2116     nonLinearCost_ = NULL;
2117     gutsOfCopy(rhs);
2118     solveType_ = 1; // say simplex based life form
[2]2119}
2120// Copy constructor from model
[393]2121ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
[1525]2122     ClpModel(rhs, scalingMode),
2123     bestPossibleImprovement_(0.0),
2124     zeroTolerance_(1.0e-13),
2125     columnPrimalSequence_(-2),
2126     rowPrimalSequence_(-2),
2127     bestObjectiveValue_(-COIN_DBL_MAX),
2128     moreSpecialOptions_(2),
2129     baseIteration_(0),
2130     primalToleranceToGetOptimal_(-1.0),
2131     largeValue_(1.0e15),
2132     largestPrimalError_(0.0),
2133     largestDualError_(0.0),
2134     alphaAccuracy_(-1.0),
2135     dualBound_(1.0e10),
2136     alpha_(0.0),
2137     theta_(0.0),
2138     lowerIn_(0.0),
2139     valueIn_(0.0),
2140     upperIn_(-COIN_DBL_MAX),
2141     dualIn_(0.0),
2142     lowerOut_(-1),
2143     valueOut_(-1),
2144     upperOut_(-1),
2145     dualOut_(-1),
2146     dualTolerance_(1.0e-7),
2147     primalTolerance_(1.0e-7),
2148     sumDualInfeasibilities_(0.0),
2149     sumPrimalInfeasibilities_(0.0),
2150     infeasibilityCost_(1.0e10),
2151     sumOfRelaxedDualInfeasibilities_(0.0),
2152     sumOfRelaxedPrimalInfeasibilities_(0.0),
2153     acceptablePivot_(1.0e-8),
2154     lower_(NULL),
2155     rowLowerWork_(NULL),
2156     columnLowerWork_(NULL),
2157     upper_(NULL),
2158     rowUpperWork_(NULL),
2159     columnUpperWork_(NULL),
2160     cost_(NULL),
2161     rowObjectiveWork_(NULL),
2162     objectiveWork_(NULL),
2163     sequenceIn_(-1),
2164     directionIn_(-1),
2165     sequenceOut_(-1),
2166     directionOut_(-1),
2167     pivotRow_(-1),
2168     lastGoodIteration_(-100),
2169     dj_(NULL),
2170     rowReducedCost_(NULL),
2171     reducedCostWork_(NULL),
2172     solution_(NULL),
2173     rowActivityWork_(NULL),
2174     columnActivityWork_(NULL),
2175     numberDualInfeasibilities_(0),
2176     numberDualInfeasibilitiesWithoutFree_(0),
2177     numberPrimalInfeasibilities_(100),
2178     numberRefinements_(0),
2179     pivotVariable_(NULL),
2180     factorization_(NULL),
2181     savedSolution_(NULL),
2182     numberTimesOptimal_(0),
2183     disasterArea_(NULL),
2184     changeMade_(1),
2185     algorithm_(0),
2186     forceFactorization_(-1),
2187     perturbation_(100),
2188     nonLinearCost_(NULL),
2189     lastBadIteration_(-999999),
2190     lastFlaggedIteration_(-999999),
2191     numberFake_(0),
2192     numberChanged_(0),
2193     progressFlag_(0),
2194     firstFree_(-1),
2195     numberExtraRows_(0),
2196     maximumBasic_(0),
2197     dontFactorizePivots_(0),
2198     incomingInfeasibility_(1.0),
2199     allowedInfeasibility_(10.0),
2200     automaticScale_(0),
2201     maximumPerturbationSize_(0),
2202     perturbationArray_(NULL),
2203     baseModel_(NULL)
[2]2204{
[1525]2205     int i;
2206     for (i = 0; i < 6; i++) {
2207          rowArray_[i] = NULL;
2208          columnArray_[i] = NULL;
2209     }
2210     for (i = 0; i < 4; i++) {
2211          spareIntArray_[i] = 0;
2212          spareDoubleArray_[i] = 0.0;
2213     }
2214     saveStatus_ = NULL;
2215     // get an empty factorization so we can set tolerances etc
2216     getEmptyFactorization();
2217     // say Steepest pricing
2218     dualRowPivot_ = new ClpDualRowSteepest();
2219     // say Steepest pricing
2220     primalColumnPivot_ = new ClpPrimalColumnSteepest();
2221     solveType_ = 1; // say simplex based life form
2222
[2]2223}
2224// Assignment operator. This copies the data
[1525]2225ClpSimplex &
[2]2226ClpSimplex::operator=(const ClpSimplex & rhs)
2227{
[1525]2228     if (this != &rhs) {
2229          gutsOfDelete(0);
2230          delete nonLinearCost_;
2231          nonLinearCost_ = NULL;
2232          ClpModel::operator=(rhs);
2233          gutsOfCopy(rhs);
2234     }
2235     return *this;
[2]2236}
[1525]2237void
[2]2238ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
2239{
[1525]2240     assert (numberRows_ == rhs.numberRows_);
2241     assert (numberColumns_ == rhs.numberColumns_);
2242     numberExtraRows_ = rhs.numberExtraRows_;
2243     maximumBasic_ = rhs.maximumBasic_;
2244     dontFactorizePivots_ = rhs.dontFactorizePivots_;
2245     int numberRows2 = numberRows_ + numberExtraRows_;
2246     moreSpecialOptions_ = rhs.moreSpecialOptions_;
2247     if ((whatsChanged_ & 1) != 0) {
2248          int numberTotal = numberColumns_ + numberRows2;
2249          if ((specialOptions_ & 65536) != 0 && maximumRows_ >= 0) {
2250               assert (maximumInternalRows_ >= numberRows2);
2251               assert (maximumInternalColumns_ >= numberColumns_);
2252               numberTotal = 2 * (maximumInternalColumns_ + maximumInternalRows_);
2253          }
2254          lower_ = ClpCopyOfArray(rhs.lower_, numberTotal);
2255          rowLowerWork_ = lower_ + numberColumns_;
2256          columnLowerWork_ = lower_;
2257          upper_ = ClpCopyOfArray(rhs.upper_, numberTotal);
2258          rowUpperWork_ = upper_ + numberColumns_;
2259          columnUpperWork_ = upper_;
2260          cost_ = ClpCopyOfArray(rhs.cost_, numberTotal);
2261          objectiveWork_ = cost_;
2262          rowObjectiveWork_ = cost_ + numberColumns_;
2263          dj_ = ClpCopyOfArray(rhs.dj_, numberTotal);
2264          if (dj_) {
2265               reducedCostWork_ = dj_;
2266               rowReducedCost_ = dj_ + numberColumns_;
2267          }
2268          solution_ = ClpCopyOfArray(rhs.solution_, numberTotal);
2269          if (solution_) {
2270               columnActivityWork_ = solution_;
2271               rowActivityWork_ = solution_ + numberColumns_;
2272          }
2273          if (rhs.pivotVariable_) {
2274               pivotVariable_ = new int[numberRows2];
2275               CoinMemcpyN ( rhs.pivotVariable_, numberRows2 , pivotVariable_);
2276          } else {
2277               pivotVariable_ = NULL;
2278          }
2279          savedSolution_ = ClpCopyOfArray(rhs.savedSolution_, numberTotal);
2280          int i;
2281          for (i = 0; i < 6; i++) {
2282               rowArray_[i] = NULL;
2283               if (rhs.rowArray_[i])
2284                    rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
2285               columnArray_[i] = NULL;
2286               if (rhs.columnArray_[i])
2287                    columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
2288          }
2289          if (rhs.saveStatus_) {
2290               saveStatus_ = ClpCopyOfArray( rhs.saveStatus_, numberTotal);
2291          }
2292     } else {
2293          lower_ = NULL;
2294          rowLowerWork_ = NULL;
2295          columnLowerWork_ = NULL;
2296          upper_ = NULL;
2297          rowUpperWork_ = NULL;
2298          columnUpperWork_ = NULL;
2299          cost_ = NULL;
2300          objectiveWork_ = NULL;
2301          rowObjectiveWork_ = NULL;
2302          dj_ = NULL;
2303          reducedCostWork_ = NULL;
2304          rowReducedCost_ = NULL;
2305          solution_ = NULL;
2306          columnActivityWork_ = NULL;
2307          rowActivityWork_ = NULL;
2308          pivotVariable_ = NULL;
2309          savedSolution_ = NULL;
2310          int i;
2311          for (i = 0; i < 6; i++) {
2312               rowArray_[i] = NULL;
2313               columnArray_[i] = NULL;
2314          }
2315          saveStatus_ = NULL;
2316     }
2317     if (rhs.factorization_) {
2318          setFactorization(*rhs.factorization_);
2319     } else {
2320          delete factorization_;
2321          factorization_ = NULL;
2322     }
2323     bestPossibleImprovement_ = rhs.bestPossibleImprovement_;
2324     columnPrimalSequence_ = rhs.columnPrimalSequence_;
2325     zeroTolerance_ = rhs.zeroTolerance_;
2326     rowPrimalSequence_ = rhs.rowPrimalSequence_;
2327     bestObjectiveValue_ = rhs.bestObjectiveValue_;
2328     baseIteration_ = rhs.baseIteration_;
2329     primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
2330     largeValue_ = rhs.largeValue_;
2331     largestPrimalError_ = rhs.largestPrimalError_;
2332     largestDualError_ = rhs.largestDualError_;
2333     alphaAccuracy_ = rhs.alphaAccuracy_;
2334     dualBound_ = rhs.dualBound_;
2335     alpha_ = rhs.alpha_;
2336     theta_ = rhs.theta_;
2337     lowerIn_ = rhs.lowerIn_;
2338     valueIn_ = rhs.valueIn_;
2339     upperIn_ = rhs.upperIn_;
2340     dualIn_ = rhs.dualIn_;
2341     sequenceIn_ = rhs.sequenceIn_;
2342     directionIn_ = rhs.directionIn_;
2343     lowerOut_ = rhs.lowerOut_;
2344     valueOut_ = rhs.valueOut_;
2345     upperOut_ = rhs.upperOut_;
2346     dualOut_ = rhs.dualOut_;
2347     sequenceOut_ = rhs.sequenceOut_;
2348     directionOut_ = rhs.directionOut_;
2349     pivotRow_ = rhs.pivotRow_;
2350     lastGoodIteration_ = rhs.lastGoodIteration_;
2351     numberRefinements_ = rhs.numberRefinements_;
2352     dualTolerance_ = rhs.dualTolerance_;
2353     primalTolerance_ = rhs.primalTolerance_;
2354     sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
2355     numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
2356     numberDualInfeasibilitiesWithoutFree_ =
2357          rhs.numberDualInfeasibilitiesWithoutFree_;
2358     sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
2359     numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
2360     dualRowPivot_ = rhs.dualRowPivot_->clone(true);
2361     dualRowPivot_->setModel(this);
2362     primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
2363     primalColumnPivot_->setModel(this);
2364     numberTimesOptimal_ = rhs.numberTimesOptimal_;
2365     disasterArea_ = NULL;
2366     changeMade_ = rhs.changeMade_;
2367     algorithm_ = rhs.algorithm_;
2368     forceFactorization_ = rhs.forceFactorization_;
2369     perturbation_ = rhs.perturbation_;
2370     infeasibilityCost_ = rhs.infeasibilityCost_;
2371     lastBadIteration_ = rhs.lastBadIteration_;
2372     lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
2373     numberFake_ = rhs.numberFake_;
2374     numberChanged_ = rhs.numberChanged_;
2375     progressFlag_ = rhs.progressFlag_;
2376     firstFree_ = rhs.firstFree_;
2377     incomingInfeasibility_ = rhs.incomingInfeasibility_;
2378     allowedInfeasibility_ = rhs.allowedInfeasibility_;
2379     automaticScale_ = rhs.automaticScale_;
2380     maximumPerturbationSize_ = rhs.maximumPerturbationSize_;
2381     if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) {
2382          perturbationArray_ = CoinCopyOfArray(rhs.perturbationArray_,
2383                                               maximumPerturbationSize_);
2384     } else {
2385          maximumPerturbationSize_ = 0;
2386          perturbationArray_ = NULL;
2387     }
2388     if (rhs.baseModel_) {
2389          baseModel_ = new ClpSimplex(*rhs.baseModel_);
2390     } else {
2391          baseModel_ = NULL;
2392     }
2393     progress_ = rhs.progress_;
2394     for (int i = 0; i < 4; i++) {
2395          spareIntArray_[i] = rhs.spareIntArray_[i];
2396          spareDoubleArray_[i] = rhs.spareDoubleArray_[i];
2397     }
2398     sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
2399     sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
2400     acceptablePivot_ = rhs.acceptablePivot_;
2401     if (rhs.nonLinearCost_ != NULL)
2402          nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
2403     else
2404          nonLinearCost_ = NULL;
2405     solveType_ = rhs.solveType_;
[2]2406}
[50]2407// type == 0 do everything, most + pivot data, 2 factorization data as well
[1525]2408void
[2]2409ClpSimplex::gutsOfDelete(int type)
2410{
[1525]2411     if (!type || (specialOptions_ & 65536) == 0) {
2412          maximumInternalColumns_ = -1;
2413          maximumInternalRows_ = -1;
2414          delete [] lower_;
2415          lower_ = NULL;
2416          rowLowerWork_ = NULL;
2417          columnLowerWork_ = NULL;
2418          delete [] upper_;
2419          upper_ = NULL;
2420          rowUpperWork_ = NULL;
2421          columnUpperWork_ = NULL;
2422          delete [] cost_;
2423          cost_ = NULL;
2424          objectiveWork_ = NULL;
2425          rowObjectiveWork_ = NULL;
2426          delete [] dj_;
2427          dj_ = NULL;
2428          reducedCostWork_ = NULL;
2429          rowReducedCost_ = NULL;
2430          delete [] solution_;
2431          solution_ = NULL;
2432          rowActivityWork_ = NULL;
2433          columnActivityWork_ = NULL;
2434          delete [] savedSolution_;
2435          savedSolution_ = NULL;
2436     }
2437     if ((specialOptions_ & 2) == 0) {
2438          delete nonLinearCost_;
2439          nonLinearCost_ = NULL;
2440     }
2441     int i;
2442     if ((specialOptions_ & 65536) == 0) {
2443          for (i = 0; i < 6; i++) {
2444               delete rowArray_[i];
2445               rowArray_[i] = NULL;
2446               delete columnArray_[i];
2447               columnArray_[i] = NULL;
2448          }
2449     }
2450     delete [] saveStatus_;
2451     saveStatus_ = NULL;
2452     if (type != 1) {
2453          delete rowCopy_;
2454          rowCopy_ = NULL;
2455     }
2456     if (!type) {
2457          // delete everything
2458          setEmptyFactorization();
2459          delete [] pivotVariable_;
2460          pivotVariable_ = NULL;
2461          delete dualRowPivot_;
2462          dualRowPivot_ = NULL;
2463          delete primalColumnPivot_;
2464          primalColumnPivot_ = NULL;
2465          delete baseModel_;
2466          baseModel_ = NULL;
2467          delete [] perturbationArray_;
2468          perturbationArray_ = NULL;
2469          maximumPerturbationSize_ = 0;
2470     } else {
2471          // delete any size information in methods
2472          if (type > 1) {
2473               //assert (factorization_);
2474               if (factorization_)
2475                    factorization_->clearArrays();
2476               delete [] pivotVariable_;
2477               pivotVariable_ = NULL;
2478          }
2479          dualRowPivot_->clearArrays();
2480          primalColumnPivot_->clearArrays();
2481     }
[2]2482}
2483// This sets largest infeasibility and most infeasible
[1525]2484void
[2]2485ClpSimplex::checkPrimalSolution(const double * rowActivities,
[1525]2486                                const double * columnActivities)
[2]2487{
[1525]2488     double * solution;
2489     int iRow, iColumn;
[2]2490
[1525]2491     objectiveValue_ = 0.0;
2492     // now look at primal solution
2493     solution = rowActivityWork_;
2494     sumPrimalInfeasibilities_ = 0.0;
2495     numberPrimalInfeasibilities_ = 0;
2496     double primalTolerance = primalTolerance_;
2497     double relaxedTolerance = primalTolerance_;
2498     // we can't really trust infeasibilities if there is primal error
2499     double error = CoinMin(1.0e-2, largestPrimalError_);
2500     // allow tolerance at least slightly bigger than standard
2501     relaxedTolerance = relaxedTolerance +  error;
2502     sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2503     for (iRow = 0; iRow < numberRows_; iRow++) {
2504          //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
2505          double infeasibility = 0.0;
2506          objectiveValue_ += solution[iRow] * rowObjectiveWork_[iRow];
2507          if (solution[iRow] > rowUpperWork_[iRow]) {
2508               infeasibility = solution[iRow] - rowUpperWork_[iRow];
2509          } else if (solution[iRow] < rowLowerWork_[iRow]) {
2510               infeasibility = rowLowerWork_[iRow] - solution[iRow];
2511          }
2512          if (infeasibility > primalTolerance) {
2513               sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2514               if (infeasibility > relaxedTolerance)
2515                    sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2516               numberPrimalInfeasibilities_ ++;
2517          }
2518          infeasibility = fabs(rowActivities[iRow] - solution[iRow]);
2519     }
2520     // Check any infeasibilities from dynamic rows
2521     matrix_->primalExpanded(this, 2);
2522     solution = columnActivityWork_;
2523     if (!matrix_->rhsOffset(this)) {
2524          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2525               //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2526               double infeasibility = 0.0;
2527               objectiveValue_ += objectiveWork_[iColumn] * solution[iColumn];
2528               if (solution[iColumn] > columnUpperWork_[iColumn]) {
2529                    infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2530               } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2531                    infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2532               }
2533               if (infeasibility > primalTolerance) {
2534                    sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2535                    if (infeasibility > relaxedTolerance)
2536                         sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2537                    numberPrimalInfeasibilities_ ++;
2538               }
2539               infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2540          }
2541     } else {
2542          // as we are using effective rhs we only check basics
2543          // But we do need to get objective
2544          objectiveValue_ += innerProduct(objectiveWork_, numberColumns_, solution);
2545          for (int j = 0; j < numberRows_; j++) {
2546               int iColumn = pivotVariable_[j];
2547               //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2548               double infeasibility = 0.0;
2549               if (solution[iColumn] > columnUpperWork_[iColumn]) {
2550                    infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2551               } else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2552                    infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2553               }
2554               if (infeasibility > primalTolerance) {
2555                    sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2556                    if (infeasibility > relaxedTolerance)
2557                         sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2558                    numberPrimalInfeasibilities_ ++;
2559               }
2560               infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2561          }
2562     }
2563     objectiveValue_ += objective_->nonlinearOffset();
2564     objectiveValue_ /= (objectiveScale_ * rhsScale_);
[2]2565}
[1525]2566void
[2]2567ClpSimplex::checkDualSolution()
2568{
2569
[1525]2570     int iRow, iColumn;
2571     sumDualInfeasibilities_ = 0.0;
2572     numberDualInfeasibilities_ = 0;
2573     numberDualInfeasibilitiesWithoutFree_ = 0;
2574     if (matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) {
2575          // pretend we found dual infeasibilities
2576          sumOfRelaxedDualInfeasibilities_ = 1.0;
2577          sumDualInfeasibilities_ = 1.0;
2578          numberDualInfeasibilities_ = 1;
2579          return;
2580     }
2581     int firstFreePrimal = -1;
2582     int firstFreeDual = -1;
2583     int numberSuperBasicWithDj = 0;
2584     bestPossibleImprovement_ = 0.0;
2585     // we can't really trust infeasibilities if there is dual error
2586     double error = CoinMin(1.0e-2, largestDualError_);
2587     // allow tolerance at least slightly bigger than standard
2588     double relaxedTolerance = dualTolerance_ +  error;
2589     // allow bigger tolerance for possible improvement
2590     double possTolerance = 5.0 * relaxedTolerance;
2591     sumOfRelaxedDualInfeasibilities_ = 0.0;
[50]2592
[1525]2593     // Check any djs from dynamic rows
2594     matrix_->dualExpanded(this, NULL, NULL, 3);
2595     numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_;
2596     objectiveValue_ = 0.0;
2597     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2598          objectiveValue_ += objectiveWork_[iColumn] * columnActivityWork_[iColumn];
2599          if (getColumnStatus(iColumn) != basic && !flagged(iColumn)) {
2600               // not basic
2601               double distanceUp = columnUpperWork_[iColumn] -
2602                                   columnActivityWork_[iColumn];
2603               double distanceDown = columnActivityWork_[iColumn] -
2604                                     columnLowerWork_[iColumn];
2605               if (distanceUp > primalTolerance_) {
2606                    double value = reducedCostWork_[iColumn];
2607                    // Check if "free"
2608                    if (distanceDown > primalTolerance_) {
2609                         if (fabs(value) > 1.0e2 * relaxedTolerance) {
2610                              numberSuperBasicWithDj++;
2611                              if (firstFreeDual < 0)
2612                                   firstFreeDual = iColumn;
2613                         }
2614                         if (firstFreePrimal < 0)
2615                              firstFreePrimal = iColumn;
2616                    }
2617                    // should not be negative
2618                    if (value < 0.0) {
2619                         value = - value;
2620                         if (value > dualTolerance_) {
2621                              if (getColumnStatus(iColumn) != isFree) {
2622                                   numberDualInfeasibilitiesWithoutFree_ ++;
2623                                   sumDualInfeasibilities_ += value - dualTolerance_;
2624                                   if (value > possTolerance)
2625                                        bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2626                                   if (value > relaxedTolerance)
2627                                        sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2628                                   numberDualInfeasibilities_ ++;
2629                              } else {
2630                                   // free so relax a lot
2631                                   value *= 0.01;
2632                                   if (value > dualTolerance_) {
2633                                        sumDualInfeasibilities_ += value - dualTolerance_;
2634                                        if (value > possTolerance)
2635                                             bestPossibleImprovement_ = 1.0e100;
2636                                        if (value > relaxedTolerance)
2637                                             sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2638                                        numberDualInfeasibilities_ ++;
2639                                   }
2640                              }
2641                         }
2642                    }
2643               }
2644               if (distanceDown > primalTolerance_) {
2645                    double value = reducedCostWork_[iColumn];
2646                    // should not be positive
2647                    if (value > 0.0) {
2648                         if (value > dualTolerance_) {
2649                              sumDualInfeasibilities_ += value - dualTolerance_;
2650                              if (value > possTolerance)
2651                                   bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2652                              if (value > relaxedTolerance)
2653                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2654                              numberDualInfeasibilities_ ++;
2655                              if (getColumnStatus(iColumn) != isFree)
2656                                   numberDualInfeasibilitiesWithoutFree_ ++;
2657                              // maybe we can make feasible by increasing tolerance
2658                         }
2659                    }
2660               }
2661          }
2662     }
2663     for (iRow = 0; iRow < numberRows_; iRow++) {
2664          objectiveValue_ += rowActivityWork_[iRow] * rowObjectiveWork_[iRow];
2665          if (getRowStatus(iRow) != basic && !flagged(iRow + numberColumns_)) {
2666               // not basic
2667               double distanceUp = rowUpperWork_[iRow] - rowActivityWork_[iRow];
2668               double distanceDown = rowActivityWork_[iRow] - rowLowerWork_[iRow];
2669               if (distanceUp > primalTolerance_) {
2670                    double value = rowReducedCost_[iRow];
2671                    // Check if "free"
2672                    if (distanceDown > primalTolerance_) {
2673                         if (fabs(value) > 1.0e2 * relaxedTolerance) {
2674                              numberSuperBasicWithDj++;
2675                              if (firstFreeDual < 0)
2676                                   firstFreeDual = iRow + numberColumns_;
2677                         }
2678                         if (firstFreePrimal < 0)
2679                              firstFreePrimal = iRow + numberColumns_;
2680                    }
2681                    // should not be negative
2682                    if (value < 0.0) {
2683                         value = - value;
2684                         if (value > dualTolerance_) {
2685                              sumDualInfeasibilities_ += value - dualTolerance_;
2686                              if (value > possTolerance)
2687                                   bestPossibleImprovement_ += value * CoinMin(distanceUp, 1.0e10);
2688                              if (value > relaxedTolerance)
2689                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2690                              numberDualInfeasibilities_ ++;
2691                              if (getRowStatus(iRow) != isFree)
2692                                   numberDualInfeasibilitiesWithoutFree_ ++;
2693                         }
2694                    }
2695               }
2696               if (distanceDown > primalTolerance_) {
2697                    double value = rowReducedCost_[iRow];
2698                    // should not be positive
2699                    if (value > 0.0) {
2700                         if (value > dualTolerance_) {
2701                              sumDualInfeasibilities_ += value - dualTolerance_;
2702                              if (value > possTolerance)
2703                                   bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2704                              if (value > relaxedTolerance)
2705                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2706                              numberDualInfeasibilities_ ++;
2707                              if (getRowStatus(iRow) != isFree)
2708                                   numberDualInfeasibilitiesWithoutFree_ ++;
2709                              // maybe we can make feasible by increasing tolerance
2710                         }
2711                    }
2712               }
2713          }
2714     }
2715     if (algorithm_ < 0 && firstFreeDual >= 0) {
2716          // dual
2717          firstFree_ = firstFreeDual;
2718     } else if (numberSuperBasicWithDj ||
2719                (progress_.lastIterationNumber(0) <= 0)) {
2720          firstFree_ = firstFreePrimal;
2721     }
2722     objectiveValue_ += objective_->nonlinearOffset();
2723     objectiveValue_ /= (objectiveScale_ * rhsScale_);
[2]2724}
[568]2725/* This sets sum and number of infeasibilities (Dual and Primal) */
[1525]2726void
[568]2727ClpSimplex::checkBothSolutions()
2728{
[1525]2729     if ((matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) ||
2730               matrix_->rhsOffset(this)) {
2731          // Say may be free or superbasic
2732          moreSpecialOptions_ &= ~8;
2733          // old way
2734          checkPrimalSolution(rowActivityWork_, columnActivityWork_);
2735          checkDualSolution();
2736          return;
2737     }
2738     int iSequence;
2739     assert (dualTolerance_ > 0.0 && dualTolerance_ < 1.0e10);
2740     assert (primalTolerance_ > 0.0 && primalTolerance_ < 1.0e10);
2741     objectiveValue_ = 0.0;
2742     sumPrimalInfeasibilities_ = 0.0;
2743     numberPrimalInfeasibilities_ = 0;
2744     double primalTolerance = primalTolerance_;
2745     double relaxedToleranceP = primalTolerance_;
2746     // we can't really trust infeasibilities if there is primal error
2747     double error = CoinMin(1.0e-2, largestPrimalError_);
2748     // allow tolerance at least slightly bigger than standard
2749     relaxedToleranceP = relaxedToleranceP +  error;
2750     sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2751     sumDualInfeasibilities_ = 0.0;
2752     numberDualInfeasibilities_ = 0;
2753     double dualTolerance = dualTolerance_;
2754     double relaxedToleranceD = dualTolerance;
2755     // we can't really trust infeasibilities if there is dual error
2756     error = CoinMin(1.0e-2, largestDualError_);
2757     // allow tolerance at least slightly bigger than standard
2758     relaxedToleranceD = relaxedToleranceD +  error;
2759     // allow bigger tolerance for possible improvement
2760     double possTolerance = 5.0 * relaxedToleranceD;
2761     sumOfRelaxedDualInfeasibilities_ = 0.0;
2762     bestPossibleImprovement_ = 0.0;
[568]2763
[1525]2764     // Check any infeasibilities from dynamic rows
2765     matrix_->primalExpanded(this, 2);
2766     // Check any djs from dynamic rows
2767     matrix_->dualExpanded(this, NULL, NULL, 3);
2768     int numberDualInfeasibilitiesFree = 0;
2769     int firstFreePrimal = -1;
2770     int firstFreeDual = -1;
2771     int numberSuperBasicWithDj = 0;
[568]2772
[1525]2773     int numberTotal = numberRows_ + numberColumns_;
2774     // Say no free or superbasic
2775     moreSpecialOptions_ |= 8;
[1696]2776     //#define PRINT_INFEAS
2777#ifdef PRINT_INFEAS
2778     int seqInf[10];
2779#endif
[1525]2780     for (iSequence = 0; iSequence < numberTotal; iSequence++) {
2781          double value = solution_[iSequence];
[618]2782#ifdef COIN_DEBUG
[1525]2783          if (fabs(value) > 1.0e20)
2784               printf("%d values %g %g %g - status %d\n", iSequence, lower_[iSequence],
2785                      solution_[iSequence], upper_[iSequence], status_[iSequence]);
[618]2786#endif
[1525]2787          objectiveValue_ += value * cost_[iSequence];
2788          double distanceUp = upper_[iSequence] - value;
2789          double distanceDown = value - lower_[iSequence];
2790          if (distanceUp < -primalTolerance) {
2791               double infeasibility = -distanceUp;
2792               sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2793               if (infeasibility > relaxedToleranceP)
2794                    sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
[1696]2795#ifdef PRINT_INFEAS
2796               if (numberPrimalInfeasibilities_<10) {
2797                 seqInf[numberPrimalInfeasibilities_]=iSequence;
2798               }
2799#endif
[1525]2800               numberPrimalInfeasibilities_ ++;
2801          } else if (distanceDown < -primalTolerance) {
2802               double infeasibility = -distanceDown;
2803               sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2804               if (infeasibility > relaxedToleranceP)
2805                    sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
[1696]2806#ifdef PRINT_INFEAS
2807               if (numberPrimalInfeasibilities_<10) {
2808                 seqInf[numberPrimalInfeasibilities_]=iSequence;
2809               }
2810#endif
[1525]2811               numberPrimalInfeasibilities_ ++;
2812          } else {
2813               // feasible (so could be free)
2814               if (getStatus(iSequence) != basic && !flagged(iSequence)) {
2815                    // not basic
2816                    double djValue = dj_[iSequence];
2817                    if (distanceDown < primalTolerance) {
2818                         if (distanceUp > primalTolerance && djValue < -dualTolerance) {
2819                              sumDualInfeasibilities_ -= djValue + dualTolerance;
2820                              if (djValue < -possTolerance)
2821                                   bestPossibleImprovement_ -= distanceUp * djValue;
2822                              if (djValue < -relaxedToleranceD)
2823                                   sumOfRelaxedDualInfeasibilities_ -= djValue + relaxedToleranceD;
2824                              numberDualInfeasibilities_ ++;
2825                         }
2826                    } else if (distanceUp < primalTolerance) {
2827                         if (djValue > dualTolerance) {
2828                              sumDualInfeasibilities_ += djValue - dualTolerance;
2829                              if (djValue > possTolerance)
2830                                   bestPossibleImprovement_ += distanceDown * djValue;
2831                              if (djValue > relaxedToleranceD)
2832                                   sumOfRelaxedDualInfeasibilities_ += djValue - relaxedToleranceD;
2833                              numberDualInfeasibilities_ ++;
2834                         }
2835                    } else {
2836                         // may be free
2837                         // Say free or superbasic
2838                         moreSpecialOptions_ &= ~8;
2839                         djValue *= 0.01;
2840                         if (fabs(djValue) > dualTolerance) {
2841                              if (getStatus(iSequence) == isFree)
2842                                   numberDualInfeasibilitiesFree++;
2843                              sumDualInfeasibilities_ += fabs(djValue) - dualTolerance;
2844                              bestPossibleImprovement_ = 1.0e100;
2845                              numberDualInfeasibilities_ ++;
2846                              if (fabs(djValue) > relaxedToleranceD) {
2847                                   sumOfRelaxedDualInfeasibilities_ += value - relaxedToleranceD;
2848                                   numberSuperBasicWithDj++;
2849                                   if (firstFreeDual < 0)
2850                                        firstFreeDual = iSequence;
2851                              }
2852                         }
2853                         if (firstFreePrimal < 0)
2854                              firstFreePrimal = iSequence;
2855                    }
2856               }
[568]2857          }
[1525]2858     }
2859     objectiveValue_ += objective_->nonlinearOffset();
2860     objectiveValue_ /= (objectiveScale_ * rhsScale_);
2861     numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ -
2862                                             numberDualInfeasibilitiesFree;
[1696]2863#ifdef PRINT_INFEAS
2864     if (numberPrimalInfeasibilities_<=10) {
2865       printf("---------------start-----------\n");
2866       if (!rowScale_) {
2867         for (int i=0;i<numberPrimalInfeasibilities_;i++) {
2868           int iSeq = seqInf[i];
2869           double infeas;
2870           if (solution_[iSeq]<lower_[iSeq])
2871             infeas = lower_[iSeq]-solution_[iSeq];
2872           else
2873             infeas = solution_[iSeq]-upper_[iSeq];
2874           if (iSeq<numberColumns_) {
2875             printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g\n",
2876                    iSeq,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas);
2877           } else {
2878             printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g\n",
2879                    iSeq-numberColumns_,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas);
2880           }
2881         }
2882       } else {
2883         for (int i=0;i<numberPrimalInfeasibilities_;i++) {
2884           int iSeq = seqInf[i];
2885           double infeas;
2886           if (solution_[iSeq]<lower_[iSeq])
2887             infeas = lower_[iSeq]-solution_[iSeq];
2888           else
2889             infeas = solution_[iSeq]-upper_[iSeq];
2890           double unscaled = infeas;
2891           if (iSeq<numberColumns_) {
2892             unscaled *= columnScale_[iSeq];
2893             printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
2894                    iSeq,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas,unscaled);
2895           } else {
2896             unscaled /= rowScale_[iSeq-numberColumns_];
2897             printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
2898                    iSeq-numberColumns_,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas,unscaled);
2899           }
2900         }
2901       }
2902     }
2903#endif
[1525]2904     if (algorithm_ < 0 && firstFreeDual >= 0) {
2905          // dual
2906          firstFree_ = firstFreeDual;
2907     } else if (numberSuperBasicWithDj ||
2908                (progress_.lastIterationNumber(0) <= 0)) {
2909          firstFree_ = firstFreePrimal;
2910     }
[568]2911}
[286]2912/* Adds multiple of a column into an array */
[1525]2913void
[286]2914ClpSimplex::add(double * array,
[1525]2915                int sequence, double multiplier) const
[286]2916{
[1525]2917     if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
2918          //slack
2919          array [sequence-numberColumns_] -= multiplier;
2920     } else {
2921          // column
2922          matrix_->add(this, array, sequence, multiplier);
2923     }
[286]2924}
[2]2925/*
[1525]2926  Unpacks one column of the matrix into indexed array
[2]2927*/
[1525]2928void
[119]2929ClpSimplex::unpack(CoinIndexedVector * rowArray) const
[2]2930{
[1525]2931     rowArray->clear();
2932     if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
2933          //slack
2934          rowArray->insert(sequenceIn_ - numberColumns_, -1.0);
2935     } else {
2936          // column
2937          matrix_->unpack(this, rowArray, sequenceIn_);
2938     }
[2]2939}
[1525]2940void
2941ClpSimplex::unpack(CoinIndexedVector * rowArray, int sequence) const
[2]2942{
[1525]2943     rowArray->clear();
2944     if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
2945          //slack
2946          rowArray->insert(sequence - numberColumns_, -1.0);
2947     } else {
2948          // column
2949          matrix_->unpack(this, rowArray, sequence);
2950     }
[2]2951}
[225]2952/*
[1525]2953  Unpacks one column of the matrix into indexed array
[225]2954*/
[1525]2955void
2956ClpSimplex::unpackPacked(CoinIndexedVector * rowArray)
[225]2957{
[1525]2958     rowArray->clear();
2959     if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
2960          //slack
2961          int * index = rowArray->getIndices();
2962          double * array = rowArray->denseVector();
2963          array[0] = -1.0;
2964          index[0] = sequenceIn_ - numberColumns_;
2965          rowArray->setNumElements(1);
2966          rowArray->setPackedMode(true);
2967     } else {
2968          // column
2969          matrix_->unpackPacked(this, rowArray, sequenceIn_);
2970     }
[225]2971}
[1525]2972void
2973ClpSimplex::unpackPacked(CoinIndexedVector * rowArray, int sequence)
[225]2974{
[1525]2975     rowArray->clear();
2976     if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
2977          //slack
2978          int * index = rowArray->getIndices();
2979          double * array = rowArray->denseVector();
2980          array[0] = -1.0;
2981          index[0] = sequence - numberColumns_;
2982          rowArray->setNumElements(1);
2983          rowArray->setPackedMode(true);
2984     } else {
2985          // column
2986          matrix_->unpackPacked(this, rowArray, sequence);
2987     }
[225]2988}
[1266]2989//static int x_gaps[4]={0,0,0,0};
[1034]2990//static int scale_times[]={0,0,0,0};
[50]2991bool
[1525]2992ClpSimplex::createRim(int what, bool makeRowCopy, int startFinishOptions)
[2]2993{
[1525]2994     bool goodMatrix = true;
2995     int saveLevel = handler_->logLevel();
2996     spareIntArray_[0] = 0;
2997     if (!matrix_->canGetRowCopy())
2998          makeRowCopy = false; // switch off row copy if can't produce
2999     // Arrays will be there and correct size unless what is 63
3000     bool newArrays = (what == 63);
3001     // We may be restarting with same size
3002     bool keepPivots = false;
3003     if (startFinishOptions == -1) {
3004          startFinishOptions = 0;
3005          keepPivots = true;
3006     }
3007     bool oldMatrix = ((startFinishOptions & 4) != 0 && (whatsChanged_ & 1) != 0);
3008     if (what == 63) {
3009          pivotRow_ = -1;
3010          if (!status_)
3011               createStatus();
3012          if (oldMatrix)
3013               newArrays = false;
3014          if (problemStatus_ == 10) {
3015               handler_->setLogLevel(0); // switch off messages
3016               if (rowArray_[0]) {
3017                    // stuff is still there
3018                    oldMatrix = true;
3019                    newArrays = false;
3020                    keepPivots = true;
3021                    for (int iRow = 0; iRow < 4; iRow++) {
3022                         rowArray_[iRow]->clear();
3023                    }
3024                    for (int iColumn = 0; iColumn < 2; iColumn++) {
3025                         columnArray_[iColumn]->clear();
3026                    }
3027               }
3028          } else if (factorization_) {
3029               // match up factorization messages
3030               if (handler_->logLevel() < 3)
3031                    factorization_->messageLevel(0);
3032               else
3033                    factorization_->messageLevel(CoinMax(3, factorization_->messageLevel()));
3034               /* Faster to keep pivots rather than re-scan matrix.  Matrix may have changed
3035                  i.e. oldMatrix false but okay as long as same number rows and status array exists
3036               */
3037               if ((startFinishOptions & 2) != 0 && factorization_->numberRows() == numberRows_ && status_)
3038                    keepPivots = true;
3039          }
3040          numberExtraRows_ = matrix_->generalExpanded(this, 2, maximumBasic_);
3041          if (numberExtraRows_ && newArrays) {
3042               // make sure status array large enough
3043               assert (status_);
3044               int numberOld = numberRows_ + numberColumns_;
3045               int numberNew = numberRows_ + numberColumns_ + numberExtraRows_;
3046               unsigned char * newStatus = new unsigned char [numberNew];
3047               memset(newStatus + numberOld, 0, numberExtraRows_);
3048               CoinMemcpyN(status_, numberOld, newStatus);
3049               delete [] status_;
3050               status_ = newStatus;
3051          }
3052     }
3053     int numberRows2 = numberRows_ + numberExtraRows_;
3054     int numberTotal = numberRows2 + numberColumns_;
3055     if ((specialOptions_ & 65536) != 0) {
3056          assert (!numberExtraRows_);
3057          if (!cost_ || numberRows2 > maximumInternalRows_ ||
3058                    numberColumns_ > maximumInternalColumns_) {
3059               newArrays = true;
3060               keepPivots = false;
3061               printf("createrim a %d rows, %d maximum rows %d maxinternal\n",
3062                      numberRows_, maximumRows_, maximumInternalRows_);
3063               int oldMaximumRows = maximumInternalRows_;
3064               int oldMaximumColumns = maximumInternalColumns_;
3065               if (cost_) {
3066                    if (numberRows2 > maximumInternalRows_)
3067                         maximumInternalRows_ = numberRows2;
3068                    if (numberColumns_ > maximumInternalColumns_)
3069                         maximumInternalColumns_ = numberColumns_;
3070               } else {
3071                    maximumInternalRows_ = numberRows2;
3072                    maximumInternalColumns_ = numberColumns_;
3073               }
3074               assert(maximumInternalRows_ == maximumRows_);
3075               assert(maximumInternalColumns_ == maximumColumns_);
3076               printf("createrim b %d rows, %d maximum rows, %d maxinternal\n",
3077                      numberRows_, maximumRows_, maximumInternalRows_);
3078               int numberTotal2 = (maximumInternalRows_ + maximumInternalColumns_) * 2;
3079               delete [] cost_;
3080               cost_ = new double[numberTotal2];
3081               delete [] lower_;
3082               delete [] upper_;
3083               lower_ = new double[numberTotal2];
3084               upper_ = new double[numberTotal2];
3085               delete [] dj_;
3086               dj_ = new double[numberTotal2];
3087               delete [] solution_;
3088               solution_ = new double[numberTotal2];
3089               // ***** should be non NULL but seems to be too much
3090               //printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_);
3091               if (savedRowScale_) {
3092                    assert (oldMaximumRows > 0);
3093                    double * temp;
3094                    temp = new double [4*maximumRows_];
3095                    CoinFillN(temp, 4 * maximumRows_, 1.0);
3096                    CoinMemcpyN(savedRowScale_, numberRows_, temp);
3097                    CoinMemcpyN(savedRowScale_ + oldMaximumRows, numberRows_, temp + maximumRows_);
3098                    CoinMemcpyN(savedRowScale_ + 2 * oldMaximumRows, numberRows_, temp + 2 * maximumRows_);
3099                    CoinMemcpyN(savedRowScale_ + 3 * oldMaximumRows, numberRows_, temp + 3 * maximumRows_);
3100                    delete [] savedRowScale_;
3101                    savedRowScale_ = temp;
3102                    temp = new double [4*maximumColumns_];
3103                    CoinFillN(temp, 4 * maximumColumns_, 1.0);
3104                    CoinMemcpyN(savedColumnScale_, numberColumns_, temp);
3105                    CoinMemcpyN(savedColumnScale_ + oldMaximumColumns, numberColumns_, temp + maximumColumns_);
3106                    CoinMemcpyN(savedColumnScale_ + 2 * oldMaximumColumns, numberColumns_, temp + 2 * maximumColumns_);
3107                    CoinMemcpyN(savedColumnScale_ + 3 * oldMaximumColumns, numberColumns_, temp + 3 * maximumColumns_);
3108                    delete [] savedColumnScale_;
3109                    savedColumnScale_ = temp;
3110               }
3111          }
3112     }
3113     int i;
3114     bool doSanityCheck = true;
3115     if (what == 63) {
3116          // We may want to switch stuff off for speed
3117          if ((specialOptions_ & 256) != 0)
3118               makeRowCopy = false; // no row copy
3119          if ((specialOptions_ & 128) != 0)
3120               doSanityCheck = false; // no sanity check
3121          //check matrix
3122          if (!matrix_)
3123               matrix_ = new ClpPackedMatrix();
3124          int checkType = (doSanityCheck) ? 15 : 14;
3125          if (oldMatrix)
3126               checkType = 14;
3127          bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
3128          if (inCbcOrOther)
3129               checkType -= 4; // don't check for duplicates
3130          if (!matrix_->allElementsInRange(this, smallElement_, 1.0e20, checkType)) {
3131               problemStatus_ = 4;
3132               secondaryStatus_ = 8;
3133               //goodMatrix= false;
3134               return false;
3135          }
3136          bool rowCopyIsScaled;
3137          if (makeRowCopy) {
3138               if(!oldMatrix || !rowCopy_) {
3139                    delete rowCopy_;
3140                    // may return NULL if can't give row copy
3141                    rowCopy_ = matrix_->reverseOrderedCopy();
3142                    rowCopyIsScaled = false;
3143               } else {
3144                    rowCopyIsScaled = true;
3145               }
3146          }
[1034]3147#if 0
[1525]3148          if (what == 63) {
3149               int k = rowScale_ ? 1 : 0;
3150               if (oldMatrix)
3151                    k += 2;
3152               scale_times[k]++;
3153               if ((scale_times[0] + scale_times[1] + scale_times[2] + scale_times[3]) % 1000 == 0)
3154                    printf("scale counts %d %d %d %d\n",
3155                           scale_times[0], scale_times[1], scale_times[2], scale_times[3]);
3156          }
[1034]3157#endif
[1525]3158          // do scaling if needed
3159          if (!oldMatrix && scalingFlag_ < 0) {
3160               if (scalingFlag_ < 0 && rowScale_) {
3161                    //if (handler_->logLevel()>0)
3162                    printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
3163                           scalingFlag_);
3164                    scalingFlag_ = 0;
3165               }
3166               delete [] rowScale_;
3167               delete [] columnScale_;
3168               rowScale_ = NULL;
3169               columnScale_ = NULL;
3170          }
3171          inverseRowScale_ = NULL;
3172          inverseColumnScale_ = NULL;
3173          if (scalingFlag_ > 0 && !rowScale_) {
3174               if ((specialOptions_ & 65536) != 0) {
3175                    assert (!rowScale_);
3176                    rowScale_ = savedRowScale_;
3177                    columnScale_ = savedColumnScale_;
3178                    // put back original
3179                    if (savedRowScale_) {
3180                         inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
3181                         inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
3182                         CoinMemcpyN(savedRowScale_ + 2 * maximumInternalRows_,
3183                                     numberRows2, savedRowScale_);
3184                         CoinMemcpyN(savedRowScale_ + 3 * maximumInternalRows_,
3185                                     numberRows2, inverseRowScale_);
3186                         CoinMemcpyN(savedColumnScale_ + 2 * maximumColumns_,
3187                                     numberColumns_, savedColumnScale_);
3188                         CoinMemcpyN(savedColumnScale_ + 3 * maximumColumns_,
3189                                     numberColumns_, inverseColumnScale_);
3190                    }
3191               }
3192               if (matrix_->scale(this))
3193                    scalingFlag_ = -scalingFlag_; // not scaled after all
3194               if (rowScale_ && automaticScale_) {
3195                    // try automatic scaling
3196                    double smallestObj = 1.0e100;
3197                    double largestObj = 0.0;
3198                    double largestRhs = 0.0;
3199                    const double * obj = objective();
3200                    for (i = 0; i < numberColumns_; i++) {
3201                         double value = fabs(obj[i]);
3202                         value *= columnScale_[i];
3203                         if (value && columnLower_[i] != columnUpper_[i]) {
3204                              smallestObj = CoinMin(smallestObj, value);
3205                              largestObj = CoinMax(largestObj, value);
3206                         }
3207                         if (columnLower_[i] > 0.0 || columnUpper_[i] < 0.0) {
3208                              double scale = 1.0 * inverseColumnScale_[i];
3209                              //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3210                              if (columnLower_[i] > 0)
3211                                   largestRhs = CoinMax(largestRhs, columnLower_[i] * scale);
3212                              if (columnUpper_[i] < 0.0)
3213                                   largestRhs = CoinMax(largestRhs, -columnUpper_[i] * scale);
3214                         }
3215                    }
3216                    for (i = 0; i < numberRows_; i++) {
3217                         if (rowLower_[i] > 0.0 || rowUpper_[i] < 0.0) {
3218                              double scale = rowScale_[i];
3219                              //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3220                              if (rowLower_[i] > 0)
3221                                   largestRhs = CoinMax(largestRhs, rowLower_[i] * scale);
3222                              if (rowUpper_[i] < 0.0)
3223                                   largestRhs = CoinMax(largestRhs, -rowUpper_[i] * scale);
3224                         }
3225                    }
3226                    printf("small obj %g, large %g - rhs %g\n", smallestObj, largestObj, largestRhs);
3227                    bool scalingDone = false;
3228                    // look at element range
3229                    double smallestNegative;
3230                    double largestNegative;
3231                    double smallestPositive;
3232                    double largestPositive;
3233                    matrix_->rangeOfElements(smallestNegative, largestNegative,
3234                                             smallestPositive, largestPositive);
3235                    smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
3236                    largestPositive = CoinMax(fabs(largestNegative), largestPositive);
3237                    if (largestObj) {
3238                         double ratio = largestObj / smallestObj;
3239                         double scale = 1.0;
3240                         if (ratio < 1.0e8) {
3241                              // reasonable
3242                              if (smallestObj < 1.0e-4) {
3243                                   // may as well scale up
3244                                   scalingDone = true;
3245                                   scale = 1.0e-3 / smallestObj;
3246                              } else if (largestObj < 1.0e6 || (algorithm_ > 0 && largestObj < 1.0e-4 * infeasibilityCost_)) {
3247                                   //done=true;
3248                              } else {
3249                                   scalingDone = true;
3250                                   if (algorithm_ < 0) {
3251                                        scale = 1.0e6 / largestObj;
3252                                   } else {
3253                                        scale = CoinMax(1.0e6, 1.0e-4 * infeasibilityCost_) / largestObj;
3254                                   }
3255                              }
3256                         } else if (ratio < 1.0e12) {
3257                              // not so good
3258                              if (smallestObj < 1.0e-7) {
3259                                   // may as well scale up
3260                                   scalingDone = true;
3261                                   scale = 1.0e-6 / smallestObj;
3262                              } else if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3263                                   //done=true;
3264                              } else {
3265                                   scalingDone = true;
3266                                   if (algorithm_ < 0) {
3267                                        scale = 1.0e7 / largestObj;
3268                                   } else {
3269                                        scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3270                                   }
3271                              }
3272                         } else {
3273                              // Really nasty problem
3274                              if (smallestObj < 1.0e-8) {
3275                                   // may as well scale up
3276                                   scalingDone = true;
3277                                   scale = 1.0e-7 / smallestObj;
3278                                   largestObj *= scale;
3279                              }
3280                              if (largestObj < 1.0e7 || (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3281                                   //done=true;
3282                              } else {
3283                                   scalingDone = true;
3284                                   if (algorithm_ < 0) {
3285                                        scale = 1.0e7 / largestObj;
3286                                   } else {
3287                                        scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3288                                   }
3289                              }
3290                         }
3291                         objectiveScale_ = scale;
3292                    }
3293                    if (largestRhs > 1.0e12) {
3294                         scalingDone = true;
3295                         rhsScale_ = 1.0e9 / largestRhs;
3296                    } else if (largestPositive > 1.0e-14 * smallestPositive && largestRhs > 1.0e6) {
3297                         scalingDone = true;
3298                         rhsScale_ = 1.0e6 / largestRhs;
3299                    } else {
3300                         rhsScale_ = 1.0;
3301                    }
3302                    if (scalingDone) {
3303                         handler_->message(CLP_RIM_SCALE, messages_)
3304                                   << objectiveScale_ << rhsScale_
3305                                   << CoinMessageEol;
3306                    }
3307               }
3308          } else if (makeRowCopy && scalingFlag_ > 0 && !rowCopyIsScaled) {
3309               matrix_->scaleRowCopy(this);
3310          }
3311          if (rowScale_ && !savedRowScale_) {
3312               inverseRowScale_ = rowScale_ + numberRows2;
3313               inverseColumnScale_ = columnScale_ + numberColumns_;
3314          }
3315          // See if we can try for faster row copy
3316          if (makeRowCopy && !oldMatrix) {
3317               ClpPackedMatrix* clpMatrix =
3318                    dynamic_cast< ClpPackedMatrix*>(matrix_);
3319               if (clpMatrix && numberThreads_)
3320                    clpMatrix->specialRowCopy(this, rowCopy_);
3321               if (clpMatrix)
3322                    clpMatrix->specialColumnCopy(this);
3323          }
3324     }
3325     if (what == 63) {
[1266]3326#if 0
[1525]3327          {
3328               x_gaps[0]++;
3329               ClpPackedMatrix* clpMatrix =
3330               dynamic_cast< ClpPackedMatrix*>(matrix_);
3331               if (clpMatrix) {
3332                    if (!clpMatrix->getPackedMatrix()->hasGaps())
3333                         x_gaps[1]++;
3334                    if ((clpMatrix->flags() & 2) == 0)
3335                         x_gaps[3]++;
3336               } else {
3337                    x_gaps[2]++;
3338               }
3339               if ((x_gaps[0] % 1000) == 0)
3340                    printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n",
3341                           x_gaps[0], x_gaps[1], x_gaps[2], x_gaps[3]);
3342          }
[1266]3343#endif
[1525]3344          if (newArrays && (specialOptions_ & 65536) == 0) {
3345               delete [] cost_;
3346               cost_ = new double[2*numberTotal];
3347               delete [] lower_;
3348               delete [] upper_;
3349               lower_ = new double[numberTotal];
3350               upper_ = new double[numberTotal];
3351               delete [] dj_;
3352               dj_ = new double[numberTotal];
3353               delete [] solution_;
3354               solution_ = new double[numberTotal];
3355          }
3356          reducedCostWork_ = dj_;
3357          rowReducedCost_ = dj_ + numberColumns_;
3358          columnActivityWork_ = solution_;
3359          rowActivityWork_ = solution_ + numberColumns_;
3360          objectiveWork_ = cost_;
3361          rowObjectiveWork_ = cost_ + numberColumns_;
3362          rowLowerWork_ = lower_ + numberColumns_;
3363          columnLowerWork_ = lower_;
3364          rowUpperWork_ = upper_ + numberColumns_;
3365          columnUpperWork_ = upper_;
3366     }
3367     if ((what & 4) != 0) {
3368          double direction = optimizationDirection_ * objectiveScale_;
3369          const double * obj = objective();
3370          const double * rowScale = rowScale_;
3371          const double * columnScale = columnScale_;
3372          // and also scale by scale factors
3373          if (rowScale) {
3374               if (rowObjective_) {
3375                    for (i = 0; i < numberRows_; i++)
3376                         rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
3377               } else {
3378                    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3379               }
3380               // If scaled then do all columns later in one loop
3381               if (what != 63) {
3382                    for (i = 0; i < numberColumns_; i++) {
3383                         CoinAssert(fabs(obj[i]) < 1.0e25);
3384                         objectiveWork_[i] = obj[i] * direction * columnScale[i];
3385                    }
3386               }
3387          } else {
3388               if (rowObjective_) {
3389                    for (i = 0; i < numberRows_; i++)
3390                         rowObjectiveWork_[i] = rowObjective_[i] * direction;
3391               } else {
3392                    memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3393               }
3394               for (i = 0; i < numberColumns_; i++) {
3395                    CoinAssert(fabs(obj[i]) < 1.0e25);
3396                    objectiveWork_[i] = obj[i] * direction;
3397               }
3398          }
3399     }
3400     if ((what & 1) != 0) {
3401          const double * rowScale = rowScale_;
3402          // clean up any mismatches on infinity
3403          // and fix any variables with tiny gaps
3404          double primalTolerance = dblParam_[ClpPrimalTolerance];
3405          if(rowScale) {
3406               // If scaled then do all columns later in one loop
3407               if (what != 63) {
3408                    const double * inverseScale = inverseColumnScale_;
3409                    for (i = 0; i < numberColumns_; i++) {
3410                         double multiplier = rhsScale_ * inverseScale[i];
3411                         double lowerValue = columnLower_[i];
3412                         double upperValue = columnUpper_[i];
3413                         if (lowerValue > -1.0e20) {
3414                              columnLowerWork_[i] = lowerValue * multiplier;
3415                              if (upperValue >= 1.0e20) {
3416                                   columnUpperWork_[i] = COIN_DBL_MAX;