Changeset 50


Ignore:
Timestamp:
Nov 5, 2002 1:04:44 AM (17 years ago)
Author:
ladanyi
Message:

devel-1 merged into HEAD

Location:
trunk
Files:
40 added
34 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpDualRowDantzig.cpp

    r2 r50  
    99#include "ClpSimplex.hpp"
    1010#include "ClpDualRowDantzig.hpp"
    11 #include "OsiIndexedVector.hpp"
     11#include "CoinIndexedVector.hpp"
    1212
    1313//#############################################################################
     
    8686*/
    8787void
    88 ClpDualRowDantzig::updatePrimalSolution(
    89                                         OsiIndexedVector * primalUpdate,
     88ClpDualRowDantzig::updatePrimalSolution(CoinIndexedVector * primalUpdate,
    9089                                        double primalRatio,
    9190                                        double & objectiveChange)
  • trunk/ClpDualRowPivot.cpp

    r2 r50  
    6666
    6767void
    68 ClpDualRowPivot::updateWeights(OsiIndexedVector * input,
    69                                OsiIndexedVector * spare,
    70                                OsiIndexedVector * updatedColumn)
     68ClpDualRowPivot::updateWeights(CoinIndexedVector * input,
     69                               CoinIndexedVector * spare,
     70                               CoinIndexedVector * updatedColumn)
    7171{
    7272}
     
    7575{
    7676}
     77// Gets rid of all arrays
     78void
     79ClpDualRowPivot::clearArrays()
     80{
     81}
  • trunk/ClpDualRowSteepest.cpp

    r2 r50  
    99#include "ClpSimplex.hpp"
    1010#include "ClpDualRowSteepest.hpp"
    11 #include "OsiIndexedVector.hpp"
     11#include "CoinIndexedVector.hpp"
    1212#include "ClpFactorization.hpp"
    1313#include "CoinHelperFunctions.hpp"
    14 #include <stdio.h>
     14#include <cstdio>
    1515
    1616//#############################################################################
     
    3030    savedWeights_(NULL)
    3131{
    32   type_=2;
     32  type_=2+64*mode;
    3333}
    3434
     
    4343  model_ = rhs.model_;
    4444  if (rhs.infeasible_) {
    45     infeasible_= new OsiIndexedVector(rhs.infeasible_);
     45    infeasible_= new CoinIndexedVector(rhs.infeasible_);
    4646  } else {
    4747    infeasible_=NULL;
     
    5151    int number = model_->numberRows();
    5252    weights_= new double[number];
    53     CoinDisjointCopyN(rhs.weights_,number,weights_);
     53    ClpDisjointCopyN(rhs.weights_,number,weights_);
    5454  } else {
    5555    weights_=NULL;
    5656  }
    5757  if (rhs.alternateWeights_) {
    58     alternateWeights_= new OsiIndexedVector(rhs.alternateWeights_);
     58    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
    5959  } else {
    6060    alternateWeights_=NULL;
    6161  }
    6262  if (rhs.savedWeights_) {
    63     savedWeights_= new OsiIndexedVector(rhs.savedWeights_);
     63    savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
    6464  } else {
    6565    savedWeights_=NULL;
     
    9494    delete savedWeights_;
    9595    if (rhs.infeasible_!=NULL) {
    96       infeasible_= new OsiIndexedVector(rhs.infeasible_);
     96      infeasible_= new CoinIndexedVector(rhs.infeasible_);
    9797    } else {
    9898      infeasible_=NULL;
     
    102102      int number = model_->numberRows();
    103103      weights_= new double[number];
    104       CoinDisjointCopyN(rhs.weights_,number,weights_);
     104      ClpDisjointCopyN(rhs.weights_,number,weights_);
    105105    } else {
    106106      weights_=NULL;
    107107    }
    108108    if (rhs.alternateWeights_!=NULL) {
    109       alternateWeights_= new OsiIndexedVector(rhs.alternateWeights_);
     109      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
    110110    } else {
    111111      alternateWeights_=NULL;
    112112    }
    113113    if (rhs.savedWeights_!=NULL) {
    114       savedWeights_= new OsiIndexedVector(rhs.savedWeights_);
     114      savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
    115115    } else {
    116116      savedWeights_=NULL;
     
    153153// Updates weights
    154154void
    155 ClpDualRowSteepest::updateWeights(OsiIndexedVector * input,
    156                                   OsiIndexedVector * spare,
    157                                   OsiIndexedVector * updatedColumn)
     155ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
     156                                  CoinIndexedVector * spare,
     157                                  CoinIndexedVector * updatedColumn)
    158158{
    159159  // clear other region
     
    172172  {
    173173    int numberRows = model_->numberRows();
    174     OsiIndexedVector * temp = new OsiIndexedVector();
     174    CoinIndexedVector * temp = new CoinIndexedVector();
    175175    temp->reserve(numberRows+
    176176                  model_->factorization()->maximumPivots());
     
    245245      }
    246246    }
     247#ifdef CLP_DEBUG
    247248    assert(work3[pivotRow]&&work[pivotRow]);
     249#endif
    248250    alternateWeights_->setNumElements(nSave);
    249251    if (norm < TRY_NORM)
     
    263265void
    264266ClpDualRowSteepest::updatePrimalSolution(
    265                                         OsiIndexedVector * primalUpdate,
     267                                        CoinIndexedVector * primalUpdate,
    266268                                        double primalRatio,
    267269                                        double & objectiveChange)
     
    275277  const int * pivotVariable = model_->pivotVariable();
    276278  double * infeas = infeasible_->denseVector();
     279  int pivotRow = model_->pivotRow();
     280  double * solution = model_->solutionRegion();
    277281  for (i=0;i<number;i++) {
    278282    int iRow=which[i];
    279283    int iPivot=pivotVariable[iRow];
    280     double & value = model_->solutionAddress(iPivot);
     284    double value = solution[iPivot];
    281285    double cost = model_->cost(iPivot);
    282286    double change = primalRatio*work[iRow];
    283287    value -= change;
     288    changeObj -= change*cost;
     289    solution[iPivot] = value;
    284290    double lower = model_->lower(iPivot);
    285291    double upper = model_->upper(iPivot);
     292    // But if pivot row then use value of incoming
     293    if (iRow==pivotRow) {
     294      iPivot = model_->sequenceIn();
     295      // make last resort choice
     296      lower = 1.0e-6*model_->lower(iPivot);
     297      upper = 1.0e-6*model_->upper(iPivot);
     298      value = 1.0e-6*model_->valueIncomingDual();
     299    }
    286300    if (value>upper+tolerance) {
    287301      // store square in list
     
    301315        infeas[iRow] = 1.0e-100;
    302316    }
    303     changeObj -= change*cost;
    304317    work[iRow]=0.0;
    305318  }
     
    332345    }
    333346    state_=1;
    334   } else if (mode==2||mode==4) {
     347  } else if (mode==2||mode==4||mode==5) {
    335348    // restore
    336     if (!weights_||state_==-1) {
     349    if (!weights_||state_==-1||mode==5) {
    337350      // initialize weights
    338351      delete [] weights_;
    339352      delete alternateWeights_;
    340353      weights_ = new double[numberRows];
    341       alternateWeights_ = new OsiIndexedVector();
     354      alternateWeights_ = new CoinIndexedVector();
    342355      // enough space so can use it for factorization
    343356      alternateWeights_->reserve(numberRows+
    344357                                 model_->factorization()->maximumPivots());
    345       if (!mode_) {
     358      if (!mode_||mode==5) {
    346359        // initialize to 1.0 (can we do better?)
    347360        for (i=0;i<numberRows;i++) {
     
    349362        }
    350363      } else {
    351         OsiIndexedVector * temp = new OsiIndexedVector();
     364        CoinIndexedVector * temp = new CoinIndexedVector();
    352365        temp->reserve(numberRows+
    353366                      model_->factorization()->maximumPivots());
     
    374387      }
    375388      // create saved weights (not really indexedvector)
    376       savedWeights_ = new OsiIndexedVector();
     389      savedWeights_ = new CoinIndexedVector();
    377390      savedWeights_->reserve(numberRows);
    378391     
     
    415428    // set up infeasibilities
    416429    if (!infeasible_) {
    417       infeasible_ = new OsiIndexedVector();
     430      infeasible_ = new CoinIndexedVector();
    418431      infeasible_->reserve(numberRows);
    419432    }
     
    468481  }
    469482}
    470 
     483// Gets rid of all arrays
     484void
     485ClpDualRowSteepest::clearArrays()
     486{
     487  delete [] weights_;
     488  weights_=NULL;
     489  delete infeasible_;
     490  infeasible_ = NULL;
     491  delete alternateWeights_;
     492  alternateWeights_ = NULL;
     493  delete savedWeights_;
     494  savedWeights_ = NULL;
     495  state_ =-1;
     496}
     497
  • trunk/ClpFactorization.cpp

    r2 r50  
    1919// Default Constructor
    2020//-------------------------------------------------------------------
    21 ClpFactorization::ClpFactorization ()
    22   : OsiFactorization()
    23 {
    24 }
     21ClpFactorization::ClpFactorization () :
     22   CoinFactorization() {}
    2523
    2624//-------------------------------------------------------------------
    2725// Copy constructor
    2826//-------------------------------------------------------------------
    29 ClpFactorization::ClpFactorization (const ClpFactorization & rhs)
    30 : OsiFactorization(rhs)
    31 
    32  
    33 }
     27ClpFactorization::ClpFactorization (const ClpFactorization & rhs) :
     28   CoinFactorization(rhs) {}
    3429
    35 ClpFactorization::ClpFactorization (const OsiFactorization & rhs)
    36 : OsiFactorization(rhs)
    37 
    38  
    39 }
     30ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
     31   CoinFactorization(rhs) {}
    4032
    4133//-------------------------------------------------------------------
    4234// Destructor
    4335//-------------------------------------------------------------------
    44 ClpFactorization::~ClpFactorization ()
    45 {
    46 }
     36ClpFactorization::~ClpFactorization () {}
    4737
    4838//----------------------------------------------------------------
     
    5343{
    5444  if (this != &rhs) {
    55     OsiFactorization::operator=(rhs);
     45    CoinFactorization::operator=(rhs);
    5646  }
    5747  return *this;
     
    7060    areaFactor_ = areaFactor;
    7161  int numberBasic = 0;
    72   OsiBigIndex numberElements=0;
     62  CoinBigIndex numberElements=0;
    7363  int numberRowBasic=0;
    7464
     
    9181  if ( numberBasic > numberRows ) {
    9282    return -2; // say too many in basis
    93   }                             
     83  }
    9484  numberElements = 3 * numberBasic + 3 * numberElements + 10000;
    9585  getAreas ( numberRows, numberBasic, numberElements,
     
    134124        // these arrays start off as copies of permute
    135125        // (and we could use permute_ instead of pivotColumn (not back though))
    136         CoinDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
    137         CoinDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
     126        ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
     127        ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
    138128      }
    139129    } else {
  • trunk/ClpMessage.cpp

    r2 r50  
    2323  {CLP_SIMPLEX_STATUS,5,1,"%d  Objective %g%? Primal infeas %g (%d)%? Dual infeas %g (%d)%? without free dual infeas (%d)"},
    2424  {CLP_DUAL_BOUNDS,5,3,"Looking optimal checking bounds with %g"},
    25   {CLP_SIMPLEX_ACCURACY,6,1,"Primal error %g, dual error %g"},
     25  {CLP_SIMPLEX_ACCURACY,6,3,"Primal error %g, dual error %g"},
    2626  {CLP_SIMPLEX_BADFACTOR,7,1,"Singular factorization of basis - status %d"},
    2727  {CLP_SIMPLEX_BOUNDTIGHTEN,8,1,"Bounds were tightened %d times"},
    2828  {CLP_SIMPLEX_INFEASIBILITIES,9,1,"%d infeasibilities"},
    29   {CLP_SIMPLEX_FLAG,10,1,"Flagging variable %c%d"},
     29  {CLP_SIMPLEX_FLAG,10,3,"Flagging variable %c%d"},
    3030  {CLP_SIMPLEX_GIVINGUP,11,2,"Stopping as close enough"},
    3131  {CLP_DUAL_CHECKB,12,3,"Looking infeasible checking bounds with %g"},
     
    3535  {CLP_PRIMAL_WEIGHT,16,2,"New infeasibility weight of %g"},
    3636  {CLP_PRIMAL_OPTIMAL,17,2,"Looking optimal with tolerance of %g"},
     37  {CLP_SINGULARITIES,18,2,"%d total structurals rejected in initial factorization"},
     38  {CLP_MODIFIEDBOUNDS,19,1,"%d variables/rows fixed as scaled bounds too close"},
     39  {CLP_RIMSTATISTICS1,20,2,"Absolute values of scaled objective range from %g to %g"},
     40  {CLP_RIMSTATISTICS2,21,2,"Absolute values of scaled bounds range from %g to %g, minimum gap %g"},
     41  {CLP_RIMSTATISTICS3,22,2,"Absolute values of scaled rhs range from %g to %g, minimum gap %g"},
     42  {CLP_POSSIBLELOOP,23,2,"Possible loop - %d matches (%x) after %d checks"},
     43  {CLP_SMALLELEMENTS,24,0,"Matrix will be packed to eliminate %d small elements"},
    3744  {CLP_SIMPLEX_HOUSE1,101,32,"dirOut %d, dirIn %d, theta %g, out %g, dj %g, alpha %g"},
    3845  {CLP_SIMPLEX_HOUSE2,102,4,"%d %g In: %c%d Out: %c%d%? dj ratio %g distance %g%? dj %g distance %g"},
     
    4249  {CLP_DUAL_CHECK,106,4,"Btran alpha %g, ftran alpha %g"},
    4350  {CLP_PRIMAL_DJ,107,4,"Btran dj %g, ftran dj %g"},
     51  {CLP_PRESOLVE_COLINFEAS,501,2,"Problem is infeasible due to column %d, %g %g"},
     52  {CLP_PRESOLVE_ROWINFEAS,502,2,"Problem is infeasible due to row %d, %g %g"},
     53  {CLP_PRESOLVE_COLUMNBOUNDA,503,2,"Problem looks unbounded above due to column %d, %g %g"},
     54  {CLP_PRESOLVE_COLUMNBOUNDB,504,2,"Problem looks unbounded below due to column %d, %g %g"},
     55  {CLP_PRESOLVE_NONOPTIMAL,505,1,"Problem not optimal, resolve after postsolve"},
     56  {CLP_PRESOLVE_STATS,506,1,"Presolve %d (%d) rows, %d (%d) columns and %d (%d) elements"},
     57  {CLP_PRESOLVE_INFEAS,507,0,"Presolve determined that the problem was infeasible with tolerance of %g"},
     58  {CLP_PRESOLVE_UNBOUND,508,0,"Presolve thinks problem is unbounded"},
     59  {CLP_PRESOLVE_INFEASUNBOUND,509,0,"Presolve thinks problem is infeasible AND unbounded???"},
     60  {CLP_PRESOLVE_INTEGERMODS,510,1,"Presolve is modifying %d integer bounds and re-presolving"},
     61  {CLP_PRESOLVE_POSTSOLVE,511,0,"After Postsolve, objective %g, infeasibilities - dual %g (%d), primal %g (%d)"},
     62  {CLP_PRESOLVE_NEEDS_CLEANING,512,0,"Presolved model was optimal, full model needs cleaning up"},
    4463  {CLP_PACKEDSCALE_INITIAL,1001,2,"Initial range of elements is %g to %g"},
    4564  {CLP_PACKEDSCALE_WHILE,1002,3,"Range of elements is %g to %g"},
     
    4867  {CLP_INITIALIZE_STEEP,1005,1,"Initializing steepest edge weights - old %g, new %g"},
    4968  {CLP_UNABLE_OPEN,6001,0,"Unable to open file %s for reading"},
     69  {CLP_BAD_BOUNDS,6002,0,"%d bad bound pairs were found - first at %c%d"},
     70  {CLP_BAD_MATRIX,6003,0,"Matrix has %d large values, first at column %d, row %d is %g"},
     71  {CLP_LOOP,6004,0,"Can't get out of loop - stopping"},
    5072  {CLP_IMPORT_RESULT,18,1,"Model was imported from %s in %g seconds"},
    5173  {CLP_IMPORT_ERRORS,3001,1," There were %d errors when importing model from %s"},
     74  {CLP_EMPTY_PROBLEM,3002,0,"Not solving empty problem - %d rows, %d columns and %d elements"},
    5275  {CLP_DUMMY_END,999999,0,""}
    5376};
     
    6083/* Constructor */
    6184ClpMessage::ClpMessage(Language language) :
    62   OsiMessages(sizeof(us_english)/sizeof(Clp_message))
     85  CoinMessages(sizeof(us_english)/sizeof(Clp_message))
    6386{
    6487  language_=language;
     
    6790
    6891  while (message->internalNumber!=CLP_DUMMY_END) {
    69       OsiOneMessage oneMessage(message->externalNumber,message->detail,
    70                 message->message);
    71     addMessage(message->internalNumber,oneMessage);
    72     message ++;
     92     CoinOneMessage oneMessage(message->externalNumber,message->detail,
     93                               message->message);
     94     addMessage(message->internalNumber,oneMessage);
     95     message ++;
    7396}
    7497
  • trunk/ClpModel.cpp

    r2 r50  
    66#endif
    77
    8 #include <math.h>
     8#include <cmath>
     9#include <cassert>
     10#include <cfloat>
     11#include <string>
     12#include <cstdio>
     13#include <iostream>
     14
     15#include <time.h>
     16#include <sys/times.h>
     17#include <sys/resource.h>
     18#include <unistd.h>
    919
    1020#include "CoinHelperFunctions.hpp"
    1121#include "ClpModel.hpp"
    1222#include "ClpPackedMatrix.hpp"
    13 #include "OsiIndexedVector.hpp"
    14 #include "OsiMpsReader.hpp"
     23#include "CoinIndexedVector.hpp"
     24#include "CoinMpsIO.hpp"
    1525#include "ClpMessage.hpp"
    16 #include <cassert>
    17 #include <cfloat>
    18 
    19 #include <string>
    20 #include <stdio.h>
    21 #include <iostream>
    22 #include  <time.h>
    23 #include <sys/times.h>
    24 #include <sys/resource.h>
    25 #include <unistd.h>
    26 // This returns a non const array filled with input from scalar
    27 // or actual array
    28 template <class T> inline T*
    29 copyOfArray( const T * array, const int size, T value)
    30 {
    31   T * arrayNew = new T[size];
    32   if (array)
    33     CoinDisjointCopyN(array,size,arrayNew);
    34   else
    35     CoinFillN ( arrayNew, size,value);
    36   return arrayNew;
    37 }
    38 
    39 // This returns a non const array filled with actual array (or NULL)
    40 template <class T> inline T*
    41 copyOfArray( const T * array, const int size)
    42 {
    43   if (array) {
    44     T * arrayNew = new T[size];
    45     CoinDisjointCopyN(array,size,arrayNew);
    46     return arrayNew;
    47   } else {
    48     return NULL;
    49   }
    50 }
     26
    5127static double cpuTime()
    5228{
     
    9167  maximumIterations_(1000000000),
    9268  defaultHandler_(true),
     69  status_(NULL),
    9370  lengthNames_(0),
    9471  rowNames_(),
    95   columnNames_()
    96 {
    97   intParam_[OsiMaxNumIteration] = 9999999;
    98   intParam_[OsiMaxNumIterationHotStart] = 9999999;
    99 
    100   dblParam_[OsiDualObjectiveLimit] = DBL_MAX;
    101   dblParam_[OsiPrimalObjectiveLimit] = DBL_MAX;
    102   dblParam_[OsiDualTolerance] = 1e-7;
    103   dblParam_[OsiPrimalTolerance] = 1e-7;
    104   dblParam_[OsiObjOffset] = 0.0;
    105 
    106   strParam_[OsiProbName] = "OsiDefaultName";
    107   handler_ = new OsiMessageHandler();
     72  columnNames_(),
     73  integerType_(NULL)
     74{
     75  intParam_[ClpMaxNumIteration] = 9999999;
     76  intParam_[ClpMaxNumIterationHotStart] = 9999999;
     77
     78  dblParam_[ClpDualObjectiveLimit] = DBL_MAX;
     79  dblParam_[ClpPrimalObjectiveLimit] = DBL_MAX;
     80  dblParam_[ClpDualTolerance] = 1e-7;
     81  dblParam_[ClpPrimalTolerance] = 1e-7;
     82  dblParam_[ClpObjOffset] = 0.0;
     83
     84  strParam_[ClpProbName] = "ClpDefaultName";
     85  handler_ = new CoinMessageHandler();
    10886  handler_->setLogLevel(2);
    10987  messages_ = ClpMessage();
     
    148126  delete [] ray_;
    149127  ray_ = NULL;
     128  delete [] integerType_;
     129  integerType_ = NULL;
     130  delete [] status_;
     131  status_=NULL;
    150132}
    151133//#############################################################################
     
    153135{
    154136  if (value>0.0&&value<1.0e10)
    155     dblParam_[OsiPrimalTolerance]=value;
     137    dblParam_[ClpPrimalTolerance]=value;
    156138}
    157139void ClpModel::setDualTolerance( double value)
    158140{
    159141  if (value>0.0&&value<1.0e10)
    160     dblParam_[OsiDualTolerance]=value;
     142    dblParam_[ClpDualTolerance]=value;
    161143}
    162144void ClpModel::setOptimizationDirection( int value)
     
    180162  reducedCost_=new double[numberColumns_];
    181163
    182   CoinFillN(dual_,numberRows_,0.0);
    183   CoinFillN(reducedCost_,numberColumns_,0.0);
     164  ClpFillN(dual_,numberRows_,0.0);
     165  ClpFillN(reducedCost_,numberColumns_,0.0);
    184166  int iRow,iColumn;
    185167
    186   rowLower_=copyOfArray(rowlb,numberRows_,-DBL_MAX);
    187   rowUpper_=copyOfArray(rowub,numberRows_,DBL_MAX);
    188   objective_=copyOfArray(obj,numberColumns_,0.0);
    189   rowObjective_=copyOfArray(rowObjective,numberRows_);
    190   columnLower_=copyOfArray(collb,numberColumns_,0.0);
    191   columnUpper_=copyOfArray(colub,numberColumns_,DBL_MAX);
     168  rowLower_=ClpCopyOfArray(rowlb,numberRows_,-DBL_MAX);
     169  rowUpper_=ClpCopyOfArray(rowub,numberRows_,DBL_MAX);
     170  objective_=ClpCopyOfArray(obj,numberColumns_,0.0);
     171  rowObjective_=ClpCopyOfArray(rowObjective,numberRows_);
     172  columnLower_=ClpCopyOfArray(collb,numberColumns_,0.0);
     173  columnUpper_=ClpCopyOfArray(colub,numberColumns_,DBL_MAX);
    192174  // set default solution
    193175  for (iRow=0;iRow<numberRows_;iRow++) {
     
    223205  } else {
    224206    // later may want to keep as unknown class
    225     OsiPackedMatrix matrix2;
     207    CoinPackedMatrix matrix2;
    226208    matrix2.reverseOrderedCopyOf(*matrix.getPackedMatrix());
    227209    matrix.releasePackedMatrix();
     
    230212}
    231213void
    232 ClpModel::loadProblem (  const OsiPackedMatrix& matrix,
     214ClpModel::loadProblem (  const CoinPackedMatrix& matrix,
    233215                     const double* collb, const double* colub,   
    234216                     const double* obj,
     
    241223    matrix_=new ClpPackedMatrix(matrix);
    242224  } else {
    243     OsiPackedMatrix matrix2;
     225    CoinPackedMatrix matrix2;
    244226    matrix2.reverseOrderedCopyOf(matrix);
    245227    matrix_=new ClpPackedMatrix(matrix2);
     
    249231ClpModel::loadProblem (
    250232                              const int numcols, const int numrows,
    251                               const int* start, const int* index,
     233                              const CoinBigIndex* start, const int* index,
    252234                              const double* value,
    253235                              const double* collb, const double* colub,   
     
    258240  gutsOfLoadModel(numrows, numcols,
    259241                  collb, colub, obj, rowlb, rowub, rowObjective);
    260   OsiPackedMatrix matrix(true,numrows,numcols,start[numcols],
     242  CoinPackedMatrix matrix(true,numrows,numcols,start[numcols],
    261243                              value,index,start,NULL);
     244  matrix_ = new ClpPackedMatrix(matrix);
     245}
     246void
     247ClpModel::loadProblem (
     248                              const int numcols, const int numrows,
     249                              const CoinBigIndex* start, const int* index,
     250                              const double* value,const int* length,
     251                              const double* collb, const double* colub,   
     252                              const double* obj,
     253                              const double* rowlb, const double* rowub,
     254                              const double * rowObjective)
     255{
     256  gutsOfLoadModel(numrows, numcols,
     257                  collb, colub, obj, rowlb, rowub, rowObjective);
     258  // Compute number of elements
     259  int numberElements = 0;
     260  int i;
     261  for (i=0;i<numcols;i++)
     262    numberElements += length[i];
     263  CoinPackedMatrix matrix(true,numrows,numcols,numberElements,
     264                              value,index,start,length);
    262265  matrix_ = new ClpPackedMatrix(matrix);
    263266}
     
    304307  defaultHandler_ = rhs.defaultHandler_;
    305308  if (defaultHandler_)
    306     handler_ = new OsiMessageHandler(*rhs.handler_);
     309    handler_ = new CoinMessageHandler(*rhs.handler_);
    307310   else
    308311    handler_ = rhs.handler_;
    309312  messages_ = rhs.messages_;
    310   intParam_[OsiMaxNumIteration] = rhs.intParam_[OsiMaxNumIteration];
    311   intParam_[OsiMaxNumIterationHotStart] =
    312     rhs.intParam_[OsiMaxNumIterationHotStart];
    313 
    314   dblParam_[OsiDualObjectiveLimit] = rhs.dblParam_[OsiDualObjectiveLimit];
    315   dblParam_[OsiPrimalObjectiveLimit] = rhs.dblParam_[OsiPrimalObjectiveLimit];
    316   dblParam_[OsiDualTolerance] = rhs.dblParam_[OsiDualTolerance];
    317   dblParam_[OsiPrimalTolerance] = rhs.dblParam_[OsiPrimalTolerance];
    318   dblParam_[OsiObjOffset] = rhs.dblParam_[OsiObjOffset];
    319 
    320   strParam_[OsiProbName] = rhs.strParam_[OsiProbName];
     313  intParam_[ClpMaxNumIteration] = rhs.intParam_[ClpMaxNumIteration];
     314  intParam_[ClpMaxNumIterationHotStart] =
     315    rhs.intParam_[ClpMaxNumIterationHotStart];
     316
     317  dblParam_[ClpDualObjectiveLimit] = rhs.dblParam_[ClpDualObjectiveLimit];
     318  dblParam_[ClpPrimalObjectiveLimit] = rhs.dblParam_[ClpPrimalObjectiveLimit];
     319  dblParam_[ClpDualTolerance] = rhs.dblParam_[ClpDualTolerance];
     320  dblParam_[ClpPrimalTolerance] = rhs.dblParam_[ClpPrimalTolerance];
     321  dblParam_[ClpObjOffset] = rhs.dblParam_[ClpObjOffset];
     322
     323  strParam_[ClpProbName] = rhs.strParam_[ClpProbName];
    321324
    322325  objectiveValue_=rhs.objectiveValue_;
     
    328331    rowNames_ = rhs.rowNames_;
    329332    columnNames_ = rhs.columnNames_;
     333    if (rhs.integerType_) {
     334      integerType_ = new char[numberColumns_];
     335      memcpy(integerType_,rhs.integerType_,numberColumns_*sizeof(char));
     336    } else {
     337      integerType_ = NULL;
     338    }
    330339    if (rhs.rowActivity_) {
    331340      rowActivity_=new double[numberRows_];
     
    333342      dual_=new double[numberRows_];
    334343      reducedCost_=new double[numberColumns_];
    335       CoinDisjointCopyN ( rhs.rowActivity_, numberRows_ ,
     344      ClpDisjointCopyN ( rhs.rowActivity_, numberRows_ ,
    336345                          rowActivity_);
    337       CoinDisjointCopyN ( rhs.columnActivity_, numberColumns_ ,
     346      ClpDisjointCopyN ( rhs.columnActivity_, numberColumns_ ,
    338347                          columnActivity_);
    339       CoinDisjointCopyN ( rhs.dual_, numberRows_ ,
     348      ClpDisjointCopyN ( rhs.dual_, numberRows_ ,
    340349                          dual_);
    341       CoinDisjointCopyN ( rhs.reducedCost_, numberColumns_ ,
     350      ClpDisjointCopyN ( rhs.reducedCost_, numberColumns_ ,
    342351                          reducedCost_);
    343352    } else {
     
    347356      reducedCost_=NULL;
    348357    }
    349     rowLower_ = copyOfArray ( rhs.rowLower_, numberRows_ );
    350     rowUpper_ = copyOfArray ( rhs.rowUpper_, numberRows_ );
    351     columnLower_ = copyOfArray ( rhs.columnLower_, numberColumns_ );
    352     columnUpper_ = copyOfArray ( rhs.columnUpper_, numberColumns_ );
    353     objective_ = copyOfArray ( rhs.objective_, numberColumns_ );
    354     rowObjective_ = copyOfArray ( rhs.rowObjective_, numberRows_ );
     358    rowLower_ = ClpCopyOfArray ( rhs.rowLower_, numberRows_ );
     359    rowUpper_ = ClpCopyOfArray ( rhs.rowUpper_, numberRows_ );
     360    columnLower_ = ClpCopyOfArray ( rhs.columnLower_, numberColumns_ );
     361    columnUpper_ = ClpCopyOfArray ( rhs.columnUpper_, numberColumns_ );
     362    objective_ = ClpCopyOfArray ( rhs.objective_, numberColumns_ );
     363    rowObjective_ = ClpCopyOfArray ( rhs.rowObjective_, numberRows_ );
     364    status_ = ClpCopyOfArray( rhs.status_,numberColumns_+numberRows_);
    355365    ray_ = NULL;
    356366    if (problemStatus_==1)
    357       ray_ = copyOfArray (rhs.ray_,numberRows_);
     367      ray_ = ClpCopyOfArray (rhs.ray_,numberRows_);
    358368    else if (problemStatus_==2)
    359       ray_ = copyOfArray (rhs.ray_,numberColumns_);
     369      ray_ = ClpCopyOfArray (rhs.ray_,numberColumns_);
    360370    if (rhs.rowCopy_) {
    361371      rowCopy_ = rhs.rowCopy_->clone();
     
    384394    rowNames_ = std::vector<std::string> ();
    385395    columnNames_ = std::vector<std::string> ();
     396    integerType_ = NULL;
     397    status_ = rhs.status_;
    386398  }
    387399}
     
    400412  numberRows_ = rhs.numberRows_;
    401413  numberColumns_ = rhs.numberColumns_;
     414  delete [] rhs.ray_;
     415  rhs.ray_=NULL;
    402416  gutsOfCopy(rhs,false);
    403417}
     
    421435  matrix_ = NULL;
    422436  rowCopy_ = NULL;
     437  delete [] otherModel.ray_;
     438  otherModel.ray_ = ray_;
    423439  ray_ = NULL;
     440  // do status
     441  if (otherModel.status_!=status_) {
     442    delete [] otherModel.status_;
     443    otherModel.status_ = status_;
     444  }
     445  status_ = NULL;
    424446}
    425447//#############################################################################
     
    428450
    429451bool
    430 ClpModel::setIntParam(OsiIntParam key, int value)
     452ClpModel::setIntParam(ClpIntParam key, int value)
    431453{
    432454  switch (key) {
    433   case OsiMaxNumIteration:
     455  case ClpMaxNumIteration:
    434456    if (value < 0)
    435457      return false;
    436458    break;
    437   case OsiMaxNumIterationHotStart:
     459  case ClpMaxNumIterationHotStart:
    438460    if (value < 0)
    439461      return false;
    440462    break;
    441   case OsiLastIntParam:
     463  case ClpLastIntParam:
    442464    return false;
    443465  }
     
    449471
    450472bool
    451 ClpModel::setDblParam(OsiDblParam key, double value)
     473ClpModel::setDblParam(ClpDblParam key, double value)
    452474{
    453475
    454476  switch (key) {
    455   case OsiDualObjectiveLimit:
     477  case ClpDualObjectiveLimit:
    456478    break;
    457479
    458   case OsiPrimalObjectiveLimit:
     480  case ClpPrimalObjectiveLimit:
    459481    break;
    460482
    461   case OsiDualTolerance:
     483  case ClpDualTolerance:
    462484    if (value<=0.0||value>1.0e10)
    463485      return false;
    464486    break;
    465487   
    466   case OsiPrimalTolerance:
     488  case ClpPrimalTolerance:
    467489    if (value<=0.0||value>1.0e10)
    468490      return false;
    469491    break;
    470492   
    471   case OsiObjOffset:
     493  case ClpObjOffset:
    472494    break;
    473495
    474   case OsiLastDblParam:
     496  case ClpLastDblParam:
    475497    return false;
    476498  }
     
    482504
    483505bool
    484 ClpModel::setStrParam(OsiStrParam key, const std::string & value)
     506ClpModel::setStrParam(ClpStrParam key, const std::string & value)
    485507{
    486508
    487509  switch (key) {
    488   case OsiProbName:
     510  case ClpProbName:
    489511    break;
    490512
    491   case OsiLastStrParam:
     513  case ClpLastStrParam:
    492514    return false;
    493515  }
     
    535557    }
    536558    delete [] array;
     559    array = newArray;
     560    delete [] deleted;
     561  }
     562  return array;
     563}
     564char * deleteChar(char * array , int size,
     565                  int number, const int * which,int & newSize,
     566                  bool ifDelete)
     567{
     568  if (array) {
     569    int i ;
     570    char * deleted = new char[size];
     571    int numberDeleted=0;
     572    memset(deleted,0,size*sizeof(char));
     573    for (i=0;i<number;i++) {
     574      int j = which[i];
     575      if (j>=0&&j<size&&!deleted[j]) {
     576        numberDeleted++;
     577        deleted[j]=1;
     578      }
     579    }
     580    newSize = size-numberDeleted;
     581    char * newArray = new char[newSize];
     582    int put=0;
     583    for (i=0;i<size;i++) {
     584      if (!deleted[i]) {
     585        newArray[put++]=array[i];
     586      }
     587    }
     588    if (ifDelete)
     589      delete [] array;
    537590    array = newArray;
    538591    delete [] deleted;
     
    578631    ray_ = NULL;
    579632  }
     633  if (status_) {
     634    unsigned char * tempC = new unsigned char [newNumberColumns+newNumberRows];
     635    unsigned char * tempR = tempC + newNumberColumns;
     636    memset(tempC,0,(newNumberColumns+newNumberRows)*sizeof(unsigned char));
     637    memcpy(tempC,status_,min(newNumberColumns,numberColumns_)*sizeof(unsigned char));
     638    memcpy(tempR,status_+numberColumns_,min(newNumberRows,numberRows_)*sizeof(unsigned char));
     639    delete [] status_;
     640    status_ = tempC;
     641  }
    580642  numberRows_ = newNumberRows;
    581643  if (newNumberColumns<numberColumns_) {
     
    586648    matrix_->deleteCols(numberColumns_-newNumberColumns,which);
    587649    delete [] which;
     650  }
     651  if (integerType_) {
     652    char * temp = new char [newNumberColumns];
     653    memset(temp,0,newNumberColumns*sizeof(char));
     654    memcpy(temp,integerType_,
     655           min(newNumberColumns,numberColumns_)*sizeof(char));
     656    delete [] integerType_;
     657    integerType_ = temp;
    588658  }
    589659  numberColumns_ = newNumberColumns;
     
    618688  rowNames_ = std::vector<std::string> ();
    619689  columnNames_ = std::vector<std::string> ();
     690  // status
     691  if (status_) {
     692    unsigned char * tempR  = (unsigned char *) deleteChar((char *)status_+numberColumns_,
     693                                        numberRows_,
     694                                        number, which, newSize,false);
     695    unsigned char * tempC = new unsigned char [numberColumns_+newSize];
     696    memcpy(tempC,status_,numberColumns_*sizeof(unsigned char));
     697    memcpy(tempC+numberColumns_,tempR,newSize*sizeof(unsigned char));
     698    delete [] tempR;
     699    delete [] status_;
     700    status_ = tempC;
     701  }
    620702}
    621703// Deletes columns
     
    644726  rowNames_ = std::vector<std::string> ();
    645727  columnNames_ = std::vector<std::string> ();
     728  integerType_ = deleteChar(integerType_,numberColumns_,
     729                            number, which, newSize,true);
     730  // status
     731  if (status_) {
     732    unsigned char * tempC  = (unsigned char *) deleteChar((char *)status_,
     733                                        numberColumns_,
     734                                        number, which, newSize,false);
     735    unsigned char * temp = new unsigned char [numberRows_+newSize];
     736    memcpy(temp,status_,newSize*sizeof(unsigned char));
     737    memcpy(temp+newSize,status_+numberColumns_,
     738           numberRows_*sizeof(unsigned char));
     739    delete [] tempC;
     740    delete [] status_;
     741    status_ = temp;
     742  }
    646743}
    647744// Infeasibility/unbounded ray (NULL returned if none/wrong)
     
    651748  double * array = NULL;
    652749  if (problemStatus_==1)
    653     array = copyOfArray(ray_,numberRows_);
     750    array = ClpCopyOfArray(ray_,numberRows_);
    654751  return array;
    655752}
     
    659756  double * array = NULL;
    660757  if (problemStatus_==2)
    661     array = copyOfArray(ray_,numberColumns_);
     758    array = ClpCopyOfArray(ray_,numberColumns_);
    662759  return array;
    663760}
     
    670767// Pass in Message handler (not deleted at end)
    671768void
    672 ClpModel::passInMessageHandler(OsiMessageHandler * handler)
     769ClpModel::passInMessageHandler(CoinMessageHandler * handler)
    673770{
    674771  defaultHandler_=false;
     
    677774// Set language
    678775void
    679 ClpModel::newLanguage(OsiMessages::Language language)
     776ClpModel::newLanguage(CoinMessages::Language language)
    680777{
    681778  messages_ = ClpMessage(language);
     
    700797    } else {
    701798      handler_->message(CLP_UNABLE_OPEN,messages_)
    702         <<fileName<<OsiMessageEol;
     799        <<fileName<<CoinMessageEol;
    703800      return -1;
    704801    }
    705802  }
    706   OsiMpsReader m;
     803  CoinMpsIO m;
    707804  double time1 = cpuTime(),time2;
    708805  int status=m.readMps(fileName,"");
     
    712809                m.getObjCoefficients(),
    713810                m.getRowLower(),m.getRowUpper());
     811    if (m.integerColumns()) {
     812      integerType_ = new char[numberColumns_];
     813      memcpy(integerType_,m.integerColumns(),numberColumns_*sizeof(char));
     814    } else {
     815      integerType_ = NULL;
     816    }
     817    // set problem name
     818    setStrParam(ClpProbName,m.getProblemName());
    714819    // do names
    715820    if (keepNames) {
     
    734839      lengthNames_=0;
    735840    }
    736     setDblParam(OsiObjOffset,m.objectiveOffset());
     841    setDblParam(ClpObjOffset,m.objectiveOffset());
    737842    time2 = cpuTime();
    738843    handler_->message(CLP_IMPORT_RESULT,messages_)
    739844      <<fileName
    740       <<time2-time1<<OsiMessageEol;
     845      <<time2-time1<<CoinMessageEol;
    741846  } else {
    742847    // errors
    743848    handler_->message(CLP_IMPORT_ERRORS,messages_)
    744       <<status<<fileName<<OsiMessageEol;
     849      <<status<<fileName<<CoinMessageEol;
    745850  }
    746851
    747852  return status;
    748853}
     854bool ClpModel::isPrimalObjectiveLimitReached() const
     855{
     856  double limit = 0.0;
     857  getDblParam(ClpPrimalObjectiveLimit, limit);
     858  if (limit > 1e30) {
     859    // was not ever set
     860    return false;
     861  }
     862   
     863  const double obj = objectiveValue();
     864  const int maxmin = optimizationDirection();
     865
     866  if (problemStatus_ == 0) // optimal
     867    return maxmin > 0 ? (obj < limit) /*minim*/ : (obj > limit) /*maxim*/;
     868  else if (problemStatus_==2)
     869    return true;
     870  else
     871    return false;
     872}
     873
     874bool ClpModel::isDualObjectiveLimitReached() const
     875{
     876
     877  double limit = 0.0;
     878  getDblParam(ClpDualObjectiveLimit, limit);
     879  if (limit > 1e30) {
     880    // was not ever set
     881    return false;
     882  }
     883   
     884  const double obj = objectiveValue();
     885  const int maxmin = optimizationDirection();
     886
     887  if (problemStatus_ == 0) // optimal
     888    return maxmin > 0 ? (obj > limit) /*minim*/ : (obj < limit) /*maxim*/;
     889  else if (problemStatus_==1)
     890    return true;
     891  else
     892    return false;
     893
     894}
     895void
     896ClpModel::copyInIntegerInformation(const char * information)
     897{
     898  delete [] integerType_;
     899  if (information) {
     900    integerType_ = new char[numberColumns_];
     901    memcpy(integerType_,information,numberColumns_*sizeof(char));
     902  } else {
     903    integerType_ = NULL;
     904  }
     905}
     906// Drops names - makes lengthnames 0 and names empty
     907void
     908ClpModel::dropNames()
     909{
     910  lengthNames_=0;
     911  rowNames_ = std::vector<std::string> ();
     912  columnNames_ = std::vector<std::string> ();
     913}
     914// Drop integer informations
     915void
     916ClpModel::deleteIntegerInformation()
     917{
     918  delete [] integerType_;
     919  integerType_ = NULL;
     920}
     921/* Return copy of status array (char[numberRows+numberColumns]),
     922   use delete [] */
     923unsigned char * 
     924ClpModel::statusCopy() const
     925{
     926  return ClpCopyOfArray(status_,numberRows_+numberColumns_);
     927}
     928// Copy in status vector
     929void
     930ClpModel::copyinStatus(const unsigned char * statusArray)
     931{
     932  delete [] status_;
     933  if (statusArray) {
     934    status_ = new unsigned char [numberRows_+numberColumns_];
     935    memcpy(status_,statusArray,(numberRows_+numberColumns_)*sizeof(unsigned char));
     936  } else {
     937    status_=NULL;
     938  }
     939}
  • trunk/ClpNonLinearCost.cpp

    r2 r50  
    88#include <iostream>
    99
     10#include "CoinIndexedVector.hpp"
     11
    1012#include "ClpNonLinearCost.hpp"
    1113#include "ClpSimplex.hpp"
    12 #include "OsiIndexedVector.hpp"
    1314
    1415//#############################################################################
     
    3132  changeCost_(0.0),
    3233  largestInfeasibility_(0.0),
     34  sumInfeasibilities_(0.0),
    3335  convex_(true)
    3436{
     
    5456  changeCost_=0.0;
    5557  double infeasibilityCost = model_->infeasibilityCost();
     58  sumInfeasibilities_=0.0;
    5659  largestInfeasibility_=0.0;
    5760
     
    115118  double infeasibilityCost = model_->infeasibilityCost();
    116119  largestInfeasibility_=0.0;
     120  sumInfeasibilities_=0.0;
    117121
    118122  int iSequence;
     
    196200  changeCost_(0.0),
    197201  largestInfeasibility_(0.0),
     202  sumInfeasibilities_(0.0),
    198203  convex_(true)
    199204
     
    215220    changeCost_ = rhs.changeCost_;
    216221    largestInfeasibility_ = rhs.largestInfeasibility_;
     222    sumInfeasibilities_ = rhs.sumInfeasibilities_;
    217223    convex_ = rhs.convex_;
    218224  }
     
    268274    changeCost_ = rhs.changeCost_;
    269275    largestInfeasibility_ = rhs.largestInfeasibility_;
     276    sumInfeasibilities_ = rhs.sumInfeasibilities_;
    270277    convex_ = rhs.convex_;
    271278  }
     
    284291  changeCost_=0.0;
    285292  largestInfeasibility_ = 0.0;
     293  sumInfeasibilities_ = 0.0;
    286294  double primalTolerance = model_->currentPrimalTolerance();
    287295 
     
    328336      // slot in here
    329337      if (value<lower[iSequence]-primalTolerance) {
    330         largestInfeasibility_ = max(largestInfeasibility_,
    331                                     lower[iSequence]-value);
     338        value = lower[iSequence]-value;
     339        sumInfeasibilities_ += value;
     340        largestInfeasibility_ = max(largestInfeasibility_,value);
    332341        changeCost_ -= lower[iSequence]*
    333342          (cost_[iRange]-cost[iSequence]);
    334343        numberInfeasibilities_++;
    335344      } else if (value>upper[iSequence]+primalTolerance) {
    336         largestInfeasibility_ = max(largestInfeasibility_,
    337                                     value-upper[iSequence]);
     345        value = value-upper[iSequence];
     346        sumInfeasibilities_ += value;
     347        largestInfeasibility_ = max(largestInfeasibility_,value);
    338348        changeCost_ -= upper[iSequence]*
    339349          (cost_[iRange]-cost[iSequence]);
     
    350360    case ClpSimplex::atUpperBound:
    351361      if (!toNearest) {
    352         assert(fabs(value-upperValue)<=primalTolerance*1.0001) ;
     362        // With increasing tolerances - we may be at wrong place
     363        if (fabs(value-upperValue)>primalTolerance*1.0001) {
     364          assert(fabs(value-lowerValue)<=primalTolerance*1.0001);
     365          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
     366        }
    353367      } else {
    354368        solution[iSequence] = upperValue;
     
    357371    case ClpSimplex::atLowerBound:
    358372      if (!toNearest) {
    359         assert(fabs(value-lowerValue)<=primalTolerance*1.0001);
     373        // With increasing tolerances - we may be at wrong place
     374        if (fabs(value-lowerValue)>primalTolerance*1.0001) {
     375          assert(fabs(value-upperValue)<=primalTolerance*1.0001);
     376          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
     377        }
    360378      } else {
    361379        solution[iSequence] = lowerValue;
     
    438456}
    439457void
    440 ClpNonLinearCost::goBackAll(const OsiIndexedVector * update)
     458ClpNonLinearCost::goBackAll(const CoinIndexedVector * update)
    441459{
    442460  assert (model_!=NULL);
     
    493511   for costs.
    494512   On input array is empty (but indices exist).  On exit just
    495    changed costs will be stored as normal OsiIndexedVector
     513   changed costs will be stored as normal CoinIndexedVector
    496514*/
    497515void
    498 ClpNonLinearCost::checkChanged(int numberInArray, OsiIndexedVector * update)
     516ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update)
    499517{
    500518  assert (model_!=NULL);
     
    573591  case ClpSimplex::atLowerBound:
    574592    // set correctly
    575     if (fabs(value-lower)<=primalTolerance*1.001) 
     593    if (fabs(value-lower)<=primalTolerance*1.001){
    576594      model_->setStatus(iPivot,ClpSimplex::atLowerBound);
    577     else if (fabs(value-upper)<=primalTolerance*1.001)
     595    } else if (fabs(value-upper)<=primalTolerance*1.001){
    578596      model_->setStatus(iPivot,ClpSimplex::atUpperBound);
    579     else
    580       assert(fabs(value-lower)<=primalTolerance*1.0001||
    581              fabs(value-upper)<=primalTolerance*1.0001);
     597    } else {
     598      // set superBasic
     599      model_->setStatus(iPivot,ClpSimplex::superBasic);
     600    }
    582601    break;
    583602  }
  • trunk/ClpPackedMatrix.cpp

    r2 r50  
    66#endif
    77
     8#include <cstdio>
     9
     10#include "CoinIndexedVector.hpp"
     11#include "CoinHelperFunctions.hpp"
     12
    813#include "ClpSimplex.hpp"
    914#include "ClpFactorization.hpp"
    10 #include "OsiIndexedVector.hpp"
    11 #include "CoinHelperFunctions.hpp"
    12 #include <stdio.h>
    1315// at end to get min/max!
    1416#include "ClpPackedMatrix.hpp"
     
    3537: ClpMatrixBase(rhs)
    3638
    37   matrix_ = new OsiPackedMatrix(*(rhs.matrix_));
     39  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
    3840 
    3941}
    4042
    41 ClpPackedMatrix::ClpPackedMatrix (const OsiPackedMatrix & rhs)
     43//-------------------------------------------------------------------
     44// assign matrix (for space reasons)
     45//-------------------------------------------------------------------
     46ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs)
    4247: ClpMatrixBase()
    4348
    44   matrix_ = new OsiPackedMatrix(rhs);
     49  matrix_ = rhs;
     50  setType(1);
     51 
     52}
     53
     54ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs)
     55: ClpMatrixBase()
     56
     57  matrix_ = new CoinPackedMatrix(rhs);
    4558  setType(1);
    4659 
     
    6477    ClpMatrixBase::operator=(rhs);
    6578    delete matrix_;
    66     matrix_ = new OsiPackedMatrix(*(rhs.matrix_));
     79    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
    6780  }
    6881  return *this;
     
    8194{
    8295  ClpPackedMatrix * copy = new ClpPackedMatrix();
    83   copy->matrix_= new OsiPackedMatrix();
     96  copy->matrix_= new CoinPackedMatrix();
    8497  copy->matrix_->reverseOrderedCopyOf(*matrix_);
    8598  copy->matrix_->removeGaps();
     
    94107  // get matrix data pointers
    95108  const int * row = matrix_->getIndices();
    96   const int * columnStart = matrix_->getVectorStarts();
     109  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    97110  const int * columnLength = matrix_->getVectorLengths();
    98111  const double * elementByColumn = matrix_->getElements();
    99112  int numberColumns = matrix_->getNumCols();
     113  //memset(y,0,matrix_->getNumRows()*sizeof(double));
    100114  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    101     int j;
     115    CoinBigIndex j;
    102116    double value = scalar*x[iColumn];
    103117    if (value) {
     
    117131  // get matrix data pointers
    118132  const int * row = matrix_->getIndices();
    119   const int * columnStart = matrix_->getVectorStarts();
     133  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    120134  const int * columnLength = matrix_->getVectorLengths();
    121135  const double * elementByColumn = matrix_->getElements();
    122136  int numberColumns = matrix_->getNumCols();
     137  //memset(y,0,numberColumns*sizeof(double));
    123138  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    124     int j;
     139    CoinBigIndex j;
    125140    double value=0.0;
    126141    for (j=columnStart[iColumn];
     
    141156  // get matrix data pointers
    142157  const int * row = matrix_->getIndices();
    143   const int * columnStart = matrix_->getVectorStarts();
     158  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    144159  const int * columnLength = matrix_->getVectorLengths();
    145160  const double * elementByColumn = matrix_->getElements();
    146161  int numberColumns = matrix_->getNumCols();
     162  //memset(y,0,matrix_->getNumRows()*sizeof(double));
    147163  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    148     int j;
     164    CoinBigIndex j;
    149165    double value = x[iColumn];
    150166    if (value) {
     
    172188  // get matrix data pointers
    173189  const int * row = matrix_->getIndices();
    174   const int * columnStart = matrix_->getVectorStarts();
     190  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    175191  const int * columnLength = matrix_->getVectorLengths();
    176192  const double * elementByColumn = matrix_->getElements();
    177193  int numberColumns = matrix_->getNumCols();
     194  //memset(y,0,numberColumns*sizeof(double));
    178195  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    179     int j;
     196    CoinBigIndex j;
    180197    double value=0.0;
    181198    // scaled
     
    192209void
    193210ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
    194                               const OsiIndexedVector * rowArray,
    195                               OsiIndexedVector * y,
    196                               OsiIndexedVector * columnArray) const
     211                              const CoinIndexedVector * rowArray,
     212                              CoinIndexedVector * y,
     213                              CoinIndexedVector * columnArray) const
    197214{
    198215  columnArray->clear();
     
    214231    // get matrix data pointers
    215232    const int * row = matrix_->getIndices();
    216     const int * columnStart = matrix_->getVectorStarts();
     233    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    217234    const int * columnLength = matrix_->getVectorLengths();
    218235    const double * elementByColumn = matrix_->getElements();
     
    225242        markVector[iColumn]=0.0;
    226243        double value2 = 0.0;
    227         int j;
     244        CoinBigIndex j;
    228245        for (j=columnStart[iColumn];
    229246             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     
    242259        double value = markVector[iColumn];
    243260        markVector[iColumn]=0.0;
    244         int j;
     261        CoinBigIndex j;
    245262        const double * columnScale = model->columnScale();
    246263        double value2 = 0.0;
     
    280297
    281298    const int * column = rowCopy->getIndices();
    282     const int * rowStart = rowCopy->getVectorStarts();
     299    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    283300    const double * element = rowCopy->getElements();
    284301    const int * whichRow = rowArray->getIndices();
     
    286303      iRow = whichRow[i];
    287304      double value = pi[iRow]*scalar;
    288       int j;
     305      CoinBigIndex j;
    289306      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    290307        int iColumn = column[j];
     
    316333
    317334    const int * column = rowCopy->getIndices();
    318     const int * rowStart = rowCopy->getVectorStarts();
     335    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    319336    const double * element = rowCopy->getElements();
    320337    const int * whichRow = rowArray->getIndices();
     
    322339    iRow = whichRow[0];
    323340    value = pi[iRow]*scalar;
    324     int j;
     341    CoinBigIndex j;
    325342    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    326343      int iColumn = column[j];
     
    357374    numberNonZero=0;
    358375    const int * column = rowCopy->getIndices();
    359     const int * rowStart = rowCopy->getVectorStarts();
     376    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    360377    const double * element = rowCopy->getElements();
    361378    double value = pi[iRow]*scalar;
    362     int j;
     379    CoinBigIndex j;
    363380    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    364381      int iColumn = column[j];
     
    378395void
    379396ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
    380                               const OsiIndexedVector * rowArray,
    381                               const OsiIndexedVector * y,
    382                               OsiIndexedVector * columnArray) const
     397                              const CoinIndexedVector * rowArray,
     398                              const CoinIndexedVector * y,
     399                              CoinIndexedVector * columnArray) const
    383400{
    384401  columnArray->clear();
     
    392409  // get matrix data pointers
    393410  const int * row = matrix_->getIndices();
    394   const int * columnStart = matrix_->getVectorStarts();
     411  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    395412  const int * columnLength = matrix_->getVectorLengths();
    396413  const double * elementByColumn = matrix_->getElements();
     
    402419      int iColumn = which[jColumn];
    403420      double value = 0.0;
    404       int j;
     421      CoinBigIndex j;
    405422      for (j=columnStart[iColumn];
    406423           j<columnStart[iColumn]+columnLength[iColumn];j++) {
     
    418435      int iColumn = which[jColumn];
    419436      double value = 0.0;
    420       int j;
     437      CoinBigIndex j;
    421438      const double * columnScale = model->columnScale();
    422439      for (j=columnStart[iColumn];
     
    435452/* Returns number of elements in basis
    436453   column is basic if entry >=0 */
    437 OsiBigIndex
     454CoinBigIndex
    438455ClpPackedMatrix::numberInBasis(const int * columnIsBasic) const
    439456{
     
    441458  int numberColumns = getNumCols();
    442459  const int * columnLength = matrix_->getVectorLengths();
    443   OsiBigIndex numberElements=0;
     460  CoinBigIndex numberElements=0;
    444461  for (i=0;i<numberColumns;i++) {
    445462    if (columnIsBasic[i]>=0) {
     
    450467}
    451468// Fills in basis (Returns number of elements and updates numberBasic)
    452 OsiBigIndex
     469CoinBigIndex
    453470ClpPackedMatrix::fillBasis(const ClpSimplex * model,
    454471                                const int * columnIsBasic, int & numberBasic,
     
    458475  const double * rowScale = model->rowScale();
    459476  const int * row = matrix_->getIndices();
    460   const int * columnStart = matrix_->getVectorStarts();
     477  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    461478  const int * columnLength = matrix_->getVectorLengths();
    462479  const double * elementByColumn = matrix_->getElements();
    463480  int i;
    464481  int numberColumns = getNumCols();
    465   OsiBigIndex numberElements=0;
     482  CoinBigIndex numberElements=0;
    466483  if (!rowScale) {
    467484    // no scaling
    468485    for (i=0;i<numberColumns;i++) {
    469486      if (columnIsBasic[i]>=0) {
    470         int j;
     487        CoinBigIndex j;
    471488        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
    472489          indexRowU[numberElements]=row[j];
     
    482499    for (i=0;i<numberColumns;i++) {
    483500      if (columnIsBasic[i]>=0) {
    484         int j;
     501        CoinBigIndex j;
    485502        double scale = columnScale[i];
    486503        for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     
    513530  int numberColumns = model->numberColumns();
    514531  const int * column = rowCopy->getIndices();
    515   const int * rowStart = rowCopy->getVectorStarts();
     532  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    516533  const double * element = rowCopy->getElements();
    517534  double * rowScale = new double [numberRows];
     
    539556  // get matrix data pointers
    540557  const int * row = matrix_->getIndices();
    541   const int * columnStart = matrix_->getVectorStarts();
     558  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    542559  const int * columnLength = matrix_->getVectorLengths();
    543560  const double * elementByColumn = matrix_->getElements();
    544561  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    545     int j;
     562    CoinBigIndex j;
    546563    char useful=0;
    547564    if (columnUpper[iColumn]>
     
    561578  model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
    562579    <<smallest<<largest
    563     <<OsiMessageEol;
     580    <<CoinMessageEol;
    564581  if (smallest>=0.5&&largest<=2.0) {
    565582    // don't bother scaling
    566583    model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
    567       <<OsiMessageEol;
     584      <<CoinMessageEol;
    568585    delete [] rowScale;
    569586    delete [] usefulRow;
     
    575592    for (iRow=0;iRow<numberRows;iRow++) {
    576593      if (usefulRow[iRow]) {
    577         int j;
     594        CoinBigIndex j;
    578595        int useful=0;
    579596        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     
    587604      }
    588605    }
    589     CoinFillN ( rowScale, numberRows,1.0);
    590     CoinFillN ( columnScale, numberColumns,1.0);
     606    ClpFillN ( rowScale, numberRows,1.0);
     607    ClpFillN ( columnScale, numberColumns,1.0);
    591608    int numberPass=3;
    592609    double overallLargest;
     
    599616      for (iRow=0;iRow<numberRows;iRow++) {
    600617        if (usefulRow[iRow]) {
    601           int j;
     618          CoinBigIndex j;
    602619          largest=1.0e-20;
    603620          smallest=1.0e50;
     
    618635        <<overallSmallest
    619636        <<overallLargest
    620         <<OsiMessageEol;
     637        <<CoinMessageEol;
    621638      // skip last column round
    622639      if (numberPass==1)
     
    625642      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    626643        if (usefulColumn[iColumn]) {
    627           int j;
     644          CoinBigIndex j;
    628645          largest=1.0e-20;
    629646          smallest=1.0e50;
     
    646663    for (iColumn=0;iColumn<numberColumns;iColumn++) {
    647664      if (usefulColumn[iColumn]) {
    648         int j;
     665        CoinBigIndex j;
    649666        largest=1.0e-20;
    650667        smallest=1.0e50;
     
    666683      <<overallSmallest
    667684      <<overallLargest
    668       <<OsiMessageEol;
     685      <<CoinMessageEol;
    669686   
    670687    delete [] usefulRow;
     
    682699        const int * columnsInThisRow = column + rowStart[iRow];
    683700        int number = rowStart[iRow+1]-rowStart[iRow];
     701        assert (number<=numberColumns);
    684702        for (j=0;j<number;j++) {
    685703          int iColumn = columnsInThisRow[j];
     
    696714  }
    697715}
    698 /* Unpacks a column into an OsiIndexedvector
     716// Creates row copy and scales if necessary
     717ClpMatrixBase *
     718ClpPackedMatrix::scaledRowCopy(ClpSimplex * model) const
     719{
     720  ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
     721  ClpPackedMatrix* rowCopy =
     722    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
     723
     724  // Make sure it is really a ClpPackedMatrix
     725  assert (rowCopy!=NULL);
     726
     727  const double * rowScale = model->rowScale();
     728  const double * columnScale = model->columnScale();
     729
     730  if (rowScale) {
     731    // scale row copy
     732    int numberRows = model->numberRows();
     733    int numberColumns = model->numberColumns();
     734    const int * column = rowCopy->getIndices();
     735    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     736    const double * element = rowCopy->getElements();
     737    int iRow;
     738    // need to replace row by row
     739    double * newElement = new double[numberColumns];
     740    // scale row copy
     741    for (iRow=0;iRow<numberRows;iRow++) {
     742      int j;
     743      double scale = rowScale[iRow];
     744      const double * elementsInThisRow = element + rowStart[iRow];
     745      const int * columnsInThisRow = column + rowStart[iRow];
     746      int number = rowStart[iRow+1]-rowStart[iRow];
     747      assert (number<=numberColumns);
     748      for (j=0;j<number;j++) {
     749        int iColumn = columnsInThisRow[j];
     750        newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
     751      }
     752      rowCopy->replaceVector(iRow,number,newElement);
     753    }
     754    delete [] newElement;
     755  }
     756  return rowCopyBase;
     757}
     758/* Unpacks a column into an CoinIndexedvector
    699759      Note that model is NOT const.  Bounds and objective could
    700760      be modified if doing column generation */
    701761void
    702 ClpPackedMatrix::unpack(ClpSimplex * model,OsiIndexedVector * rowArray,
     762ClpPackedMatrix::unpack(ClpSimplex * model,CoinIndexedVector * rowArray,
    703763                   int iColumn) const
    704764{
    705765  const double * rowScale = model->rowScale();
    706766  const int * row = matrix_->getIndices();
    707   const int * columnStart = matrix_->getVectorStarts();
    708   const int * columnLength = matrix_->getVectorLengths();
    709   const double * elementByColumn = matrix_->getElements();
    710   int i;
     767  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     768  const int * columnLength = matrix_->getVectorLengths();
     769  const double * elementByColumn = matrix_->getElements();
     770  CoinBigIndex i;
    711771  if (!rowScale) {
    712772    for (i=columnStart[iColumn];
     
    724784  }
    725785}
    726 /* Adds multiple of a column into an OsiIndexedvector
     786/* Adds multiple of a column into an CoinIndexedvector
    727787      You can use quickAdd to add to vector */
    728788void
    729 ClpPackedMatrix::add(const ClpSimplex * model,OsiIndexedVector * rowArray,
     789ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
    730790                   int iColumn, double multiplier) const
    731791{
    732792  const double * rowScale = model->rowScale();
    733793  const int * row = matrix_->getIndices();
    734   const int * columnStart = matrix_->getVectorStarts();
    735   const int * columnLength = matrix_->getVectorLengths();
    736   const double * elementByColumn = matrix_->getElements();
    737   int i;
     794  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     795  const int * columnLength = matrix_->getVectorLengths();
     796  const double * elementByColumn = matrix_->getElements();
     797  CoinBigIndex i;
    738798  if (!rowScale) {
    739799    for (i=columnStart[iColumn];
     
    752812  }
    753813}
    754 
     814/* Checks if all elements are in valid range.  Can just
     815   return true if you are not paranoid.  For Clp I will
     816   probably expect no zeros.  Code can modify matrix to get rid of
     817   small elements.
     818*/
     819bool
     820ClpPackedMatrix::allElementsInRange(ClpSimplex * model,
     821                                    double smallest, double largest)
     822{
     823  int iColumn;
     824  CoinBigIndex numberLarge=0;;
     825  CoinBigIndex numberSmall=0;;
     826  int firstBadColumn=-1;
     827  int firstBadRow=-1;
     828  double firstBadElement=0.0;
     829  // get matrix data pointers
     830  const int * row = matrix_->getIndices();
     831  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     832  const int * columnLength = matrix_->getVectorLengths();
     833  const double * elementByColumn = matrix_->getElements();
     834  int numberColumns = matrix_->getNumCols();
     835  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     836    CoinBigIndex j;
     837    for (j=columnStart[iColumn];
     838         j<columnStart[iColumn]+columnLength[iColumn];j++) {
     839      double value = fabs(elementByColumn[j]);
     840      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     841      if (value<smallest) {
     842        numberSmall++;
     843      } else if (value>largest) {
     844        numberLarge++;
     845        if (firstBadColumn<0) {
     846          firstBadColumn=iColumn;
     847          firstBadRow=row[j];
     848          firstBadElement=elementByColumn[j];
     849        }
     850      }
     851    }
     852  }
     853  if (numberLarge) {
     854    model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
     855      <<numberLarge
     856      <<firstBadColumn<<firstBadRow<<firstBadElement
     857      <<CoinMessageEol;
     858    return false;
     859  }
     860  if (numberSmall) {
     861    model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
     862      <<numberSmall
     863      <<CoinMessageEol;
     864    matrix_->compress(smallest);
     865  }
     866  return true;
     867}
     868
     869
  • trunk/ClpPrimalColumnDantzig.cpp

    r2 r50  
    66#endif
    77
     8#include <cstdio>
     9
     10#include "CoinIndexedVector.hpp"
    811
    912#include "ClpSimplex.hpp"
    1013#include "ClpPrimalColumnDantzig.hpp"
    11 #include "OsiIndexedVector.hpp"
    1214#include "ClpFactorization.hpp"
    1315#include "ClpPackedMatrix.hpp"
    14 #include <stdio.h>
     16
    1517//#############################################################################
    1618// Constructors / Destructor / Assignment
     
    5759// Returns pivot column, -1 if none
    5860int
    59 ClpPrimalColumnDantzig::pivotColumn(OsiIndexedVector * updates,
    60                                     OsiIndexedVector * spareRow1,
    61                                     OsiIndexedVector * spareRow2,
    62                                     OsiIndexedVector * spareColumn1,
    63                                     OsiIndexedVector * spareColumn2)
     61ClpPrimalColumnDantzig::pivotColumn(CoinIndexedVector * updates,
     62                                    CoinIndexedVector * spareRow1,
     63                                    CoinIndexedVector * spareRow2,
     64                                    CoinIndexedVector * spareColumn1,
     65                                    CoinIndexedVector * spareColumn2)
    6466{
    6567  assert(model_);
     
    145147     
    146148      switch(status) {
    147        
     149
    148150      case ClpSimplex::basic:
    149151        break;
  • trunk/ClpPrimalColumnPivot.cpp

    r2 r50  
    6262
    6363void
    64 ClpPrimalColumnPivot::updateWeights(OsiIndexedVector * input)
     64ClpPrimalColumnPivot::updateWeights(CoinIndexedVector * input)
    6565{
    6666}
    6767
     68// Gets rid of all arrays
     69void
     70ClpPrimalColumnPivot::clearArrays()
     71{
     72}
  • trunk/ClpPrimalColumnSteepest.cpp

    r5 r50  
    99#include "ClpSimplex.hpp"
    1010#include "ClpPrimalColumnSteepest.hpp"
    11 #include "OsiIndexedVector.hpp"
     11#include "CoinIndexedVector.hpp"
    1212#include "ClpFactorization.hpp"
    1313#include "ClpMessage.hpp"
     
    3535    devex_(0.0)
    3636{
    37   type_=2;
     37  type_=2+64*mode;
    3838}
    3939
     
    5252  devex_ = rhs.devex_;
    5353  if (rhs.infeasible_) {
    54     infeasible_= new OsiIndexedVector(rhs.infeasible_);
     54    infeasible_= new CoinIndexedVector(rhs.infeasible_);
    5555  } else {
    5656    infeasible_=NULL;
     
    6161    int number = model_->numberRows()+model_->numberColumns();
    6262    weights_= new double[number];
    63     CoinDisjointCopyN(rhs.weights_,number,weights_);
     63    ClpDisjointCopyN(rhs.weights_,number,weights_);
    6464    savedWeights_= new double[number];
    65     CoinDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
     65    ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
    6666    if (!mode_) {
    6767      reference_ = new unsigned int[(number+31)>>5];
     
    7373  }
    7474  if (rhs.alternateWeights_) {
    75     alternateWeights_= new OsiIndexedVector(rhs.alternateWeights_);
     75    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
    7676  } else {
    7777    alternateWeights_=NULL;
     
    114114    savedWeights_ = NULL;
    115115    if (rhs.infeasible_!=NULL) {
    116       infeasible_= new OsiIndexedVector(rhs.infeasible_);
     116      infeasible_= new CoinIndexedVector(rhs.infeasible_);
    117117    } else {
    118118      infeasible_=NULL;
     
    122122      int number = model_->numberRows()+model_->numberColumns();
    123123      weights_= new double[number];
    124       CoinDisjointCopyN(rhs.weights_,number,weights_);
     124      ClpDisjointCopyN(rhs.weights_,number,weights_);
    125125      savedWeights_= new double[number];
    126       CoinDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
     126      ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
    127127      if (!mode_) {
    128128        reference_ = new unsigned int[(number+31)>>5];
     
    134134    }
    135135    if (rhs.alternateWeights_!=NULL) {
    136       alternateWeights_= new OsiIndexedVector(rhs.alternateWeights_);
     136      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
    137137    } else {
    138138      alternateWeights_=NULL;
     
    146146// Returns pivot column, -1 if none
    147147int
    148 ClpPrimalColumnSteepest::pivotColumn(OsiIndexedVector * updates,
    149                                     OsiIndexedVector * spareRow1,
    150                                     OsiIndexedVector * spareRow2,
    151                                     OsiIndexedVector * spareColumn1,
    152                                     OsiIndexedVector * spareColumn2)
     148ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
     149                                    CoinIndexedVector * spareRow1,
     150                                    CoinIndexedVector * spareRow2,
     151                                    CoinIndexedVector * spareColumn1,
     152                                    CoinIndexedVector * spareColumn2)
    153153{
    154154  assert(model_);
     
    225225          continue;
    226226        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    227        
     227
    228228        switch(status) {
    229229         
     
    288288     
    289289      switch(status) {
    290        
     290
    291291      case ClpSimplex::basic:
    292292        break;
     
    352352     
    353353      switch(status) {
    354        
     354
    355355      case ClpSimplex::basic:
    356356        infeasible_->zero(iSequence+addSequence);
     
    457457      }
    458458      weight[iSequence] = thisWeight;
    459     }                   
     459    }
    460460   
    461461    // columns
     
    495495      }
    496496      weight[iSequence] = thisWeight;
    497     }                   
     497    }
    498498    // restore outgoing weight
    499499    weights_[sequenceOut]=outgoingWeight;
     
    534534  index = infeasible_->getIndices();
    535535  number = infeasible_->getNumElements();
    536   // we can't really trust infeasibilities if there is dual error
    537   double checkTolerance = 1.0e-8;
    538   if (!model_->factorization()->pivots())
    539     checkTolerance = 1.0e-6;
    540   if (model_->largestDualError()>checkTolerance)
    541     tolerance *= model_->largestDualError()/checkTolerance;
    542 
     536  if(model_->numberIterations()<model_->lastBadIteration()+200) {
     537    // we can't really trust infeasibilities if there is dual error
     538    double checkTolerance = 1.0e-8;
     539    if (!model_->factorization()->pivots())
     540      checkTolerance = 1.0e-6;
     541    if (model_->largestDualError()>checkTolerance)
     542      tolerance *= model_->largestDualError()/checkTolerance;
     543  }
     544  // stop last one coming immediately
     545  double saveOutInfeasibility=0.0;
     546  if (sequenceOut>=0) {
     547    saveOutInfeasibility = infeas[sequenceOut];
     548    infeas[sequenceOut]=0.0;
     549  }
    543550  tolerance *= tolerance; // as we are using squares
    544551  for (i=0;i<number;i++) {
     
    559566    }
    560567  }
     568  if (sequenceOut>=0) {
     569    infeas[sequenceOut]=saveOutInfeasibility;
     570  }
    561571  /*if (model_->numberIterations()%100==0)
    562572    printf("%d best %g\n",bestSequence,bestDj);*/
     
    605615      delete alternateWeights_;
    606616      weights_ = new double[numberRows+numberColumns];
    607       alternateWeights_ = new OsiIndexedVector();
     617      alternateWeights_ = new CoinIndexedVector();
    608618      // enough space so can use it for factorization
    609619      alternateWeights_->reserve(numberRows+
     
    628638        memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
    629639               sizeof(double));
    630         pivotSequence_= savedPivotSequence_;
    631         model_->setSequenceOut(savedSequenceOut_);
     640        // was - but I think should not be
     641        //pivotSequence_= savedPivotSequence_;
     642        //model_->setSequenceOut(savedSequenceOut_);
     643        pivotSequence_= -1;
     644        model_->setSequenceOut(-1);
    632645      }
    633646    }
     
    635648    // set up infeasibilities
    636649    if (!infeasible_) {
    637       infeasible_ = new OsiIndexedVector();
     650      infeasible_ = new CoinIndexedVector();
    638651      infeasible_->reserve(numberColumns+numberRows);
    639652    }
     
    664677      }
    665678      alternateWeights_->setNumElements(number);
     679#ifdef CLP_DEBUG
     680      // Can happen but I should clean up
    666681      assert(found>=0);
     682#endif
    667683      pivotSequence_ = found;
    668684      delete [] temp;
     
    682698     
    683699      switch(status) {
    684        
     700
    685701      case ClpSimplex::basic:
    686702        break;
     
    734750}
    735751void
    736 ClpPrimalColumnSteepest::updateWeights(OsiIndexedVector * input)
     752ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
    737753{
    738754  int number=input->getNumElements();
     
    824840                                        *model_->messagesPointer())
    825841                                          <<oldDevex<<devex_
    826                                           <<OsiMessageEol;
     842                                          <<CoinMessageEol;
    827843      initializeWeights();
    828     }                           
     844    }
    829845  }
    830846  if (pivotRow>=0) {
     
    837853ClpPrimalColumnSteepest::checkAccuracy(int sequence,
    838854                                       double relativeTolerance,
    839                                        OsiIndexedVector * rowArray1,
    840                                        OsiIndexedVector * rowArray2)
     855                                       CoinIndexedVector * rowArray1,
     856                                       CoinIndexedVector * rowArray2)
    841857{
    842858  model_->unpack(rowArray1,sequence);
     
    908924    }
    909925  } else {
    910     OsiIndexedVector * temp = new OsiIndexedVector();
     926    CoinIndexedVector * temp = new CoinIndexedVector();
    911927    temp->reserve(numberRows+
    912928                  model_->factorization()->maximumPivots());
     
    936952  }
    937953}
     954// Gets rid of all arrays
     955void
     956ClpPrimalColumnSteepest::clearArrays()
     957{
     958  delete [] weights_;
     959  weights_=NULL;
     960  delete infeasible_;
     961  infeasible_ = NULL;
     962  delete alternateWeights_;
     963  alternateWeights_ = NULL;
     964  delete [] savedWeights_;
     965  savedWeights_ = NULL;
     966  delete [] reference_;
     967  reference_ = NULL;
     968  pivotSequence_=-1;
     969  state_ = -1;
     970  savedPivotSequence_ = -1;
     971  savedSequenceOut_ = -1; 
     972  devex_ = 0.0;
     973}
  • trunk/ClpSimplex.cpp

    r2 r50  
    1414#include "ClpSimplex.hpp"
    1515#include "ClpFactorization.hpp"
    16 #include "OsiPackedMatrix.hpp"
    17 #include "OsiIndexedVector.hpp"
    18 #include "OsiWarmStartBasis.hpp"
     16#include "ClpPackedMatrix.hpp"
     17#include "CoinIndexedVector.hpp"
     18#include "ClpDualRowDantzig.hpp"
    1919#include "ClpDualRowSteepest.hpp"
    2020#include "ClpPrimalColumnDantzig.hpp"
     
    2727#include <stdio.h>
    2828#include <iostream>
    29 // This returns a non const array filled with input from scalar
    30 // or actual array
    31 template <class T> inline T*
    32 copyOfArray( const T * array, const int size, T value)
    33 {
    34   T * arrayNew = new T[size];
    35   if (array)
    36     CoinDisjointCopyN(array,size,arrayNew);
    37   else
    38     CoinFillN ( arrayNew, size,value);
    39   return arrayNew;
    40 }
    41 
    42 // This returns a non const array filled with actual array (or NULL)
    43 template <class T> inline T*
    44 copyOfArray( const T * array, const int size)
    45 {
    46   if (array) {
    47     T * arrayNew = new T[size];
    48     CoinDisjointCopyN(array,size,arrayNew);
    49     return arrayNew;
    50   } else {
    51     return NULL;
    52   }
    53 }
    54 
    5529//#############################################################################
    5630
     
    7246  largestDualError_(0.0),
    7347  largestSolutionError_(0.0),
    74   dualBound_(1.0e7),
     48  dualBound_(1.0e6),
    7549  lower_(NULL),
    7650  rowLowerWork_(NULL),
     
    9771  directionOut_(-1),
    9872  pivotRow_(-1),
    99   status_(NULL),
    10073  dj_(NULL),
    10174  rowReducedCost_(NULL),
     
    11790  savedSolution_(NULL),
    11891  columnScale_(NULL),
    119   scalingFlag_(0),
     92  scalingFlag_(1),
    12093  numberTimesOptimal_(0),
    12194  changeMade_(1),
     
    12598  infeasibilityCost_(1.0e7),
    12699  nonLinearCost_(NULL),
    127   specialOptions_(0)
     100  specialOptions_(0),
     101  lastBadIteration_(-999999),
     102  numberFake_(0),
     103  progressFlag_(0),
     104  sumOfRelaxedDualInfeasibilities_(0.0),
     105  sumOfRelaxedPrimalInfeasibilities_(0.0)
    128106
    129107{
     
    179157    for (iRow=0;iRow<numberRows_;iRow++) {
    180158      int iPivot=pivotVariable_[iRow];
    181       if (iPivot>=numberColumns_) {
    182         // slack
    183         save[iRow] = rowActivityWork_[iPivot-numberColumns_];
    184       } else {
    185         // column
    186         save[iRow] = columnActivityWork_[iPivot];
    187       }
     159      save[iRow] = solution_[iPivot];
    188160    }
    189161  }
     
    200172        <<nonLinearCost_->changeInCost()
    201173        <<nonLinearCost_->numberInfeasibilities()
    202         <<OsiMessageEol;
     174        <<CoinMessageEol;
    203175  }
    204176  if (valuesPass) {
     177#ifdef CLP_DEBUG
    205178    std::cout<<"Largest given infeasibility "<<oldValue
    206179             <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl;
     180#endif
    207181    int numberOut=0;
    208182    if (oldValue<1.0
     
    212186      int iRow;
    213187      int * sort = new int[numberRows_];
     188      // first put back solution and store difference
    214189      for (iRow=0;iRow<numberRows_;iRow++) {
    215190        int iPivot=pivotVariable_[iRow];
     191        double difference = fabs(solution_[iPivot]=save[iRow]);
     192        solution_[iPivot]=save[iRow];
     193        save[iRow]=difference;
     194      }
     195      for (iRow=0;iRow<numberRows_;iRow++) {
     196        int iPivot=pivotVariable_[iRow];
    216197
    217198        if (iPivot<numberColumns_) {
    218199          // column
    219           double difference= fabs(save[iRow] - columnActivityWork_[iPivot]);
     200          double difference= save[iRow];
    220201          if (difference>1.0e-4) {
    221202            sort[numberOut]=iPivot;
     
    229210      for (iRow=0;iRow<numberOut;iRow++) {
    230211        int iColumn=sort[iRow];
    231         setColumnStatus(iColumn,ClpSimplex::superBasic);
    232        
     212        setColumnStatus(iColumn,superBasic);
     213
    233214      }
    234215      delete [] sort;
     
    245226  objectiveValue_ += objectiveModification;
    246227  checkDualSolution();
    247   if (handler_->logLevel()>3||(largestPrimalError_>1.0e-4||
    248                                largestDualError_>1.0e-4))
     228  if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
     229                               largestDualError_>1.0e-2))
    249230    handler_->message(CLP_SIMPLEX_ACCURACY,messages_)
    250231      <<largestPrimalError_
    251232      <<largestDualError_
    252       <<OsiMessageEol;
     233      <<CoinMessageEol;
    253234  return 0;
    254235}
     
    259240
    260241  //work space
    261   OsiIndexedVector  * workSpace = rowArray_[0];
     242  CoinIndexedVector  * workSpace = rowArray_[0];
    262243
    263244  double * array = new double [numberRows_];
     
    266247
    267248  // accumulate non basic stuff
    268   CoinFillN(array,numberRows_,0.0);
     249  ClpFillN(array,numberRows_,0.0);
    269250
    270251  int iRow;
     
    273254  // So zero out basic
    274255  if (columnActivities!=columnActivityWork_)
    275     CoinDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
     256    ClpDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
    276257  if (rowActivities!=rowActivityWork_)
    277     CoinDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
     258    ClpDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
    278259  for (iRow=0;iRow<numberRows_;iRow++) {
    279260    int iPivot=pivotVariable_[iRow];
    280     if (iPivot>=numberColumns_) {
    281       // slack
    282       rowActivityWork_[iPivot-numberColumns_] = 0.0;
    283     } else {
    284       // column
    285       columnActivityWork_[iPivot] = 0.0;
    286     }
     261    solution_[iPivot] = 0.0;
    287262  }
    288263  times(-1.0,columnActivityWork_,array);
     
    300275  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
    301276    // save in case we want to restore
    302     CoinDisjointCopyN ( array, numberRows_ , save);
     277    ClpDisjointCopyN ( array, numberRows_ , save);
    303278
    304279    // put solution in correct place
    305280    for (iRow=0;iRow<numberRows_;iRow++) {
    306281      int iPivot=pivotVariable_[iRow];
    307       if (iPivot>=numberColumns_) {
    308         // slack
    309         rowActivityWork_[iPivot-numberColumns_] = array[iRow];
    310       } else {
    311         // column
    312         columnActivityWork_[iPivot] = array[iRow];
    313       }
     282      solution_[iPivot] = array[iRow];
    314283    }
    315284    // check Ax == b  (for all)
     
    355324
    356325  // solution as accurate as we are going to get
    357   CoinFillN(work,numberRows_,0.0);
     326  ClpFillN(work,numberRows_,0.0);
    358327  // put solution in correct place
    359328  for (iRow=0;iRow<numberRows_;iRow++) {
    360329    int iPivot=pivotVariable_[iRow];
    361     if (iPivot>=numberColumns_) {
    362       // slack
    363       rowActivityWork_[iPivot-numberColumns_] = array[iRow];
    364     } else {
    365       // column
    366       columnActivityWork_[iPivot] = array[iRow];
    367      }
     330    solution_[iPivot] = array[iRow];
    368331  }
    369332  delete [] previous;
     
    377340  double slackValue = factorization_->slackValue();
    378341  //work space
    379   OsiIndexedVector  * workSpace = rowArray_[0];
     342  CoinIndexedVector  * workSpace = rowArray_[0];
    380343
    381344  double * array = new double [numberRows_];
     
    389352  for (iRow=0;iRow<numberRows_;iRow++) {
    390353    int iPivot=pivotVariable_[iRow];
    391     if (iPivot>=numberColumns_) {
    392       // slack
    393       array[iRow] = rowObjectiveWork_[iPivot-numberColumns_];
    394     } else {
    395       // column
    396       array[iRow]=objectiveWork_[iPivot];
    397     }
    398   }
    399   CoinDisjointCopyN ( array, numberRows_ , save);
     354    array[iRow]=cost_[iPivot];
     355  }
     356  ClpDisjointCopyN ( array, numberRows_ , save);
    400357
    401358  // Btran basic costs and get as accurate as possible
     
    408365    largestDualError_=0.0;
    409366    // would be faster to do just for basic but this reduces code
    410     CoinDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
     367    ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    411368    transposeTimes(-1.0,array,reducedCostWork_);
    412369    for (iRow=0;iRow<numberRows_;iRow++) {
     
    453410    }
    454411  }
    455   CoinFillN(work,numberRows_,0.0);
     412  ClpFillN(work,numberRows_,0.0);
    456413  // now look at dual solution
    457414  for (iRow=0;iRow<numberRows_;iRow++) {
     
    462419    rowReducedCost_[iRow]=value;
    463420  }
    464   CoinDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
     421  ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    465422  transposeTimes(-1.0,dual_,reducedCostWork_);
    466423  delete [] previous;
     
    480437    gutsOfSolution ( rowActivities, columnActivities);
    481438    // release extra memory
    482     deleteRim();
     439    deleteRim(false);
    483440  }
    484441  return factorization_->status();
     
    491448  double * rowActivities = new double[numberRows_];
    492449  double * columnActivities = new double[numberColumns_];
    493   CoinDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
    494   CoinDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
     450  ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
     451  ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
    495452  int status = getSolution( rowActivities, columnActivities);
    496453  delete [] rowActivities;
     
    507464  int status = internalFactorize(-1);
    508465  // release extra memory
    509   deleteRim();
     466  deleteRim(false);
    510467
    511468  return status;
     
    516473      If 10 added then in primal values pass
    517474*/
    518 // Return codes are as from ClpFactorization
     475/* Return codes are as from ClpFactorization unless initial factorization
     476   when total number of singularities is returned */
    519477int ClpSimplex::internalFactorize ( int solveType)
    520478{
    521479
    522480  int iRow,iColumn;
     481  int totalSlacks=numberRows_;
    523482
    524483  bool valuesPass=false;
     
    530489    if (!valuesPass) {
    531490      // not values pass so set to bounds
     491      bool allSlack=true;
    532492      if (status_) {
     493        for (iRow=0;iRow<numberRows_;iRow++) {
     494          if (getRowStatus(iRow)!=basic) {
     495            allSlack=false;
     496            break;
     497          }
     498        }
     499      }
     500      if (!allSlack) {
    533501        // set values from warm start (if sensible)
    534502        int numberBasic=0;
     
    536504          switch(getRowStatus(iRow)) {
    537505           
    538           case ClpSimplex::basic:
     506          case basic:
    539507            numberBasic++;
    540508            break;
    541           case ClpSimplex::isFree:
     509          case isFree:
    542510            assert(rowLowerWork_[iRow]<-largeValue_);
    543511            assert(rowUpperWork_[iRow]>largeValue_);
    544512            rowActivityWork_[iRow]=0.0;
    545513            break;
    546           case ClpSimplex::atUpperBound:
     514          case atUpperBound:
    547515            rowActivityWork_[iRow]=rowUpperWork_[iRow];
    548516            if (rowActivityWork_[iRow]>largeValue_) {
    549               assert(rowLowerWork_[iRow]>-largeValue_);
    550               rowActivityWork_[iRow]=rowLowerWork_[iRow];
    551               setRowStatus(iRow,ClpSimplex::atLowerBound);
     517              if (rowLowerWork_[iRow]>-largeValue_) {
     518                rowActivityWork_[iRow]=rowLowerWork_[iRow];
     519                setRowStatus(iRow,atLowerBound);
     520              } else {
     521                // say free
     522                setRowStatus(iRow,isFree);
     523                rowActivityWork_[iRow]=0.0;
     524              }
    552525            }
    553526            break;
    554           case ClpSimplex::atLowerBound:
     527          case atLowerBound:
    555528            rowActivityWork_[iRow]=rowLowerWork_[iRow];
    556529            if (rowActivityWork_[iRow]<-largeValue_) {
    557               assert(rowUpperWork_[iRow]<largeValue_);
    558               rowActivityWork_[iRow]=rowUpperWork_[iRow];
    559               setRowStatus(iRow,ClpSimplex::atUpperBound);
     530              if (rowUpperWork_[iRow]<largeValue_) {
     531                rowActivityWork_[iRow]=rowUpperWork_[iRow];
     532                setRowStatus(iRow,atUpperBound);
     533              } else {
     534                // say free
     535                setRowStatus(iRow,isFree);
     536                rowActivityWork_[iRow]=0.0;
     537              }
    560538            }
    561539            break;
    562           case ClpSimplex::superBasic:
    563             abort();
     540          case superBasic:
     541            if (rowUpperWork_[iRow]>largeValue_) {
     542              if (rowLowerWork_[iRow]>-largeValue_) {
     543                rowActivityWork_[iRow]=rowLowerWork_[iRow];
     544                setRowStatus(iRow,atLowerBound);
     545              } else {
     546                // say free
     547                setRowStatus(iRow,isFree);
     548                rowActivityWork_[iRow]=0.0;
     549              }
     550            } else {
     551              if (rowLowerWork_[iRow]>-largeValue_) {
     552                // set to nearest
     553                if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
     554                    <fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])) {
     555                  rowActivityWork_[iRow]=rowLowerWork_[iRow];
     556                  setRowStatus(iRow,atLowerBound);
     557                } else {
     558                  rowActivityWork_[iRow]=rowUpperWork_[iRow];
     559                  setRowStatus(iRow,atUpperBound);
     560                }
     561              } else {
     562                rowActivityWork_[iRow]=rowUpperWork_[iRow];
     563                setRowStatus(iRow,atUpperBound);
     564              }
     565            }
    564566            break;
    565567          }
    566568        }
     569        totalSlacks=numberBasic;
    567570        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    568571          switch(getColumnStatus(iColumn)) {
    569572           
    570           case ClpSimplex::basic:
     573          case basic:
    571574            if (numberBasic==numberRows_) {
    572575              // take out of basis
     
    575578                    columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
    576579                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    577                   setColumnStatus(iColumn,ClpSimplex::atLowerBound);
     580                  setColumnStatus(iColumn,atLowerBound);
    578581                } else {
    579582                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    580                   setColumnStatus(iColumn,ClpSimplex::atUpperBound);
     583                  setColumnStatus(iColumn,atUpperBound);
    581584                }
    582585              } else if (columnUpperWork_[iColumn]<largeValue_) {
    583586                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    584                 setColumnStatus(iColumn,ClpSimplex::atUpperBound);
     587                setColumnStatus(iColumn,atUpperBound);
    585588              } else {
    586589                columnActivityWork_[iColumn]=0.0;
    587                 setColumnStatus(iColumn,ClpSimplex::isFree);
     590                setColumnStatus(iColumn,isFree);
    588591              }
    589592            } else {
     
    591594            }
    592595            break;
    593           case ClpSimplex::isFree:
     596          case isFree:
    594597            columnActivityWork_[iColumn]=0.0;
    595598            break;
    596           case ClpSimplex::atUpperBound:
     599          case atUpperBound:
    597600            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    598601            if (columnActivityWork_[iColumn]>largeValue_) {
    599602              assert(columnLowerWork_[iColumn]>-largeValue_);
    600603              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    601               setColumnStatus(iColumn,ClpSimplex::atLowerBound);
     604              setColumnStatus(iColumn,atLowerBound);
    602605            }
    603606            break;
    604           case ClpSimplex::atLowerBound:
     607          case atLowerBound:
    605608            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    606609            if (columnActivityWork_[iColumn]<-largeValue_) {
    607610              assert(columnUpperWork_[iColumn]<largeValue_);
    608611              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    609               setColumnStatus(iColumn,ClpSimplex::atUpperBound);
     612              setColumnStatus(iColumn,atUpperBound);
    610613            }
    611614            break;
    612           case ClpSimplex::superBasic:
    613             abort();
     615          case superBasic:
     616            if (columnUpperWork_[iColumn]>largeValue_) {
     617              if (columnLowerWork_[iColumn]>-largeValue_) {
     618                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     619                setColumnStatus(iColumn,atLowerBound);
     620              } else {
     621                // say free
     622                setColumnStatus(iColumn,isFree);
     623                columnActivityWork_[iColumn]=0.0;
     624              }
     625            } else {
     626              if (columnLowerWork_[iColumn]>-largeValue_) {
     627                // set to nearest
     628                if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
     629                    <fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])) {
     630                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     631                  setColumnStatus(iColumn,atLowerBound);
     632                } else {
     633                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     634                  setColumnStatus(iColumn,atUpperBound);
     635                }
     636              } else {
     637                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     638                setColumnStatus(iColumn,atUpperBound);
     639              }
     640            }
    614641            break;
    615642          }
     
    621648        // changed to put free variables in basis
    622649        if (!status_) {
    623           status_ = new unsigned char [numberColumns_+numberRows_];
    624           memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
     650          createStatus();
    625651        }
    626652        for (iRow=0;iRow<numberRows_;iRow++) {
     
    629655          if (lower>-largeValue_||upper<largeValue_) {
    630656            if (fabs(lower)<=fabs(upper)) {
    631               setRowStatus(iRow,ClpSimplex::atLowerBound);
     657              setRowStatus(iRow,atLowerBound);
    632658              rowActivityWork_[iRow]=lower;
    633659            } else {
    634               setRowStatus(iRow,ClpSimplex::atUpperBound);
     660              setRowStatus(iRow,atUpperBound);
    635661              rowActivityWork_[iRow]=upper;
    636662            }
    637663          } else {
    638             setRowStatus(iRow,ClpSimplex::isFree);
     664            setRowStatus(iRow,isFree);
    639665            rowActivityWork_[iRow]=0.0;
    640666          }
    641667#ifdef TESTFREE
    642           setRowStatus(iRow,ClpSimplex::basic);
     668          setRowStatus(iRow,basic);
    643669          numberBasic++;
    644670#endif
     
    650676          if (lower>-big_bound||upper<big_bound) {
    651677            if (fabs(lower)<=fabs(upper)) {
    652               setColumnStatus(iColumn,ClpSimplex::atLowerBound);
     678              setColumnStatus(iColumn,atLowerBound);
    653679              columnActivityWork_[iColumn]=lower;
    654680            } else {
    655               setColumnStatus(iColumn,ClpSimplex::atUpperBound);
     681              setColumnStatus(iColumn,atUpperBound);
    656682              columnActivityWork_[iColumn]=upper;
    657683            }
     
    659685#ifndef TESTFREE
    660686            numberBasic++;
    661             setColumnStatus(iColumn,ClpSimplex::basic);
     687            setColumnStatus(iColumn,basic);
    662688#else
    663             setColumnStatus(iColumn,ClpSimplex::isFree);
     689            setColumnStatus(iColumn,isFree);
    664690#endif
    665691            columnActivityWork_[iColumn]=0.0;
     
    670696          // might as well do all slack basis
    671697          for (iRow=0;iRow<numberRows_;iRow++) {
    672             setRowStatus(iRow,ClpSimplex::basic);
     698            setRowStatus(iRow,basic);
    673699          }
    674700        }
     
    677703      // values pass has less coding
    678704      // make row activities correct
    679       times(-1.0,columnActivityWork_,rowActivityWork_);
     705      memset(rowActivityWork_,0,numberRows_*sizeof(double));
     706      times(1.0,columnActivityWork_,rowActivityWork_);
    680707      if (status_) {
    681708        // set values from warm start (if sensible)
    682709        int numberBasic=0;
    683710        for (iRow=0;iRow<numberRows_;iRow++) {
    684           if (getRowStatus(iRow)==ClpSimplex::basic)
     711          if (getRowStatus(iRow)==basic)
    685712            numberBasic++;
    686           else
    687             setRowStatus(iRow,ClpSimplex::superBasic);
     713          else {
     714            setRowStatus(iRow,superBasic);
     715            // but put to bound if close
     716            if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
     717                <=primalTolerance_)
     718              rowActivityWork_[iRow]=rowLowerWork_[iRow];
     719            else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
     720                <=primalTolerance_)
     721              rowActivityWork_[iRow]=rowUpperWork_[iRow];
     722          }
     723         
    688724        }
     725        totalSlacks=numberBasic;
    689726        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    690           if (getColumnStatus(iColumn)==ClpSimplex::basic) {
    691             if (numberBasic==numberRows_)
     727          if (getColumnStatus(iColumn)==basic) {
     728            if (numberBasic==numberRows_) {
    692729              // take out of basis
    693               setColumnStatus(iColumn,ClpSimplex::superBasic);
    694             else
     730              setColumnStatus(iColumn,superBasic);
     731              // but put to bound if close
     732              if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
     733                  <=primalTolerance_)
     734                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     735              else if (fabs(columnActivityWork_[iColumn]
     736                            -columnUpperWork_[iColumn])
     737                       <=primalTolerance_)
     738                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     739            } else
    695740              numberBasic++;
    696741          } else {
    697             setColumnStatus(iColumn,ClpSimplex::superBasic);
     742            setColumnStatus(iColumn,superBasic);
     743            // but put to bound if close
     744            if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
     745                <=primalTolerance_)
     746              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     747            else if (fabs(columnActivityWork_[iColumn]
     748                          -columnUpperWork_[iColumn])
     749                     <=primalTolerance_)
     750              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    698751          }
    699752        }
     
    702755        int numberBasic=0;
    703756        if (!status_) {
    704           status_ = new unsigned char [numberColumns_+numberRows_];
    705           memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
     757          createStatus();
    706758        }
    707759        for (iRow=0;iRow<numberRows_;iRow++) {
    708           setRowStatus(iRow,ClpSimplex::basic);
     760          setRowStatus(iRow,basic);
    709761          numberBasic++;
    710762        }
    711763        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    712           setColumnStatus(iColumn,ClpSimplex::superBasic);
     764          setColumnStatus(iColumn,superBasic);
     765          // but put to bound if close
     766          if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
     767              <=primalTolerance_)
     768            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
     769          else if (fabs(columnActivityWork_[iColumn]
     770                        -columnUpperWork_[iColumn])
     771                   <=primalTolerance_)
     772            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    713773        }
    714774      }
     
    725785    int numberBasic=0;
    726786    for (i=0;i<numberRows_;i++) {
    727       if (getRowStatus(i) == ClpSimplex::basic) {
     787      if (getRowStatus(i) == basic) {
    728788        rowIsBasic[i]=1;
    729789        numberBasic++;
     
    733793    }
    734794    for (i=0;i<numberColumns_;i++) {
    735       if (getColumnStatus(i) == ClpSimplex::basic) {
     795      if (getColumnStatus(i) == basic) {
    736796        columnIsBasic[i]=1;
    737797        numberBasic++;
     
    754814      // do pivot information
    755815      for (i=0;i<numberRows_;i++) {
    756         if (getRowStatus(i) == ClpSimplex::basic) {
     816        if (getRowStatus(i) == basic) {
    757817          pivotVariable_[rowIsBasic[i]]=i+numberColumns_;
    758818        }
    759819      }
    760820      for (i=0;i<numberColumns_;i++) {
    761         if (getColumnStatus(i) == ClpSimplex::basic) {
     821        if (getColumnStatus(i) == basic) {
    762822          pivotVariable_[columnIsBasic[i]]=i;
    763823        }
     
    769829      }
    770830      for (i=0;i<numberRows_;i++) {
    771         if (getRowStatus(i) == ClpSimplex::basic) {
     831        if (getRowStatus(i) == basic) {
    772832          int iPivot = rowIsBasic[i];
    773833          if (iPivot>=0)
     
    776836      }
    777837      for (i=0;i<numberColumns_;i++) {
    778         if (getColumnStatus(i) == ClpSimplex::basic) {
     838        if (getColumnStatus(i) == basic) {
    779839          int iPivot = columnIsBasic[i];
    780840          if (iPivot>=0)
     
    788848        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    789849          if (getColumnStatus(iColumn)==
    790               ClpSimplex::basic) {
     850              basic) {
    791851            // take out
    792             double lower=columnLowerWork_[iColumn];
    793             double upper=columnUpperWork_[iColumn];
    794             double value=columnActivityWork_[iColumn];
    795             if (lower>-largeValue_||upper<largeValue_) {
    796               if (fabs(value-lower)<fabs(value-upper)) {
    797                 setColumnStatus(iColumn,ClpSimplex::atLowerBound);
    798                 columnActivityWork_[iColumn]=lower;
     852            if (!valuesPass) {
     853              double lower=columnLowerWork_[iColumn];
     854              double upper=columnUpperWork_[iColumn];
     855              double value=columnActivityWork_[iColumn];
     856              if (lower>-largeValue_||upper<largeValue_) {
     857                if (fabs(value-lower)<fabs(value-upper)) {
     858                  setColumnStatus(iColumn,atLowerBound);
     859                  columnActivityWork_[iColumn]=lower;
     860                } else {
     861                  setColumnStatus(iColumn,atUpperBound);
     862                  columnActivityWork_[iColumn]=upper;
     863                }
    799864              } else {
    800                 setColumnStatus(iColumn,ClpSimplex::atUpperBound);
    801                 columnActivityWork_[iColumn]=upper;
     865                setColumnStatus(iColumn,isFree);
    802866              }
    803867            } else {
    804               setColumnStatus(iColumn,ClpSimplex::isFree);
     868              setColumnStatus(iColumn,superBasic);
    805869            }
    806870          }
     
    815879            } else {
    816880              // put back structural
    817               setColumnStatus(iSequence,ClpSimplex::basic);
     881              setColumnStatus(iSequence,basic);
    818882            }
    819883          } else {
    820884            // put in slack
    821             setRowStatus(iRow,ClpSimplex::basic);
     885            setRowStatus(iRow,basic);
    822886          }
    823887        }
     
    832896    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
    833897      <<status
    834       <<OsiMessageEol;
     898      <<CoinMessageEol;
    835899    return -1;
     900  } else if (!solveType) {
     901    // Initial basis - return number of singularities
     902    int numberSlacks=0;
     903    for (iRow=0;iRow<numberRows_;iRow++) {
     904      if (getRowStatus(iRow) == basic)
     905        numberSlacks++;
     906    }
     907    status= max(numberSlacks-totalSlacks,0);
    836908  }
    837909
     
    860932    <<directionIn_<<theta_
    861933    <<dualOut_<<dualIn_<<alpha_
    862     <<OsiMessageEol;
    863   if (getStatus(sequenceIn_)==ClpSimplex::isFree) {
     934    <<CoinMessageEol;
     935  if (getStatus(sequenceIn_)==isFree) {
    864936    handler_->message(CLP_SIMPLEX_FREEIN,messages_)
    865937      <<sequenceIn_
    866       <<OsiMessageEol;
     938      <<CoinMessageEol;
    867939  }
    868940  // change of incoming
     
    880952    if (sequenceIn_!=sequenceOut_) {
    881953      // variable becoming basic
    882       setStatus(sequenceIn_,ClpSimplex::basic);
     954      setStatus(sequenceIn_,basic);
    883955      if (algorithm_<0) {
    884956        value = upperIn_+dualOut_;
     
    889961    } else {
    890962      value=lowerIn_;
    891       setStatus(sequenceIn_, ClpSimplex::atLowerBound);
     963      setStatus(sequenceIn_, atLowerBound);
    892964    }
    893965  } else {
     
    895967    if (sequenceIn_!=sequenceOut_) {
    896968      // variable becoming basic
    897       setStatus(sequenceIn_,ClpSimplex::basic);
     969      setStatus(sequenceIn_,basic);
    898970      if (algorithm_<0) {
    899971        value = lowerIn_+dualOut_;
     
    904976    } else {
    905977      value=upperIn_;
    906       setStatus(sequenceIn_, ClpSimplex::atUpperBound);
     978      setStatus(sequenceIn_, atUpperBound);
    907979    }
    908980  }
     
    915987  // outgoing
    916988  if (sequenceIn_!=sequenceOut_) {
    917     assert( getStatus(sequenceOut_)== ClpSimplex::basic);
     989    assert( getStatus(sequenceOut_)== basic);
     990    if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
     991      progressFlag_ |= 1; // making real progress
    918992    if (algorithm_<0) {
    919993      if (directionOut_>0) {
    920994        value = lowerOut_;
    921         setStatus(sequenceOut_,ClpSimplex::atLowerBound);
     995        setStatus(sequenceOut_,atLowerBound);
    922996        dj_[sequenceOut_]=theta_;
    923997      } else {
    924998        value = upperOut_;
    925         setStatus(sequenceOut_,ClpSimplex::atUpperBound);
     999        setStatus(sequenceOut_,atUpperBound);
    9261000        dj_[sequenceOut_]=-theta_;
    9271001      }
     
    9331007        value = upperOut_;
    9341008      }
    935       nonLinearCost_->setOne(sequenceOut_,value);
    9361009      double lowerValue = lower_[sequenceOut_];
    9371010      double upperValue = upper_[sequenceOut_];
     
    9401013      // may not be exactly at bound and bounds may have changed
    9411014      if (value<=lowerValue+primalTolerance_) {
    942         setStatus(sequenceOut_,ClpSimplex::atLowerBound);
     1015        setStatus(sequenceOut_,atLowerBound);
    9431016      } else if (value>=upperValue-primalTolerance_) {
    944         setStatus(sequenceOut_,ClpSimplex::atUpperBound);
     1017        setStatus(sequenceOut_,atUpperBound);
    9451018      } else {
    9461019        printf("*** variable wandered off bound %g %g %g!\n",
    9471020               lowerValue,value,upperValue);
    9481021        if (value-lowerValue<=upperValue-value) {
    949           setStatus(sequenceOut_,ClpSimplex::atLowerBound);
     1022          setStatus(sequenceOut_,atLowerBound);
    9501023          value=lowerValue;
    9511024        } else {
    952           setStatus(sequenceOut_,ClpSimplex::atUpperBound);
     1025          setStatus(sequenceOut_,atUpperBound);
    9531026          value=upperValue;
    9541027        }
    9551028      }
    9561029      solution_[sequenceOut_]=value;
     1030      nonLinearCost_->setOne(sequenceOut_,value);
    9571031    }
    9581032  }
     
    9671041  handler_->printing(algorithm_<0)<<theta_<<dualOut_;
    9681042  handler_->printing(algorithm_>0)<<dualIn_<<theta_;
    969   handler_->message()<<OsiMessageEol;
    970 
     1043  handler_->message()<<CoinMessageEol;
    9711044  if (numberIterations_>=maximumIterations_)
    9721045    return 2;
     
    10311104  directionOut_(-1),
    10321105  pivotRow_(-1),
    1033   status_(NULL),
    10341106  dj_(NULL),
    10351107  rowReducedCost_(NULL),
     
    10591131  infeasibilityCost_(1.0e7),
    10601132  nonLinearCost_(NULL),
    1061   specialOptions_(0)
     1133  specialOptions_(0),
     1134  lastBadIteration_(-999999),
     1135  numberFake_(0),
     1136  progressFlag_(0),
     1137  sumOfRelaxedDualInfeasibilities_(0.0),
     1138  sumOfRelaxedPrimalInfeasibilities_(0.0)
    10621139{
    10631140  int i;
     
    11151192  directionOut_(-1),
    11161193  pivotRow_(-1),
    1117   status_(NULL),
    11181194  dj_(NULL),
    11191195  rowReducedCost_(NULL),
     
    11431219  infeasibilityCost_(1.0e7),
    11441220  nonLinearCost_(NULL),
    1145   specialOptions_(0)
     1221  specialOptions_(0),
     1222  lastBadIteration_(-999999),
     1223  numberFake_(0),
     1224  progressFlag_(0),
     1225  sumOfRelaxedDualInfeasibilities_(0.0),
     1226  sumOfRelaxedPrimalInfeasibilities_(0.0)
    11461227{
    11471228  int i;
     
    11731254ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
    11741255{
    1175   lower_ = copyOfArray(rhs.lower_,numberColumns_+numberRows_);
     1256  lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows_);
    11761257  rowLowerWork_ = lower_+numberColumns_;
    11771258  columnLowerWork_ = lower_;
    1178   upper_ = copyOfArray(rhs.upper_,numberColumns_+numberRows_);
     1259  upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows_);
    11791260  rowUpperWork_ = upper_+numberColumns_;
    11801261  columnUpperWork_ = upper_;
    1181   cost_ = copyOfArray(rhs.cost_,(numberColumns_+numberRows_));
     1262  cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows_));
    11821263  objectiveWork_ = cost_;
    11831264  rowObjectiveWork_ = cost_+numberColumns_;
    1184   dj_ = copyOfArray(rhs.dj_,numberRows_+numberColumns_);
     1265  dj_ = ClpCopyOfArray(rhs.dj_,numberRows_+numberColumns_);
    11851266  if (dj_) {
    11861267    reducedCostWork_ = dj_;
    11871268    rowReducedCost_ = dj_+numberColumns_;
    11881269  }
    1189   solution_ = copyOfArray(rhs.solution_,numberRows_+numberColumns_);
     1270  solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
    11901271  if (solution_) {
    11911272    columnActivityWork_ = solution_;
     
    11941275  if (rhs.pivotVariable_) {
    11951276    pivotVariable_ = new int[numberRows_];
    1196     CoinDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
     1277    ClpDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
    11971278  } else {
    11981279    pivotVariable_=NULL;
     
    12031284    factorization_=NULL;
    12041285  }
    1205   rowScale_ = copyOfArray(rhs.rowScale_,numberRows_);
    1206   savedSolution_ = copyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
    1207   columnScale_ = copyOfArray(rhs.columnScale_,numberColumns_);
     1286  rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_);
     1287  savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
     1288  columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_);
    12081289  int i;
    12091290  for (i=0;i<6;i++) {
    12101291    rowArray_[i]=NULL;
    12111292    if (rhs.rowArray_[i])
    1212       rowArray_[i] = new OsiIndexedVector(*rhs.rowArray_[i]);
     1293      rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
    12131294    columnArray_[i]=NULL;
    12141295    if (rhs.columnArray_[i])
    1215       columnArray_[i] = new OsiIndexedVector(*rhs.columnArray_[i]);
    1216   }
    1217   if (rhs.status_) {
    1218     status_ = copyOfArray( rhs.status_,numberColumns_+numberRows_);
     1296      columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
    12191297  }
    12201298  if (rhs.saveStatus_) {
    1221     saveStatus_ = copyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
     1299    saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
    12221300  }
    12231301  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
     
    12701348  infeasibilityCost_ = rhs.infeasibilityCost_;
    12711349  specialOptions_ = rhs.specialOptions_;
     1350  lastBadIteration_ = rhs.lastBadIteration_;
     1351  numberFake_ = rhs.numberFake_;
     1352  progressFlag_ = rhs.progressFlag_;
     1353  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
     1354  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    12721355  if (rhs.nonLinearCost_!=NULL)
    12731356    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
    12741357}
    1275 // type == 0 do everything
     1358// type == 0 do everything, most + pivot data, 2 factorization data as well
    12761359void
    12771360ClpSimplex::gutsOfDelete(int type)
     
    13141397  delete rowCopy_;
    13151398  rowCopy_=NULL;
     1399  delete [] saveStatus_;
     1400  saveStatus_=NULL;
    13161401  if (!type) {
    13171402    // delete everything
     1403    delete factorization_;
     1404    factorization_ = NULL;
    13181405    delete [] pivotVariable_;
    13191406    pivotVariable_=NULL;
    1320     delete factorization_;
    1321     factorization_ = NULL;
    13221407    delete dualRowPivot_;
    13231408    dualRowPivot_ = NULL;
    13241409    delete primalColumnPivot_;
    13251410    primalColumnPivot_ = NULL;
    1326     delete [] saveStatus_;
    1327     saveStatus_=NULL;
    1328     delete status_;
    1329     status_=NULL;
     1411  } else {
     1412    // delete any size information in methods
     1413    if (type>1) {
     1414      factorization_->clearArrays();
     1415      delete [] pivotVariable_;
     1416      pivotVariable_=NULL;
     1417    }
     1418    dualRowPivot_->clearArrays();
     1419    primalColumnPivot_->clearArrays();
    13301420  }
    13311421}
     
    13481438  sumPrimalInfeasibilities_=0.0;
    13491439  numberPrimalInfeasibilities_=0;
     1440  double primalTolerance = primalTolerance_;
     1441  double relaxedTolerance=dualTolerance_;
     1442  // we can't really trust infeasibilities if there is primal error
     1443  double error = min(1.0e-3,largestPrimalError_);
     1444  relaxedTolerance = max(relaxedTolerance, error);
     1445  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
     1446
    13501447  for (iRow=0;iRow<numberRows_;iRow++) {
    13511448    double infeasibility=0.0;
     
    13561453      infeasibility=rowLowerWork_[iRow]-solution[iRow];
    13571454    }
    1358     if (infeasibility>primalTolerance_) {
     1455    if (infeasibility>primalTolerance) {
    13591456      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
     1457      if (infeasibility>relaxedTolerance)
     1458        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
    13601459      numberPrimalInfeasibilities_ ++;
    13611460    }
     
    13851484      columnPrimalSequence_=iColumn;
    13861485    }
    1387     if (infeasibility>primalTolerance_) {
     1486    if (infeasibility>primalTolerance) {
    13881487      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
     1488      if (infeasibility>relaxedTolerance)
     1489        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
    13891490      numberPrimalInfeasibilities_ ++;
    13901491    }
     
    13921493    if (infeasibility>largestSolutionError_)
    13931494      largestSolutionError_=infeasibility;
    1394     if (columnLowerWork_[iColumn]!=columnUpperWork_[iColumn])
     1495    if (columnUpperWork_[iColumn]-columnLowerWork_[iColumn]>primalTolerance_)
    13951496      clearFixed(iColumn);
    13961497    else
     
    14151516  remainingDualInfeasibility_=0.0;
    14161517  solution = rowActivityWork_;
    1417   double dualTolerance=dualTolerance_;
    1418   if (algorithm_>0) {
    1419     // primal
    1420     // we can't really trust infeasibilities if there is dual error
    1421     if (largestDualError_>1.0e-6)
    1422       dualTolerance *= largestDualError_/1.0e-6;
    1423   }
     1518  double relaxedTolerance=dualTolerance_;
     1519  // we can't really trust infeasibilities if there is dual error
     1520  double error = min(1.0e-3,largestDualError_);
     1521  relaxedTolerance = max(relaxedTolerance, error);
     1522  sumOfRelaxedDualInfeasibilities_ = 0.0;
     1523
    14241524  for (iRow=0;iRow<numberRows_;iRow++) {
    1425     if (getRowStatus(iRow) != ClpSimplex::basic) {
     1525    if (getRowStatus(iRow) != basic) {
    14261526      // not basic
    14271527      double value = rowReducedCost_[iRow];
     
    14361536            rowDualSequence_=iRow;
    14371537          }
    1438           if (value>dualTolerance) {
    1439             sumDualInfeasibilities_ += value-dualTolerance;
     1538          if (value>dualTolerance_) {
     1539            sumDualInfeasibilities_ += value-dualTolerance_;
     1540            if (value>relaxedTolerance)
     1541              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    14401542            numberDualInfeasibilities_ ++;
    1441             if (getRowStatus(iRow) != ClpSimplex::isFree)
     1543            if (getRowStatus(iRow) != isFree)
    14421544              numberDualInfeasibilitiesWithoutFree_ ++;
    14431545            // maybe we can make feasible by increasing tolerance
     
    14611563            rowDualSequence_=iRow;
    14621564          }
    1463           if (value>dualTolerance) {
    1464             sumDualInfeasibilities_ += value-dualTolerance;
     1565          if (value>dualTolerance_) {
     1566            sumDualInfeasibilities_ += value-dualTolerance_;
     1567            if (value>relaxedTolerance)
     1568              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    14651569            numberDualInfeasibilities_ ++;
    1466             if (getRowStatus(iRow) != ClpSimplex::isFree)
     1570            if (getRowStatus(iRow) != isFree)
    14671571              numberDualInfeasibilitiesWithoutFree_ ++;
    14681572            // maybe we can make feasible by increasing tolerance
     
    14771581  solution = columnActivityWork_;
    14781582  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1479     if (getColumnStatus(iColumn) != ClpSimplex::basic) {
     1583    if (getColumnStatus(iColumn) != basic) {
    14801584      // not basic
    14811585      double value = reducedCostWork_[iColumn];
     
    14901594            columnDualSequence_=iColumn;
    14911595          }
    1492           if (value>dualTolerance) {
    1493             sumDualInfeasibilities_ += value-dualTolerance;
     1596          if (value>dualTolerance_) {
     1597            sumDualInfeasibilities_ += value-dualTolerance_;
     1598            if (value>relaxedTolerance)
     1599              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    14941600            numberDualInfeasibilities_ ++;
    1495             if (getColumnStatus(iColumn) != ClpSimplex::isFree)
     1601            if (getColumnStatus(iColumn) != isFree)
    14961602              numberDualInfeasibilitiesWithoutFree_ ++;
    14971603            // maybe we can make feasible by increasing tolerance
     
    15151621            columnDualSequence_=iColumn;
    15161622          }
    1517           if (value>dualTolerance) {
    1518             sumDualInfeasibilities_ += value-dualTolerance;
     1623          if (value>dualTolerance_) {
     1624            sumDualInfeasibilities_ += value-dualTolerance_;
     1625            if (value>relaxedTolerance)
     1626              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    15191627            numberDualInfeasibilities_ ++;
    1520             if (getColumnStatus(iColumn) != ClpSimplex::isFree)
     1628            if (getColumnStatus(iColumn) != isFree)
    15211629              numberDualInfeasibilitiesWithoutFree_ ++;
    15221630            // maybe we can make feasible by increasing tolerance
     
    15341642*/
    15351643void
    1536 ClpSimplex::unpack(OsiIndexedVector * rowArray)
     1644ClpSimplex::unpack(CoinIndexedVector * rowArray)
    15371645{
    15381646  rowArray->clear();
     
    15461654}
    15471655void
    1548 ClpSimplex::unpack(OsiIndexedVector * rowArray,int sequence)
     1656ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence)
    15491657{
    15501658  rowArray->clear();
     
    15571665  }
    15581666}
    1559 void
     1667bool
    15601668ClpSimplex::createRim(int what,bool makeRowCopy)
    15611669{
     1670  bool goodMatrix=true;
     1671  int i;
    15621672  if ((what&16)!=0) {
    15631673    // move information to work arrays
     1674    if (optimizationDirection_<0.0) {
     1675      // reverse all dual signs
     1676      for (i=0;i<numberColumns_;i++)
     1677        reducedCost_[i] = -reducedCost_[i];
     1678      for (i=0;i<numberRows_;i++)
     1679        dual_[i] = -dual_[i];
     1680    }
    15641681    // row reduced costs
    15651682    if (!dj_) {
     
    15681685      rowReducedCost_ = dj_+numberColumns_;
    15691686      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
     1687      memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double));
    15701688    }
    15711689    if (!solution_) {
     
    15781696             numberRows_*sizeof(double));
    15791697    }
    1580    
     1698    //check matrix
     1699    if (!matrix_->allElementsInRange(this,1.0e-20,1.0e20)) {
     1700      problemStatus_=4;
     1701      goodMatrix= false;
     1702    }
    15811703    if (makeRowCopy) {
    15821704      delete rowCopy_;
     
    15851707    }
    15861708  }
    1587   int i;
    15881709  if ((what&4)!=0) {
    15891710    delete [] cost_;
     
    16751796      rowActivityWork_[i] *= rowScale_[i];
    16761797      dual_[i] /= rowScale_[i];
     1798      rowReducedCost_[i] = dual_[i];
    16771799    }
    16781800  }
    16791801 
     1802  if ((what&16)!=0) {
     1803    // check rim of problem okay
     1804    if (!sanityCheck())
     1805      goodMatrix=false;
     1806  }
    16801807  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
    16811808  // maybe we need to move scales to SimplexModel for factorization?
     
    16931820    for (iRow=0;iRow<4;iRow++) {
    16941821      if (!rowArray_[iRow]) {
    1695         rowArray_[iRow]=new OsiIndexedVector();
     1822        rowArray_[iRow]=new CoinIndexedVector();
    16961823        int length =numberRows_+factorization_->maximumPivots();
    16971824        if (iRow==3)
     
    17031830    for (iColumn=0;iColumn<2;iColumn++) {
    17041831      if (!columnArray_[iColumn]) {
    1705         columnArray_[iColumn]=new OsiIndexedVector();
     1832        columnArray_[iColumn]=new CoinIndexedVector();
    17061833        columnArray_[iColumn]->reserve(numberColumns_);
    17071834      }
    17081835    }   
    17091836  }
    1710 
     1837  return goodMatrix;
    17111838}
    17121839void
    1713 ClpSimplex::deleteRim()
     1840ClpSimplex::deleteRim(bool getRidOfFactorizationData)
    17141841{
    17151842  int i;
     
    17461873  }
    17471874  // direction may have been modified by scaling - clean up
    1748   if (optimizationDirection_>0.0)
     1875  if (optimizationDirection_>0.0) {
    17491876    optimizationDirection_ = 1.0;
    1750   else if (optimizationDirection_<0.0)
     1877  }  else if (optimizationDirection_<0.0) {
    17511878    optimizationDirection_ = -1.0;
     1879    // and reverse all dual signs
     1880    for (i=0;i<numberColumns_;i++)
     1881      reducedCost_[i] = -reducedCost_[i];
     1882    for (i=0;i<numberRows_;i++)
     1883      dual_[i] = -dual_[i];
     1884  }
    17521885  // scaling may have been turned off
    17531886  scalingFlag_ = abs(scalingFlag_);
    1754   gutsOfDelete(1);
     1887  if(getRidOfFactorizationData)
     1888    gutsOfDelete(2);
     1889  else
     1890    gutsOfDelete(1);
    17551891}
    17561892void
     
    17991935  }
    18001936}
    1801 // Sets up basis
    1802 void
    1803 ClpSimplex::setBasis ( const OsiWarmStartBasis & basis)
    1804 {
    1805   if (basis.getNumStructural()) {
    1806     // transform basis to status arrays
    1807     int iRow,iColumn;
    1808     if (!status_) {
    1809       /*
    1810         get status arrays
    1811         OsiWarmStartBasis would seem to have overheads and we will need
    1812         extra bits anyway.
    1813       */
    1814       status_ = new unsigned char [numberColumns_+numberRows_];
    1815       memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
    1816     }
    1817     OsiWarmStartBasis basis2 = basis;
    1818     // resize if necessary
    1819     basis2.resize(numberRows_,numberColumns_);
    1820     // move status
    1821     for (iRow=0;iRow<numberRows_;iRow++) {
    1822       setRowStatus(iRow,
    1823                    (Status) basis2.getArtifStatus(iRow));
    1824     }
    1825     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1826       setColumnStatus(iColumn,
    1827                       (Status) basis2.getStructStatus(iColumn));
    1828     }
    1829   } else {
    1830     // null basis so just delete any status arrays
    1831     delete [] status_;
    1832     status_=NULL;
    1833   }
    1834 }
    18351937// Passes in factorization
    18361938void
     
    18391941  delete factorization_;
    18401942  factorization_= new ClpFactorization(factorization);
    1841 }
    1842 // Warm start
    1843 OsiWarmStartBasis 
    1844 ClpSimplex::getBasis() const
    1845 {
    1846   int iRow,iColumn;
    1847   OsiWarmStartBasis basis;
    1848   basis.setSize(numberColumns_,numberRows_);
    1849 
    1850   if(status_) {
    1851     for (iRow=0;iRow<numberRows_;iRow++) {
    1852       basis.setArtifStatus(iRow,
    1853                            (OsiWarmStartBasis::Status) getRowStatus(iRow));
    1854     }
    1855     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1856       basis.setStructStatus(iColumn,
    1857                        (OsiWarmStartBasis::Status) getColumnStatus(iColumn));
    1858     }
    1859   }
    1860   return basis;
    18611943}
    18621944void
     
    19182000 
    19192001  // Get a row copy in standard format
    1920   OsiPackedMatrix copy;
     2002  CoinPackedMatrix copy;
    19212003  copy.reverseOrderedCopyOf(*matrix());
    19222004  // get matrix data pointers
    19232005  const int * column = copy.getIndices();
    1924   const int * rowStart = copy.getVectorStarts();
     2006  const CoinBigIndex * rowStart = copy.getVectorStarts();
    19252007  const int * rowLength = copy.getVectorLengths();
    19262008  const double * element = copy.getElements();
     
    19632045        double maximumDown = 0.0;
    19642046        double newBound;
    1965         int rStart = rowStart[iRow];
    1966         int rEnd = rowStart[iRow]+rowLength[iRow];
    1967         int j;
     2047        CoinBigIndex rStart = rowStart[iRow];
     2048        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
     2049        CoinBigIndex j;
    19682050
    19692051        // Compute possible lower and upper ranges
     
    21062188    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
    21072189      <<totalTightened
    2108       <<OsiMessageEol;
     2190      <<CoinMessageEol;
    21092191    // Set bounds slightly loose
    21102192    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     
    21312213    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
    21322214      <<numberInfeasible
    2133       <<OsiMessageEol;
     2215      <<CoinMessageEol;
    21342216    // restore column bounds
    21352217    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
     
    21522234  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
    21532235}
     2236/* For strong branching.  On input lower and upper are new bounds
     2237   while on output they are objective function values (>1.0e50 infeasible).
     2238   Return code is 0 if nothing interesting, -1 if infeasible both
     2239   ways and +1 if infeasible one way (check values to see which one(s))
     2240*/
     2241int ClpSimplex::strongBranching(int numberVariables,const int * variables,
     2242                                double * newLower, double * newUpper,
     2243                                bool stopOnFirstInfeasible,
     2244                                bool alwaysFinish)
     2245{
     2246  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
     2247                                                    newLower,  newUpper,
     2248                                                    stopOnFirstInfeasible,
     2249                                                    alwaysFinish);
     2250}
     2251/* Borrow model.  This is so we dont have to copy large amounts
     2252   of data around.  It assumes a derived class wants to overwrite
     2253   an empty model with a real one - while it does an algorithm.
     2254   This is same as ClpModel one, but sets scaling on etc. */
     2255void
     2256ClpSimplex::borrowModel(ClpModel & otherModel)
     2257{
     2258  ClpModel::borrowModel(otherModel);
     2259  createStatus();
     2260  scaling();
     2261  ClpDualRowSteepest steep1;
     2262  setDualRowPivotAlgorithm(steep1);
     2263  ClpPrimalColumnSteepest steep2;
     2264  setPrimalColumnPivotAlgorithm(steep2);
     2265}
     2266typedef struct {
     2267  double optimizationDirection;
     2268  double dblParam[ClpLastDblParam];
     2269  double objectiveValue;
     2270  double dualBound;
     2271  double dualTolerance;
     2272  double primalTolerance;
     2273  double sumDualInfeasibilities;
     2274  double sumPrimalInfeasibilities;
     2275  double infeasibilityCost;
     2276  int numberRows;
     2277  int numberColumns;
     2278  int intParam[ClpLastIntParam];
     2279  int numberIterations;
     2280  int problemStatus;
     2281  int maximumIterations;
     2282  int lengthNames;
     2283  int numberDualInfeasibilities;
     2284  int numberDualInfeasibilitiesWithoutFree;
     2285  int numberPrimalInfeasibilities;
     2286  int numberRefinements;
     2287  int scalingFlag;
     2288  int algorithm;
     2289  unsigned int specialOptions;
     2290  int dualPivotChoice;
     2291  int primalPivotChoice;
     2292  int matrixStorageChoice;
     2293} Clp_scalars;
     2294
     2295int outDoubleArray(double * array, int length, FILE * fp)
     2296{
     2297  int numberWritten;
     2298  if (array&&length) {
     2299    numberWritten = fwrite(&length,sizeof(int),1,fp);
     2300    if (numberWritten!=1)
     2301      return 1;
     2302    numberWritten = fwrite(array,sizeof(double),length,fp);
     2303    if (numberWritten!=length)
     2304      return 1;
     2305  } else {
     2306    length = 0;
     2307    numberWritten = fwrite(&length,sizeof(int),1,fp);
     2308    if (numberWritten!=1)
     2309      return 1;
     2310  }
     2311  return 0;
     2312}
     2313// Save model to file, returns 0 if success
     2314int
     2315ClpSimplex::saveModel(const char * fileName)
     2316{
     2317  FILE * fp = fopen(fileName,"wb");
     2318  if (fp) {
     2319    Clp_scalars scalars;
     2320    int i;
     2321    CoinBigIndex numberWritten;
     2322    // Fill in scalars
     2323    scalars.optimizationDirection = optimizationDirection_;
     2324    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
     2325    scalars.objectiveValue = objectiveValue_;
     2326    scalars.dualBound = dualBound_;
     2327    scalars.dualTolerance = dualTolerance_;
     2328    scalars.primalTolerance = primalTolerance_;
     2329    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
     2330    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
     2331    scalars.infeasibilityCost = infeasibilityCost_;
     2332    scalars.numberRows = numberRows_;
     2333    scalars.numberColumns = numberColumns_;
     2334    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
     2335    scalars.numberIterations = numberIterations_;
     2336    scalars.problemStatus = problemStatus_;
     2337    scalars.maximumIterations = maximumIterations_;
     2338    scalars.lengthNames = lengthNames_;
     2339    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
     2340    scalars.numberDualInfeasibilitiesWithoutFree
     2341      = numberDualInfeasibilitiesWithoutFree_;
     2342    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
     2343    scalars.numberRefinements = numberRefinements_;
     2344    scalars.scalingFlag = scalingFlag_;
     2345    scalars.algorithm = algorithm_;
     2346    scalars.specialOptions = specialOptions_;
     2347    scalars.dualPivotChoice = dualRowPivot_->type();
     2348    scalars.primalPivotChoice = primalColumnPivot_->type();
     2349    scalars.matrixStorageChoice = matrix_->type();
     2350
     2351    // put out scalars
     2352    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
     2353    if (numberWritten!=1)
     2354      return 1;
     2355    // strings
     2356    CoinBigIndex length;
     2357    for (i=0;i<ClpLastStrParam;i++) {
     2358      length = strParam_[i].size();
     2359      numberWritten = fwrite(&length,sizeof(int),1,fp);
     2360      if (numberWritten!=1)
     2361        return 1;
     2362      numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
     2363      if (numberWritten!=1)
     2364        return 1;
     2365    }
     2366    // arrays - in no particular order
     2367    if (outDoubleArray(rowActivity_,numberRows_,fp))
     2368        return 1;
     2369    if (outDoubleArray(columnActivity_,numberColumns_,fp))
     2370        return 1;
     2371    if (outDoubleArray(dual_,numberRows_,fp))
     2372        return 1;
     2373    if (outDoubleArray(reducedCost_,numberColumns_,fp))
     2374        return 1;
     2375    if (outDoubleArray(rowLower_,numberRows_,fp))
     2376        return 1;
     2377    if (outDoubleArray(rowUpper_,numberRows_,fp))
     2378        return 1;
     2379    if (outDoubleArray(objective_,numberColumns_,fp))
     2380        return 1;
     2381    if (outDoubleArray(rowObjective_,numberRows_,fp))
     2382        return 1;
     2383    if (outDoubleArray(columnLower_,numberColumns_,fp))
     2384        return 1;
     2385    if (outDoubleArray(columnUpper_,numberColumns_,fp))
     2386        return 1;
     2387    if (ray_) {
     2388      if (problemStatus_==1)
     2389        if (outDoubleArray(ray_,numberRows_,fp))
     2390          return 1;
     2391      else if (problemStatus_==2)
     2392        if (outDoubleArray(ray_,numberColumns_,fp))
     2393          return 1;
     2394      else
     2395        if (outDoubleArray(NULL,0,fp))
     2396          return 1;
     2397    } else {
     2398      if (outDoubleArray(NULL,0,fp))
     2399        return 1;
     2400    }
     2401    if (status_&&(numberRows_+numberColumns_)>0) {
     2402      length = numberRows_+numberColumns_;
     2403      numberWritten = fwrite(&length,sizeof(int),1,fp);
     2404      if (numberWritten!=1)
     2405        return 1;
     2406      numberWritten = fwrite(status_,sizeof(char),length, fp);
     2407      if (numberWritten!=length)
     2408        return 1;
     2409    } else {
     2410      length = 0;
     2411      numberWritten = fwrite(&length,sizeof(int),1,fp);
     2412      if (numberWritten!=1)
     2413        return 1;
     2414    }
     2415    if (lengthNames_) {
     2416      assert (numberRows_ == (int) rowNames_.size());
     2417      for (i=0;i<numberRows_;i++) {
     2418        length = rowNames_[i].size();
     2419        numberWritten = fwrite(&length,sizeof(int),1,fp);
     2420        if (numberWritten!=1)
     2421          return 1;
     2422        numberWritten = fwrite(rowNames_[i].c_str(),length,1,fp);
     2423        if (numberWritten!=1)
     2424          return 1;
     2425      }
     2426      assert (numberColumns_ == (int) columnNames_.size());
     2427      for (i=0;i<numberColumns_;i++) {
     2428        length = columnNames_[i].size();
     2429        numberWritten = fwrite(&length,sizeof(int),1,fp);
     2430        if (numberWritten!=1)
     2431          return 1;
     2432        numberWritten = fwrite(columnNames_[i].c_str(),length,1,fp);
     2433        if (numberWritten!=1)
     2434          return 1;
     2435      }
     2436    }
     2437    // just standard type at present
     2438    assert (matrix_->type()==1);
     2439    assert (matrix_->getNumCols() == numberColumns_);
     2440    assert (matrix_->getNumRows() == numberRows_);
     2441    // we are going to save with gaps
     2442    length = matrix_->getVectorStarts()[numberColumns_-1]
     2443      + matrix_->getVectorLengths()[numberColumns_-1];
     2444    numberWritten = fwrite(&length,sizeof(int),1,fp);
     2445    if (numberWritten!=1)
     2446      return 1;
     2447    numberWritten = fwrite(matrix_->getElements(),
     2448                           sizeof(double),length,fp);
     2449    if (numberWritten!=length)
     2450      return 1;
     2451    numberWritten = fwrite(matrix_->getIndices(),
     2452                           sizeof(int),length,fp);
     2453    if (numberWritten!=length)
     2454      return 1;
     2455    numberWritten = fwrite(matrix_->getVectorStarts(),
     2456                           sizeof(int),numberColumns_,fp);
     2457    if (numberWritten!=numberColumns_)
     2458      return 1;
     2459    numberWritten = fwrite(matrix_->getVectorLengths(),
     2460                           sizeof(int),numberColumns_,fp);
     2461    if (numberWritten!=numberColumns_)
     2462      return 1;
     2463    // finished
     2464    fclose(fp);
     2465    return 0;
     2466  } else {
     2467    return -1;
     2468  }
     2469}
     2470
     2471int inDoubleArray(double * &array, int length, FILE * fp)
     2472{
     2473  int numberRead;
     2474  int length2;
     2475  numberRead = fread(&length2,sizeof(int),1,fp);
     2476  if (numberRead!=1)
     2477    return 1;
     2478  if (length2) {
     2479    // lengths must match
     2480    if (length!=length2)
     2481      return 2;
     2482    array = new double[length];
     2483    numberRead = fread(array,sizeof(double),length,fp);
     2484    if (numberRead!=length)
     2485      return 1;
     2486  }
     2487  return 0;
     2488}
     2489/* Restore model from file, returns 0 if success,
     2490   deletes current model */
     2491int
     2492ClpSimplex::restoreModel(const char * fileName)
     2493{
     2494  FILE * fp = fopen(fileName,"rb");
     2495  if (fp) {
     2496    // Get rid of current model
     2497    ClpModel::gutsOfDelete();
     2498    gutsOfDelete(0);
     2499    int i;
     2500    for (i=0;i<6;i++) {
     2501      rowArray_[i]=NULL;
     2502      columnArray_[i]=NULL;
     2503    }
     2504    // get an empty factorization so we can set tolerances etc
     2505    factorization_ = new ClpFactorization();
     2506    // Say sparse
     2507    factorization_->sparseThreshold(1);
     2508    Clp_scalars scalars;
     2509    CoinBigIndex numberRead;
     2510
     2511    // get scalars
     2512    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
     2513    if (numberRead!=1)
     2514      return 1;
     2515    // Fill in scalars
     2516    optimizationDirection_ = scalars.optimizationDirection;
     2517    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
     2518    objectiveValue_ = scalars.objectiveValue;
     2519    dualBound_ = scalars.dualBound;
     2520    dualTolerance_ = scalars.dualTolerance;
     2521    primalTolerance_ = scalars.primalTolerance;
     2522    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
     2523    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
     2524    infeasibilityCost_ = scalars.infeasibilityCost;
     2525    numberRows_ = scalars.numberRows;
     2526    numberColumns_ = scalars.numberColumns;
     2527    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
     2528    numberIterations_ = scalars.numberIterations;
     2529    problemStatus_ = scalars.problemStatus;
     2530    maximumIterations_ = scalars.maximumIterations;
     2531    lengthNames_ = scalars.lengthNames;
     2532    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
     2533    numberDualInfeasibilitiesWithoutFree_
     2534      = scalars.numberDualInfeasibilitiesWithoutFree;
     2535    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
     2536    numberRefinements_ = scalars.numberRefinements;
     2537    scalingFlag_ = scalars.scalingFlag;
     2538    algorithm_ = scalars.algorithm;
     2539    specialOptions_ = scalars.specialOptions;
     2540    // strings
     2541    CoinBigIndex length;
     2542    for (i=0;i<ClpLastStrParam;i++) {
     2543      numberRead = fread(&length,sizeof(int),1,fp);
     2544      if (numberRead!=1)
     2545        return 1;
     2546      char * array = new char[length+1];
     2547      numberRead = fread(array,length,1,fp);
     2548      if (numberRead!=1)
     2549        return 1;
     2550      array[length]='\0';
     2551      strParam_[i]=array;
     2552      delete [] array;
     2553    }
     2554    // arrays - in no particular order
     2555    if (inDoubleArray(rowActivity_,numberRows_,fp))
     2556        return 1;
     2557    if (inDoubleArray(columnActivity_,numberColumns_,fp))
     2558        return 1;
     2559    if (inDoubleArray(dual_,numberRows_,fp))
     2560        return 1;
     2561    if (inDoubleArray(reducedCost_,numberColumns_,fp))
     2562        return 1;
     2563    if (inDoubleArray(rowLower_,numberRows_,fp))
     2564        return 1;
     2565    if (inDoubleArray(rowUpper_,numberRows_,fp))
     2566        return 1;
     2567    if (inDoubleArray(objective_,numberColumns_,fp))
     2568        return 1;
     2569    if (inDoubleArray(rowObjective_,numberRows_,fp))
     2570        return 1;
     2571    if (inDoubleArray(columnLower_,numberColumns_,fp))
     2572        return 1;
     2573    if (inDoubleArray(columnUpper_,numberColumns_,fp))
     2574        return 1;
     2575    if (problemStatus_==1) {
     2576      if (inDoubleArray(ray_,numberRows_,fp))
     2577        return 1;
     2578    } else if (problemStatus_==2) {
     2579      if (inDoubleArray(ray_,numberColumns_,fp))
     2580        return 1;
     2581    } else {
     2582      // ray should be null
     2583      numberRead = fread(&length,sizeof(int),1,fp);
     2584      if (numberRead!=1)
     2585        return 1;
     2586      if (length)
     2587        return 2;
     2588    }
     2589    delete [] status_;
     2590    status_=NULL;
     2591    // status region
     2592    numberRead = fread(&length,sizeof(int),1,fp);
     2593    if (numberRead!=1)
     2594        return 1;
     2595    if (length) {
     2596      if (length!=numberRows_+numberColumns_)
     2597        return 1;
     2598      status_ = new char unsigned[length];
     2599      numberRead = fread(status_,sizeof(char),length, fp);
     2600      if (numberRead!=length)
     2601        return 1;
     2602    }
     2603    if (lengthNames_) {
     2604      char * array = new char[lengthNames_+1];
     2605      rowNames_.resize(0);
     2606      for (i=0;i<numberRows_;i++) {
     2607        numberRead = fread(&length,sizeof(int),1,fp);
     2608        if (numberRead!=1)
     2609          return 1;
     2610        numberRead = fread(array,length,1,fp);
     2611        if (numberRead!=1)
     2612          return 1;
     2613        rowNames_[i]=array;
     2614      }
     2615      columnNames_.resize(0);
     2616      for (i=0;i<numberColumns_;i++) {
     2617        numberRead = fread(&length,sizeof(int),1,fp);
     2618        if (numberRead!=1)
     2619          return 1;
     2620        numberRead = fread(array,length,1,fp);
     2621        if (numberRead!=1)
     2622          return 1;
     2623        columnNames_[i]=array;
     2624      }
     2625      delete [] array;
     2626    }
     2627    // Pivot choices
     2628    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
     2629    delete dualRowPivot_;
     2630    switch ((scalars.dualPivotChoice&63)) {
     2631    case 1:
     2632      // Dantzig
     2633      dualRowPivot_ = new ClpDualRowDantzig();
     2634      break;
     2635    case 2:
     2636      // Steepest - use mode
     2637      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
     2638      break;
     2639    default:
     2640      abort();
     2641    }
     2642    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
     2643    delete primalColumnPivot_;
     2644    switch ((scalars.primalPivotChoice&63)) {
     2645    case 1:
     2646      // Dantzig
     2647      primalColumnPivot_ = new ClpPrimalColumnDantzig();
     2648      break;
     2649    case 2:
     2650      // Steepest - use mode
     2651      primalColumnPivot_
     2652        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
     2653      break;
     2654    default:
     2655      abort();
     2656    }
     2657    assert(scalars.matrixStorageChoice==1);
     2658    delete matrix_;
     2659    // get arrays
     2660    numberRead = fread(&length,sizeof(int),1,fp);
     2661    if (numberRead!=1)
     2662      return 1;
     2663    double * elements = new double[length];
     2664    int * indices = new int[length];
     2665    CoinBigIndex * starts = new CoinBigIndex[numberColumns_];
     2666    int * lengths = new int[numberColumns_];
     2667    numberRead = fread(elements, sizeof(double),length,fp);
     2668    if (numberRead!=length)
     2669      return 1;
     2670    numberRead = fread(indices, sizeof(int),length,fp);
     2671    if (numberRead!=length)
     2672      return 1;
     2673    numberRead = fread(starts, sizeof(int),numberColumns_,fp);
     2674    if (numberRead!=numberColumns_)
     2675      return 1;
     2676    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
     2677    if (numberRead!=numberColumns_)
     2678      return 1;
     2679    // assign matrix
     2680    CoinPackedMatrix * matrix = new CoinPackedMatrix();
     2681    matrix->assignMatrix(true, numberRows_, numberColumns_,
     2682                         length, elements, indices, starts, lengths);
     2683    // and transfer to Clp
     2684    matrix_ = new ClpPackedMatrix(matrix);
     2685    // finished
     2686    fclose(fp);
     2687    return 0;
     2688  } else {
     2689    return -1;
     2690  }
     2691  return 0;
     2692}
     2693// value of incoming variable (in Dual)
     2694double
     2695ClpSimplex::valueIncomingDual() const
     2696{
     2697  // Need value of incoming for list of infeasibilities as may be infeasible
     2698  double valueIncoming = (dualOut_/alpha_)*directionOut_;
     2699  if (directionIn_==-1)
     2700    valueIncoming = upperIn_-valueIncoming;
     2701  else
     2702    valueIncoming = lowerIn_-valueIncoming;
     2703  return valueIncoming;
     2704}
     2705//#############################################################################
     2706
     2707ClpSimplexProgress::ClpSimplexProgress ()
     2708{
     2709  int i;
     2710  for (i=0;i<CLP_PROGRESS;i++) {
     2711    objective_[i] = 0.0;
     2712    infeasibility_[i] = -1.0; // set to an impossible value
     2713    numberInfeasibilities_[i]=-1;
     2714  }
     2715  numberTimes_ = 0;
     2716  numberBadTimes_ = 0;
     2717  model_ = NULL;
     2718}
     2719
     2720
     2721//-----------------------------------------------------------------------------
     2722
     2723ClpSimplexProgress::~ClpSimplexProgress ()
     2724{
     2725}
     2726// Copy constructor.
     2727ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
     2728{
     2729  int i;
     2730  for (i=0;i<CLP_PROGRESS;i++) {
     2731    objective_[i] = rhs.objective_[i];
     2732    infeasibility_[i] = rhs.infeasibility_[i];
     2733    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
     2734  }
     2735  numberTimes_ = rhs.numberTimes_;
     2736  numberBadTimes_ = rhs.numberBadTimes_;
     2737  model_ = rhs.model_;
     2738}
     2739// Copy constructor.from model
     2740ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model)
     2741{
     2742  int i;
     2743  for (i=0;i<CLP_PROGRESS;i++) {
     2744    objective_[i] = 0.0;
     2745    infeasibility_[i] = -1.0; // set to an impossible value
     2746    numberInfeasibilities_[i]=-1;
     2747  }
     2748  numberTimes_ = 0;
     2749  numberBadTimes_ = 0;
     2750  model_ = model;
     2751}
     2752// Assignment operator. This copies the data
     2753ClpSimplexProgress &
     2754ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
     2755{
     2756  if (this != &rhs) {
     2757    int i;
     2758    for (i=0;i<CLP_PROGRESS;i++) {
     2759      objective_[i] = rhs.objective_[i];
     2760      infeasibility_[i] = rhs.infeasibility_[i];
     2761      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
     2762    }
     2763    numberTimes_ = rhs.numberTimes_;
     2764    numberBadTimes_ = rhs.numberBadTimes_;
     2765    model_ = rhs.model_;
     2766  }
     2767  return *this;
     2768}
     2769// Seems to be something odd about exact comparison of doubles on linux
     2770static bool equalDouble(double value1, double value2)
     2771{
     2772  assert(sizeof(unsigned int)*2==sizeof(double));
     2773  unsigned int *i1 = (unsigned int *) &value1;
     2774  unsigned int *i2 = (unsigned int *) &value2;
     2775  return (i1[0]==i2[0]&&i1[1]==i2[1]);
     2776}
     2777int
     2778ClpSimplexProgress::looping()
     2779{
     2780  assert(model_);
     2781  double objective = model_->objectiveValue();
     2782  double infeasibility;
     2783  int numberInfeasibilities;
     2784  if (model_->algorithm()<0) {
     2785    // dual
     2786    infeasibility = model_->sumPrimalInfeasibilities();
     2787    numberInfeasibilities = model_->numberPrimalInfeasibilities();
     2788  } else {
     2789    //primal
     2790    infeasibility = model_->sumDualInfeasibilities();
     2791    numberInfeasibilities = model_->numberDualInfeasibilities();
     2792  }
     2793  int i;
     2794  int numberMatched=0;
     2795  int matched=0;
     2796
     2797  for (i=0;i<CLP_PROGRESS;i++) {
     2798    bool matchedOnObjective = equalDouble(objective,objective_[i]);
     2799    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
     2800    bool matchedOnInfeasibilities =
     2801      (numberInfeasibilities==numberInfeasibilities_[i]);
     2802
     2803    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
     2804      matched |= (1<<i);
     2805      numberMatched++;
     2806      // here mainly to get over compiler bug?
     2807      if (model_->messageHandler()->logLevel()>10)
     2808        printf("%d %d %d %d %d loop check\n",i,numberMatched,
     2809               matchedOnObjective, matchedOnInfeasibility,
     2810               matchedOnInfeasibilities);
     2811    }
     2812    if (i) {
     2813      objective_[i-1] = objective_[i];
     2814      infeasibility_[i-1] = infeasibility_[i];
     2815      numberInfeasibilities_[i-1]=numberInfeasibilities_[i];
     2816    }
     2817  }
     2818  objective_[CLP_PROGRESS-1] = objective;
     2819  infeasibility_[CLP_PROGRESS-1] = infeasibility;
     2820  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
     2821  if (model_->progressFlag())
     2822    numberMatched=0;
     2823  numberTimes_++;
     2824  if (numberTimes_<10)
     2825    numberMatched=0;
     2826  // skip if just last time as may be checking something
     2827  if (matched==(1<<(CLP_PROGRESS-1)))
     2828    numberMatched=0;
     2829  if (numberMatched) {
     2830    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
     2831      <<numberMatched
     2832      <<matched
     2833      <<numberTimes_
     2834      <<CoinMessageEol;
     2835    numberBadTimes_++;
     2836    if (numberBadTimes_<10) {
     2837      if (model_->algorithm()<0) {
     2838        // dual - change tolerance
     2839        model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
     2840        // if infeasible increase dual bound
     2841        if (model_->dualBound()<1.0e19) {
     2842          model_->setDualBound(model_->dualBound()*1.1);
     2843        }
     2844      } else {
     2845        // primal - change tolerance
     2846        model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
     2847        // if infeasible increase infeasibility cost
     2848        if (model_->infeasibilityCost()<1.0e19) {
     2849          model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
     2850        }
     2851      }
     2852    } else {
     2853      model_->messageHandler()->message(CLP_LOOP,model_->messages())
     2854        <<CoinMessageEol;
     2855      // look at solution and maybe declare victory
     2856      if (infeasibility<1.0e-4)
     2857        return 0;
     2858      else
     2859        return 3;
     2860    }
     2861  }
     2862  return -1;
     2863}
     2864// Sanity check on input data - returns true if okay
     2865bool
     2866ClpSimplex::sanityCheck()
     2867{
     2868  // bad if empty
     2869  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
     2870    handler_->message(CLP_EMPTY_PROBLEM,messages_)
     2871      <<numberRows_
     2872      <<numberColumns_
     2873      <<matrix_->getNumElements()
     2874      <<CoinMessageEol;
     2875    problemStatus_=4;
     2876    return false;
     2877  }
     2878  int numberBad , numberBadBounds;
     2879  double largestBound, smallestBound, minimumGap;
     2880  double smallestObj, largestObj;
     2881  int firstBad;
     2882  int modifiedBounds=0;
     2883  int i;
     2884  numberBad=0;
     2885  numberBadBounds=0;
     2886  firstBad=-1;
     2887  minimumGap=1.0e100;
     2888  smallestBound=1.0e100;
     2889  largestBound=0.0;
     2890  smallestObj=1.0e100;
     2891  largestObj=0.0;
     2892  // If bounds are too close - fix
     2893  double fixTolerance = 1.1*primalTolerance_;
     2894  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
     2895    double value;
     2896    value = fabs(cost_[i]);
     2897    if (value>1.0e50) {
     2898      numberBad++;
     2899      if (firstBad<0)
     2900        firstBad=i;
     2901    } else if (value) {
     2902      if (value>largestObj)
     2903        largestObj=value;
     2904      if (value<smallestObj)
     2905        smallestObj=value;
     2906    }
     2907    value=upper_[i]-lower_[i];
     2908    if (value<-primalTolerance_) {
     2909      numberBadBounds++;
     2910      if (firstBad<0)
     2911        firstBad=i;
     2912    } else if (value<=fixTolerance) {
     2913      if (value) {
     2914        // modify
     2915        upper_[i] = lower_[i];
     2916        modifiedBounds++;
     2917      }
     2918    } else {
     2919      if (value<minimumGap)
     2920        minimumGap=value;
     2921    }
     2922    if (lower_[i]>-1.0e100&&lower_[i]) {
     2923      value = fabs(lower_[i]);
     2924      if (value>largestBound)
     2925        largestBound=value;
     2926      if (value<smallestBound)
     2927        smallestBound=value;
     2928    }
     2929    if (upper_[i]<1.0e100&&upper_[i]) {
     2930      value = fabs(upper_[i]);
     2931      if (value>largestBound)
     2932        largestBound=value;
     2933      if (value<smallestBound)
     2934        smallestBound=value;
     2935    }
     2936  }
     2937  if (largestBound)
     2938    handler_->message(CLP_RIMSTATISTICS3,messages_)
     2939      <<smallestBound
     2940      <<largestBound
     2941      <<minimumGap
     2942      <<CoinMessageEol;
     2943  minimumGap=1.0e100;
     2944  smallestBound=1.0e100;
     2945  largestBound=0.0;
     2946  for (i=0;i<numberColumns_;i++) {
     2947    double value;
     2948    value = fabs(cost_[i]);
     2949    if (value>1.0e50) {
     2950      numberBad++;
     2951      if (firstBad<0)
     2952        firstBad=i;
     2953    } else if (value) {
     2954      if (value>largestObj)
     2955        largestObj=value;
     2956      if (value<smallestObj)
     2957        smallestObj=value;
     2958    }
     2959    value=upper_[i]-lower_[i];
     2960    if (value<-primalTolerance_) {
     2961      numberBadBounds++;
     2962      if (firstBad<0)
     2963        firstBad=i;
     2964    } else if (value<=fixTolerance) {
     2965      if (value) {
     2966        // modify
     2967        upper_[i] = lower_[i];
     2968        modifiedBounds++;
     2969      }
     2970    } else {
     2971      if (value<minimumGap)
     2972        minimumGap=value;
     2973    }
     2974    if (lower_[i]>-1.0e100&&lower_[i]) {
     2975      value = fabs(lower_[i]);
     2976      if (value>largestBound)
     2977        largestBound=value;
     2978      if (value<smallestBound)
     2979        smallestBound=value;
     2980    }
     2981    if (upper_[i]<1.0e100&&upper_[i]) {
     2982      value = fabs(upper_[i]);
     2983      if (value>largestBound)
     2984        largestBound=value;
     2985      if (value<smallestBound)
     2986        smallestBound=value;
     2987    }
     2988  }
     2989  char rowcol[]={'R','C'};
     2990  if (numberBad) {
     2991    handler_->message(CLP_BAD_BOUNDS,messages_)
     2992      <<numberBad
     2993      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
     2994      <<CoinMessageEol;
     2995    problemStatus_=4;
     2996    return false;
     2997  }
     2998  if (modifiedBounds)
     2999    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
     3000      <<modifiedBounds
     3001      <<CoinMessageEol;
     3002  handler_->message(CLP_RIMSTATISTICS1,messages_)
     3003    <<smallestObj
     3004    <<largestObj
     3005    <<CoinMessageEol;
     3006  if (largestBound)
     3007    handler_->message(CLP_RIMSTATISTICS2,messages_)
     3008      <<smallestBound
     3009      <<largestBound
     3010      <<minimumGap
     3011      <<CoinMessageEol;
     3012  return true;
     3013}
     3014// Set up status array (for OsiClp)
     3015void
     3016ClpSimplex::createStatus()
     3017{
     3018  if(!status_)
     3019    status_ = new unsigned char [numberColumns_+numberRows_];
     3020  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
     3021  int i;
     3022  // set column status to one nearest zero
     3023  for (i=0;i<numberColumns_;i++) {
     3024#if 0
     3025    if (columnLower_[i]>=0.0) {
     3026      setColumnStatus(i,atLowerBound);
     3027    } else if (columnUpper_[i]<=0.0) {
     3028      setColumnStatus(i,atUpperBound);
     3029    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
     3030      // free
     3031      setColumnStatus(i,isFree);
     3032    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
     3033      setColumnStatus(i,atLowerBound);
     3034    } else {
     3035      setColumnStatus(i,atUpperBound);
     3036    }
     3037#else
     3038    setColumnStatus(i,atLowerBound);
     3039#endif
     3040  }
     3041  for (i=0;i<numberRows_;i++) {
     3042    setRowStatus(i,basic);
     3043  }
     3044}
     3045/* Loads a problem (the constraints on the
     3046   rows are given by lower and upper bounds). If a pointer is 0 then the
     3047   following values are the default:
     3048   <ul>
     3049   <li> <code>colub</code>: all columns have upper bound infinity
     3050   <li> <code>collb</code>: all columns have lower bound 0
     3051   <li> <code>rowub</code>: all rows have upper bound infinity
     3052   <li> <code>rowlb</code>: all rows have lower bound -infinity
     3053   <li> <code>obj</code>: all variables have 0 objective coefficient
     3054   </ul>
     3055*/
     3056void
     3057ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
     3058                    const double* collb, const double* colub,   
     3059                    const double* obj,
     3060                    const double* rowlb, const double* rowub,
     3061                    const double * rowObjective)
     3062{
     3063  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
     3064                        rowObjective);
     3065  createStatus();
     3066}
     3067void
     3068ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
     3069                    const double* collb, const double* colub,   
     3070                    const double* obj,
     3071                    const double* rowlb, const double* rowub,
     3072                    const double * rowObjective)
     3073{
     3074  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
     3075                        rowObjective);
     3076  createStatus();
     3077}
     3078
     3079/* Just like the other loadProblem() method except that the matrix is
     3080   given in a standard column major ordered format (without gaps). */
     3081void
     3082ClpSimplex::loadProblem (  const int numcols, const int numrows,
     3083                    const CoinBigIndex* start, const int* index,
     3084                    const double* value,
     3085                    const double* collb, const double* colub,   
     3086                    const double* obj,
     3087                    const double* rowlb, const double* rowub,
     3088                    const double * rowObjective)
     3089{
     3090  ClpModel::loadProblem(numcols, numrows, start, index, value,
     3091                          collb, colub, obj, rowlb, rowub,
     3092                          rowObjective);
     3093  createStatus();
     3094}
     3095void
     3096ClpSimplex::loadProblem (  const int numcols, const int numrows,
     3097                           const CoinBigIndex* start, const int* index,
     3098                           const double* value,const int * length,
     3099                           const double* collb, const double* colub,   
     3100                           const double* obj,
     3101                           const double* rowlb, const double* rowub,
     3102                           const double * rowObjective)
     3103{
     3104  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
     3105                          collb, colub, obj, rowlb, rowub,
     3106                          rowObjective);
     3107  createStatus();
     3108}
     3109// Read an mps file from the given filename
     3110int
     3111ClpSimplex::readMps(const char *filename,
     3112            bool keepNames,
     3113            bool ignoreErrors)
     3114{
     3115  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
     3116  createStatus();
     3117  return status;
     3118}
     3119// Just check solution (for external use)
     3120void
     3121ClpSimplex::checkSolution()
     3122{
     3123  // put in standard form
     3124  createRim(7+8+16);
     3125  dualTolerance_=dblParam_[ClpDualTolerance];
     3126  primalTolerance_=dblParam_[ClpPrimalTolerance];
     3127  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
     3128  checkDualSolution();
     3129#ifdef CLP_DEBUG
     3130  int i;
     3131  double value=0.0;
     3132  for (i=0;i<numberRows_+numberColumns_;i++)
     3133    value += dj_[i]*solution_[i];
     3134  printf("dual value %g, primal %g\n",value,objectiveValue_);
     3135#endif
     3136  // release extra memory
     3137  deleteRim(false);
     3138}
  • trunk/ClpSimplexDual.cpp

    r6 r50  
    9797#include "ClpSimplexDual.hpp"
    9898#include "ClpFactorization.hpp"
    99 #include "OsiPackedMatrix.hpp"
    100 #include "OsiIndexedVector.hpp"
    101 #include "OsiWarmStartBasis.hpp"
     99#include "CoinPackedMatrix.hpp"
     100#include "CoinIndexedVector.hpp"
     101#include "CoinWarmStartBasis.hpp"
    102102#include "ClpDualRowDantzig.hpp"
    103103#include "ClpMessage.hpp"
     
    107107#include <stdio.h>
    108108#include <iostream>
    109 // This returns a non const array filled with input from scalar
    110 // or actual array
    111 template <class T> inline T*
    112 copyOfArray( const T * array, const int size, T value)
    113 {
    114   T * arrayNew = new T[size];
    115   if (array)
    116     CoinDisjointCopyN(array,size,arrayNew);
    117   else
    118     CoinFillN ( arrayNew, size,value);
    119   return arrayNew;
    120 }
    121 
    122 // This returns a non const array filled with actual array (or NULL)
    123 template <class T> inline T*
    124 copyOfArray( const T * array, const int size)
    125 {
    126   if (array) {
    127     T * arrayNew = new T[size];
    128     CoinDisjointCopyN(array,size,arrayNew);
    129     return arrayNew;
    130   } else {
    131     return NULL;
    132   }
    133 }
     109#define CHECK_DJ
    134110// dual
    135111int ClpSimplexDual::dual ( )
     
    232208  assert(rowUpper_);
    233209
     210
     211  algorithm_ = -1;
     212  dualTolerance_=dblParam_[ClpDualTolerance];
     213  primalTolerance_=dblParam_[ClpPrimalTolerance];
     214
     215  // put in standard form (and make row copy)
     216  // create modifiable copies of model rim and do optional scaling
     217  bool goodMatrix = createRim(7+8+16,true);
     218
     219  // save dual bound
     220  double saveDualBound_ = dualBound_;
     221  // save perturbation
     222  int savePerturbation = perturbation_;
     223  // save if sparse factorization wanted
     224  int saveSparse = factorization_->sparseThreshold();
     225
     226  if (goodMatrix) {
     227    // Problem looks okay
     228    int iRow,iColumn;
     229    // Do initial factorization
     230    // and set certain stuff
     231    // We can either set increasing rows so ...IsBasic gives pivot row
     232    // or we can just increment iBasic one by one
     233    // for now let ...iBasic give pivot row
     234    factorization_->increasingRows(2);
     235    // row activities have negative sign
     236    factorization_->slackValue(-1.0);
     237    factorization_->zeroTolerance(1.0e-13);
     238   
     239    int factorizationStatus = internalFactorize(0);
     240    if (factorizationStatus<0)
     241      return 1; // some error
     242    else if (factorizationStatus)
     243      handler_->message(CLP_SINGULARITIES,messages_)
     244        <<factorizationStatus
     245        <<CoinMessageEol;
     246   
     247    // If user asked for perturbation - do it
     248   
     249    if (perturbation_<100)
     250      perturb();
     251   
     252    double objectiveChange;
     253    // for dual we will change bounds using dualBound_
     254    // for this we need clean basis so it is after factorize
     255    gutsOfSolution(rowActivityWork_,columnActivityWork_);
     256   
     257    numberFake_ =0; // Number of variables at fake bounds
     258    changeBounds(true,NULL,objectiveChange);
     259   
     260    problemStatus_ = -1;
     261    numberIterations_=0;
     262   
     263    int lastCleaned=0; // last time objective or bounds cleaned up
     264   
     265    // number of times we have declared optimality
     266    numberTimesOptimal_=0;
     267   
     268    // Progress indicator
     269    ClpSimplexProgress progress(this);
     270   
     271    // This says whether to restore things etc
     272    int factorType=0;
     273    /*
     274      Status of problem:
     275      0 - optimal
     276      1 - infeasible
     277      2 - unbounded
     278      -1 - iterating
     279      -2 - factorization wanted
     280      -3 - redo checking without factorization
     281      -4 - looks infeasible
     282    */
     283    while (problemStatus_<0) {
     284      // clear
     285      for (iRow=0;iRow<4;iRow++) {
     286        rowArray_[iRow]->clear();
     287      }   
     288     
     289      for (iColumn=0;iColumn<2;iColumn++) {
     290        columnArray_[iColumn]->clear();
     291      }   
     292     
     293      // give matrix (and model costs and bounds a chance to be
     294      // refreshed (normally null)
     295      matrix_->refresh(this);
     296      // If getting nowhere - why not give it a kick
     297#if 0
     298      // does not seem to work too well - do some more work
     299      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
     300        perturb();
     301#endif
     302      // may factorize, checks if problem finished
     303      statusOfProblemInDual(lastCleaned,factorType,progress);
     304     
     305      // Say good factorization
     306      factorType=1;
     307      if (saveSparse) {
     308        // use default at present
     309        factorization_->sparseThreshold(0);
     310        factorization_->goSparse();
     311      }
     312     
     313      // Do iterations
     314      whileIterating();
     315    }
     316  }
     317
     318  //assert(!numberFake_||problemStatus_); // all bounds should be okay
     319  factorization_->sparseThreshold(saveSparse);
     320  // Get rid of some arrays and empty factorization
     321  deleteRim();
     322  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
     323    <<objectiveValue()
     324    <<CoinMessageEol;
     325  // Restore any saved stuff
     326  perturbation_ = savePerturbation;
     327  dualBound_ = saveDualBound_;
     328  return problemStatus_;
     329}
     330/* Reasons to come out:
     331   -1 iterations etc
     332   -2 inaccuracy
     333   -3 slight inaccuracy (and done iterations)
     334   +0 looks optimal (might be unbounded - but we will investigate)
     335   +1 looks infeasible
     336   +3 max iterations
     337 */
     338int
     339ClpSimplexDual::whileIterating()
     340{
    234341#ifdef CLP_DEBUG
    235342  int debugIteration=-1;
    236343#endif
    237 
    238   algorithm_ = -1;
    239   dualTolerance_=dblParam_[OsiDualTolerance];
    240   primalTolerance_=dblParam_[OsiPrimalTolerance];
    241 
    242   // put in standard form (and make row copy)
    243   // create modifiable copies of model rim and do optional scaling
    244   createRim(7+8+16,true);
    245 
    246   // save dual bound
    247   double saveDualBound_ = dualBound_;
    248 
    249   int iRow,iColumn;
    250   // Do initial factorization
    251   // and set certain stuff
    252   // We can either set increasing rows so ...IsBasic gives pivot row
    253   // or we can just increment iBasic one by one
    254   // for now let ...iBasic give pivot row
    255   factorization_->increasingRows(2);
    256   // row activities have negative sign
    257   factorization_->slackValue(-1.0);
    258   factorization_->zeroTolerance(1.0e-13);
    259   if (internalFactorize(0))
    260     return 1; // some error
    261 
    262   // If user asked for perturbation - do it
    263   int savePerturbation = perturbation_;
    264 
    265   if (perturbation_<100)
    266     perturb();
    267 
    268   double objectiveChange;
    269   // for dual we will change bounds using dualBound_
    270   // for this we need clean basis so it is after factorize
    271   gutsOfSolution(rowActivityWork_,columnActivityWork_);
    272   changeBounds(true,NULL,objectiveChange);
    273 
    274   problemStatus_ = -1;
    275   numberIterations_=0;
    276 
    277   int lastCleaned=0; // last time objective or bounds cleaned up
    278 
    279   // number of times we have declared optimality
    280   numberTimesOptimal_=0;
    281 
    282   // This says whether to restore things etc
    283   int factorType=0;
    284   /*
    285     Status of problem:
    286     0 - optimal
    287     1 - infeasible
    288     2 - unbounded
    289     -1 - iterating
    290     -2 - factorization wanted
    291     -3 - redo checking without factorization
    292     -4 - looks infeasible
    293   */
    294   while (problemStatus_<0) {
    295     // clear
    296     for (iRow=0;iRow<4;iRow++) {
    297       rowArray_[iRow]->clear();
    298     }   
    299    
    300     for (iColumn=0;iColumn<2;iColumn++) {
    301       columnArray_[iColumn]->clear();
    302     }   
    303 
    304     // give matrix (and model costs and bounds a chance to be
    305     // refreshed (normally null)
    306     matrix_->refresh(this);
    307     // If getting nowhere - why not give it a kick
    308 #if 0
    309     // does not seem to work too well - do some more work
    310     if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
    311       perturb();
    312 #endif
    313     // may factorize, checks if problem finished
    314     statusOfProblemInDual(lastCleaned,factorType);
    315 
    316     // Say good factorization
    317     factorType=1;
    318 
    319     // status stays at -1 while iterating, >=0 finished, -2 to invert
    320     // status -3 to go to top without an invert
    321     while (problemStatus_==-1) {
     344  // status stays at -1 while iterating, >=0 finished, -2 to invert
     345  // status -3 to go to top without an invert
     346  int returnCode = -1;
     347  while (problemStatus_==-1) {
    322348#ifdef CLP_DEBUG
    323       {
    324         int i;
    325         for (i=0;i<4;i++) {
    326           rowArray_[i]->checkClear();
    327         }   
    328         for (i=0;i<2;i++) {
    329           columnArray_[i]->checkClear();
    330         }   
    331       }     
     349    {
     350      int i;
     351      for (i=0;i<4;i++) {
     352        rowArray_[i]->checkClear();
     353      }   
     354      for (i=0;i<2;i++) {
     355        columnArray_[i]->checkClear();
     356      }   
     357    }     
    332358#endif
    333359#if CLP_DEBUG>2
    334       // very expensive
    335       if (numberIterations_>0&&numberIterations_<-801) {
    336         handler_->setLogLevel(63);
    337         double saveValue = objectiveValue_;
    338         double * saveRow1 = new double[numberRows_];
    339         double * saveRow2 = new double[numberRows_];
    340         memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
    341         memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
    342         double * saveColumn1 = new double[numberColumns_];
    343         double * saveColumn2 = new double[numberColumns_];
    344         memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
    345         memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
    346         gutsOfSolution(rowActivityWork_,columnActivityWork_);
    347         printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
    348                numberIterations_,
    349                saveValue,objectiveValue_,sumDualInfeasibilities_);
    350         if (saveValue>objectiveValue_+1.0e-2)
    351           printf("**bad**\n");
    352         memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
    353         memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
    354         memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
    355         memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
    356         delete [] saveRow1;
    357         delete [] saveRow2;
    358         delete [] saveColumn1;
    359         delete [] saveColumn2;
    360         objectiveValue_=saveValue;
     360    // very expensive
     361    if (numberIterations_>0&&numberIterations_<-801) {
     362      handler_->setLogLevel(63);
     363      double saveValue = objectiveValue_;
     364      double * saveRow1 = new double[numberRows_];
     365      double * saveRow2 = new double[numberRows_];
     366      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
     367      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
     368      double * saveColumn1 = new double[numberColumns_];
     369      double * saveColumn2 = new double[numberColumns_];
     370      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
     371      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
     372      gutsOfSolution(rowActivityWork_,columnActivityWork_);
     373      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
     374             numberIterations_,
     375             saveValue,objectiveValue_,sumDualInfeasibilities_);
     376      if (saveValue>objectiveValue_+1.0e-2)
     377        printf("**bad**\n");
     378      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
     379      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
     380      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
     381      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
     382      delete [] saveRow1;
     383      delete [] saveRow2;
     384      delete [] saveColumn1;
     385      delete [] saveColumn2;
     386      objectiveValue_=saveValue;
     387    }
     388#endif
     389#ifdef CLP_DEBUG
     390    {
     391      int iSequence, number=numberRows_+numberColumns_;
     392      for (iSequence=0;iSequence<number;iSequence++) {
     393        double lowerValue=lower_[iSequence];
     394        double upperValue=upper_[iSequence];
     395        double value=solution_[iSequence];
     396        if(getStatus(iSequence)!=basic) {
     397          assert(lowerValue>-1.0e20);
     398          assert(upperValue<1.0e20);
     399        }
     400        switch(getStatus(iSequence)) {
     401         
     402        case basic:
     403          break;
     404        case isFree:
     405        case superBasic:
     406          break;
     407        case atUpperBound:
     408          assert (fabs(value-upperValue)<=primalTolerance_) ;
     409          break;
     410        case atLowerBound:
     411          assert (fabs(value-lowerValue)<=primalTolerance_) ;
     412          break;
     413        }
    361414      }
    362 #endif
     415    }
     416    if(numberIterations_==debugIteration) {
     417      printf("dodgy iteration coming up\n");
     418    }
     419#endif
     420    // choose row to go out
     421    // dualRow will go to virtual row pivot choice algorithm
     422    dualRow();
     423    if (pivotRow_>=0) {
     424      // we found a pivot row
     425      handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
     426        <<pivotRow_
     427        <<CoinMessageEol;
     428      // check accuracy of weights
     429      dualRowPivot_->checkAccuracy();
     430      // get sign for finding row of tableau
     431      rowArray_[0]->insert(pivotRow_,directionOut_);
     432      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
     433      // put row of tableau in rowArray[0] and columnArray[0]
     434      matrix_->transposeTimes(this,-1.0,
     435                              rowArray_[0],columnArray_[1],columnArray_[0]);
     436      // rowArray has pi equivalent
     437      // do ratio test
     438      dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
     439                 rowArray_[3]);
     440      if (sequenceIn_>=0) {
     441        // normal iteration
     442        // update the incoming column
     443        unpack(rowArray_[1]);
     444        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
     445        // and update dual weights (can do in parallel - with extra array)
     446        dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
     447                                     rowArray_[1]);
     448        // see if update stable
     449        double btranAlpha = -alpha_*directionOut_; // for check
     450        alpha_=(*rowArray_[1])[pivotRow_];
    363451#ifdef CLP_DEBUG
    364       {
    365         int iSequence, number=numberRows_+numberColumns_;
    366         for (iSequence=0;iSequence<number;iSequence++) {
    367           double lowerValue=lower_[iSequence];
    368           double upperValue=upper_[iSequence];
    369           double value=solution_[iSequence];
    370           if(getStatus(iSequence)!=ClpSimplex::basic) {
    371             assert(lowerValue>-1.0e20);
    372             assert(upperValue<1.0e20);
    373           }
    374           switch(getStatus(iSequence)) {
    375            
    376           case ClpSimplex::basic:
     452        if ((handler_->logLevel()&32))
     453          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
     454#endif
     455        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
     456            fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
     457          handler_->message(CLP_DUAL_CHECK,messages_)
     458            <<btranAlpha
     459            <<alpha_
     460            <<CoinMessageEol;
     461          dualRowPivot_->unrollWeights();
     462          if (factorization_->pivots()) {
     463            problemStatus_=-2; // factorize now
     464            rowArray_[0]->clear();
     465            rowArray_[1]->clear();
     466            columnArray_[0]->clear();
     467            returnCode=-2;
    377468            break;
    378           case ClpSimplex::isFree:
    379           case ClpSimplex::superBasic:
    380             break;
    381           case ClpSimplex::atUpperBound:
    382             assert (fabs(value-upperValue)<=primalTolerance_) ;
    383             break;
    384           case ClpSimplex::atLowerBound:
    385             assert (fabs(value-lowerValue)<=primalTolerance_) ;
    386               break;
    387           }
    388         }
    389       }
    390       if(numberIterations_==debugIteration) {
    391         printf("dodgy iteration coming up\n");
    392       }
    393 #endif
    394       // choose row to go out
    395       dualRow();
    396       if (pivotRow_>=0) {
    397         // we found a pivot row
    398         handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
    399           <<pivotRow_
    400           <<OsiMessageEol;
    401         // check accuracy of weights
    402         dualRowPivot_->checkAccuracy();
    403         // get sign for finding row of tableau
    404         rowArray_[0]->insert(pivotRow_,directionOut_);
    405         factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
    406         // put row of tableau in rowArray[0] and columnArray[0]
    407         matrix_->transposeTimes(this,-1.0,
    408                                 rowArray_[0],columnArray_[1],columnArray_[0]);
    409         // rowArray has pi equivalent
    410         // do ratio test
    411         dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
    412                    rowArray_[3]);
    413         if (sequenceIn_>=0) {
    414           // normal iteration
    415           // update the incoming column
    416           unpack(rowArray_[1]);
    417           factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
    418           // and update dual weights (can do in parallel - with extra array)
    419           dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
    420                                        rowArray_[1]);
    421           // see if update stable
    422           double btranAlpha = -alpha_*directionOut_; // for check
    423           alpha_=(*rowArray_[1])[pivotRow_];
    424 #ifdef CLP_DEBUG
    425           if ((handler_->logLevel()&32))
    426             printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
    427 #endif
    428           if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    429               fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
    430             handler_->message(CLP_DUAL_CHECK,messages_)
    431               <<btranAlpha
    432               <<alpha_
    433               <<OsiMessageEol;
    434             dualRowPivot_->unrollWeights();
    435             if (factorization_->pivots()) {
    436               problemStatus_=-2; // factorize now
    437               rowArray_[0]->clear();
    438               rowArray_[1]->clear();
    439               columnArray_[0]->clear();
    440               break;
    441             } else {
    442               // take on more relaxed criterion
    443               if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    444                   fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
    445                 // need to reject something
    446                 char x = isColumn(sequenceOut_) ? 'C' :'R';
    447                 handler_->message(CLP_SIMPLEX_FLAG,messages_)
    448                   <<x<<sequenceWithin(sequenceOut_)
    449                   <<OsiMessageEol;
    450                 setFlagged(sequenceOut_);
    451                 rowArray_[0]->clear();
    452                 rowArray_[1]->clear();
    453                 columnArray_[0]->clear();
    454                 continue;
    455               }
    456             }
    457           }
    458           // update duals BEFORE replaceColumn so can do updateColumn
    459           objectiveChange=0.0;
    460           // do duals first as variables may flip bounds
    461           // rowArray_[0] and columnArray_[0] may have flips
    462           // so use rowArray_[3] for work array from here on
    463           int nswapped =
    464             updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
    465                         objectiveChange);
    466           // which will change basic solution
    467           if (nswapped) {
    468 #ifdef CLP_DEBUG
    469           if ((handler_->logLevel()&16))
    470             printf("old dualOut_ %g, v %g, l %g, u %g - new ",
    471                    dualOut_,valueOut_,lowerOut_,upperOut_);
    472             double oldOut=dualOut_;
    473 #endif
    474             factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
    475             dualRowPivot_->updatePrimalSolution(rowArray_[2],
    476                                                 1.0,objectiveChange);
    477 
    478             // recompute dualOut_
    479             valueOut_ = solution_[sequenceOut_];
    480             if (directionOut_<0) {
    481               dualOut_ = valueOut_ - upperOut_;
    482             } else {
    483               dualOut_ = lowerOut_ - valueOut_;
    484             }
    485 #ifdef CLP_DEBUG
    486             if ((handler_->logLevel()&16))
    487               printf("%g\n",dualOut_);
    488             assert(dualOut_<=oldOut);
    489 #endif
    490             if(dualOut_<0.0&&factorization_->pivots()) {
    491               // going backwards - factorize
    492               dualRowPivot_->unrollWeights();
    493               problemStatus_=-2; // factorize now
    494               break;
    495             }
    496           }
    497           // amount primal will move
    498           double movement = -dualOut_*directionOut_/alpha_;
    499           // so objective should increase by fabs(dj)*movement
    500           // but we already have objective change - so check will be good
    501           if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
    502 #ifdef CLP_DEBUG
    503             if (handler_->logLevel()&32)
    504               printf("movement %g, swap change %g, rest %g  * %g\n",
    505                      objectiveChange+fabs(movement*dualIn_),
    506                      objectiveChange,movement,dualIn_);
    507 #endif
    508             if(factorization_->pivots()>5) {
    509               // going backwards - factorize
    510               dualRowPivot_->unrollWeights();
    511               problemStatus_=-2; // factorize now
    512               break;
    513             }
    514           }
    515           // if stable replace in basis
    516           int updateStatus = factorization_->replaceColumn(rowArray_[2],
    517                                                            pivotRow_,
    518                                                            alpha_);
    519           if (updateStatus==1) {
    520             // slight error
    521             if (factorization_->pivots()>5)
    522               problemStatus_=-2; // factorize now
    523           } else if (updateStatus==2) {
    524             // major error
    525             dualRowPivot_->unrollWeights();
    526             // later we may need to unwind more e.g. fake bounds
    527             if (factorization_->pivots()) {
    528               problemStatus_=-2; // factorize now
    529               break;
    530             } else {
     469          } else {
     470            // take on more relaxed criterion
     471            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
     472                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
    531473              // need to reject something
    532474              char x = isColumn(sequenceOut_) ? 'C' :'R';
    533475              handler_->message(CLP_SIMPLEX_FLAG,messages_)
    534476                <<x<<sequenceWithin(sequenceOut_)
    535                 <<OsiMessageEol;
     477                <<CoinMessageEol;
    536478              setFlagged(sequenceOut_);
     479              lastBadIteration_ = numberIterations_; // say be more cautious
    537480              rowArray_[0]->clear();
    538481              rowArray_[1]->clear();
     
    540483              continue;
    541484            }
    542           } else if (updateStatus==3) {
    543             // out of memory
    544             // increase space if not many iterations
    545             if (factorization_->pivots()<
    546                 0.5*factorization_->maximumPivots()&&
    547                 factorization_->pivots()<200)
    548               factorization_->areaFactor(
    549                                          factorization_->areaFactor() * 1.1);
     485          }
     486        }
     487        // update duals BEFORE replaceColumn so can do updateColumn
     488        double objectiveChange=0.0;
     489        // do duals first as variables may flip bounds
     490        // rowArray_[0] and columnArray_[0] may have flips
     491        // so use rowArray_[3] for work array from here on
     492        int nswapped =
     493          updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
     494                            objectiveChange);
     495        // which will change basic solution
     496        if (nswapped) {
     497#ifdef CLP_DEBUG
     498          if ((handler_->logLevel()&16))
     499            printf("old dualOut_ %g, v %g, l %g, u %g - new ",
     500                   dualOut_,valueOut_,lowerOut_,upperOut_);
     501          double oldOut=dualOut_;
     502#endif
     503          factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
     504          dualRowPivot_->updatePrimalSolution(rowArray_[2],
     505                                              1.0,objectiveChange);
     506         
     507          // recompute dualOut_
     508          valueOut_ = solution_[sequenceOut_];
     509          if (directionOut_<0) {
     510            dualOut_ = valueOut_ - upperOut_;
     511          } else {
     512            dualOut_ = lowerOut_ - valueOut_;
     513          }
     514#ifdef CLP_DEBUG
     515          if ((handler_->logLevel()&16))
     516            printf("%g\n",dualOut_);
     517          assert(dualOut_<=oldOut);
     518#endif
     519          if(dualOut_<0.0&&factorization_->pivots()) {
     520            // going backwards - factorize
     521            dualRowPivot_->unrollWeights();
    550522            problemStatus_=-2; // factorize now
    551           }
    552           // update primal solution
    553           if (theta_<0.0) {
    554 #ifdef CLP_DEBUG
    555             if (handler_->logLevel()&32)
    556               printf("negative theta %g\n",theta_);
    557 #endif
    558             theta_=0.0;
    559           }
    560           // do actual flips
    561           flipBounds(rowArray_[0],columnArray_[0],theta_);
    562           dualRowPivot_->updatePrimalSolution(rowArray_[1],
    563                                               movement,
    564                                               objectiveChange);
    565 #ifdef CLP_DEBUG
    566           double oldobj=objectiveValue_;
    567 #endif
    568           int whatNext=housekeeping(objectiveChange);
    569           // and set bounds correctly
    570           originalBound(sequenceIn_);
    571           changeBound(sequenceOut_);
    572 #ifdef CLP_DEBUG
    573           if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
    574             printf("obj backwards %g %g\n",objectiveValue_,oldobj);
    575 #endif
    576           if (whatNext==1) {
    577             problemStatus_ =-2; // refactorize
    578           } else if (whatNext==2) {
    579             // maximum iterations or equivalent
    580             problemStatus_= 3;
     523            returnCode=-2;
    581524            break;
    582525          }
    583         } else {
    584           // no incoming column is valid
     526        }
     527        // amount primal will move
     528        double movement = -dualOut_*directionOut_/alpha_;
     529        // so objective should increase by fabs(dj)*movement
     530        // but we already have objective change - so check will be good
     531        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
    585532#ifdef CLP_DEBUG
    586533          if (handler_->logLevel()&32)
    587             printf("** no column pivot\n");
    588 #endif
     534            printf("movement %g, swap change %g, rest %g  * %g\n",
     535                   objectiveChange+fabs(movement*dualIn_),
     536                   objectiveChange,movement,dualIn_);
     537#endif
     538          if(factorization_->pivots()>5) {
     539            // going backwards - factorize
     540            dualRowPivot_->unrollWeights();
     541            problemStatus_=-2; // factorize now
     542            returnCode=-2;
     543            break;
     544          }
     545        }
     546        // if stable replace in basis
     547        int updateStatus = factorization_->replaceColumn(rowArray_[2],
     548                                                         pivotRow_,
     549                                                         alpha_);
     550        if (updateStatus==1) {
     551          // slight error
     552          if (factorization_->pivots()>5) {
     553            problemStatus_=-2; // factorize now
     554            returnCode=-3;
     555          }
     556        } else if (updateStatus==2) {
     557          // major error
     558          dualRowPivot_->unrollWeights();
     559          // later we may need to unwind more e.g. fake bounds
     560          if (factorization_->pivots()) {
     561            problemStatus_=-2; // factorize now
     562            returnCode=-2;
     563            break;
     564          } else {
     565            // need to reject something
     566            char x = isColumn(sequenceOut_) ? 'C' :'R';
     567            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     568              <<x<<sequenceWithin(sequenceOut_)
     569              <<CoinMessageEol;
     570            setFlagged(sequenceOut_);
     571            lastBadIteration_ = numberIterations_; // say be more cautious
     572            rowArray_[0]->clear();
     573            rowArray_[1]->clear();
     574            columnArray_[0]->clear();
     575            // make sure dual feasible
     576            // look at all rows and columns
     577            CoinIotaN(rowArray_[0]->getIndices(),numberRows_,0);
     578            rowArray_[0]->setNumElements(numberRows_);
     579            CoinIotaN(columnArray_[0]->getIndices(),numberColumns_,0);
     580            columnArray_[0]->setNumElements(numberColumns_);
     581            double objectiveChange=0.0;
     582            updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
     583                              0.0,objectiveChange);
     584            continue;
     585          }
     586        } else if (updateStatus==3) {
     587          // out of memory
     588          // increase space if not many iterations
     589          if (factorization_->pivots()<
     590              0.5*factorization_->maximumPivots()&&
     591              factorization_->pivots()<200)
     592            factorization_->areaFactor(
     593                                       factorization_->areaFactor() * 1.1);
     594          problemStatus_=-2; // factorize now
     595        }
     596        // update primal solution
     597        if (theta_<0.0) {
     598#ifdef CLP_DEBUG
     599          if (handler_->logLevel()&32)
     600            printf("negative theta %g\n",theta_);
     601#endif
     602          theta_=0.0;
     603        }
     604        // do actual flips
     605        flipBounds(rowArray_[0],columnArray_[0],theta_);
     606        dualRowPivot_->updatePrimalSolution(rowArray_[1],
     607                                            movement,
     608                                            objectiveChange);
     609#ifdef CLP_DEBUG
     610        double oldobj=objectiveValue_;
     611#endif
     612        int whatNext=housekeeping(objectiveChange);
     613        // and set bounds correctly
     614        originalBound(sequenceIn_);
     615        changeBound(sequenceOut_);
     616#ifdef CLP_DEBUG
     617        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
     618          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
     619#endif
     620
     621#if 0 // *MERGE*
    589622          if (factorization_->pivots()<5) {
    590623            problemStatus_=-4; //say looks infeasible
     
    596629          rowArray_[0]->clear();
    597630          columnArray_[0]->clear();
     631#else // devel-1
     632        if (whatNext==1) {
     633          problemStatus_ =-2; // refactorize
     634        } else if (whatNext==2) {
     635          // maximum iterations or equivalent
     636          problemStatus_= 3;
     637          returnCode=3;
     638#endif
    598639          break;
    599640        }
    600641      } else {
    601         // no pivot row
     642        // no incoming column is valid
    602643#ifdef CLP_DEBUG
    603644        if (handler_->logLevel()&32)
    604           printf("** no row pivot\n");
    605 #endif
     645          printf("** no column pivot\n");
     646#endif
     647#if 0 // *MERGE*
    606648        if (!factorization_->pivots()) {
    607649          // may have crept through - so may be optimal
     
    625667            CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
    626668          }
    627         }
     669#else // devel-1
     670        if (factorization_->pivots()<5) {
     671          problemStatus_=-4; //say looks infeasible
     672          // create ray anyway
     673          delete [] ray_;
     674          ray_ = new double [ numberRows_];
     675          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
     676#endif
     677        }
     678        rowArray_[0]->clear();
     679        columnArray_[0]->clear();
     680        returnCode=1;
    628681        break;
    629682      }
    630     }
    631   }
    632 
    633   // at present we are leaving factorization around
    634   // maybe we should empty it
    635   deleteRim();
    636   handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
    637     <<objectiveValue()
    638     <<OsiMessageEol;
    639   // Restore any saved stuff
    640   perturbation_ = savePerturbation;
    641   dualBound_ = saveDualBound_;
    642   return problemStatus_;
     683    } else {
     684      // no pivot row
     685#ifdef CLP_DEBUG
     686      if (handler_->logLevel()&32)
     687        printf("** no row pivot\n");
     688#endif
     689      if (!factorization_->pivots()) {
     690        // may have crept through - so may be optimal
     691        //problemStatus_=-5; //say looks unbounded
     692        problemStatus_=0;
     693        // check any flagged variables
     694        int iRow;
     695        for (iRow=0;iRow<numberRows_;iRow++) {
     696          int iPivot=pivotVariable_[iRow];
     697          if (flagged(iPivot))
     698            break;
     699        }
     700        if (iRow<numberRows_) {
     701#ifdef CLP_DEBUG
     702          std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
     703#endif
     704          problemStatus_=-4; //say looks infeasible
     705          // create ray anyway
     706          delete [] ray_;
     707          ray_ = new double [ numberRows_];
     708          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
     709        }
     710      }
     711      returnCode=0;
     712      break;
     713    }
     714  }
     715  return returnCode;
    643716}
    644717/* The duals are updated by the given arrays.
    645     Returns number of infeasibilities.
    646     rowArray and columnarray will have flipped
    647     The output vector has movement (row length array) */
     718   Returns number of infeasibilities.
     719   rowArray and columnarray will have flipped
     720   The output vector has movement (row length array) */
    648721int
    649 ClpSimplexDual::updateDualsInDual(OsiIndexedVector * rowArray,
    650                                 OsiIndexedVector * columnArray,
    651                                 OsiIndexedVector * outputArray,
    652                                 double theta,
    653                                 double & objectiveChange)
     722ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
     723                                  CoinIndexedVector * columnArray,
     724                                  CoinIndexedVector * outputArray,
     725                                  double theta,
     726                                  double & objectiveChange)
    654727{
    655 
     728 
    656729  outputArray->clear();
    657730 
     
    667740                       columnArray->getNumElements()==numberColumns_);
    668741  int numberAtFake=0;
    669 
     742 
    670743  // use a tighter tolerance except for all being okay
    671744  double tolerance = dualTolerance_;
    672 
     745 
    673746  double changeObj=0.0;
    674 
     747 
    675748  int iSection;
    676749 
     
    710783        FakeBound bound = getFakeBound(iSequence+addSequence);
    711784        Status status = getStatus(iSequence+addSequence);
    712        
     785
    713786        switch(status) {
    714787         
    715         case ClpSimplex::basic:
    716         case ClpSimplex::superBasic:
     788        case basic:
     789        case superBasic:
    717790          break;
    718         case ClpSimplex::isFree:
     791        case isFree:
    719792          if (fabs(value)>tolerance) {
    720793#ifdef CLP_DEBUG
     
    725798          }
    726799          break;
    727         case ClpSimplex::atUpperBound:
     800        case atUpperBound:
    728801          if (value>tolerance) {
    729802            // to lower bound (if swap)
     
    749822                  value>=-tolerance) {
    750823                movement = lower[iSequence]-upper[iSequence];
    751                 setStatus(iSequence+addSequence,ClpSimplex::atLowerBound);
     824                setStatus(iSequence+addSequence,atLowerBound);
    752825                solution[iSequence] = lower[iSequence];
    753826                changeObj += movement*cost[iSequence];
     
    758831          }
    759832          break;
    760         case ClpSimplex::atLowerBound:
     833        case atLowerBound:
    761834          if (value<-tolerance) {
    762835            // to upper bound
     
    782855                  value<=tolerance) {
    783856                movement = upper[iSequence] - lower[iSequence];
    784                 setStatus(iSequence+addSequence,ClpSimplex::atUpperBound);
     857                setStatus(iSequence+addSequence,atUpperBound);
    785858                solution[iSequence] = upper[iSequence];
    786859                changeObj += movement*cost[iSequence];
     
    816889    // do actual flips
    817890    flipBounds(rowArray,columnArray,theta);
     891    numberFake_ = numberAtFake;
    818892  }
    819893  objectiveChange += changeObj;
     
    870944int
    871945ClpSimplexDual::changeBounds(bool initialize,
    872                                  OsiIndexedVector * outputArray,
     946                                 CoinIndexedVector * outputArray,
    873947                                 double & changeCost)
    874948{
     
    889963      setFakeBound(iSequence,ClpSimplexDual::noFake);
    890964      switch(getStatus(iSequence)) {
    891        
    892       case ClpSimplex::basic:
     965
     966      case basic:
    893967        break;
    894       case ClpSimplex::isFree:
    895       case ClpSimplex::superBasic:
     968      case isFree:
     969      case superBasic:
    896970        break;
    897       case ClpSimplex::atUpperBound:
     971      case atUpperBound:
    898972        if (fabs(value-upperValue)>primalTolerance_)
    899973          numberInfeasibilities++;
    900974        break;
    901       case ClpSimplex::atLowerBound:
     975      case atLowerBound:
    902976        if (fabs(value-lowerValue)>primalTolerance_)
    903977          numberInfeasibilities++;
     
    913987        double newUpperValue;
    914988        Status status = getStatus(iSequence);
    915         if (status==ClpSimplex::atUpperBound||
    916             status==ClpSimplex::atLowerBound) {
     989        if (status==atUpperBound||
     990            status==atLowerBound) {
    917991          double value = solution_[iSequence];
    918992          if (value-lowerValue<=upperValue-value) {
     
    9341008              setFakeBound(iSequence,ClpSimplexDual::upperFake);
    9351009          }
    936           if (status==ClpSimplex::atUpperBound)
     1010          if (status==atUpperBound)
    9371011            solution_[iSequence] = newUpperValue;
    9381012          else
     
    9621036    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    9631037      Status status = getStatus(iSequence);
    964       if (status==ClpSimplex::atUpperBound||
    965           status==ClpSimplex::atLowerBound) {
     1038      if (status==atUpperBound||
     1039          status==atLowerBound) {
    9661040        double lowerValue=lower_[iSequence];
    9671041        double upperValue=upper_[iSequence];
     
    10011075*/
    10021076void
    1003 ClpSimplexDual::dualColumn(OsiIndexedVector * rowArray,
    1004                        OsiIndexedVector * columnArray,
    1005                        OsiIndexedVector * spareArray,
    1006                        OsiIndexedVector * spareArray2)
     1077ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
     1078                       CoinIndexedVector * columnArray,
     1079                       CoinIndexedVector * spareArray,
     1080                       CoinIndexedVector * spareArray2)
    10071081{
    10081082  double * work;
     
    11191193      double value = oldValue-tentativeTheta*alpha;
    11201194      int keep = 0;
    1121        
     1195
    11221196      switch(getStatus(iSequence+addSequence)) {
    11231197         
    1124       case ClpSimplex::basic:
     1198      case basic:
    11251199        break;
    1126       case ClpSimplex::isFree:
    1127       case ClpSimplex::superBasic:
     1200      case isFree:
     1201      case superBasic:
    11281202        if (oldValue>dualTolerance_) {
    11291203          if (value<-newTolerance)
     
    11391213        }
    11401214        break;
    1141       case ClpSimplex::atUpperBound:
     1215      case atUpperBound:
     1216#ifdef CHECK_DJ
     1217        // For Debug so we can find out where problem is
     1218        perturbation_ = iSequence+addSequence;
     1219#endif
    11421220        assert (oldValue<=dualTolerance_*1.0001);
    11431221        if (value>newTolerance) {
     
    11481226        }
    11491227        break;
    1150       case ClpSimplex::atLowerBound:
     1228      case atLowerBound:
     1229#ifdef CHECK_DJ
     1230        // For Debug so we can find out where problem is
     1231        perturbation_ = iSequence+addSequence;
     1232#endif
    11511233        assert (oldValue>=-dualTolerance_*1.0001);
    11521234        if (value<-newTolerance) {
     
    14811563      oldValue = dj_[sequenceIn_];
    14821564      theta_ = oldValue/alpha;
     1565#if 0
     1566      if (numberIterations_>2000)
     1567        exit(99);
     1568      if (numberIterations_>2000-20)
     1569        handler_->setLogLevel(63);
     1570      if (numberIterations_>2000-20)
     1571        printf("theta %g oldValue %g tol %g %g\n",theta_,oldValue,dualTolerance_,
     1572               newTolerance);
     1573#endif
    14831574      if (theta_<MINIMUMTHETA) {
    14841575        // can't pivot to zero
     
    15071598            alpha=workColumn[iSequence];
    15081599          double value = dj_[iSequence]-theta_*alpha;
     1600#if 0
     1601          if (numberIterations_>2000-20)
     1602            printf("%d alpha %g value %g\n",iSequence,alpha,value);
     1603#endif
    15091604           
    15101605          // can't be free here
     
    15231618              dj_[iSequence] += modification;
    15241619              cost_[iSequence] +=  modification;
     1620#if 0
     1621              if (numberIterations_>2000-20)
     1622                printf("%d acost %g mod %g\n",iSequence,cost_[iSequence],
     1623                       modification);
     1624#endif
    15251625#endif
    15261626            }
     
    15371637              dj_[iSequence] += modification;
    15381638              cost_[iSequence] +=  modification;
     1639#if 0
     1640              if (numberIterations_>2000-20)
     1641                printf("%d bcost %g mod %g\n",iSequence,cost_[iSequence],
     1642                       modification);
     1643#endif
    15391644#endif
    15401645            }
     
    15881693  }
    15891694  if (thisIncrease) {
    1590     newTolerance = dualTolerance_+1.0e-4*dblParam_[OsiDualTolerance];
     1695    newTolerance = dualTolerance_+1.0e-4*dblParam_[ClpDualTolerance];
    15911696  }
    15921697
     
    16021707   Returns -3 if not, 2 if is unbounded */
    16031708int
    1604 ClpSimplexDual::checkUnbounded(OsiIndexedVector * ray,
    1605                                    OsiIndexedVector * spare,
     1709ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
     1710                                   CoinIndexedVector * spare,
    16061711                                   double changeCost)
    16071712{
     
    16491754    delete [] ray_;
    16501755    ray_ = new double [numberColumns_];
    1651     CoinFillN(ray_,numberColumns_,0.0);
     1756    ClpFillN(ray_,numberColumns_,0.0);
    16521757    for (i=0;i<number;i++) {
    16531758      int iRow=index[i];
     
    16631768/* Checks if finished.  Updates status */
    16641769void
    1665 ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type)
     1770ClpSimplexDual::statusOfProblemInDual(int & lastCleaned,int type,
     1771                               ClpSimplexProgress &progress)
    16661772{
    16671773  if (type==2) {
     
    16841790    // save dual weights
    16851791    dualRowPivot_->saveWeights(this,1);
    1686     // is factorization okay?
    1687     if (internalFactorize(1)) {
    1688       // no - restore previous basis
    1689       assert (type==1);
    1690       changeMade_++; // say something changed
    1691       memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    1692       memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
    1693              numberRows_*sizeof(double));
    1694       memcpy(columnActivityWork_,savedSolution_ ,
    1695              numberColumns_*sizeof(double));
    1696       // get correct bounds on all variables
    1697       double dummyChangeCost=0.0;
    1698       changeBounds(true,rowArray_[2],dummyChangeCost);
    1699       // throw away change
    1700       rowArray_[2]->clear();
    1701       forceFactorization_=1; // a bit drastic but ..
    1702       type = 2;
    1703       assert (internalFactorize(1)==0);
     1792    if (type) {
     1793      // is factorization okay?
     1794      if (internalFactorize(1)) {
     1795        // no - restore previous basis
     1796        assert (type==1);
     1797        changeMade_++; // say something changed
     1798        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
     1799        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
     1800               numberRows_*sizeof(double));
     1801        memcpy(columnActivityWork_,savedSolution_ ,
     1802               numberColumns_*sizeof(double));
     1803        // get correct bounds on all variables
     1804        double dummyChangeCost=0.0;
     1805        changeBounds(true,rowArray_[2],dummyChangeCost);
     1806        // throw away change
     1807        rowArray_[2]->clear();
     1808        forceFactorization_=1; // a bit drastic but ..
     1809        type = 2;
     1810        assert (internalFactorize(1)==0);
     1811      }
    17041812    }
    17051813    problemStatus_=-3;
     
    17081816  // get primal and dual solutions
    17091817  gutsOfSolution(rowActivityWork_,columnActivityWork_);
     1818  // Check if looping
     1819  int loop = progress.looping();
     1820  bool situationChanged=false;
     1821  if (loop>=0) {
     1822    problemStatus_ = loop; //exit if in loop
     1823    return;
     1824  } else if (loop<-1) {
     1825    // something may have changed
     1826    gutsOfSolution(rowActivityWork_,columnActivityWork_);
     1827    situationChanged=true;
     1828  }
     1829  progressFlag_ = 0; //reset progress flag
    17101830#ifdef CLP_DEBUG
    17111831  if (!rowScale_&&(handler_->logLevel()&32)) {
    17121832    double * objectiveSimplex
    1713       = copyOfArray(objective_,numberColumns_,0.0);
     1833      = ClpCopyOfArray(objective_,numberColumns_,0.0);
    17141834    double * rowObjectiveSimplex
    1715       = copyOfArray(rowObjective_,numberRows_,0.0);
     1835      = ClpCopyOfArray(rowObjective_,numberRows_,0.0);
    17161836    int i;
    17171837    double largest;
     
    17451865                       <<numberDualInfeasibilities_-
    17461866    numberDualInfeasibilitiesWithoutFree_;
    1747   handler_->message()<<OsiMessageEol;
     1867  handler_->message()<<CoinMessageEol;
     1868
    17481869  while (problemStatus_<=-3) {
    1749     bool cleanDuals=false;
     1870    bool cleanDuals=situationChanged;
    17501871    int numberChangedBounds=0;
    17511872    int doOriginalTolerance=0;
     
    17531874      doOriginalTolerance=1;
    17541875    // check optimal
     1876    // give code benefit of doubt
     1877    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
     1878        sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     1879      // say optimal (with these bounds etc)
     1880      numberDualInfeasibilities_ = 0;
     1881      sumDualInfeasibilities_ = 0.0;
     1882      numberPrimalInfeasibilities_ = 0;
     1883      sumPrimalInfeasibilities_ = 0.0;
     1884    }
    17551885    if (dualFeasible()||problemStatus_==-4) {
    17561886      if (primalFeasible()) {
     
    17581888        handler_->message(CLP_DUAL_BOUNDS,messages_)
    17591889          <<dualBound_
    1760           <<OsiMessageEol;
     1890          <<CoinMessageEol;
    17611891        // save solution in case unbounded
    1762         CoinDisjointCopyN(columnActivityWork_,numberColumns_,
     1892        ClpDisjointCopyN(columnActivityWork_,numberColumns_,
    17631893                          columnArray_[0]->denseVector());
    1764         CoinDisjointCopyN(rowActivityWork_,numberRows_,
     1894        ClpDisjointCopyN(rowActivityWork_,numberRows_,