Changeset 213


Ignore:
Timestamp:
Oct 2, 2003 2:06:08 PM (16 years ago)
Author:
forrest
Message:

ClpInterior? compiles

Location:
branches/pre
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpInterior.cpp

    r212 r213  
    1111#include "CoinHelperFunctions.hpp"
    1212#include "ClpInterior.hpp"
    13 #include "ClpFactorization.hpp"
    1413#include "ClpPackedMatrix.hpp"
    1514#include "CoinIndexedVector.hpp"
    16 #include "ClpDualRowDantzig.hpp"
    17 #include "ClpDualRowSteepest.hpp"
    18 #include "ClpPrimalColumnDantzig.hpp"
    19 #include "ClpPrimalColumnSteepest.hpp"
    20 #include "ClpNonLinearCost.hpp"
    2115#include "ClpMessage.hpp"
    2216#include "ClpLinearObjective.hpp"
     
    3125
    3226  ClpModel(),
    33   columnPrimalInfeasibility_(0.0),
    34   rowPrimalInfeasibility_(0.0),
    35   columnPrimalSequence_(-2),
    36   rowPrimalSequence_(-2),
    37   columnDualInfeasibility_(0.0),
    38   rowDualInfeasibility_(0.0),
    39   columnDualSequence_(-2),
    40   rowDualSequence_(-2),
    41   primalToleranceToGetOptimal_(-1.0),
    42   remainingDualInfeasibility_(0.0),
    43   largeValue_(1.0e15),
    4427  largestPrimalError_(0.0),
    4528  largestDualError_(0.0),
    46   largestSolutionError_(0.0),
    47   dualBound_(1.0e10),
    48   alpha_(0.0),
    49   theta_(0.0),
    50   lowerIn_(0.0),
    51   valueIn_(0.0),
    52   upperIn_(0.0),
    53   dualIn_(0.0),
    54   lowerOut_(-1),
    55   valueOut_(-1),
    56   upperOut_(-1),
    57   dualOut_(-1),
    58   dualTolerance_(0.0),
    59   primalTolerance_(0.0),
    6029  sumDualInfeasibilities_(0.0),
    6130  sumPrimalInfeasibilities_(0.0),
    62   infeasibilityCost_(1.0e10),
    63   sumOfRelaxedDualInfeasibilities_(0.0),
    64   sumOfRelaxedPrimalInfeasibilities_(0.0),
    6531  lower_(NULL),
    6632  rowLowerWork_(NULL),
     
    6935  rowUpperWork_(NULL),
    7036  columnUpperWork_(NULL),
    71   cost_(NULL),
    72   rowObjectiveWork_(NULL),
    73   objectiveWork_(NULL),
    74   sequenceIn_(-1),
    75   directionIn_(-1),
    76   sequenceOut_(-1),
    77   directionOut_(-1),
    78   pivotRow_(-1),
    79   lastGoodIteration_(-100),
    80   dj_(NULL),
    81   rowReducedCost_(NULL),
    82   reducedCostWork_(NULL),
    83   solution_(NULL),
    84   rowActivityWork_(NULL),
    85   columnActivityWork_(NULL),
    86   numberDualInfeasibilities_(0),
    87   numberDualInfeasibilitiesWithoutFree_(0),
    88   numberPrimalInfeasibilities_(100),
    89   numberRefinements_(0),
    90   pivotVariable_(NULL),
    91   factorization_(NULL),
    92   rowScale_(NULL),
    93   savedSolution_(NULL),
    94   columnScale_(NULL),
    95   scalingFlag_(3),
    96   numberTimesOptimal_(0),
    97   changeMade_(1),
    98   algorithm_(0),
    99   forceFactorization_(-1),
    100   perturbation_(100),
    101   nonLinearCost_(NULL),
    102   specialOptions_(0),
    103   lastBadIteration_(-999999),
    104   numberFake_(0),
    105   progressFlag_(0),
    106   firstFree_(-1),
    107   incomingInfeasibility_(1.0),
    108   allowedInfeasibility_(10.0),
    109   progress_(NULL)
    110 {
    111   int i;
    112   for (i=0;i<6;i++) {
    113     rowArray_[i]=NULL;
    114     columnArray_[i]=NULL;
    115   }
    116   saveStatus_=NULL;
    117   // get an empty factorization so we can set tolerances etc
    118   factorization_ = new ClpFactorization();
    119   // Say sparse
    120   factorization_->sparseThreshold(1);
    121   // say Steepest pricing
    122   dualRowPivot_ = new ClpDualRowSteepest();
    123   // say Steepest pricing
    124   primalColumnPivot_ = new ClpPrimalColumnSteepest();
    125   solveType_=1; // say simplex based life form
    126  
     37  cost_(NULL)
     38{
     39  solveType_=2; // say interior based life form
    12740}
    12841
     
    13447  : ClpModel(rhs, numberRows, whichRow,
    13548             numberColumns,whichColumn,dropNames,dropIntegers),
    136   columnPrimalInfeasibility_(0.0),
    137   rowPrimalInfeasibility_(0.0),
    138   columnPrimalSequence_(-2),
    139   rowPrimalSequence_(-2),
    140   columnDualInfeasibility_(0.0),
    141   rowDualInfeasibility_(0.0),
    142   columnDualSequence_(-2),
    143   rowDualSequence_(-2),
    144   primalToleranceToGetOptimal_(-1.0),
    145   remainingDualInfeasibility_(0.0),
    146   largeValue_(1.0e15),
    14749  largestPrimalError_(0.0),
    14850  largestDualError_(0.0),
    149   largestSolutionError_(0.0),
    150   dualBound_(1.0e10),
    151   alpha_(0.0),
    152   theta_(0.0),
    153   lowerIn_(0.0),
    154   valueIn_(0.0),
    155   upperIn_(0.0),
    156   dualIn_(0.0),
    157   lowerOut_(-1),
    158   valueOut_(-1),
    159   upperOut_(-1),
    160   dualOut_(-1),
    161   dualTolerance_(0.0),
    162   primalTolerance_(0.0),
    16351  sumDualInfeasibilities_(0.0),
    16452  sumPrimalInfeasibilities_(0.0),
    165   infeasibilityCost_(1.0e10),
    166   sumOfRelaxedDualInfeasibilities_(0.0),
    167   sumOfRelaxedPrimalInfeasibilities_(0.0),
    16853  lower_(NULL),
    16954  rowLowerWork_(NULL),
     
    17257  rowUpperWork_(NULL),
    17358  columnUpperWork_(NULL),
    174   cost_(NULL),
    175   rowObjectiveWork_(NULL),
    176   objectiveWork_(NULL),
    177   sequenceIn_(-1),
    178   directionIn_(-1),
    179   sequenceOut_(-1),
    180   directionOut_(-1),
    181   pivotRow_(-1),
    182   lastGoodIteration_(-100),
    183   dj_(NULL),
    184   rowReducedCost_(NULL),
    185   reducedCostWork_(NULL),
    186   solution_(NULL),
    187   rowActivityWork_(NULL),
    188   columnActivityWork_(NULL),
    189   numberDualInfeasibilities_(0),
    190   numberDualInfeasibilitiesWithoutFree_(0),
    191   numberPrimalInfeasibilities_(100),
    192   numberRefinements_(0),
    193   pivotVariable_(NULL),
    194   factorization_(NULL),
    195   rowScale_(NULL),
    196   savedSolution_(NULL),
    197   columnScale_(NULL),
    198   scalingFlag_(3),
    199   numberTimesOptimal_(0),
    200   changeMade_(1),
    201   algorithm_(0),
    202   forceFactorization_(-1),
    203   perturbation_(100),
    204   nonLinearCost_(NULL),
    205   specialOptions_(0),
    206   lastBadIteration_(-999999),
    207   numberFake_(0),
    208   progressFlag_(0),
    209   firstFree_(-1),
    210   incomingInfeasibility_(1.0),
    211   allowedInfeasibility_(10.0),
    212   progress_(NULL)
    213 {
    214   int i;
    215   for (i=0;i<6;i++) {
    216     rowArray_[i]=NULL;
    217     columnArray_[i]=NULL;
    218   }
    219   saveStatus_=NULL;
    220   // get an empty factorization so we can set tolerances etc
    221   factorization_ = new ClpFactorization();
    222   // say Steepest pricing
    223   dualRowPivot_ = new ClpDualRowSteepest();
    224   // say Steepest pricing
    225   primalColumnPivot_ = new ClpPrimalColumnSteepest();
    226   solveType_=1; // say simplex based life form
     59  cost_(NULL)
     60{
     61  solveType_=2; // say interior based life form
    22762 
    22863}
     
    23267ClpInterior::~ClpInterior ()
    23368{
    234   gutsOfDelete(0);
    235   delete nonLinearCost_;
     69  gutsOfDelete();
    23670}
    23771//#############################################################################
    238 void ClpInterior::setLargeValue( double value)
    239 {
    240   if (value>0.0&&value<COIN_DBL_MAX)
    241     largeValue_=value;
    242 }
    243 int
    244 ClpInterior::gutsOfSolution ( double * givenDuals,
    245                              const double * givenPrimals,
    246                              bool valuesPass)
    247 {
    248 
    249 
    250   // if values pass, save values of basic variables
    251   double * save = NULL;
    252   double oldValue=0.0;
    253   if (valuesPass) {
    254     assert(algorithm_>0); // only primal at present
    255     assert(nonLinearCost_);
    256     int iRow;
    257     checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    258     // get correct bounds on all variables
    259     nonLinearCost_->checkInfeasibilities(primalTolerance_);
    260     oldValue = nonLinearCost_->largestInfeasibility();
    261     save = new double[numberRows_];
    262     for (iRow=0;iRow<numberRows_;iRow++) {
    263       int iPivot=pivotVariable_[iRow];
    264       save[iRow] = solution_[iPivot];
    265     }
    266   }
    267   // do work
    268   computePrimals(rowActivityWork_, columnActivityWork_);
    269   // If necessary - override results
    270   if (givenPrimals) {
    271     memcpy(columnActivityWork_,givenPrimals,numberColumns_*sizeof(double));
    272     memset(rowActivityWork_,0,numberRows_*sizeof(double));
    273     times(-1.0,columnActivityWork_,rowActivityWork_);
    274   }
    275   double objectiveModification = 0.0;
    276   if (algorithm_>0&&nonLinearCost_!=NULL) {
    277     // primal algorithm
    278     // get correct bounds on all variables
    279     // If  4 bit set - Force outgoing variables to exact bound (primal)
    280     if ((specialOptions_&4)==0)
    281       nonLinearCost_->checkInfeasibilities(primalTolerance_);
    282     else
    283       nonLinearCost_->checkInfeasibilities(0.0);
    284     objectiveModification += nonLinearCost_->changeInCost();
    285     if (nonLinearCost_->numberInfeasibilities())
    286       handler_->message(CLP_SIMPLEX_NONLINEAR,messages_)
    287         <<nonLinearCost_->changeInCost()
    288         <<nonLinearCost_->numberInfeasibilities()
    289         <<CoinMessageEol;
    290   }
    291   if (valuesPass) {
    292 #ifdef CLP_DEBUG
    293     std::cout<<"Largest given infeasibility "<<oldValue
    294              <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl;
    295 #endif
    296     int numberOut=0;
    297     if (oldValue<incomingInfeasibility_
    298         &&nonLinearCost_->largestInfeasibility()>
    299         max(incomingInfeasibility_,allowedInfeasibility_)||
    300         largestPrimalError_>1.0e-3) {
    301       printf("Original largest infeas %g, now %g, primalError %g\n",
    302              oldValue,nonLinearCost_->largestInfeasibility(),
    303              largestPrimalError_);
    304       // throw out up to 1000 structurals
    305       int iRow;
    306       int * sort = new int[numberRows_];
    307       // first put back solution and store difference
    308       for (iRow=0;iRow<numberRows_;iRow++) {
    309         int iPivot=pivotVariable_[iRow];
    310         double difference = fabs(solution_[iPivot]-save[iRow]);
    311         solution_[iPivot]=save[iRow];
    312         save[iRow]=difference;
    313       }
    314       for (iRow=0;iRow<numberRows_;iRow++) {
    315         int iPivot=pivotVariable_[iRow];
    316 
    317         if (iPivot<numberColumns_) {
    318           // column
    319           double difference= save[iRow];
    320           if (difference>1.0e-4) {
    321             sort[numberOut]=iPivot;
    322             save[numberOut++]=difference;
    323           }
    324         }
    325       }
    326       CoinSort_2(save, save + numberOut, sort,
    327                  CoinFirstGreater_2<double, int>());
    328       numberOut = min(1000,numberOut);
    329       for (iRow=0;iRow<numberOut;iRow++) {
    330         int iColumn=sort[iRow];
    331         setColumnStatus(iColumn,superBasic);
    332 
    333       }
    334       delete [] sort;
    335     }
    336     delete [] save;
    337     if (numberOut)
    338       return numberOut;
    339   }
    340 
    341   computeDuals(givenDuals);
    342 
    343   // now check solutions
    344   checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    345   objectiveValue_ += objectiveModification;
    346   checkDualSolution();
    347   if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
    348                                largestDualError_>1.0e-2))
    349     handler_->message(CLP_SIMPLEX_ACCURACY,messages_)
    350       <<largestPrimalError_
    351       <<largestDualError_
    352       <<CoinMessageEol;
    353   // Switch off false values pass indicator
    354   if (!valuesPass&&algorithm_>0)
    355     firstFree_ = -1;
    356   return 0;
    357 }
    358 void
    359 ClpInterior::computePrimals ( const double * rowActivities,
    360                                      const double * columnActivities)
    361 {
    362 
    363   //work space
    364   CoinIndexedVector  * workSpace = rowArray_[0];
    365 
    366   CoinIndexedVector arrayVector;
    367   arrayVector.reserve(numberRows_+1);
    368   CoinIndexedVector previousVector;
    369   previousVector.reserve(numberRows_+1);
    370 
    371   // accumulate non basic stuff
    372 
    373   int iRow;
    374   // order is this way for scaling
    375   // Use whole matrix every time to make it easier for ClpMatrixBase
    376   // So zero out basic
    377   if (columnActivities!=columnActivityWork_)
    378     ClpDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
    379   if (rowActivities!=rowActivityWork_)
    380     ClpDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
    381   for (iRow=0;iRow<numberRows_;iRow++) {
    382     int iPivot=pivotVariable_[iRow];
    383     solution_[iPivot] = 0.0;
    384   }
    385   double * array = arrayVector.denseVector();
    386   times(-1.0,columnActivityWork_,array);
    387 
    388   int * index = arrayVector.getIndices();
    389   int number=0;
    390   for (iRow=0;iRow<numberRows_;iRow++) {
    391     double value = array[iRow] + rowActivityWork_[iRow];
    392     if (value) {
    393       array[iRow]=value;
    394       index[number++]=iRow;
    395     } else {
    396       array[iRow]=0.0;
    397     }
    398   }
    399   arrayVector.setNumElements(number);
    400 
    401   // Ftran adjusted RHS and iterate to improve accuracy
    402   double lastError=COIN_DBL_MAX;
    403   int iRefine;
    404   double * work = workSpace->denseVector();
    405   CoinIndexedVector * thisVector = &arrayVector;
    406   CoinIndexedVector * lastVector = &previousVector;
    407   factorization_->updateColumn(workSpace,thisVector);
    408   bool goodSolution=true;
    409   for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
    410 
    411     int numberIn = thisVector->getNumElements();
    412     int * indexIn = thisVector->getIndices();
    413     double * arrayIn = thisVector->denseVector();
    414     // put solution in correct place
    415     int j;
    416     for (j=0;j<numberIn;j++) {
    417       iRow = indexIn[j];
    418       int iPivot=pivotVariable_[iRow];
    419       solution_[iPivot] = arrayIn[iRow];
    420     }
    421     // check Ax == b  (for all)
    422     times(-1.0,columnActivityWork_,work);
    423     largestPrimalError_=0.0;
    424     double multiplier = 131072.0;
    425     for (iRow=0;iRow<numberRows_;iRow++) {
    426       double value = work[iRow] + rowActivityWork_[iRow];
    427       work[iRow] = value*multiplier;
    428       if (fabs(value)>largestPrimalError_) {
    429         largestPrimalError_=fabs(value);
    430       }
    431     }
    432     if (largestPrimalError_>=lastError) {
    433       // restore
    434       CoinIndexedVector * temp = thisVector;
    435       thisVector = lastVector;
    436       lastVector=temp;
    437       goodSolution=false;
    438       break;
    439     }
    440     if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-10) {
    441       // try and make better
    442       // save this
    443       CoinIndexedVector * temp = thisVector;
    444       thisVector = lastVector;
    445       lastVector=temp;
    446       int * indexOut = thisVector->getIndices();
    447       int number=0;
    448       array = thisVector->denseVector();
    449       thisVector->clear();
    450       for (iRow=0;iRow<numberRows_;iRow++) {
    451         double value = work[iRow];
    452         if (value) {
    453           array[iRow]=value;
    454           indexOut[number++]=iRow;
    455           work[iRow]=0.0;
    456         }
    457       }
    458       thisVector->setNumElements(number);
    459       lastError=largestPrimalError_;
    460       factorization_->updateColumn(workSpace,thisVector);
    461       multiplier = 1.0/multiplier;
    462       double * previous = lastVector->denseVector();
    463       number=0;
    464       for (iRow=0;iRow<numberRows_;iRow++) {
    465         double value = previous[iRow] + multiplier*array[iRow];
    466         if (value) {
    467           array[iRow]=value;
    468           indexOut[number++]=iRow;
    469         } else {
    470           array[iRow]=0.0;
    471         }
    472       }
    473       thisVector->setNumElements(number);
    474     } else {
    475       break;
    476     }
    477   }
    478 
    479   // solution as accurate as we are going to get
    480   ClpFillN(work,numberRows_,0.0);
    481   if (!goodSolution) {
    482     array = thisVector->denseVector();
    483     // put solution in correct place
    484     for (iRow=0;iRow<numberRows_;iRow++) {
    485       int iPivot=pivotVariable_[iRow];
    486       solution_[iPivot] = array[iRow];
    487     }
    488   }
    489 }
    490 // now dual side
    491 void
    492 ClpInterior::computeDuals(double * givenDjs)
    493 {
    494   //work space
    495   CoinIndexedVector  * workSpace = rowArray_[0];
    496 
    497   CoinIndexedVector arrayVector;
    498   arrayVector.reserve(numberRows_+1);
    499   CoinIndexedVector previousVector;
    500   previousVector.reserve(numberRows_+1);
    501 
    502 
    503   int iRow;
    504 #ifdef CLP_DEBUG
    505   workSpace->checkClear();
    506 #endif
    507   double * array = arrayVector.denseVector();
    508   int * index = arrayVector.getIndices();
    509   int number=0;
    510   if (!givenDjs) {
    511     for (iRow=0;iRow<numberRows_;iRow++) {
    512       int iPivot=pivotVariable_[iRow];
    513       double value = cost_[iPivot];
    514       if (value) {
    515         array[iRow]=value;
    516         index[number++]=iRow;
    517       }
    518     }
    519   } else {
    520     // dual values pass - djs may not be zero
    521     for (iRow=0;iRow<numberRows_;iRow++) {
    522       int iPivot=pivotVariable_[iRow];
    523       // make sure zero if done
    524       if (!pivoted(iPivot))
    525         givenDjs[iPivot]=0.0;
    526       double value =cost_[iPivot]-givenDjs[iPivot];
    527       if (value) {
    528         array[iRow]=value;
    529         index[number++]=iRow;
    530       }
    531     }
    532   }
    533   arrayVector.setNumElements(number);
    534 
    535   // Btran basic costs and get as accurate as possible
    536   double lastError=COIN_DBL_MAX;
    537   int iRefine;
    538   double * work = workSpace->denseVector();
    539   CoinIndexedVector * thisVector = &arrayVector;
    540   CoinIndexedVector * lastVector = &previousVector;
    541   factorization_->updateColumnTranspose(workSpace,thisVector);
    542 
    543   for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
    544     // check basic reduced costs zero
    545     largestDualError_=0.0;
    546     // would be faster to do just for basic but this reduces code
    547     ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    548     transposeTimes(-1.0,array,reducedCostWork_);
    549     if (!givenDjs) {
    550       for (iRow=0;iRow<numberRows_;iRow++) {
    551         int iPivot=pivotVariable_[iRow];
    552         double value;
    553         if (iPivot>=numberColumns_) {
    554           // slack
    555           value = rowObjectiveWork_[iPivot-numberColumns_]
    556             + array[iPivot-numberColumns_];
    557         } else {
    558           // column
    559           value = reducedCostWork_[iPivot];
    560         }
    561         work[iRow]=value;
    562         if (fabs(value)>largestDualError_) {
    563           largestDualError_=fabs(value);
    564         }
    565       }
    566     } else {
    567       for (iRow=0;iRow<numberRows_;iRow++) {
    568         int iPivot=pivotVariable_[iRow];
    569         if (iPivot>=numberColumns_) {
    570           // slack
    571           work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
    572             + array[iPivot-numberColumns_]-givenDjs[iPivot];
    573         } else {
    574           // column
    575           work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot];
    576         }
    577         if (fabs(work[iRow])>largestDualError_) {
    578           largestDualError_=fabs(work[iRow]);
    579           //assert (largestDualError_<1.0e-7);
    580           //if (largestDualError_>1.0e-7)
    581           //printf("large dual error %g\n",largestDualError_);
    582         }
    583       }
    584     }
    585     if (largestDualError_>=lastError) {
    586       // restore
    587       CoinIndexedVector * temp = thisVector;
    588       thisVector = lastVector;
    589       lastVector=temp;
    590       break;
    591     }
    592     if (iRefine<numberRefinements_&&largestDualError_>1.0e-10
    593         &&!givenDjs) {
    594       // try and make better
    595       // save this
    596       CoinIndexedVector * temp = thisVector;
    597       thisVector = lastVector;
    598       lastVector=temp;
    599       int * indexOut = thisVector->getIndices();
    600       int number=0;
    601       array = thisVector->denseVector();
    602       thisVector->clear();
    603       double multiplier = 131072.0;
    604       for (iRow=0;iRow<numberRows_;iRow++) {
    605         double value = multiplier*work[iRow];
    606         if (value) {
    607           array[iRow]=value;
    608           indexOut[number++]=iRow;
    609           work[iRow]=0.0;
    610         }
    611         work[iRow]=0.0;
    612       }
    613       thisVector->setNumElements(number);
    614       lastError=largestDualError_;
    615       factorization_->updateColumnTranspose(workSpace,thisVector);
    616       multiplier = 1.0/multiplier;
    617       double * previous = lastVector->denseVector();
    618       number=0;
    619       for (iRow=0;iRow<numberRows_;iRow++) {
    620         double value = previous[iRow] + multiplier*array[iRow];
    621         if (value) {
    622           array[iRow]=value;
    623           indexOut[number++]=iRow;
    624         } else {
    625           array[iRow]=0.0;
    626         }
    627       }
    628       thisVector->setNumElements(number);
    629     } else {
    630       break;
    631     }
    632   }
    633   ClpFillN(work,numberRows_,0.0);
    634   // now look at dual solution
    635   array = thisVector->denseVector();
    636   for (iRow=0;iRow<numberRows_;iRow++) {
    637     // slack
    638     double value = array[iRow];
    639     dual_[iRow]=value;
    640     value += rowObjectiveWork_[iRow];
    641     rowReducedCost_[iRow]=value;
    642   }
    643   ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
    644   transposeTimes(-1.0,dual_,reducedCostWork_);
    645   // If necessary - override results
    646   if (givenDjs) {
    647     // restore accurate duals
    648     memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double));
    649   }
    650 
    651 }
    652 /* Given an existing factorization computes and checks
    653    primal and dual solutions.  Uses input arrays for variables at
    654    bounds.  Returns feasibility states */
    655 int ClpInterior::getSolution ( const double * rowActivities,
    656                                const double * columnActivities)
    657 {
    658   if (!factorization_->status()) {
    659     // put in standard form
    660     createRim(7+8+16+32);
    661     // do work
    662     gutsOfSolution ( NULL,NULL);
    663     // release extra memory
    664     deleteRim(0);
    665   }
    666   return factorization_->status();
    667 }
    668 /* Given an existing factorization computes and checks
    669    primal and dual solutions.  Uses current problem arrays for
    670    bounds.  Returns feasibility states */
    671 int ClpInterior::getSolution ( )
    672 {
    673   double * rowActivities = new double[numberRows_];
    674   double * columnActivities = new double[numberColumns_];
    675   ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
    676   ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
    677   int status = getSolution( rowActivities, columnActivities);
    678   delete [] rowActivities;
    679   delete [] columnActivities;
    680   return status;
    681 }
    682 // Factorizes using current basis.  This is for external use
    683 // Return codes are as from ClpFactorization
    684 int ClpInterior::factorize ()
    685 {
    686   // put in standard form
    687   createRim(7+8+16+32,false);
    688   // do work
    689   int status = internalFactorize(-1);
    690   // release extra memory
    691   deleteRim(0);
    692 
    693   return status;
    694 }
    695 // Clean up status
    696 void
    697 ClpInterior::cleanStatus()
    698 {
    699   int iRow,iColumn;
    700   int numberBasic=0;
    701   // make row activities correct
    702   memset(rowActivityWork_,0,numberRows_*sizeof(double));
    703   times(1.0,columnActivityWork_,rowActivityWork_);
    704   if (!status_)
    705     createStatus();
    706   for (iRow=0;iRow<numberRows_;iRow++) {
    707     if (getRowStatus(iRow)==basic)
    708       numberBasic++;
    709     else {
    710       setRowStatus(iRow,superBasic);
    711       // but put to bound if close
    712       if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
    713           <=primalTolerance_) {
    714         rowActivityWork_[iRow]=rowLowerWork_[iRow];
    715         setRowStatus(iRow,atLowerBound);
    716       } else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
    717                  <=primalTolerance_) {
    718         rowActivityWork_[iRow]=rowUpperWork_[iRow];
    719         setRowStatus(iRow,atUpperBound);
    720       }
    721     }
    722   }
    723   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    724     if (getColumnStatus(iColumn)==basic) {
    725       if (numberBasic==numberRows_) {
    726         // take out of basis
    727         setColumnStatus(iColumn,superBasic);
    728         // but put to bound if close
    729         if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    730             <=primalTolerance_) {
    731           columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    732           setColumnStatus(iColumn,atLowerBound);
    733         } else if (fabs(columnActivityWork_[iColumn]
    734                         -columnUpperWork_[iColumn])
    735                    <=primalTolerance_) {
    736           columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    737           setColumnStatus(iColumn,atUpperBound);
    738         }
    739       } else
    740         numberBasic++;
    741     } else {
    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         setColumnStatus(iColumn,atLowerBound);
    748       } else if (fabs(columnActivityWork_[iColumn]
    749                       -columnUpperWork_[iColumn])
    750                  <=primalTolerance_) {
    751         columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    752         setColumnStatus(iColumn,atUpperBound);
    753       }
    754     }
    755   }
    756 }
    757 
    758 /* Factorizes using current basis. 
    759       solveType - 1 iterating, 0 initial, -1 external
    760       - 2 then iterating but can throw out of basis
    761       If 10 added then in primal values pass
    762 */
    763 /* Return codes are as from ClpFactorization unless initial factorization
    764    when total number of singularities is returned */
    765 int ClpInterior::internalFactorize ( int solveType)
    766 {
    767 
    768   int iRow,iColumn;
    769   int totalSlacks=numberRows_;
    770   if (!status_)
    771     createStatus();
    772 
    773   bool valuesPass=false;
    774   if (solveType>=10) {
    775     valuesPass=true;
    776     solveType -= 10;
    777   }
    778 #ifdef CLP_DEBUG
    779   if (solveType>0) {
    780     int numberFreeIn=0,numberFreeOut=0;
    781     double biggestDj=0.0;
    782     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    783       switch(getColumnStatus(iColumn)) {
    784 
    785       case basic:
    786         if (columnLower_[iColumn]<-largeValue_
    787             &&columnUpper_[iColumn]>largeValue_)
    788           numberFreeIn++;
    789         break;
    790       default:
    791         if (columnLower_[iColumn]<-largeValue_
    792             &&columnUpper_[iColumn]>largeValue_) {
    793           numberFreeOut++;
    794           biggestDj = max(fabs(dj_[iColumn]),biggestDj);
    795         }
    796         break;
    797       }
    798     }
    799     if (numberFreeIn+numberFreeOut)
    800       printf("%d in basis, %d out - largest dj %g\n",
    801              numberFreeIn,numberFreeOut,biggestDj);
    802   }
    803 #endif
    804   if (solveType<=0) {
    805     // Make sure everything is clean
    806     for (iRow=0;iRow<numberRows_;iRow++) {
    807       if(getRowStatus(iRow)==isFixed) {
    808         // double check fixed
    809         if (rowUpperWork_[iRow]>rowLowerWork_[iRow])
    810           setRowStatus(iRow,atLowerBound);
    811       } else if (getRowStatus(iRow)==isFree) {
    812         // may not be free after all
    813         if (rowLowerWork_[iRow]>-largeValue_||rowUpperWork_[iRow]<largeValue_)
    814           setRowStatus(iRow,superBasic);
    815       }
    816     }
    817     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    818       if(getColumnStatus(iColumn)==isFixed) {
    819         // double check fixed
    820         if (columnUpperWork_[iColumn]>columnLowerWork_[iColumn])
    821           setColumnStatus(iColumn,atLowerBound);
    822       } else if (getColumnStatus(iColumn)==isFree) {
    823         // may not be free after all
    824         if (columnLowerWork_[iColumn]>-largeValue_||columnUpperWork_[iColumn]<largeValue_)
    825           setColumnStatus(iColumn,superBasic);
    826       }
    827     }
    828     if (!valuesPass) {
    829       // not values pass so set to bounds
    830       bool allSlack=true;
    831       if (status_) {
    832         for (iRow=0;iRow<numberRows_;iRow++) {
    833           if (getRowStatus(iRow)!=basic) {
    834             allSlack=false;
    835             break;
    836           }
    837         }
    838       }
    839       if (!allSlack) {
    840         // set values from warm start (if sensible)
    841         int numberBasic=0;
    842         for (iRow=0;iRow<numberRows_;iRow++) {
    843           switch(getRowStatus(iRow)) {
    844            
    845           case basic:
    846             numberBasic++;
    847             break;
    848           case atUpperBound:
    849             rowActivityWork_[iRow]=rowUpperWork_[iRow];
    850             if (rowActivityWork_[iRow]>largeValue_) {
    851               if (rowLowerWork_[iRow]>-largeValue_) {
    852                 rowActivityWork_[iRow]=rowLowerWork_[iRow];
    853                 setRowStatus(iRow,atLowerBound);
    854               } else {
    855                 // say free
    856                 setRowStatus(iRow,isFree);
    857                 rowActivityWork_[iRow]=0.0;
    858               }
    859             }
    860             break;
    861           case ClpInterior::isFixed:
    862           case atLowerBound:
    863             rowActivityWork_[iRow]=rowLowerWork_[iRow];
    864             if (rowActivityWork_[iRow]<-largeValue_) {
    865               if (rowUpperWork_[iRow]<largeValue_) {
    866                 rowActivityWork_[iRow]=rowUpperWork_[iRow];
    867                 setRowStatus(iRow,atUpperBound);
    868               } else {
    869                 // say free
    870                 setRowStatus(iRow,isFree);
    871                 rowActivityWork_[iRow]=0.0;
    872               }
    873             }
    874             break;
    875           case isFree:
    876               break;
    877             // not really free - fall through to superbasic
    878           case superBasic:
    879             if (rowUpperWork_[iRow]>largeValue_) {
    880               if (rowLowerWork_[iRow]>-largeValue_) {
    881                 rowActivityWork_[iRow]=rowLowerWork_[iRow];
    882                 setRowStatus(iRow,atLowerBound);
    883               } else {
    884                 // say free
    885                 setRowStatus(iRow,isFree);
    886                 rowActivityWork_[iRow]=0.0;
    887               }
    888             } else {
    889               if (rowLowerWork_[iRow]>-largeValue_) {
    890                 // set to nearest
    891                 if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
    892                     <fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])) {
    893                   rowActivityWork_[iRow]=rowLowerWork_[iRow];
    894                   setRowStatus(iRow,atLowerBound);
    895                 } else {
    896                   rowActivityWork_[iRow]=rowUpperWork_[iRow];
    897                   setRowStatus(iRow,atUpperBound);
    898                 }
    899               } else {
    900                 rowActivityWork_[iRow]=rowUpperWork_[iRow];
    901                 setRowStatus(iRow,atUpperBound);
    902               }
    903             }
    904             break;
    905           }
    906         }
    907         totalSlacks=numberBasic;
    908 
    909         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    910           switch(getColumnStatus(iColumn)) {
    911            
    912           case basic:
    913             if (numberBasic==numberRows_) {
    914               // take out of basis
    915               if (columnLowerWork_[iColumn]>-largeValue_) {
    916                 if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]<
    917                     columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
    918                   columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    919                   setColumnStatus(iColumn,atLowerBound);
    920                 } else {
    921                   columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    922                   setColumnStatus(iColumn,atUpperBound);
    923                 }
    924               } else if (columnUpperWork_[iColumn]<largeValue_) {
    925                 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    926                 setColumnStatus(iColumn,atUpperBound);
    927               } else {
    928                 columnActivityWork_[iColumn]=0.0;
    929                 setColumnStatus(iColumn,isFree);
    930               }
    931             } else {
    932               numberBasic++;
    933             }
    934             break;
    935           case atUpperBound:
    936             columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    937             if (columnActivityWork_[iColumn]>largeValue_) {
    938               if (columnLowerWork_[iColumn]<-largeValue_) {
    939                 columnActivityWork_[iColumn]=0.0;
    940                 setColumnStatus(iColumn,isFree);
    941               } else {
    942                 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    943                 setColumnStatus(iColumn,atLowerBound);
    944               }
    945             }
    946             break;
    947           case isFixed:
    948           case atLowerBound:
    949             columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    950             if (columnActivityWork_[iColumn]<-largeValue_) {
    951               if (columnUpperWork_[iColumn]>largeValue_) {
    952                 columnActivityWork_[iColumn]=0.0;
    953                 setColumnStatus(iColumn,isFree);
    954               } else {
    955                 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    956                 setColumnStatus(iColumn,atUpperBound);
    957               }
    958             }
    959             break;
    960           case isFree:
    961               break;
    962             // not really free - fall through to superbasic
    963           case superBasic:
    964             if (columnUpperWork_[iColumn]>largeValue_) {
    965               if (columnLowerWork_[iColumn]>-largeValue_) {
    966                 columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    967                 setColumnStatus(iColumn,atLowerBound);
    968               } else {
    969                 // say free
    970                 setColumnStatus(iColumn,isFree);
    971                 columnActivityWork_[iColumn]=0.0;
    972               }
    973             } else {
    974               if (columnLowerWork_[iColumn]>-largeValue_) {
    975                 // set to nearest
    976                 if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    977                     <fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])) {
    978                   columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    979                   setColumnStatus(iColumn,atLowerBound);
    980                 } else {
    981                   columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    982                   setColumnStatus(iColumn,atUpperBound);
    983                 }
    984               } else {
    985                 columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    986                 setColumnStatus(iColumn,atUpperBound);
    987               }
    988             }
    989             break;
    990           }
    991         }
    992       } else {
    993         // all slack basis
    994         int numberBasic=0;
    995         if (!status_) {
    996           createStatus();
    997         }
    998         for (iRow=0;iRow<numberRows_;iRow++) {
    999           double lower=rowLowerWork_[iRow];
    1000           double upper=rowUpperWork_[iRow];
    1001           if (lower>-largeValue_||upper<largeValue_) {
    1002             if (fabs(lower)<=fabs(upper)) {
    1003               rowActivityWork_[iRow]=lower;
    1004             } else {
    1005               rowActivityWork_[iRow]=upper;
    1006             }
    1007           } else {
    1008             rowActivityWork_[iRow]=0.0;
    1009           }
    1010           setRowStatus(iRow,basic);
    1011           numberBasic++;
    1012         }
    1013         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1014           double lower=columnLowerWork_[iColumn];
    1015           double upper=columnUpperWork_[iColumn];
    1016           double big_bound = largeValue_;
    1017           if (lower>-big_bound||upper<big_bound) {
    1018             if ((getColumnStatus(iColumn)==atLowerBound&&
    1019                  columnActivityWork_[iColumn]==lower)||
    1020                 (getColumnStatus(iColumn)==atUpperBound&&
    1021                  columnActivityWork_[iColumn]==upper)) {
    1022               // status looks plausible
    1023             } else {
    1024               // set to sensible
    1025               if (fabs(lower)<=fabs(upper)) {
    1026                 setColumnStatus(iColumn,atLowerBound);
    1027                 columnActivityWork_[iColumn]=lower;
    1028               } else {
    1029                 setColumnStatus(iColumn,atUpperBound);
    1030                 columnActivityWork_[iColumn]=upper;
    1031               }
    1032             }
    1033           } else {
    1034             setColumnStatus(iColumn,isFree);
    1035             columnActivityWork_[iColumn]=0.0;
    1036           }
    1037         }
    1038       }
    1039     } else {
    1040       // values pass has less coding
    1041       // make row activities correct and clean basis a bit
    1042       cleanStatus();
    1043       if (status_) {
    1044         int numberBasic=0;
    1045         for (iRow=0;iRow<numberRows_;iRow++) {
    1046           if (getRowStatus(iRow)==basic)
    1047             numberBasic++;
    1048         }
    1049         totalSlacks=numberBasic;
    1050         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1051           if (getColumnStatus(iColumn)==basic)
    1052             numberBasic++;
    1053         }
    1054       } else {
    1055         // all slack basis
    1056         int numberBasic=0;
    1057         if (!status_) {
    1058           createStatus();
    1059         }
    1060         for (iRow=0;iRow<numberRows_;iRow++) {
    1061           setRowStatus(iRow,basic);
    1062           numberBasic++;
    1063         }
    1064         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1065           setColumnStatus(iColumn,superBasic);
    1066           // but put to bound if close
    1067           if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    1068               <=primalTolerance_) {
    1069             columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    1070             setColumnStatus(iColumn,atLowerBound);
    1071           } else if (fabs(columnActivityWork_[iColumn]
    1072                         -columnUpperWork_[iColumn])
    1073                      <=primalTolerance_) {
    1074             columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
    1075             setColumnStatus(iColumn,atUpperBound);
    1076           }
    1077         }
    1078       }
    1079     }
    1080     numberRefinements_=1;
    1081     // set fixed if they are
    1082     for (iRow=0;iRow<numberRows_;iRow++) {
    1083       if (getRowStatus(iRow)!=basic ) {
    1084         if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {
    1085           rowActivityWork_[iRow]=rowLowerWork_[iRow];
    1086           setRowStatus(iRow,isFixed);
    1087         }
    1088       }
    1089     }
    1090     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1091       if (getColumnStatus(iColumn)!=basic ) {
    1092         if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {
    1093           columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    1094           setColumnStatus(iColumn,isFixed);
    1095         }
    1096       }
    1097     }
    1098   }
    1099   if (0)  {
    1100     static int k=0;
    1101     printf("start basis\n");
    1102     int i;
    1103     for (i=0;i<numberRows_;i++)
    1104       printf ("xx %d %d\n",i,pivotVariable_[i]);
    1105     for (i=0;i<numberRows_+numberColumns_;i++)
    1106       if (getColumnStatus(i)==basic)
    1107         printf ("yy %d basic\n",i);
    1108     if (k>20)
    1109       exit(0);
    1110     k++;
    1111   }
    1112   int status = factorization_->factorize(this, solveType,valuesPass);
    1113   if (status) {
    1114     handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
    1115       <<status
    1116       <<CoinMessageEol;
    1117     return -1;
    1118   } else if (!solveType) {
    1119     // Initial basis - return number of singularities
    1120     int numberSlacks=0;
    1121     for (iRow=0;iRow<numberRows_;iRow++) {
    1122       if (getRowStatus(iRow) == basic)
    1123         numberSlacks++;
    1124     }
    1125     status= max(numberSlacks-totalSlacks,0);
    1126   }
    1127 
    1128   // sparse methods
    1129   //if (factorization_->sparseThreshold()) {
    1130     // get default value
    1131     factorization_->sparseThreshold(0);
    1132     factorization_->goSparse();
    1133     //}
    1134  
    1135   return status;
    1136 }
    113772/*
    1138    This does basis housekeeping and does values for in/out variables.
    1139    Can also decide to re-factorize
     73   This does housekeeping
    114074*/
    114175int
    1142 ClpInterior::housekeeping(double objectiveChange)
     76ClpInterior::housekeeping()
    114377{
    114478  numberIterations_++;
    1145   changeMade_++; // something has happened
    1146   // incoming variable
    1147   handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
    1148     <<directionOut_
    1149     <<directionIn_<<theta_
    1150     <<dualOut_<<dualIn_<<alpha_
    1151     <<CoinMessageEol;
    1152   if (getStatus(sequenceIn_)==isFree) {
    1153     handler_->message(CLP_SIMPLEX_FREEIN,messages_)
    1154       <<sequenceIn_
    1155       <<CoinMessageEol;
    1156   }
    1157   // change of incoming
    1158   char rowcol[]={'R','C'};
    1159   if (pivotRow_>=0)
    1160     pivotVariable_[pivotRow_]=sequenceIn();
    1161   if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<-1.0e20)
    1162     progressFlag_ |= 2; // making real progress
    1163   solution_[sequenceIn_]=valueIn_;
    1164   if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
    1165     progressFlag_ |= 1; // making real progress
    1166   if (sequenceIn_!=sequenceOut_) {
    1167     //assert( getStatus(sequenceOut_)== basic);
    1168     setStatus(sequenceIn_,basic);
    1169     if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) {
    1170       // As Nonlinear costs may have moved bounds (to more feasible)
    1171       // Redo using value
    1172       if (fabs(valueOut_-lower_[sequenceOut_])<fabs(valueOut_-upper_[sequenceOut_])) {
    1173         // going to lower
    1174         setStatus(sequenceOut_,atLowerBound);
    1175       } else {
    1176         // going to upper
    1177         setStatus(sequenceOut_,atUpperBound);
    1178       }
    1179     } else {
    1180       // fixed
    1181       setStatus(sequenceOut_,isFixed);
    1182     }
    1183     solution_[sequenceOut_]=valueOut_;
    1184   } else {
    1185     // flip from bound to bound
    1186     // As Nonlinear costs may have moved bounds (to more feasible)
    1187     // Redo using value
    1188     if (fabs(valueIn_-lower_[sequenceIn_])<fabs(valueIn_-upper_[sequenceIn_])) {
    1189       // as if from upper bound
    1190       setStatus(sequenceIn_, atLowerBound);
    1191     } else {
    1192       // as if from lower bound
    1193       setStatus(sequenceIn_, atUpperBound);
    1194     }
    1195   }
    1196   objectiveValue_ += objectiveChange;
    1197   handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
    1198     <<numberIterations_<<objectiveValue()
    1199     <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
    1200     <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
    1201   handler_->printing(algorithm_<0)<<theta_<<dualOut_;
    1202   handler_->printing(algorithm_>0)<<dualIn_<<theta_;
    1203   handler_->message()<<CoinMessageEol;
    1204   if (hitMaximumIterations())
    1205     return 2;
    1206 #if 0
    1207   // check for small cycles
    1208   int cycle=progress_->cycle(sequenceIn_,sequenceOut_,
    1209                             directionIn_,directionOut_);
    1210   if (cycle>0) {
    1211     printf("Cycle of %d\n",cycle);
    1212     if (factorization_->pivots()>cycle)
    1213       forceFactorization_=cycle-1;
    1214     else
    1215       forceFactorization_=1;
    1216     return 1;
    1217   }
    1218 #endif
    1219   // only time to re-factorize if one before real time
    1220   // this is so user won't be surprised that maximumPivots has exact meaning
    1221   if (factorization_->pivots()==factorization_->maximumPivots()) {
    1222     return 1;
    1223   } else {
    1224     if (forceFactorization_>0&&
    1225         factorization_->pivots()==forceFactorization_) {
    1226       // relax
    1227       forceFactorization_ = (3+5*forceFactorization_)/4;
    1228       if (forceFactorization_>factorization_->maximumPivots())
    1229         forceFactorization_ = -1; //off
    1230       return 1;
    1231     } else {
    1232       // carry on iterating
    1233       return 0;
    1234     }
    1235   }
     79  return 0;
    123680}
    123781// Copy constructor.
    123882ClpInterior::ClpInterior(const ClpInterior &rhs) :
    123983  ClpModel(rhs),
    1240   columnPrimalInfeasibility_(0.0),
    1241   rowPrimalInfeasibility_(0.0),
    1242   columnPrimalSequence_(-2),
    1243   rowPrimalSequence_(-2),
    1244   columnDualInfeasibility_(0.0),
    1245   rowDualInfeasibility_(0.0),
    1246   columnDualSequence_(-2),
    1247   rowDualSequence_(-2),
    1248   primalToleranceToGetOptimal_(-1.0),
    1249   remainingDualInfeasibility_(0.0),
    1250   largeValue_(1.0e15),
    125184  largestPrimalError_(0.0),
    125285  largestDualError_(0.0),
    1253   largestSolutionError_(0.0),
    1254   dualBound_(1.0e10),
    1255   alpha_(0.0),
    1256   theta_(0.0),
    1257   lowerIn_(0.0),
    1258   valueIn_(0.0),
    1259   upperIn_(0.0),
    1260   dualIn_(0.0),
    1261   lowerOut_(-1),
    1262   valueOut_(-1),
    1263   upperOut_(-1),
    1264   dualOut_(-1),
    1265   dualTolerance_(0.0),
    1266   primalTolerance_(0.0),
    126786  sumDualInfeasibilities_(0.0),
    126887  sumPrimalInfeasibilities_(0.0),
    1269   infeasibilityCost_(1.0e10),
    1270   sumOfRelaxedDualInfeasibilities_(0.0),
    1271   sumOfRelaxedPrimalInfeasibilities_(0.0),
    127288  lower_(NULL),
    127389  rowLowerWork_(NULL),
     
    127692  rowUpperWork_(NULL),
    127793  columnUpperWork_(NULL),
    1278   cost_(NULL),
    1279   rowObjectiveWork_(NULL),
    1280   objectiveWork_(NULL),
    1281   sequenceIn_(-1),
    1282   directionIn_(-1),
    1283   sequenceOut_(-1),
    1284   directionOut_(-1),
    1285   pivotRow_(-1),
    1286   lastGoodIteration_(-100),
    1287   dj_(NULL),
    1288   rowReducedCost_(NULL),
    1289   reducedCostWork_(NULL),
    1290   solution_(NULL),
    1291   rowActivityWork_(NULL),
    1292   columnActivityWork_(NULL),
    1293   numberDualInfeasibilities_(0),
    1294   numberDualInfeasibilitiesWithoutFree_(0),
    1295   numberPrimalInfeasibilities_(100),
    1296   numberRefinements_(0),
    1297   pivotVariable_(NULL),
    1298   factorization_(NULL),
    1299   rowScale_(NULL),
    1300   savedSolution_(NULL),
    1301   columnScale_(NULL),
    1302   scalingFlag_(3),
    1303   numberTimesOptimal_(0),
    1304   changeMade_(1),
    1305   algorithm_(0),
    1306   forceFactorization_(-1),
    1307   perturbation_(100),
    1308   nonLinearCost_(NULL),
    1309   specialOptions_(0),
    1310   lastBadIteration_(-999999),
    1311   numberFake_(0),
    1312   progressFlag_(0),
    1313   firstFree_(-1),
    1314   incomingInfeasibility_(1.0),
    1315   allowedInfeasibility_(10.0),
    1316   progress_(NULL)
    1317 {
    1318   int i;
    1319   for (i=0;i<6;i++) {
    1320     rowArray_[i]=NULL;
    1321     columnArray_[i]=NULL;
    1322   }
    1323   saveStatus_=NULL;
    1324   factorization_ = NULL;
    1325   dualRowPivot_ = NULL;
    1326   primalColumnPivot_ = NULL;
    1327   gutsOfDelete(0);
    1328   specialOptions_ =0;
    1329   delete nonLinearCost_;
    1330   nonLinearCost_ = NULL;
     94  cost_(NULL)
     95{
     96  gutsOfDelete();
    133197  gutsOfCopy(rhs);
    1332   solveType_=1; // say simplex based life form
     98  solveType_=2; // say interior based life form
    133399}
    1334100// Copy constructor from model
    1335101ClpInterior::ClpInterior(const ClpModel &rhs) :
    1336102  ClpModel(rhs),
    1337   columnPrimalInfeasibility_(0.0),
    1338   rowPrimalInfeasibility_(0.0),
    1339   columnPrimalSequence_(-2),
    1340   rowPrimalSequence_(-2),
    1341   columnDualInfeasibility_(0.0),
    1342   rowDualInfeasibility_(0.0),
    1343   columnDualSequence_(-2),
    1344   rowDualSequence_(-2),
    1345   primalToleranceToGetOptimal_(-1.0),
    1346   remainingDualInfeasibility_(0.0),
    1347   largeValue_(1.0e15),
    1348103  largestPrimalError_(0.0),
    1349104  largestDualError_(0.0),
    1350   largestSolutionError_(0.0),
    1351   dualBound_(1.0e10),
    1352   alpha_(0.0),
    1353   theta_(0.0),
    1354   lowerIn_(0.0),
    1355   valueIn_(0.0),
    1356   upperIn_(0.0),
    1357   dualIn_(0.0),
    1358   lowerOut_(-1),
    1359   valueOut_(-1),
    1360   upperOut_(-1),
    1361   dualOut_(-1),
    1362   dualTolerance_(0.0),
    1363   primalTolerance_(0.0),
    1364105  sumDualInfeasibilities_(0.0),
    1365106  sumPrimalInfeasibilities_(0.0),
    1366   infeasibilityCost_(1.0e10),
    1367   sumOfRelaxedDualInfeasibilities_(0.0),
    1368   sumOfRelaxedPrimalInfeasibilities_(0.0),
    1369107  lower_(NULL),
    1370108  rowLowerWork_(NULL),
     
    1373111  rowUpperWork_(NULL),
    1374112  columnUpperWork_(NULL),
    1375   cost_(NULL),
    1376   rowObjectiveWork_(NULL),
    1377   objectiveWork_(NULL),
    1378   sequenceIn_(-1),
    1379   directionIn_(-1),
    1380   sequenceOut_(-1),
    1381   directionOut_(-1),
    1382   pivotRow_(-1),
    1383   lastGoodIteration_(-100),
    1384   dj_(NULL),
    1385   rowReducedCost_(NULL),
    1386   reducedCostWork_(NULL),
    1387   solution_(NULL),
    1388   rowActivityWork_(NULL),
    1389   columnActivityWork_(NULL),
    1390   numberDualInfeasibilities_(0),
    1391   numberDualInfeasibilitiesWithoutFree_(0),
    1392   numberPrimalInfeasibilities_(100),
    1393   numberRefinements_(0),
    1394   pivotVariable_(NULL),
    1395   factorization_(NULL),
    1396   rowScale_(NULL),
    1397   savedSolution_(NULL),
    1398   columnScale_(NULL),
    1399   scalingFlag_(3),
    1400   numberTimesOptimal_(0),
    1401   changeMade_(1),
    1402   algorithm_(0),
    1403   forceFactorization_(-1),
    1404   perturbation_(100),
    1405   nonLinearCost_(NULL),
    1406   specialOptions_(0),
    1407   lastBadIteration_(-999999),
    1408   numberFake_(0),
    1409   progressFlag_(0),
    1410   firstFree_(-1),
    1411   incomingInfeasibility_(1.0),
    1412   allowedInfeasibility_(10.0),
    1413   progress_(NULL)
    1414 {
    1415   int i;
    1416   for (i=0;i<6;i++) {
    1417     rowArray_[i]=NULL;
    1418     columnArray_[i]=NULL;
    1419   }
    1420   saveStatus_=NULL;
    1421   // get an empty factorization so we can set tolerances etc
    1422   factorization_ = new ClpFactorization();
    1423   // say Steepest pricing
    1424   dualRowPivot_ = new ClpDualRowSteepest();
    1425   // say Steepest pricing
    1426   primalColumnPivot_ = new ClpPrimalColumnSteepest();
    1427   solveType_=1; // say simplex based life form
    1428  
     113  cost_(NULL)
     114{
     115  solveType_=2; // say interior based life form
    1429116}
    1430117// Assignment operator. This copies the data
     
    1433120{
    1434121  if (this != &rhs) {
    1435     gutsOfDelete(0);
    1436     specialOptions_=0;
    1437     delete nonLinearCost_;
    1438     nonLinearCost_ = NULL;
     122    gutsOfDelete();
    1439123    ClpModel::operator=(rhs);
    1440124    gutsOfCopy(rhs);
     
    1452136  columnUpperWork_ = upper_;
    1453137  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
    1454   cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows_));
    1455   objectiveWork_ = cost_;
    1456   rowObjectiveWork_ = cost_+numberColumns_;
    1457   dj_ = ClpCopyOfArray(rhs.dj_,numberRows_+numberColumns_);
    1458   if (dj_) {
    1459     reducedCostWork_ = dj_;
    1460     rowReducedCost_ = dj_+numberColumns_;
    1461   }
    1462   solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
    1463   if (solution_) {
    1464     columnActivityWork_ = solution_;
    1465     rowActivityWork_ = solution_+numberColumns_;
    1466   }
    1467   if (rhs.pivotVariable_) {
    1468     pivotVariable_ = new int[numberRows_];
    1469     ClpDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
    1470   } else {
    1471     pivotVariable_=NULL;
    1472   }
    1473   if (rhs.factorization_) {
    1474     factorization_ = new ClpFactorization(*rhs.factorization_);
    1475   } else {
    1476     factorization_=NULL;
    1477   }
    1478   rowScale_ = ClpCopyOfArray(rhs.rowScale_,numberRows_);
    1479   savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
    1480   columnScale_ = ClpCopyOfArray(rhs.columnScale_,numberColumns_);
    1481   int i;
    1482   for (i=0;i<6;i++) {
    1483     rowArray_[i]=NULL;
    1484     if (rhs.rowArray_[i])
    1485       rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
    1486     columnArray_[i]=NULL;
    1487     if (rhs.columnArray_[i])
    1488       columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
    1489   }
    1490   if (rhs.saveStatus_) {
    1491     saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
    1492   }
    1493   columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
    1494   columnPrimalSequence_ = rhs.columnPrimalSequence_;
    1495   rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
    1496   rowPrimalSequence_ = rhs.rowPrimalSequence_;
    1497   columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
    1498   columnDualSequence_ = rhs.columnDualSequence_;
    1499   rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
    1500   rowDualSequence_ = rhs.rowDualSequence_;
    1501   primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
    1502   remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
    1503   largeValue_ = rhs.largeValue_;
     138  cost_ = ClpCopyOfArray(rhs.cost_,numberColumns_);
    1504139  largestPrimalError_ = rhs.largestPrimalError_;
    1505140  largestDualError_ = rhs.largestDualError_;
    1506   largestSolutionError_ = rhs.largestSolutionError_;
    1507   dualBound_ = rhs.dualBound_;
    1508   alpha_ = rhs.alpha_;
    1509   theta_ = rhs.theta_;
    1510   lowerIn_ = rhs.lowerIn_;
    1511   valueIn_ = rhs.valueIn_;
    1512   upperIn_ = rhs.upperIn_;
    1513   dualIn_ = rhs.dualIn_;
    1514   sequenceIn_ = rhs.sequenceIn_;
    1515   directionIn_ = rhs.directionIn_;
    1516   lowerOut_ = rhs.lowerOut_;
    1517   valueOut_ = rhs.valueOut_;
    1518   upperOut_ = rhs.upperOut_;
    1519   dualOut_ = rhs.dualOut_;
    1520   sequenceOut_ = rhs.sequenceOut_;
    1521   directionOut_ = rhs.directionOut_;
    1522   pivotRow_ = rhs.pivotRow_;
    1523   lastGoodIteration_ = rhs.lastGoodIteration_;
    1524   numberRefinements_ = rhs.numberRefinements_;
    1525   dualTolerance_ = rhs.dualTolerance_;
    1526   primalTolerance_ = rhs.primalTolerance_;
    1527141  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
    1528   numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    1529   numberDualInfeasibilitiesWithoutFree_ =
    1530     rhs.numberDualInfeasibilitiesWithoutFree_;
    1531142  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    1532   numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    1533   dualRowPivot_ = rhs.dualRowPivot_->clone(true);
    1534   primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
    1535   scalingFlag_ = rhs.scalingFlag_;
    1536   numberTimesOptimal_ = rhs.numberTimesOptimal_;
    1537   changeMade_ = rhs.changeMade_;
    1538   algorithm_ = rhs.algorithm_;
    1539   forceFactorization_ = rhs.forceFactorization_;
    1540   perturbation_ = rhs.perturbation_;
    1541   infeasibilityCost_ = rhs.infeasibilityCost_;
    1542   specialOptions_ = rhs.specialOptions_;
    1543   lastBadIteration_ = rhs.lastBadIteration_;
    1544   numberFake_ = rhs.numberFake_;
    1545   progressFlag_ = rhs.progressFlag_;
    1546   firstFree_ = rhs.firstFree_;
    1547   incomingInfeasibility_ = rhs.incomingInfeasibility_;
    1548   allowedInfeasibility_ = rhs.allowedInfeasibility_;
    1549   if (rhs.progress_)
    1550     progress_ = new ClpInteriorProgress(*rhs.progress_);
    1551   else
    1552     progress_=NULL;
    1553   sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    1554   sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    1555   if (rhs.nonLinearCost_!=NULL)
    1556     nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
    1557   else
    1558     nonLinearCost_=NULL;
    1559143  solveType_=rhs.solveType_;
    1560144}
    1561145// type == 0 do everything, most + pivot data, 2 factorization data as well
    1562146void
    1563 ClpInterior::gutsOfDelete(int type)
     147ClpInterior::gutsOfDelete()
    1564148{
    1565149  delete [] lower_;
     
    1573157  delete [] cost_;
    1574158  cost_=NULL;
    1575   objectiveWork_=NULL;
    1576   rowObjectiveWork_=NULL;
    1577   delete [] dj_;
    1578   dj_=NULL;
    1579   reducedCostWork_=NULL;
    1580   rowReducedCost_=NULL;
    1581   delete [] solution_;
    1582   solution_=NULL;
    1583   rowActivityWork_=NULL;
    1584   columnActivityWork_=NULL;
    1585   delete [] rowScale_;
    1586   rowScale_ = NULL;
    1587   delete [] savedSolution_;
    1588   savedSolution_ = NULL;
    1589   delete [] columnScale_;
    1590   columnScale_ = NULL;
    1591   if ((specialOptions_&2)==0) {
    1592     delete nonLinearCost_;
    1593     nonLinearCost_ = NULL;
    1594   }
    1595   int i;
    1596   for (i=0;i<6;i++) {
    1597     delete rowArray_[i];
    1598     rowArray_[i]=NULL;
    1599     delete columnArray_[i];
    1600     columnArray_[i]=NULL;
    1601   }
    1602   delete rowCopy_;
    1603   rowCopy_=NULL;
    1604   delete [] saveStatus_;
    1605   saveStatus_=NULL;
    1606   if (!type) {
    1607     // delete everything
    1608     delete factorization_;
    1609     factorization_ = NULL;
    1610     delete [] pivotVariable_;
    1611     pivotVariable_=NULL;
    1612     delete dualRowPivot_;
    1613     dualRowPivot_ = NULL;
    1614     delete primalColumnPivot_;
    1615     primalColumnPivot_ = NULL;
    1616     delete progress_;
    1617     progress_=NULL;
    1618   } else {
    1619     // delete any size information in methods
    1620     if (type>1) {
    1621       factorization_->clearArrays();
    1622       delete [] pivotVariable_;
    1623       pivotVariable_=NULL;
    1624     }
    1625     dualRowPivot_->clearArrays();
    1626     primalColumnPivot_->clearArrays();
    1627   }
    1628 }
    1629 // This sets largest infeasibility and most infeasible
    1630 void
    1631 ClpInterior::checkPrimalSolution(const double * rowActivities,
    1632                                         const double * columnActivities)
    1633 {
    1634   double * solution;
    1635   int iRow,iColumn;
    1636 
    1637   objectiveValue_ = 0.0;
    1638   // now look at primal solution
    1639   columnPrimalInfeasibility_=0.0;
    1640   columnPrimalSequence_=-1;
    1641   rowPrimalInfeasibility_=0.0;
    1642   rowPrimalSequence_=-1;
    1643   largestSolutionError_=0.0;
    1644   solution = rowActivityWork_;
    1645   sumPrimalInfeasibilities_=0.0;
    1646   numberPrimalInfeasibilities_=0;
    1647   double primalTolerance = primalTolerance_;
    1648   double relaxedTolerance=primalTolerance_;
    1649   // we can't really trust infeasibilities if there is primal error
    1650   double error = min(1.0e-3,largestPrimalError_);
    1651   // allow tolerance at least slightly bigger than standard
    1652   relaxedTolerance = relaxedTolerance +  error;
    1653   sumOfRelaxedPrimalInfeasibilities_ = 0.0;
    1654 
    1655   for (iRow=0;iRow<numberRows_;iRow++) {
    1656     //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
    1657     double infeasibility=0.0;
    1658     objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
    1659     if (solution[iRow]>rowUpperWork_[iRow]) {
    1660       infeasibility=solution[iRow]-rowUpperWork_[iRow];
    1661     } else if (solution[iRow]<rowLowerWork_[iRow]) {
    1662       infeasibility=rowLowerWork_[iRow]-solution[iRow];
    1663     }
    1664     if (infeasibility>primalTolerance) {
    1665       sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
    1666       if (infeasibility>relaxedTolerance)
    1667         sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
    1668       numberPrimalInfeasibilities_ ++;
    1669     }
    1670     if (infeasibility>rowPrimalInfeasibility_) {
    1671       rowPrimalInfeasibility_=infeasibility;
    1672       rowPrimalSequence_=iRow;
    1673     }
    1674     infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
    1675     if (infeasibility>largestSolutionError_)
    1676       largestSolutionError_=infeasibility;
    1677   }
    1678   solution = columnActivityWork_;
    1679   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1680     //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
    1681     double infeasibility=0.0;
    1682     objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
    1683     if (solution[iColumn]>columnUpperWork_[iColumn]) {
    1684       infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
    1685     } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
    1686       infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
    1687     }
    1688     if (infeasibility>columnPrimalInfeasibility_) {
    1689       columnPrimalInfeasibility_=infeasibility;
    1690       columnPrimalSequence_=iColumn;
    1691     }
    1692     if (infeasibility>primalTolerance) {
    1693       sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
    1694       if (infeasibility>relaxedTolerance)
    1695         sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
    1696       numberPrimalInfeasibilities_ ++;
    1697     }
    1698     infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
    1699     if (infeasibility>largestSolutionError_)
    1700       largestSolutionError_=infeasibility;
    1701   }
    1702 }
    1703 void
    1704 ClpInterior::checkDualSolution()
    1705 {
    1706 
    1707   int iRow,iColumn;
    1708   sumDualInfeasibilities_=0.0;
    1709   numberDualInfeasibilities_=0;
    1710   numberDualInfeasibilitiesWithoutFree_=0;
    1711   columnDualInfeasibility_=0.0;
    1712   columnDualSequence_=-1;
    1713   rowDualInfeasibility_=0.0;
    1714   rowDualSequence_=-1;
    1715   int firstFreePrimal = -1;
    1716   int firstFreeDual = -1;
    1717   int numberSuperBasicWithDj=0;
    1718   primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
    1719                                    columnPrimalInfeasibility_);
    1720   remainingDualInfeasibility_=0.0;
    1721   double relaxedTolerance=dualTolerance_;
    1722   // we can't really trust infeasibilities if there is dual error
    1723   double error = min(1.0e-3,largestDualError_);
    1724   // allow tolerance at least slightly bigger than standard
    1725   relaxedTolerance = relaxedTolerance +  error;
    1726   sumOfRelaxedDualInfeasibilities_ = 0.0;
    1727 
    1728   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1729     if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) {
    1730       // not basic
    1731       double distanceUp = columnUpperWork_[iColumn]-
    1732         columnActivityWork_[iColumn];
    1733       double distanceDown = columnActivityWork_[iColumn] -
    1734         columnLowerWork_[iColumn];
    1735       if (distanceUp>primalTolerance_) {
    1736         double value = reducedCostWork_[iColumn];
    1737         // Check if "free"
    1738         if (distanceDown>primalTolerance_) {
    1739           if (fabs(value)>1.0e2*relaxedTolerance) {
    1740             numberSuperBasicWithDj++;
    1741             if (firstFreeDual<0)
    1742               firstFreeDual = iColumn;
    1743           }
    1744           if (firstFreePrimal<0)
    1745             firstFreePrimal = iColumn;
    1746         }
    1747         // should not be negative
    1748         if (value<0.0) {
    1749           value = - value;
    1750           if (value>columnDualInfeasibility_) {
    1751             columnDualInfeasibility_=value;
    1752             columnDualSequence_=iColumn;
    1753           }
    1754           if (value>dualTolerance_) {
    1755             if (getColumnStatus(iColumn) != isFree) {
    1756               numberDualInfeasibilitiesWithoutFree_ ++;
    1757               sumDualInfeasibilities_ += value-dualTolerance_;
    1758               if (value>relaxedTolerance)
    1759                 sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    1760               numberDualInfeasibilities_ ++;
    1761             } else {
    1762               // free so relax a lot
    1763               value *= 0.01;
    1764               if (value>dualTolerance_) {
    1765                 sumDualInfeasibilities_ += value-dualTolerance_;
    1766                 if (value>relaxedTolerance)
    1767                   sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    1768                 numberDualInfeasibilities_ ++;
    1769               }
    1770             }
    1771             // maybe we can make feasible by increasing tolerance
    1772             if (distanceUp<largeValue_) {
    1773               if (distanceUp>primalToleranceToGetOptimal_)
    1774                 primalToleranceToGetOptimal_=distanceUp;
    1775             } else {
    1776               //gap too big for any tolerance
    1777               remainingDualInfeasibility_=
    1778                 max(remainingDualInfeasibility_,value);
    1779             }
    1780           }
    1781         }
    1782       }
    1783       if (distanceDown>primalTolerance_) {
    1784         double value = reducedCostWork_[iColumn];
    1785         // should not be positive
    1786         if (value>0.0) {
    1787           if (value>columnDualInfeasibility_) {
    1788             columnDualInfeasibility_=value;
    1789             columnDualSequence_=iColumn;
    1790           }
    1791           if (value>dualTolerance_) {
    1792             sumDualInfeasibilities_ += value-dualTolerance_;
    1793             if (value>relaxedTolerance)
    1794               sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    1795             numberDualInfeasibilities_ ++;
    1796             if (getColumnStatus(iColumn) != isFree)
    1797               numberDualInfeasibilitiesWithoutFree_ ++;
    1798             // maybe we can make feasible by increasing tolerance
    1799             if (distanceDown<largeValue_&&
    1800                 distanceDown>primalToleranceToGetOptimal_)
    1801               primalToleranceToGetOptimal_=distanceDown;
    1802           }
    1803         }
    1804       }
    1805     }
    1806   }
    1807   for (iRow=0;iRow<numberRows_;iRow++) {
    1808     if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) {
    1809       // not basic
    1810       double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
    1811       double distanceDown = rowActivityWork_[iRow] -rowLowerWork_[iRow];
    1812       if (distanceUp>primalTolerance_) {
    1813         double value = rowReducedCost_[iRow];
    1814         // Check if "free"
    1815         if (distanceDown>primalTolerance_) {
    1816           if (fabs(value)>1.0e2*relaxedTolerance) {
    1817             numberSuperBasicWithDj++;
    1818             if (firstFreeDual<0)
    1819               firstFreeDual = iRow+numberColumns_;
    1820           }
    1821           if (firstFreePrimal<0)
    1822             firstFreePrimal = iRow+numberColumns_;
    1823         }
    1824         // should not be negative
    1825         if (value<0.0) {
    1826           value = - value;
    1827           if (value>rowDualInfeasibility_) {
    1828             rowDualInfeasibility_=value;
    1829             rowDualSequence_=iRow;
    1830           }
    1831           if (value>dualTolerance_) {
    1832             sumDualInfeasibilities_ += value-dualTolerance_;
    1833             if (value>relaxedTolerance)
    1834               sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    1835             numberDualInfeasibilities_ ++;
    1836             if (getRowStatus(iRow) != isFree)
    1837               numberDualInfeasibilitiesWithoutFree_ ++;
    1838             // maybe we can make feasible by increasing tolerance
    1839             if (distanceUp<largeValue_) {
    1840               if (distanceUp>primalToleranceToGetOptimal_)
    1841                 primalToleranceToGetOptimal_=distanceUp;
    1842             } else {
    1843               //gap too big for any tolerance
    1844               remainingDualInfeasibility_=
    1845                 max(remainingDualInfeasibility_,value);
    1846             }
    1847           }
    1848         }
    1849       }
    1850       if (distanceDown>primalTolerance_) {
    1851         double value = rowReducedCost_[iRow];
    1852         // should not be positive
    1853         if (value>0.0) {
    1854           if (value>rowDualInfeasibility_) {
    1855             rowDualInfeasibility_=value;
    1856             rowDualSequence_=iRow;
    1857           }
    1858           if (value>dualTolerance_) {
    1859             sumDualInfeasibilities_ += value-dualTolerance_;
    1860             if (value>relaxedTolerance)
    1861               sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
    1862             numberDualInfeasibilities_ ++;
    1863             if (getRowStatus(iRow) != isFree)
    1864               numberDualInfeasibilitiesWithoutFree_ ++;
    1865             // maybe we can make feasible by increasing tolerance
    1866             if (distanceDown<largeValue_&&
    1867                 distanceDown>primalToleranceToGetOptimal_)
    1868               primalToleranceToGetOptimal_=distanceDown;
    1869           }
    1870         }
    1871       }
    1872     }
    1873   }
    1874   if (algorithm_<0&&firstFreeDual>=0) {
    1875     // dual
    1876     firstFree_ = firstFreeDual;
    1877   } else if (numberSuperBasicWithDj||
    1878              (progress_&&progress_->lastIterationNumber(0)<=0)) {
    1879     firstFree_=firstFreePrimal;
    1880   }
    1881 }
    1882 /*
    1883   Unpacks one column of the matrix into indexed array
    1884 */
    1885 void
    1886 ClpInterior::unpack(CoinIndexedVector * rowArray) const
    1887 {
    1888   rowArray->clear();
    1889   if (sequenceIn_>=numberColumns_) {
    1890     //slack
    1891     rowArray->insert(sequenceIn_-numberColumns_,-1.0);
    1892   } else {
    1893     // column
    1894     matrix_->unpack(this,rowArray,sequenceIn_);
    1895   }
    1896 }
    1897 void
    1898 ClpInterior::unpack(CoinIndexedVector * rowArray,int sequence) const
    1899 {
    1900   rowArray->clear();
    1901   if (sequence>=numberColumns_) {
    1902     //slack
    1903     rowArray->insert(sequence-numberColumns_,-1.0);
    1904   } else {
    1905     // column
    1906     matrix_->unpack(this,rowArray,sequence);
    1907   }
    1908 }
    1909 /*
    1910   Unpacks one column of the matrix into indexed array
    1911 */
    1912 void
    1913 ClpInterior::unpackPacked(CoinIndexedVector * rowArray)
    1914 {
    1915   rowArray->clear();
    1916   if (sequenceIn_>=numberColumns_) {
    1917     //slack
    1918     int * index = rowArray->getIndices();
    1919     double * array = rowArray->denseVector();
    1920     array[0]=-1.0;
    1921     index[0]=sequenceIn_-numberColumns_;
    1922     rowArray->setNumElements(1);
    1923     rowArray->setPackedMode(true);
    1924   } else {
    1925     // column
    1926     matrix_->unpackPacked(this,rowArray,sequenceIn_);
    1927   }
    1928 }
    1929 void
    1930 ClpInterior::unpackPacked(CoinIndexedVector * rowArray,int sequence)
    1931 {
    1932   rowArray->clear();
    1933   if (sequence>=numberColumns_) {
    1934     //slack
    1935     int * index = rowArray->getIndices();
    1936     double * array = rowArray->denseVector();
    1937     array[0]=-1.0;
    1938     index[0]=sequence-numberColumns_;
    1939     rowArray->setNumElements(1);
    1940     rowArray->setPackedMode(true);
    1941   } else {
    1942     // column
    1943     matrix_->unpackPacked(this,rowArray,sequence);
    1944   }
    1945159}
    1946160bool
    1947 ClpInterior::createRim(int what,bool makeRowCopy)
     161ClpInterior::createWorkingData()
    1948162{
    1949163  bool goodMatrix=true;
    1950   int saveLevel=handler_->logLevel();
    1951   if (problemStatus_==10) {
    1952     handler_->setLogLevel(0); // switch off messages
    1953   } else if (factorization_) {
    1954     // match up factorization messages
    1955     if (handler_->logLevel()<3)
    1956       factorization_->messageLevel(0);
    1957     else
    1958       factorization_->messageLevel(3);
    1959   }
    1960   int i;
    1961   if ((what&(16+32))!=0) {
    1962     // move information to work arrays
    1963     double direction = optimizationDirection_;
    1964     // direction is actually scale out not scale in
    1965     if (direction)
    1966       direction = 1.0/direction;
    1967     if (direction!=1.0) {
    1968       // reverse all dual signs
    1969       for (i=0;i<numberColumns_;i++)
    1970         reducedCost_[i] *= direction;
    1971       for (i=0;i<numberRows_;i++)
    1972         dual_[i] *= direction;
    1973     }
    1974     // row reduced costs
    1975     if (!dj_) {
    1976       dj_ = new double[numberRows_+numberColumns_];
    1977       reducedCostWork_ = dj_;
    1978       rowReducedCost_ = dj_+numberColumns_;
    1979       memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
    1980       memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double));
    1981     }
    1982     if (!solution_||(what&32)!=0) {
    1983       if (!solution_)
    1984         solution_ = new double[numberRows_+numberColumns_];
    1985       columnActivityWork_ = solution_;
    1986       rowActivityWork_ = solution_+numberColumns_;
    1987       memcpy(columnActivityWork_,columnActivity_,
    1988              numberColumns_*sizeof(double));
    1989       memcpy(rowActivityWork_,rowActivity_,
    1990              numberRows_*sizeof(double));
    1991     }
    1992   }
    1993   if ((what&16)!=0) {
    1994     //check matrix
    1995     if (!matrix_->allElementsInRange(this,smallElement_,1.0e20)) {
    1996       problemStatus_=4;
    1997       goodMatrix= false;
    1998     }
    1999     if (makeRowCopy) {
    2000       delete rowCopy_;
    2001       // may return NULL if can't give row copy
    2002       rowCopy_ = matrix_->reverseOrderedCopy();
    2003 #ifdef TAKEOUT
    2004       {
    2005 
    2006         ClpPackedMatrix* rowCopy =
    2007           dynamic_cast< ClpPackedMatrix*>(rowCopy_);
    2008         const int * column = rowCopy->getIndices();
    2009         const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
    2010         const double * element = rowCopy->getElements();
    2011         int i;
    2012         for (i=133;i<numberRows_;i++) {
    2013           if (rowStart[i+1]-rowStart[i]==10||rowStart[i+1]-rowStart[i]==15)
    2014             printf("Row %d has %d elements\n",i,rowStart[i+1]-rowStart[i]);
    2015         }
    2016       } 
    2017 #endif
    2018     }
    2019   }
    2020   if ((what&4)!=0) {
    2021     delete [] cost_;
    2022     // extra copy with original costs
    2023     int nTotal = numberRows_+numberColumns_;
    2024     //cost_ = new double[2*nTotal];
    2025     cost_ = new double[nTotal];
    2026     objectiveWork_ = cost_;
    2027     rowObjectiveWork_ = cost_+numberColumns_;
    2028     memcpy(objectiveWork_,objective(),numberColumns_*sizeof(double));
    2029     if (rowObjective_)
    2030       memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
    2031     else
    2032       memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
    2033     // and initialize changes to zero
    2034     //memset(cost_+nTotal,0,nTotal*sizeof(double));
    2035   }
    2036   if ((what&1)!=0) {
    2037     delete [] lower_;
    2038     delete [] upper_;
    2039     lower_ = new double[numberColumns_+numberRows_];
    2040     upper_ = new double[numberColumns_+numberRows_];
    2041     rowLowerWork_ = lower_+numberColumns_;
    2042     columnLowerWork_ = lower_;
    2043     rowUpperWork_ = upper_+numberColumns_;
    2044     columnUpperWork_ = upper_;
    2045     memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
    2046     memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
    2047     memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
    2048     memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
    2049     // clean up any mismatches on infinity
    2050     for (i=0;i<numberColumns_;i++) {
    2051       if (columnLowerWork_[i]<-1.0e30)
    2052         columnLowerWork_[i] = -COIN_DBL_MAX;
    2053       if (columnUpperWork_[i]>1.0e30)
    2054         columnUpperWork_[i] = COIN_DBL_MAX;
    2055     }
    2056     // clean up any mismatches on infinity
    2057     for (i=0;i<numberRows_;i++) {
    2058       if (rowLowerWork_[i]<-1.0e30)
    2059         rowLowerWork_[i] = -COIN_DBL_MAX;
    2060       if (rowUpperWork_[i]>1.0e30)
    2061         rowUpperWork_[i] = COIN_DBL_MAX;
    2062     }
    2063   }
    2064   // do scaling if needed
    2065   if (scalingFlag_>0&&!rowScale_) {
    2066     if (matrix_->scale(this))
    2067       scalingFlag_=-scalingFlag_; // not scaled after all
    2068   }
    2069   if ((what&4)!=0) {
    2070     double direction = optimizationDirection_;
    2071     // direction is actually scale out not scale in
    2072     if (direction)
    2073       direction = 1.0/direction;
    2074     // but also scale by scale factors
    2075     // not really sure about row scaling
    2076     if (!rowScale_) {
    2077       if (direction!=1.0) {
    2078         for (i=0;i<numberRows_;i++)
    2079           rowObjectiveWork_[i] *= direction;
    2080         for (i=0;i<numberColumns_;i++)
    2081           objectiveWork_[i] *= direction;
    2082       }
    2083     } else {
    2084       for (i=0;i<numberRows_;i++)
    2085         rowObjectiveWork_[i] *= direction/rowScale_[i];
    2086       for (i=0;i<numberColumns_;i++)
    2087         objectiveWork_[i] *= direction*columnScale_[i];
    2088     }
    2089   }
    2090   if ((what&(1+32))!=0&&rowScale_) {
    2091     for (i=0;i<numberColumns_;i++) {
    2092       double multiplier = 1.0/columnScale_[i];
    2093       if (columnLowerWork_[i]>-1.0e50)
    2094         columnLowerWork_[i] *= multiplier;
    2095       if (columnUpperWork_[i]<1.0e50)
    2096         columnUpperWork_[i] *= multiplier;
    2097      
    2098     }
    2099     for (i=0;i<numberRows_;i++) {
    2100       if (rowLowerWork_[i]>-1.0e50)
    2101         rowLowerWork_[i] *= rowScale_[i];
    2102       if (rowUpperWork_[i]<1.0e50)
    2103         rowUpperWork_[i] *= rowScale_[i];
    2104     }
    2105   }
    2106   if ((what&(8+32))!=0&&rowScale_) {
    2107     // on entry
    2108     for (i=0;i<numberColumns_;i++) {
    2109       columnActivityWork_[i] /= columnScale_[i];
    2110       reducedCostWork_[i] *= columnScale_[i];
    2111     }
    2112     for (i=0;i<numberRows_;i++) {
    2113       rowActivityWork_[i] *= rowScale_[i];
    2114       dual_[i] /= rowScale_[i];
    2115       rowReducedCost_[i] = dual_[i];
    2116     }
    2117   }
    2118  
    2119   if ((what&16)!=0) {
    2120     // check rim of problem okay
    2121     if (!sanityCheck())
    2122       goodMatrix=false;
    2123   }
    2124   // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
    2125   // maybe we need to move scales to SimplexModel for factorization?
    2126   if ((what&8)!=0&&!pivotVariable_) {
    2127     pivotVariable_=new int[numberRows_];
    2128   }
    2129 
    2130   if ((what&16)!=0&&!rowArray_[2]) {
    2131     // get some arrays
    2132     int iRow,iColumn;
    2133     // these are "indexed" arrays so we always know where nonzeros are
    2134     /**********************************************************
    2135       rowArray_[3] is long enough for rows+columns
    2136     *********************************************************/
    2137     for (iRow=0;iRow<4;iRow++) {
    2138       if (!rowArray_[iRow]) {
    2139         rowArray_[iRow]=new CoinIndexedVector();
    2140         int length =numberRows_+factorization_->maximumPivots();
    2141         if (iRow==3)
    2142           length += numberColumns_;
    2143         rowArray_[iRow]->reserve(length);
    2144       }
    2145     }
    2146    
    2147     for (iColumn=0;iColumn<2;iColumn++) {
    2148       if (!columnArray_[iColumn]) {
    2149         columnArray_[iColumn]=new CoinIndexedVector();
    2150         columnArray_[iColumn]->reserve(numberColumns_);
    2151       }
    2152     }   
    2153   }
    2154   double primalTolerance=dblParam_[ClpPrimalTolerance];
    2155   if ((what&1)!=0) {
    2156     // fix any variables with tiny gaps
    2157     for (i=0;i<numberColumns_;i++) {
    2158       if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
    2159         if (columnLowerWork_[i]>=0.0) {
    2160           columnUpperWork_[i] = columnLowerWork_[i];
    2161         } else if (columnUpperWork_[i]<=0.0) {
    2162           columnLowerWork_[i] = columnUpperWork_[i];
    2163         } else {
    2164           columnUpperWork_[i] = 0.0;
    2165           columnLowerWork_[i] = 0.0;
    2166         }
    2167       }
    2168     }
    2169     for (i=0;i<numberRows_;i++) {
    2170       if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
    2171         if (rowLowerWork_[i]>=0.0) {
    2172           rowUpperWork_[i] = rowLowerWork_[i];
    2173         } else if (rowUpperWork_[i]<=0.0) {
    2174           rowLowerWork_[i] = rowUpperWork_[i];
    2175         } else {
    2176           rowUpperWork_[i] = 0.0;
    2177           rowLowerWork_[i] = 0.0;
    2178         }
    2179       }
    2180     }
    2181   }
    2182   if (problemStatus_==10) {
    2183     problemStatus_=-1;
    2184     handler_->setLogLevel(saveLevel); // switch back messages
    2185   }
     164  //check matrix
     165  if (!matrix_->allElementsInRange(this,1.0e-12,1.0e20)) {
     166    problemStatus_=4;
     167    goodMatrix= false;
     168  }
     169  // check rim of problem okay
     170  if (!sanityCheck())
     171    goodMatrix=false;
    2186172  return goodMatrix;
    2187173}
    2188174void
    2189 ClpInterior::deleteRim(int getRidOfFactorizationData)
    2190 {
    2191   int i;
    2192   if (problemStatus_!=1&&problemStatus_!=2) {
    2193     delete [] ray_;
    2194     ray_=NULL;
    2195   }
    2196   // ray may be null if in branch and bound
    2197   if (rowScale_) {
    2198     for (i=0;i<numberColumns_;i++) {
    2199       columnActivity_[i] = columnActivityWork_[i]*columnScale_[i];
    2200       reducedCost_[i] = reducedCostWork_[i]/columnScale_[i];
    2201     }
    2202     for (i=0;i<numberRows_;i++) {
    2203       rowActivity_[i] = rowActivityWork_[i]/rowScale_[i];
    2204       dual_[i] *= rowScale_[i];
    2205     }
    2206     if (problemStatus_==2) {
    2207       for (i=0;i<numberColumns_;i++) {
    2208         ray_[i] *= columnScale_[i];
    2209       }
    2210     } else if (problemStatus_==1&&ray_) {
    2211       for (i=0;i<numberRows_;i++) {
    2212         ray_[i] *= rowScale_[i];
    2213       }
    2214     }
    2215   } else {
    2216     for (i=0;i<numberColumns_;i++) {
    2217       columnActivity_[i] = columnActivityWork_[i];
    2218       reducedCost_[i] = reducedCostWork_[i];
    2219     }
    2220     for (i=0;i<numberRows_;i++) {
    2221       rowActivity_[i] = rowActivityWork_[i];
    2222     }
    2223   }
    2224   if (optimizationDirection_!=1.0) {
    2225     // and modify all dual signs
    2226     for (i=0;i<numberColumns_;i++)
    2227       reducedCost_[i] *= optimizationDirection_;
    2228     for (i=0;i<numberRows_;i++)
    2229       dual_[i] *= optimizationDirection_;
    2230   }
    2231   // scaling may have been turned off
    2232   scalingFlag_ = abs(scalingFlag_);
    2233   if(getRidOfFactorizationData>=0)
    2234     gutsOfDelete(getRidOfFactorizationData+1);
    2235 }
    2236 void
    2237 ClpInterior::setDualBound(double value)
    2238 {
    2239   if (value>0.0)
    2240     dualBound_=value;
    2241 }
    2242 void
    2243 ClpInterior::setInfeasibilityCost(double value)
    2244 {
    2245   if (value>0.0)
    2246     infeasibilityCost_=value;
    2247 }
    2248 void ClpInterior::setNumberRefinements( int value)
    2249 {
    2250   if (value>=0&&value<10)
    2251     numberRefinements_=value;
    2252 }
    2253 // Sets row pivot choice algorithm in dual
    2254 void
    2255 ClpInterior::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
    2256 {
    2257   delete dualRowPivot_;
    2258   dualRowPivot_ = choice.clone(true);
    2259 }
    2260 // Sets row pivot choice algorithm in dual
    2261 void
    2262 ClpInterior::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
    2263 {
    2264   delete primalColumnPivot_;
    2265   primalColumnPivot_ = choice.clone(true);
    2266 }
    2267 // Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
    2268 void
    2269 ClpInterior::scaling(int mode)
    2270 {
    2271   if (mode>0&&mode<4) {
    2272     scalingFlag_=mode;
    2273   } else if (!mode) {
    2274     scalingFlag_=0;
    2275     delete [] rowScale_;
    2276     rowScale_ = NULL;
    2277     delete [] columnScale_;
    2278     columnScale_ = NULL;
    2279   }
    2280 }
    2281 // Passes in factorization
    2282 void
    2283 ClpInterior::setFactorization( ClpFactorization & factorization)
    2284 {
    2285   delete factorization_;
    2286   factorization_= new ClpFactorization(factorization);
    2287 }
    2288 void
    2289 ClpInterior::times(double scalar,
    2290                   const double * x, double * y) const
    2291 {
    2292   if (rowScale_)
    2293     matrix_->times(scalar,x,y,rowScale_,columnScale_);
    2294   else
    2295     matrix_->times(scalar,x,y);
    2296 }
    2297 void
    2298 ClpInterior::transposeTimes(double scalar,
    2299                            const double * x, double * y) const
    2300 {
    2301   if (rowScale_)
    2302     matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
    2303   else
    2304     matrix_->transposeTimes(scalar,x,y);
    2305 }
    2306 /* Perturbation:
    2307    -50 to +50 - perturb by this power of ten (-6 sounds good)
    2308    100 - auto perturb if takes too long (1.0e-6 largest nonzero)
    2309    101 - we are perturbed
    2310    102 - don't try perturbing again
    2311    default is 100
    2312 */
    2313 void
    2314 ClpInterior::setPerturbation(int value)
    2315 {
    2316   if(value<=100&&value >=-1000) {
    2317     perturbation_=value;
    2318   }
    2319 }
    2320 // Sparsity on or off
    2321 bool
    2322 ClpInterior::sparseFactorization() const
    2323 {
    2324   return factorization_->sparseThreshold()!=0;
    2325 }
    2326 void
    2327 ClpInterior::setSparseFactorization(bool value)
    2328 {
    2329   if (value) {
    2330     if (!factorization_->sparseThreshold())
    2331       factorization_->goSparse();
    2332   } else {
    2333     factorization_->sparseThreshold(0);
    2334   }
    2335 }
    2336 /* Tightens primal bounds to make dual faster.  Unless
    2337    fixed, bounds are slightly looser than they could be.
    2338    This is to make dual go faster and is probably not needed
    2339    with a presolve.  Returns non-zero if problem infeasible
    2340 
    2341    Fudge for branch and bound - put bounds on columns of factor *
    2342    largest value (at continuous) - should improve stability
    2343    in branch and bound on infeasible branches (0.0 is off)
    2344 */
    2345 int
    2346 ClpInterior::tightenPrimalBounds(double factor)
    2347 {
    2348  
    2349   // Get a row copy in standard format
    2350   CoinPackedMatrix copy;
    2351   copy.reverseOrderedCopyOf(*matrix());
    2352   // get matrix data pointers
    2353   const int * column = copy.getIndices();
    2354   const CoinBigIndex * rowStart = copy.getVectorStarts();
    2355   const int * rowLength = copy.getVectorLengths();
    2356   const double * element = copy.getElements();
    2357   int numberChanged=1,iPass=0;
    2358   double large = largeValue(); // treat bounds > this as infinite
    2359   int numberInfeasible=0;
    2360   int totalTightened = 0;
    2361 
    2362   double tolerance = primalTolerance();
    2363 
    2364 
    2365   // Save column bounds
    2366   double * saveLower = new double [numberColumns_];
    2367   memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
    2368   double * saveUpper = new double [numberColumns_];
    2369   memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
    2370 
    2371   int iRow, iColumn;
    2372 
    2373   // If wanted - tighten column bounds using solution
    2374   if (factor) {
    2375     assert (factor>1.0);
    2376     double largest=0.0;
    2377     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2378       if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
    2379         largest = max(largest,fabs(columnActivity_[iColumn]));
    2380       }
    2381     }
    2382     largest *= factor;
    2383     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2384       if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
    2385         columnUpper_[iColumn] = min(columnUpper_[iColumn],largest);
    2386         columnLower_[iColumn] = max(columnLower_[iColumn],-largest);
    2387       }
    2388     }
    2389   }
    2390 #define MAXPASS 10
    2391 
    2392   // Loop round seeing if we can tighten bounds
    2393   // Would be faster to have a stack of possible rows
    2394   // and we put altered rows back on stack
    2395   int numberCheck=-1;
    2396   while(numberChanged>numberCheck) {
    2397 
    2398     numberChanged = 0; // Bounds tightened this pass
    2399    
    2400     if (iPass==MAXPASS) break;
    2401     iPass++;
    2402    
    2403     for (iRow = 0; iRow < numberRows_; iRow++) {
    2404 
    2405       if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
    2406 
    2407         // possible row
    2408         int infiniteUpper = 0;
    2409         int infiniteLower = 0;
    2410         double maximumUp = 0.0;
    2411         double maximumDown = 0.0;
    2412         double newBound;
    2413         CoinBigIndex rStart = rowStart[iRow];
    2414         CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
    2415         CoinBigIndex j;
    2416         // Compute possible lower and upper ranges
    2417      
    2418         for (j = rStart; j < rEnd; ++j) {
    2419           double value=element[j];
    2420           iColumn = column[j];
    2421           if (value > 0.0) {
    2422             if (columnUpper_[iColumn] >= large) {
    2423               ++infiniteUpper;
    2424             } else {
    2425               maximumUp += columnUpper_[iColumn] * value;
    2426             }
    2427             if (columnLower_[iColumn] <= -large) {
    2428               ++infiniteLower;
    2429             } else {
    2430               maximumDown += columnLower_[iColumn] * value;
    2431             }
    2432           } else if (value<0.0) {
    2433             if (columnUpper_[iColumn] >= large) {
    2434               ++infiniteLower;
    2435             } else {
    2436               maximumDown += columnUpper_[iColumn] * value;
    2437             }
    2438             if (columnLower_[iColumn] <= -large) {
    2439               ++infiniteUpper;
    2440             } else {
    2441               maximumUp += columnLower_[iColumn] * value;
    2442             }
    2443           }
    2444         }
    2445         double maxUp = maximumUp+infiniteUpper*1.0e31;
    2446         double maxDown = maximumDown-infiniteLower*1.0e31;
    2447         if (maxUp <= rowUpper_[iRow] + tolerance &&
    2448             maxDown >= rowLower_[iRow] - tolerance) {
    2449          
    2450           // Row is redundant - make totally free
    2451           // NO - can't do this for postsolve
    2452           // rowLower_[iRow]=-COIN_DBL_MAX;
    2453           // rowUpper_[iRow]=COIN_DBL_MAX;
    2454           //printf("Redundant row in presolveX %d\n",iRow);
    2455 
    2456         } else {
    2457           if (maxUp < rowLower_[iRow] -tolerance ||
    2458               maxDown > rowUpper_[iRow]+tolerance) {
    2459             // problem is infeasible - exit at once
    2460             numberInfeasible++;
    2461             break;
    2462           }
    2463           double lower = rowLower_[iRow];
    2464           double upper = rowUpper_[iRow];
    2465           for (j = rStart; j < rEnd; ++j) {
    2466             double value=element[j];
    2467             iColumn = column[j];
    2468             double nowLower = columnLower_[iColumn];
    2469             double nowUpper = columnUpper_[iColumn];
    2470             if (value > 0.0) {
    2471               // positive value
    2472               if (lower>-large) {
    2473                 if (!infiniteUpper) {
    2474                   assert(nowUpper < large);
    2475                   newBound = nowUpper +
    2476                     (lower - maximumUp) / value;
    2477                 } else if (infiniteUpper==1&&nowUpper>large) {
    2478                   newBound = (lower -maximumUp) / value;
    2479                 } else {
    2480                   newBound = -COIN_DBL_MAX;
    2481                 }
    2482                 if (newBound > nowLower + 1.0e-12) {
    2483                   // Tighten the lower bound
    2484                   columnLower_[iColumn] = newBound;
    2485                   numberChanged++;
    2486                   // check infeasible (relaxed)
    2487                   if (nowUpper - newBound <
    2488                       -100.0*tolerance) {
    2489                     numberInfeasible++;
    2490                   }
    2491                   // adjust
    2492                   double now;
    2493                   if (nowLower<-large) {
    2494                     now=0.0;
    2495                     infiniteLower--;
    2496                   } else {
    2497                     now = nowLower;
    2498                   }
    2499                   maximumDown += (newBound-now) * value;
    2500                   nowLower = newBound;
    2501                 }
    2502               }
    2503               if (upper <large) {
    2504                 if (!infiniteLower) {
    2505                   assert(nowLower >- large);
    2506                   newBound = nowLower +
    2507                     (upper - maximumDown) / value;
    2508                 } else if (infiniteLower==1&&nowLower<-large) {
    2509                   newBound =   (upper - maximumDown) / value;
    2510                 } else {
    2511                   newBound = COIN_DBL_MAX;
    2512                 }
    2513                 if (newBound < nowUpper - 1.0e-12) {
    2514                   // Tighten the upper bound
    2515                   columnUpper_[iColumn] = newBound;
    2516                   numberChanged++;
    2517                   // check infeasible (relaxed)
    2518                   if (newBound - nowLower <
    2519                       -100.0*tolerance) {
    2520                     numberInfeasible++;
    2521                   }
    2522                   // adjust
    2523                   double now;
    2524                   if (nowUpper>large) {
    2525                     now=0.0;
    2526                     infiniteUpper--;
    2527                   } else {
    2528                     now = nowUpper;
    2529                   }
    2530                   maximumUp += (newBound-now) * value;
    2531                   nowUpper = newBound;
    2532                 }
    2533               }
    2534             } else {
    2535               // negative value
    2536               if (lower>-large) {
    2537                 if (!infiniteUpper) {
    2538                   assert(nowLower >- large);
    2539                   newBound = nowLower +
    2540                     (lower - maximumUp) / value;
    2541                 } else if (infiniteUpper==1&&nowLower<-large) {
    2542                   newBound = (lower -maximumUp) / value;
    2543                 } else {
    2544                   newBound = COIN_DBL_MAX;
    2545                 }
    2546                 if (newBound < nowUpper - 1.0e-12) {
    2547                   // Tighten the upper bound
    2548                   columnUpper_[iColumn] = newBound;
    2549                   numberChanged++;
    2550                   // check infeasible (relaxed)
    2551                   if (newBound - nowLower <
    2552                       -100.0*tolerance) {
    2553                     numberInfeasible++;
    2554                   }
    2555                   // adjust
    2556                   double now;
    2557                   if (nowUpper>large) {
    2558                     now=0.0;
    2559                     infiniteLower--;
    2560                   } else {
    2561                     now = nowUpper;
    2562                   }
    2563                   maximumDown += (newBound-now) * value;
    2564                   nowUpper = newBound;
    2565                 }
    2566               }
    2567               if (upper <large) {
    2568                 if (!infiniteLower) {
    2569                   assert(nowUpper < large);
    2570                   newBound = nowUpper +
    2571                     (upper - maximumDown) / value;
    2572                 } else if (infiniteLower==1&&nowUpper>large) {
    2573                   newBound =   (upper - maximumDown) / value;
    2574                 } else {
    2575                   newBound = -COIN_DBL_MAX;
    2576                 }
    2577                 if (newBound > nowLower + 1.0e-12) {
    2578                   // Tighten the lower bound
    2579                   columnLower_[iColumn] = newBound;
    2580                   numberChanged++;
    2581                   // check infeasible (relaxed)
    2582                   if (nowUpper - newBound <
    2583                       -100.0*tolerance) {
    2584                     numberInfeasible++;
    2585                   }
    2586                   // adjust
    2587                   double now;
    2588                   if (nowLower<-large) {
    2589                     now=0.0;
    2590                     infiniteUpper--;
    2591                   } else {
    2592                     now = nowLower;
    2593                   }
    2594                   maximumUp += (newBound-now) * value;
    2595                   nowLower = newBound;
    2596                 }
    2597               }
    2598             }
    2599           }
    2600         }
    2601       }
    2602     }
    2603     totalTightened += numberChanged;
    2604     if (iPass==1)
    2605       numberCheck=numberChanged>>4;
    2606     if (numberInfeasible) break;
    2607   }
    2608   if (!numberInfeasible) {
    2609     handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
    2610       <<totalTightened
    2611       <<CoinMessageEol;
    2612     // Set bounds slightly loose
    2613     double useTolerance = 1.0e-3;
    2614     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    2615       if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
    2616         if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e8) {
    2617           // relax enough so will have correct dj
    2618 #if 1
    2619           columnLower_[iColumn]=max(saveLower[iColumn],
    2620                                     columnLower_[iColumn]-100.0*useTolerance);
    2621           columnUpper_[iColumn]=min(saveUpper[iColumn],
    2622                                     columnUpper_[iColumn]+100.0*useTolerance);
    2623 #else
    2624           if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
    2625             if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
    2626               columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
    2627             } else {
    2628               columnLower_[iColumn]=saveLower[iColumn];
    2629               columnUpper_[iColumn]=min(saveUpper[iColumn],
    2630                                         saveLower[iColumn]+100.0*useTolerance);
    2631             }
    2632           } else {
    2633             if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
    2634               columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
    2635             } else {
    2636               columnUpper_[iColumn]=saveUpper[iColumn];
    2637               columnLower_[iColumn]=max(saveLower[iColumn],
    2638                                         saveUpper[iColumn]-100.0*useTolerance);
    2639             }
    2640           }
    2641 #endif
    2642         } else {
    2643           if (columnUpper_[iColumn]<saveUpper[iColumn]) {
    2644             // relax a bit
    2645             columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*useTolerance,
    2646                                         saveUpper[iColumn]);
    2647           }
    2648           if (columnLower_[iColumn]>saveLower[iColumn]) {
    2649             // relax a bit
    2650             columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*useTolerance,
    2651                                         saveLower[iColumn]);
    2652           }
    2653         }
    2654       }
    2655     }
    2656   } else {
    2657     handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
    2658       <<numberInfeasible
    2659       <<CoinMessageEol;
    2660     // restore column bounds
    2661     memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
    2662     memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
    2663   }
    2664   delete [] saveLower;
    2665   delete [] saveUpper;
    2666   return (numberInfeasible);
    2667 }
    2668 // dual
    2669 #include "ClpInteriorDual.hpp"
    2670 #include "ClpInteriorPrimal.hpp"
    2671 int ClpInterior::dual (int ifValuesPass )
    2672 {
    2673   int returnCode = ((ClpInteriorDual *) this)->dual(ifValuesPass);
    2674   if (problemStatus_==10) {
    2675     //printf("Cleaning up with primal\n");
    2676     int savePerturbation = perturbation_;
    2677     perturbation_=100;
    2678     bool denseFactorization = initialDenseFactorization();
    2679     // It will be safe to allow dense
    2680     setInitialDenseFactorization(true);
    2681     returnCode = ((ClpInteriorPrimal *) this)->primal(1);
    2682     setInitialDenseFactorization(denseFactorization);
    2683     perturbation_=savePerturbation;
    2684   }
    2685   return returnCode;
    2686 }
    2687 // primal
    2688 int ClpInterior::primal (int ifValuesPass )
    2689 {
    2690   int returnCode = ((ClpInteriorPrimal *) this)->primal(ifValuesPass);
    2691   if (problemStatus_==10) {
    2692     //printf("Cleaning up with dual\n");
    2693     int savePerturbation = perturbation_;
    2694     perturbation_=100;
    2695     bool denseFactorization = initialDenseFactorization();
    2696     // It will be safe to allow dense
    2697     setInitialDenseFactorization(true);
    2698     returnCode = ((ClpInteriorDual *) this)->dual(0);
    2699     setInitialDenseFactorization(denseFactorization);
    2700     perturbation_=savePerturbation;
    2701   }
    2702   return returnCode;
    2703 }
    2704 #include "ClpInteriorPrimalQuadratic.hpp"
    2705 /* Solves quadratic problem using SLP - may be used as crash
    2706    for other algorithms when number of iterations small
    2707 */
    2708 int
    2709 ClpInterior::quadraticSLP(int numberPasses, double deltaTolerance)
    2710 {
    2711   return ((ClpInteriorPrimalQuadratic *) this)->primalSLP(numberPasses,deltaTolerance);
    2712 }
    2713 // Solves quadratic using Dantzig's algorithm - primal
    2714 int
    2715 ClpInterior::quadraticPrimal(int phase)
    2716 {
    2717   return ((ClpInteriorPrimalQuadratic *) this)->primalQuadratic(phase);
    2718 }
    2719 /* For strong branching.  On input lower and upper are new bounds
    2720    while on output they are objective function values (>1.0e50 infeasible).
    2721    Return code is 0 if nothing interesting, -1 if infeasible both
    2722    ways and +1 if infeasible one way (check values to see which one(s))
    2723 */
    2724 int ClpInterior::strongBranching(int numberVariables,const int * variables,
    2725                                 double * newLower, double * newUpper,
    2726                                 double ** outputSolution,
    2727                                 int * outputStatus, int * outputIterations,
    2728                                 bool stopOnFirstInfeasible,
    2729                                 bool alwaysFinish)
    2730 {
    2731   return ((ClpInteriorDual *) this)->strongBranching(numberVariables,variables,
    2732                                                     newLower,  newUpper,outputSolution,
    2733                                                     outputStatus, outputIterations,
    2734                                                     stopOnFirstInfeasible,
    2735                                                     alwaysFinish);
    2736 }
    2737 /* Borrow model.  This is so we dont have to copy large amounts
    2738    of data around.  It assumes a derived class wants to overwrite
    2739    an empty model with a real one - while it does an algorithm.
    2740    This is same as ClpModel one, but sets scaling on etc. */
    2741 void
    2742 ClpInterior::borrowModel(ClpModel & otherModel)
    2743 {
    2744   ClpModel::borrowModel(otherModel);
    2745   createStatus();
    2746   ClpDualRowSteepest steep1;
    2747   setDualRowPivotAlgorithm(steep1);
    2748   ClpPrimalColumnSteepest steep2;
    2749   setPrimalColumnPivotAlgorithm(steep2);
    2750 }
    2751 typedef struct {
    2752   double optimizationDirection;
    2753   double dblParam[ClpLastDblParam];
    2754   double objectiveValue;
    2755   double dualBound;
    2756   double dualTolerance;
    2757   double primalTolerance;
    2758   double sumDualInfeasibilities;
    2759   double sumPrimalInfeasibilities;
    2760   double infeasibilityCost;
    2761   int numberRows;
    2762   int numberColumns;
    2763   int intParam[ClpLastIntParam];
    2764   int numberIterations;
    2765   int problemStatus;
    2766   int maximumIterations;
    2767   int lengthNames;
    2768   int numberDualInfeasibilities;
    2769   int numberDualInfeasibilitiesWithoutFree;
    2770   int numberPrimalInfeasibilities;
    2771   int numberRefinements;
    2772   int scalingFlag;
    2773   int algorithm;
    2774   unsigned int specialOptions;
    2775   int dualPivotChoice;
    2776   int primalPivotChoice;
    2777   int matrixStorageChoice;
    2778 } Clp_scalars;
    2779 
    2780 int outDoubleArray(double * array, int length, FILE * fp)
    2781 {
    2782   int numberWritten;
    2783   if (array&&length) {
    2784     numberWritten = fwrite(&length,sizeof(int),1,fp);
    2785     if (numberWritten!=1)
    2786       return 1;
    2787     numberWritten = fwrite(array,sizeof(double),length,fp);
    2788     if (numberWritten!=length)
    2789       return 1;
    2790   } else {
    2791     length = 0;
    2792     numberWritten = fwrite(&length,sizeof(int),1,fp);
    2793     if (numberWritten!=1)
    2794       return 1;
    2795   }
    2796   return 0;
    2797 }
    2798 // Save model to file, returns 0 if success
    2799 int
    2800 ClpInterior::saveModel(const char * fileName)
    2801 {
    2802   FILE * fp = fopen(fileName,"wb");
    2803   if (fp) {
    2804     Clp_scalars scalars;
    2805     int i;
    2806     CoinBigIndex numberWritten;
    2807     // Fill in scalars
    2808     scalars.optimizationDirection = optimizationDirection_;
    2809     memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
    2810     scalars.objectiveValue = objectiveValue_;
    2811     scalars.dualBound = dualBound_;
    2812     scalars.dualTolerance = dualTolerance_;
    2813     scalars.primalTolerance = primalTolerance_;
    2814     scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
    2815     scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
    2816     scalars.infeasibilityCost = infeasibilityCost_;
    2817     scalars.numberRows = numberRows_;
    2818     scalars.numberColumns = numberColumns_;
    2819     memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
    2820     scalars.numberIterations = numberIterations_;
    2821     scalars.problemStatus = problemStatus_;
    2822     scalars.maximumIterations = maximumIterations();
    2823     scalars.lengthNames = lengthNames_;
    2824     scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
    2825     scalars.numberDualInfeasibilitiesWithoutFree
    2826       = numberDualInfeasibilitiesWithoutFree_;
    2827     scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
    2828     scalars.numberRefinements = numberRefinements_;
    2829     scalars.scalingFlag = scalingFlag_;
    2830     scalars.algorithm = algorithm_;
    2831     scalars.specialOptions = specialOptions_;
    2832     scalars.dualPivotChoice = dualRowPivot_->type();
    2833     scalars.primalPivotChoice = primalColumnPivot_->type();
    2834     scalars.matrixStorageChoice = matrix_->type();
    2835 
    2836     // put out scalars
    2837     numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
    2838     if (numberWritten!=1)
    2839       return 1;
    2840     // strings
    2841     CoinBigIndex length;
    2842     for (i=0;i<ClpLastStrParam;i++) {
    2843       length = strParam_[i].size();
    2844       numberWritten = fwrite(&length,sizeof(int),1,fp);
    2845       if (numberWritten!=1)
    2846         return 1;
    2847       if (length) {
    2848         numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
    2849         if (numberWritten!=1)
    2850           return 1;
    2851       }
    2852     }
    2853     // arrays - in no particular order
    2854     if (outDoubleArray(rowActivity_,numberRows_,fp))
    2855         return 1;
    2856     if (outDoubleArray(columnActivity_,numberColumns_,fp))
    2857         return 1;
    2858     if (outDoubleArray(dual_,numberRows_,fp))
    2859         return 1;
    2860     if (outDoubleArray(reducedCost_,numberColumns_,fp))
    2861         return 1;
    2862     if (outDoubleArray(rowLower_,numberRows_,fp))
    2863         return 1;
    2864     if (outDoubleArray(rowUpper_,numberRows_,fp))
    2865         return 1;
    2866     if (outDoubleArray(objective(),numberColumns_,fp))
    2867         return 1;
    2868     if (outDoubleArray(rowObjective_,numberRows_,fp))
    2869         return 1;
    2870     if (outDoubleArray(columnLower_,numberColumns_,fp))
    2871         return 1;
    2872     if (outDoubleArray(columnUpper_,numberColumns_,fp))
    2873         return 1;
    2874     if (ray_) {
    2875       if (problemStatus_==1)
    2876         if (outDoubleArray(ray_,numberRows_,fp))
    2877           return 1;
    2878       else if (problemStatus_==2)
    2879         if (outDoubleArray(ray_,numberColumns_,fp))
    2880           return 1;
    2881       else
    2882         if (outDoubleArray(NULL,0,fp))
    2883           return 1;
    2884     } else {
    2885       if (outDoubleArray(NULL,0,fp))
    2886         return 1;
    2887     }
    2888     if (status_&&(numberRows_+numberColumns_)>0) {
    2889       length = numberRows_+numberColumns_;
    2890       numberWritten = fwrite(&length,sizeof(int),1,fp);
    2891       if (numberWritten!=1)
    2892         return 1;
    2893       numberWritten = fwrite(status_,sizeof(char),length, fp);
    2894       if (numberWritten!=length)
    2895         return 1;
    2896     } else {
    2897       length = 0;
    2898       numberWritten = fwrite(&length,sizeof(int),1,fp);
    2899       if (numberWritten!=1)
    2900         return 1;
    2901     }
    2902     if (lengthNames_) {
    2903       char * array =
    2904         new char[max(numberRows_,numberColumns_)*(lengthNames_+1)];
    2905       char * put = array;
    2906       assert (numberRows_ == (int) rowNames_.size());
    2907       for (i=0;i<numberRows_;i++) {
    2908         assert((int)rowNames_[i].size()<=lengthNames_);
    2909         strcpy(put,rowNames_[i].c_str());
    2910         put += lengthNames_+1;
    2911       }
    2912       numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
    2913       if (numberWritten!=numberRows_)
    2914         return 1;
    2915       put=array;
    2916       assert (numberColumns_ == (int) columnNames_.size());
    2917       for (i=0;i<numberColumns_;i++) {
    2918         assert((int) columnNames_[i].size()<=lengthNames_);
    2919         strcpy(put,columnNames_[i].c_str());
    2920         put += lengthNames_+1;
    2921       }
    2922       numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
    2923       if (numberWritten!=numberColumns_)
    2924         return 1;
    2925       delete [] array;
    2926     }
    2927     // just standard type at present
    2928     assert (matrix_->type()==1);
    2929     assert (matrix_->getNumCols() == numberColumns_);
    2930     assert (matrix_->getNumRows() == numberRows_);
    2931     // we are going to save with gaps
    2932     length = matrix_->getVectorStarts()[numberColumns_-1]
    2933       + matrix_->getVectorLengths()[numberColumns_-1];
    2934     numberWritten = fwrite(&length,sizeof(int),1,fp);
    2935     if (numberWritten!=1)
    2936       return 1;
    2937     numberWritten = fwrite(matrix_->getElements(),
    2938                            sizeof(double),length,fp);
    2939     if (numberWritten!=length)
    2940       return 1;
    2941     numberWritten = fwrite(matrix_->getIndices(),
    2942                            sizeof(int),length,fp);
    2943     if (numberWritten!=length)
    2944       return 1;
    2945     numberWritten = fwrite(matrix_->getVectorStarts(),
    2946                            sizeof(int),numberColumns_+1,fp);
    2947     if (numberWritten!=numberColumns_+1)
    2948       return 1;
    2949     numberWritten = fwrite(matrix_->getVectorLengths(),
    2950                            sizeof(int),numberColumns_,fp);
    2951     if (numberWritten!=numberColumns_)
    2952       return 1;
    2953     // finished
    2954     fclose(fp);
    2955     return 0;
    2956   } else {
    2957     return -1;
    2958   }
    2959 }
    2960 
    2961 int inDoubleArray(double * &array, int length, FILE * fp)
    2962 {
    2963   int numberRead;
    2964   int length2;
    2965   numberRead = fread(&length2,sizeof(int),1,fp);
    2966   if (numberRead!=1)
    2967     return 1;
    2968   if (length2) {
    2969     // lengths must match
    2970     if (length!=length2)
    2971       return 2;
    2972     array = new double[length];
    2973     numberRead = fread(array,sizeof(double),length,fp);
    2974     if (numberRead!=length)
    2975       return 1;
    2976   }
    2977   return 0;
    2978 }
    2979 /* Restore model from file, returns 0 if success,
    2980    deletes current model */
    2981 int
    2982 ClpInterior::restoreModel(const char * fileName)
    2983 {
    2984   FILE * fp = fopen(fileName,"rb");
    2985   if (fp) {
    2986     // Get rid of current model
    2987     ClpModel::gutsOfDelete();
    2988     gutsOfDelete(0);
    2989     int i;
    2990     for (i=0;i<6;i++) {
    2991       rowArray_[i]=NULL;
    2992       columnArray_[i]=NULL;
    2993     }
    2994     // get an empty factorization so we can set tolerances etc
    2995     factorization_ = new ClpFactorization();
    2996     // Say sparse
    2997     factorization_->sparseThreshold(1);
    2998     Clp_scalars scalars;
    2999     CoinBigIndex numberRead;
    3000 
    3001     // get scalars
    3002     numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
    3003     if (numberRead!=1)
    3004       return 1;
    3005     // Fill in scalars
    3006     optimizationDirection_ = scalars.optimizationDirection;
    3007     memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
    3008     objectiveValue_ = scalars.objectiveValue;
    3009     dualBound_ = scalars.dualBound;
    3010     dualTolerance_ = scalars.dualTolerance;
    3011     primalTolerance_ = scalars.primalTolerance;
    3012     sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
    3013     sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
    3014     infeasibilityCost_ = scalars.infeasibilityCost;
    3015     numberRows_ = scalars.numberRows;
    3016     numberColumns_ = scalars.numberColumns;
    3017     memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
    3018     numberIterations_ = scalars.numberIterations;
    3019     problemStatus_ = scalars.problemStatus;
    3020     setMaximumIterations(scalars.maximumIterations);
    3021     lengthNames_ = scalars.lengthNames;
    3022     numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
    3023     numberDualInfeasibilitiesWithoutFree_
    3024       = scalars.numberDualInfeasibilitiesWithoutFree;
    3025     numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
    3026     numberRefinements_ = scalars.numberRefinements;
    3027     scalingFlag_ = scalars.scalingFlag;
    3028     algorithm_ = scalars.algorithm;
    3029     specialOptions_ = scalars.specialOptions;
    3030     // strings
    3031     CoinBigIndex length;
    3032     for (i=0;i<ClpLastStrParam;i++) {
    3033       numberRead = fread(&length,sizeof(int),1,fp);
    3034       if (numberRead!=1)
    3035         return 1;
    3036       if (length) {
    3037         char * array = new char[length+1];
    3038         numberRead = fread(array,length,1,fp);
    3039         if (numberRead!=1)
    3040           return 1;
    3041         array[length]='\0';
    3042         strParam_[i]=array;
    3043         delete [] array;
    3044       }
    3045     }
    3046     // arrays - in no particular order
    3047     if (inDoubleArray(rowActivity_,numberRows_,fp))
    3048         return 1;
    3049     if (inDoubleArray(columnActivity_,numberColumns_,fp))
    3050         return 1;
    3051     if (inDoubleArray(dual_,numberRows_,fp))
    3052         return 1;
    3053     if (inDoubleArray(reducedCost_,numberColumns_,fp))
    3054         return 1;
    3055     if (inDoubleArray(rowLower_,numberRows_,fp))
    3056         return 1;
    3057     if (inDoubleArray(rowUpper_,numberRows_,fp))
    3058         return 1;
    3059     double * objective;
    3060     if (inDoubleArray(objective,numberColumns_,fp))
    3061         return 1;
    3062     delete objective_;
    3063     objective_ = new ClpLinearObjective(objective,numberColumns_);
    3064     delete [] objective;
    3065     if (inDoubleArray(rowObjective_,numberRows_,fp))
    3066         return 1;
    3067     if (inDoubleArray(columnLower_,numberColumns_,fp))
    3068         return 1;
    3069     if (inDoubleArray(columnUpper_,numberColumns_,fp))
    3070         return 1;
    3071     if (problemStatus_==1) {
    3072       if (inDoubleArray(ray_,numberRows_,fp))
    3073         return 1;
    3074     } else if (problemStatus_==2) {
    3075       if (inDoubleArray(ray_,numberColumns_,fp))
    3076         return 1;
    3077     } else {
    3078       // ray should be null
    3079       numberRead = fread(&length,sizeof(int),1,fp);
    3080       if (numberRead!=1)
    3081         return 1;
    3082       if (length)
    3083         return 2;
    3084     }
    3085     delete [] status_;
    3086     status_=NULL;
    3087     // status region
    3088     numberRead = fread(&length,sizeof(int),1,fp);
    3089     if (numberRead!=1)
    3090         return 1;
    3091     if (length) {
    3092       if (length!=numberRows_+numberColumns_)
    3093         return 1;
    3094       status_ = new char unsigned[length];
    3095       numberRead = fread(status_,sizeof(char),length, fp);
    3096       if (numberRead!=length)
    3097         return 1;
    3098     }
    3099     if (lengthNames_) {
    3100       char * array =
    3101         new char[max(numberRows_,numberColumns_)*(lengthNames_+1)];
    3102       char * get = array;
    3103       numberRead = fread(array,lengthNames_+1,numberRows_,fp);
    3104       if (numberRead!=numberRows_)
    3105         return 1;
    3106       rowNames_ = std::vector<std::string> ();
    3107       rowNames_.resize(numberRows_);
    3108       for (i=0;i<numberRows_;i++) {
    3109         rowNames_[i]=get;
    3110         get += lengthNames_+1;
    3111       }
    3112       get = array;
    3113       numberRead = fread(array,lengthNames_+1,numberColumns_,fp);
    3114       if (numberRead!=numberColumns_)
    3115         return 1;
    3116       columnNames_ = std::vector<std::string> ();
    3117       columnNames_.resize(numberColumns_);
    3118       for (i=0;i<numberColumns_;i++) {
    3119         columnNames_[i]=get;
    3120         get += lengthNames_+1;
    3121       }
    3122       delete [] array;
    3123     }
    3124     // Pivot choices
    3125     assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
    3126     delete dualRowPivot_;
    3127     switch ((scalars.dualPivotChoice&63)) {
    3128     default:
    3129       printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63);
    3130     case 1:
    3131       // Dantzig
    3132       dualRowPivot_ = new ClpDualRowDantzig();
    3133       break;
    3134     case 2:
    3135       // Steepest - use mode
    3136       dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
    3137       break;
    3138     }
    3139     assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
    3140     delete primalColumnPivot_;
    3141     switch ((scalars.primalPivotChoice&63)) {
    3142     default:
    3143       printf("Need another primalPivot case %d\n",
    3144              scalars.primalPivotChoice&63);
    3145     case 1:
    3146       // Dantzig
    3147       primalColumnPivot_ = new ClpPrimalColumnDantzig();
    3148       break;
    3149     case 2:
    3150       // Steepest - use mode
    3151       primalColumnPivot_
    3152         = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
    3153       break;
    3154     }
    3155     assert(scalars.matrixStorageChoice==1);
    3156     delete matrix_;
    3157     // get arrays
    3158     numberRead = fread(&length,sizeof(int),1,fp);
    3159     if (numberRead!=1)
    3160       return 1;
    3161     double * elements = new double[length];
    3162     int * indices = new int[length];
    3163     CoinBigIndex * starts = new CoinBigIndex[numberColumns_];
    3164     int * lengths = new int[numberColumns_];
    3165     numberRead = fread(elements, sizeof(double),length,fp);
    3166     if (numberRead!=length)
    3167       return 1;
    3168     numberRead = fread(indices, sizeof(int),length,fp);
    3169     if (numberRead!=length)
    3170       return 1;
    3171     numberRead = fread(starts, sizeof(int),numberColumns_+1,fp);
    3172     if (numberRead!=numberColumns_+1)
    3173       return 1;
    3174     numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
    3175     if (numberRead!=numberColumns_)
    3176       return 1;
    3177     // assign matrix
    3178     CoinPackedMatrix * matrix = new CoinPackedMatrix();
    3179     matrix->assignMatrix(true, numberRows_, numberColumns_,
    3180                          length, elements, indices, starts, lengths);
    3181     // and transfer to Clp
    3182     matrix_ = new ClpPackedMatrix(matrix);
    3183     // finished
    3184     fclose(fp);
    3185     return 0;
    3186   } else {
    3187     return -1;
    3188   }
    3189   return 0;
    3190 }
    3191 // value of incoming variable (in Dual)
    3192 double
    3193 ClpInterior::valueIncomingDual() const
    3194 {
    3195   // Need value of incoming for list of infeasibilities as may be infeasible
    3196   double valueIncoming = (dualOut_/alpha_)*directionOut_;
    3197   if (directionIn_==-1)
    3198     valueIncoming = upperIn_-valueIncoming;
    3199   else
    3200     valueIncoming = lowerIn_-valueIncoming;
    3201   return valueIncoming;
     175ClpInterior::deleteWorkingData()
     176{
    3202177}
    3203178// Sanity check on input data - returns true if okay
     
    3229204  largestObj=0.0;
    3230205  // If bounds are too close - fix
    3231   double fixTolerance = 1.1*primalTolerance_;
     206  double fixTolerance = 1.1*primalTolerance();
    3232207  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
    3233208    double value;
     
    3244219    }
    3245220    value=upper_[i]-lower_[i];
    3246     if (value<-primalTolerance_) {
     221    if (value<-primalTolerance()) {
    3247222      numberBad++;
    3248223      if (firstBad<0)
     
    3296271    }
    3297272    value=upper_[i]-lower_[i];
    3298     if (value<-primalTolerance_) {
     273    if (value<-primalTolerance()) {
    3299274      numberBad++;
    3300275      if (firstBad<0)
     
    3349324  return true;
    3350325}
    3351 // Set up status array (for OsiClp)
    3352 void
    3353 ClpInterior::createStatus()
    3354 {
    3355   if(!status_)
    3356     status_ = new unsigned char [numberColumns_+numberRows_];
    3357   memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
    3358   int i;
    3359   // set column status to one nearest zero
    3360   for (i=0;i<numberColumns_;i++) {
    3361 #if 0
    3362     if (columnLower_[i]>=0.0) {
    3363       setColumnStatus(i,atLowerBound);
    3364     } else if (columnUpper_[i]<=0.0) {
    3365       setColumnStatus(i,atUpperBound);
    3366     } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
    3367       // free
    3368       setColumnStatus(i,isFree);
    3369     } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
    3370       setColumnStatus(i,atLowerBound);
    3371     } else {
    3372       setColumnStatus(i,atUpperBound);
    3373     }
    3374 #else
    3375     setColumnStatus(i,atLowerBound);
    3376 #endif
    3377   }
    3378   for (i=0;i<numberRows_;i++) {
    3379     setRowStatus(i,basic);
    3380   }
    3381 }
    3382 /* Loads a problem (the constraints on the
    3383    rows are given by lower and upper bounds). If a pointer is 0 then the
    3384    following values are the default:
    3385    <ul>
    3386    <li> <code>colub</code>: all columns have upper bound infinity
    3387    <li> <code>collb</code>: all columns have lower bound 0
    3388    <li> <code>rowub</code>: all rows have upper bound infinity
    3389    <li> <code>rowlb</code>: all rows have lower bound -infinity
    3390    <li> <code>obj</code>: all variables have 0 objective coefficient
    3391    </ul>
    3392 */
    3393 void
    3394 ClpInterior::loadProblem (  const ClpMatrixBase& matrix,
    3395                     const double* collb, const double* colub,   
    3396                     const double* obj,
    3397                     const double* rowlb, const double* rowub,
    3398                     const double * rowObjective)
    3399 {
    3400   ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
    3401                         rowObjective);
    3402   createStatus();
    3403 }
    3404 void
    3405 ClpInterior::loadProblem (  const CoinPackedMatrix& matrix,
    3406                     const double* collb, const double* colub,   
    3407                     const double* obj,
    3408                     const double* rowlb, const double* rowub,
    3409                     const double * rowObjective)
    3410 {
    3411   ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
    3412                         rowObjective);
    3413   createStatus();
    3414 }
    3415 
    3416 /* Just like the other loadProblem() method except that the matrix is
    3417    given in a standard column major ordered format (without gaps). */
    3418 void
    3419 ClpInterior::loadProblem (  const int numcols, const int numrows,
    3420                     const CoinBigIndex* start, const int* index,
    3421                     const double* value,
    3422                     const double* collb, const double* colub,   
    3423                     const double* obj,
    3424                     const double* rowlb, const double* rowub,
    3425                     const double * rowObjective)
    3426 {
    3427   ClpModel::loadProblem(numcols, numrows, start, index, value,
    3428                           collb, colub, obj, rowlb, rowub,
    3429                           rowObjective);
    3430   createStatus();
    3431 }
    3432 void
    3433 ClpInterior::loadProblem (  const int numcols, const int numrows,
    3434                            const CoinBigIndex* start, const int* index,
    3435                            const double* value,const int * length,
    3436                            const double* collb, const double* colub,   
    3437                            const double* obj,
    3438                            const double* rowlb, const double* rowub,
    3439                            const double * rowObjective)
    3440 {
    3441   ClpModel::loadProblem(numcols, numrows, start, index, value, length,
    3442                           collb, colub, obj, rowlb, rowub,
    3443                           rowObjective);
    3444   createStatus();
    3445 }
    3446 // Read an mps file from the given filename
    3447 int
    3448 ClpInterior::readMps(const char *filename,
    3449             bool keepNames,
    3450             bool ignoreErrors)
    3451 {
    3452   int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
    3453   createStatus();
    3454   return status;
    3455 }
    3456 // Just check solution (for external use)
    3457 void
    3458 ClpInterior::checkSolution()
    3459 {
    3460   // put in standard form
    3461   createRim(7+8+16);
    3462   dualTolerance_=dblParam_[ClpDualTolerance];
    3463   primalTolerance_=dblParam_[ClpPrimalTolerance];
    3464   checkPrimalSolution( rowActivityWork_, columnActivityWork_);
    3465   checkDualSolution();
    3466   if (!numberDualInfeasibilities_&&
    3467       !numberPrimalInfeasibilities_)
    3468     problemStatus_=0;
    3469   else
    3470     problemStatus_=-1;
    3471 #ifdef CLP_DEBUG
    3472   int i;
    3473   double value=0.0;
    3474   for (i=0;i<numberRows_+numberColumns_;i++)
    3475     value += dj_[i]*solution_[i];
    3476   printf("dual value %g, primal %g\n",value,objectiveValue());
    3477 #endif
    3478   // release extra memory
    3479   deleteRim(0);
    3480 }
    3481 /* Crash - at present just aimed at dual, returns
    3482    -2 if dual preferred and crash basis created
    3483    -1 if dual preferred and all slack basis preferred
    3484    0 if basis going in was not all slack
    3485    1 if primal preferred and all slack basis preferred
    3486    2 if primal preferred and crash basis created.
    3487    
    3488    if gap between bounds <="gap" variables can be flipped
    3489    
    3490    If "pivot" is
    3491    0 No pivoting (so will just be choice of algorithm)
    3492    1 Simple pivoting e.g. gub
    3493    2 Mini iterations
    3494 */
    3495 int
    3496 ClpInterior::crash(double gap,int pivot)
    3497 {
    3498   assert(!rowObjective_); // not coded
    3499   int iColumn;
    3500   int numberBad=0;
    3501   int numberBasic=0;
    3502   double dualTolerance=dblParam_[ClpDualTolerance];
    3503   //double primalTolerance=dblParam_[ClpPrimalTolerance];
    3504 
    3505   // If no basis then make all slack one
    3506   if (!status_)
    3507     createStatus();
    3508  
    3509   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3510     if (getColumnStatus(iColumn)==basic)
    3511       numberBasic++;
    3512   }
    3513   if (numberBasic) {
    3514     // not all slack
    3515     return 0;
    3516   } else {
    3517     double * dj = new double [numberColumns_];
    3518     double * solution = columnActivity_;
    3519     const double * linearObjective = objective();
    3520     //double objectiveValue=0.0;
    3521     int iColumn;
    3522     double direction = optimizationDirection_;
    3523     // direction is actually scale out not scale in
    3524     if (direction)
    3525       direction = 1.0/direction;
    3526     for (iColumn=0;iColumn<numberColumns_;iColumn++)
    3527       dj[iColumn] = direction*linearObjective[iColumn];
    3528     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3529       // assume natural place is closest to zero
    3530       double lowerBound = columnLower_[iColumn];
    3531       double upperBound = columnUpper_[iColumn];
    3532       if (lowerBound>-1.0e20||upperBound<1.0e20) {
    3533         bool atLower;
    3534         if (fabs(upperBound)<fabs(lowerBound)) {
    3535           atLower=false;
    3536           setColumnStatus(iColumn,atUpperBound);
    3537           solution[iColumn]=upperBound;
    3538         } else {
    3539           atLower=true;
    3540           setColumnStatus(iColumn,atLowerBound);
    3541           solution[iColumn]=lowerBound;
    3542         }
    3543         if (dj[iColumn]<0.0) {
    3544           // should be at upper bound
    3545           if (atLower) {
    3546             // can we flip
    3547             if (upperBound-lowerBound<=gap) {
    3548               columnActivity_[iColumn]=upperBound;
    3549               setColumnStatus(iColumn,atUpperBound);
    3550             } else if (dj[iColumn]<-dualTolerance) {
    3551               numberBad++;
    3552             }
    3553           }
    3554         } else if (dj[iColumn]>0.0) {
    3555           // should be at lower bound
    3556           if (!atLower) {
    3557             // can we flip
    3558             if (upperBound-lowerBound<=gap) {
    3559               columnActivity_[iColumn]=lowerBound;
    3560               setColumnStatus(iColumn,atLowerBound);
    3561             } else if (dj[iColumn]>dualTolerance) {
    3562               numberBad++;
    3563             }
    3564           }
    3565         }
    3566       } else {
    3567         // free
    3568         setColumnStatus(iColumn,isFree);
    3569         if (fabs(dj[iColumn])>dualTolerance)
    3570           numberBad++;
    3571       }
    3572     }
    3573     if (numberBad||pivot) {
    3574       if (!pivot) {
    3575         delete [] dj;
    3576         return 1;
    3577       } else {
    3578         // see if can be made dual feasible with gubs etc
    3579         double * pi = new double[numberRows_];
    3580         memset (pi,0,numberRows_*sizeof(double));
    3581         int * way = new int[numberColumns_];
    3582         int numberIn = 0;
    3583 
    3584         // Get column copy
    3585         CoinPackedMatrix * columnCopy = matrix();
    3586         // Get a row copy in standard format
    3587         CoinPackedMatrix copy;
    3588         copy.reverseOrderedCopyOf(*columnCopy);
    3589         // get matrix data pointers
    3590         const int * column = copy.getIndices();
    3591         const CoinBigIndex * rowStart = copy.getVectorStarts();
    3592         const int * rowLength = copy.getVectorLengths();
    3593         const double * elementByRow = copy.getElements();
    3594         //const int * row = columnCopy->getIndices();
    3595         //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
    3596         //const int * columnLength = columnCopy->getVectorLengths();
    3597         //const double * element = columnCopy->getElements();
    3598 
    3599 
    3600         // if equality row and bounds mean artificial in basis bad
    3601         // then do anyway
    3602 
    3603         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3604           // - if we want to reduce dj, + if we want to increase
    3605           int thisWay = 100;
    3606           double lowerBound = columnLower_[iColumn];
    3607           double upperBound = columnUpper_[iColumn];
    3608           if (upperBound>lowerBound) {
    3609             switch(getColumnStatus(iColumn)) {
    3610              
    3611             case basic:
    3612               thisWay=0;
    3613             case ClpInterior::isFixed:
    3614               break;
    3615             case isFree:
    3616             case superBasic:
    3617               if (dj[iColumn]<-dualTolerance)
    3618                 thisWay = 1;
    3619               else if (dj[iColumn]>dualTolerance)
    3620                 thisWay = -1;
    3621               else
    3622                 thisWay =0;
    3623               break;
    3624             case atUpperBound:
    3625               if (dj[iColumn]>dualTolerance)
    3626                 thisWay = -1;
    3627               else if (dj[iColumn]<-dualTolerance)
    3628                 thisWay = -3;
    3629               else
    3630                 thisWay = -2;
    3631               break;
    3632             case atLowerBound:
    3633               if (dj[iColumn]<-dualTolerance)
    3634                 thisWay = 1;
    3635               else if (dj[iColumn]>dualTolerance)
    3636                 thisWay = 3;
    3637               else
    3638                 thisWay = 2;
    3639               break;
    3640             }
    3641           }
    3642           way[iColumn] = thisWay;
    3643         }
    3644         /*if (!numberBad)
    3645           printf("Was dual feasible before passes - rows %d\n",
    3646           numberRows_);*/
    3647         int lastNumberIn = -100000;
    3648         int numberPasses=5;
    3649         while (numberIn>lastNumberIn+numberRows_/100) {
    3650           lastNumberIn = numberIn;
    3651           // we need to maximize chance of doing good
    3652           int iRow;
    3653           for (iRow=0;iRow<numberRows_;iRow++) {
    3654             double lowerBound = rowLower_[iRow];
    3655             double upperBound = rowUpper_[iRow];
    3656             if (getRowStatus(iRow)==basic) {
    3657               // see if we can find a column to pivot on
    3658               int j;
    3659               // down is amount pi can go down
    3660               double maximumDown = COIN_DBL_MAX;
    3661               double maximumUp = COIN_DBL_MAX;
    3662               double minimumDown =0.0;
    3663               double minimumUp =0.0;
    3664               int iUp=-1;
    3665               int iDown=-1;
    3666               int iUpB=-1;
    3667               int iDownB=-1;
    3668               if (lowerBound<-1.0e20)
    3669                 maximumUp = -1.0;
    3670               if (upperBound>1.0e20)
    3671                 maximumDown = -1.0;
    3672               for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
    3673                 int iColumn = column[j];
    3674                 double value = elementByRow[j];
    3675                 double djValue = dj[iColumn];
    3676                 /* way -
    3677                    -3 - okay at upper bound with negative dj
    3678                    -2 - marginal at upper bound with zero dj - can only decrease
    3679                    -1 - bad at upper bound
    3680                    0 - we can never pivot on this row
    3681                    1 - bad at lower bound
    3682                    2 - marginal at lower bound with zero dj - can only increase
    3683                    3 - okay at lower bound with positive dj
    3684                    100 - fine we can just ignore
    3685                 */
    3686                 if (way[iColumn]!=100) {
    3687                   switch(way[iColumn]) {
    3688                    
    3689                   case -3:
    3690                     if (value>0.0) {
    3691                       if (maximumDown*value>-djValue) {
    3692                         maximumDown = -djValue/value;
    3693                         iDown = iColumn;
    3694                       }
    3695                     } else {
    3696                       if (-maximumUp*value>-djValue) {
    3697                         maximumUp = djValue/value;
    3698                         iUp = iColumn;
    3699                       }
    3700                     }
    3701                     break;
    3702                   case -2:
    3703                     if (value>0.0) {
    3704                       maximumDown = 0.0;
    3705                     } else {
    3706                       maximumUp = 0.0;
    3707                     }
    3708                     break;
    3709                   case -1:
    3710                     // see if could be satisfied
    3711                     // dj value > 0
    3712                     if (value>0.0) {
    3713                       maximumDown=0.0;
    3714                       if (maximumUp*value<djValue-dualTolerance) {
    3715                         maximumUp = 0.0; // would improve but not enough
    3716                       } else {
    3717                         if (minimumUp*value<djValue) {
    3718                           minimumUp = djValue/value;
    3719                           iUpB = iColumn;
    3720                         }
    3721                       }
    3722                     } else {
    3723                       maximumUp=0.0;
    3724                       if (-maximumDown*value<djValue-dualTolerance) {
    3725                         maximumDown = 0.0; // would improve but not enough
    3726                       } else {
    3727                         if (-minimumDown*value<djValue) {
    3728                           minimumDown = -djValue/value;
    3729                           iDownB = iColumn;
    3730                         }
    3731                       }
    3732                     }
    3733                    
    3734                     break;
    3735                   case 0:
    3736                     maximumDown = -1.0;
    3737                     maximumUp=-1.0;
    3738                     break;
    3739                   case 1:
    3740                     // see if could be satisfied
    3741                     // dj value < 0
    3742                     if (value>0.0) {
    3743                       maximumUp=0.0;
    3744                       if (maximumDown*value<-djValue-dualTolerance) {
    3745                         maximumDown = 0.0; // would improve but not enough
    3746                       } else {
    3747                         if (minimumDown*value<-djValue) {
    3748                           minimumDown = -djValue/value;
    3749                           iDownB = iColumn;
    3750                         }
    3751                       }
    3752                     } else {
    3753                       maximumDown=0.0;
    3754                       if (-maximumUp*value<-djValue-dualTolerance) {
    3755                         maximumUp = 0.0; // would improve but not enough
    3756                       } else {
    3757                         if (-minimumUp*value<-djValue) {
    3758                           minimumUp = djValue/value;
    3759                           iUpB = iColumn;
    3760                         }
    3761                       }
    3762                     }
    3763                    
    3764                     break;
    3765                   case 2:
    3766                     if (value>0.0) {
    3767                       maximumUp = 0.0;
    3768                     } else {
    3769                       maximumDown = 0.0;
    3770                     }
    3771                    
    3772                     break;
    3773                   case 3:
    3774                     if (value>0.0) {
    3775                       if (maximumUp*value>djValue) {
    3776                         maximumUp = djValue/value;
    3777                         iUp = iColumn;
    3778                       }
    3779                     } else {
    3780                       if (-maximumDown*value>djValue) {
    3781                         maximumDown = -djValue/value;
    3782                         iDown = iColumn;
    3783                       }
    3784                     }
    3785                    
    3786                     break;
    3787                   default:
    3788                     break;
    3789                   }
    3790                 }
    3791               }
    3792               if (iUpB>=0)
    3793                 iUp=iUpB;
    3794               if (maximumUp<=dualTolerance||maximumUp<minimumUp)
    3795                 iUp=-1;
    3796               if (iDownB>=0)
    3797                 iDown=iDownB;
    3798               if (maximumDown<=dualTolerance||maximumDown<minimumDown)
    3799                 iDown=-1;
    3800               if (iUp>=0||iDown>=0) {
    3801                 // do something
    3802                 if (iUp>=0&&iDown>=0) {
    3803                   if (maximumDown>maximumUp)
    3804                     iUp=-1;
    3805                 }
    3806                 double change;
    3807                 int kColumn;
    3808                 if (iUp>=0) {
    3809                   kColumn=iUp;
    3810                   change=maximumUp;
    3811                   // just do minimum if was dual infeasible
    3812                   // ? only if maximum large?
    3813                   if (minimumUp>0.0)
    3814                     change=minimumUp;
    3815                   setRowStatus(iRow,atUpperBound);
    3816                 } else {
    3817                   kColumn=iDown;
    3818                   change=-maximumDown;
    3819                   // just do minimum if was dual infeasible
    3820                   // ? only if maximum large?
    3821                   if (minimumDown>0.0)
    3822                     change=-minimumDown;
    3823                   setRowStatus(iRow,atLowerBound);
    3824                 }
    3825                 assert (fabs(change)<1.0e20);
    3826                 setColumnStatus(kColumn,basic);
    3827                 numberIn++;
    3828                 pi[iRow]=change;
    3829                 for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
    3830                   int iColumn = column[j];
    3831                   double value = elementByRow[j];
    3832                   double djValue = dj[iColumn]-change*value;
    3833                   dj[iColumn]=djValue;
    3834                   if (abs(way[iColumn])==1) {
    3835                     numberBad--;
    3836                     /*if (!numberBad)
    3837                       printf("Became dual feasible at row %d out of %d\n",
    3838                       iRow, numberRows_);*/
    3839                     lastNumberIn=-1000000;
    3840                   }
    3841                   int thisWay = 100;
    3842                   double lowerBound = columnLower_[iColumn];
    3843                   double upperBound = columnUpper_[iColumn];
    3844                   if (upperBound>lowerBound) {
    3845                     switch(getColumnStatus(iColumn)) {
    3846                      
    3847                     case basic:
    3848                       thisWay=0;
    3849                     case isFixed:
    3850                       break;
    3851                     case isFree:
    3852                     case superBasic:
    3853                       if (djValue<-dualTolerance)
    3854                         thisWay = 1;
    3855                       else if (djValue>dualTolerance)
    3856                         thisWay = -1;
    3857                       else
    3858                         { thisWay =0; abort();}
    3859                       break;
    3860                     case atUpperBound:
    3861                       if (djValue>dualTolerance)
    3862                         { thisWay =-1; abort();}
    3863                       else if (djValue<-dualTolerance)
    3864                         thisWay = -3;
    3865                       else
    3866                         thisWay = -2;
    3867                       break;
    3868                     case atLowerBound:
    3869                       if (djValue<-dualTolerance)
    3870                         { thisWay =1; abort();}
    3871                       else if (djValue>dualTolerance)
    3872                         thisWay = 3;
    3873                       else
    3874                         thisWay = 2;
    3875                       break;
    3876                     }
    3877                   }
    3878                   way[iColumn] = thisWay;
    3879                 }
    3880               }
    3881             }
    3882           }
    3883           if (numberIn==lastNumberIn||numberBad||pivot<2)
    3884             break;
    3885           if (!(--numberPasses))
    3886             break;
    3887           //printf("%d put in so far\n",numberIn);
    3888         }
    3889         // last attempt to flip
    3890         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3891           double lowerBound = columnLower_[iColumn];
    3892           double upperBound = columnUpper_[iColumn];
    3893           if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
    3894             double djValue=dj[iColumn];
    3895             switch(getColumnStatus(iColumn)) {
    3896              
    3897             case basic:
    3898             case ClpInterior::isFixed:
    3899               break;
    3900             case isFree:
    3901             case superBasic:
    3902               break;
    3903             case atUpperBound:
    3904               if (djValue>dualTolerance) {
    3905                 setColumnStatus(iColumn,atUpperBound);
    3906                 solution[iColumn]=upperBound;
    3907               }
    3908               break;
    3909             case atLowerBound:
    3910               if (djValue<-dualTolerance) {
    3911                 setColumnStatus(iColumn,atUpperBound);
    3912                 solution[iColumn]=upperBound;
    3913               }
    3914               break;
    3915             }
    3916           }
    3917         }
    3918         delete [] pi;
    3919         delete [] dj;
    3920         delete [] way;
    3921         handler_->message(CLP_CRASH,messages_)
    3922           <<numberIn
    3923           <<numberBad
    3924           <<CoinMessageEol;
    3925         return -1;
    3926       }
    3927     } else {
    3928       delete [] dj;
    3929       return -1;
    3930     }
    3931   }
    3932 }
    3933 /* Pivot in a variable and out a variable.  Returns 0 if okay,
    3934    1 if inaccuracy forced re-factorization, -1 if would be singular.
    3935    Also updates primal/dual infeasibilities.
    3936    Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
    3937 */
    3938 int ClpInterior::pivot()
    3939 {
    3940   // scaling not allowed
    3941   assert (!scalingFlag_);
    3942   // assume In_ and Out_ are correct and directionOut_ set
    3943   // (or In_ if flip
    3944   lowerIn_ = lower_[sequenceIn_];
    3945   valueIn_ = solution_[sequenceIn_];
    3946   upperIn_ = upper_[sequenceIn_];
    3947   dualIn_ = dj_[sequenceIn_];
    3948   if (sequenceOut_>=0&&sequenceIn_!=sequenceIn_) {
    3949     assert (pivotRow_>=0&&pivotRow_<numberRows_);
    3950     assert (pivotVariable_[pivotRow_]==sequenceOut_);
    3951     lowerOut_ = lower_[sequenceOut_];
    3952     valueOut_ = solution_[sequenceOut_];
    3953     upperOut_ = upper_[sequenceOut_];
    3954     // for now assume primal is feasible (or in dual)
    3955     dualOut_ = dj_[sequenceOut_];
    3956     assert(fabs(dualOut_)<1.0e-6);
    3957   } else {
    3958     assert (pivotRow_<0);
    3959   }
    3960   bool roundAgain = true;
    3961   int returnCode=0;
    3962   while (roundAgain) {
    3963     roundAgain=false;
    3964     unpack(rowArray_[1]);
    3965     factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
    3966     // we are going to subtract movement from current basic
    3967     double movement;
    3968     // see where incoming will go to
    3969     if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
    3970       // flip so go to bound
    3971       movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
    3972     } else {
    3973       // get where outgoing needs to get to
    3974       double outValue = (directionOut_>0) ? upperOut_ : lowerOut_;
    3975       // solutionOut_ - movement*alpha_ == outValue
    3976       movement = (outValue-valueOut_)/alpha_;
    3977       // set directionIn_ correctly
    3978       directionIn_ = (movement>0) ? 1 :-1;
    3979     }
    3980     // update primal solution
    3981     {
    3982       int i;
    3983       int * index = rowArray_[1]->getIndices();
    3984       int number = rowArray_[1]->getNumElements();
    3985       double * element = rowArray_[1]->denseVector();
    3986       for (i=0;i<number;i++) {
    3987         int ii = index[i];
    3988         // get column
    3989         ii = pivotVariable_[ii];
    3990         solution_[ii] -= movement*element[i];
    3991         element[i]=0.0;
    3992       }
    3993       // see where something went to
    3994       if (sequenceOut_<0) {
    3995         if (directionIn_<0) {
    3996           assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
    3997           solution_[sequenceIn_]=upperIn_;
    3998         } else {
    3999           assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
    4000           solution_[sequenceIn_]=lowerIn_;
    4001         }
    4002       } else {
    4003         if (directionOut_<0) {
    4004           assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
    4005           solution_[sequenceOut_]=upperOut_;
    4006         } else {
    4007           assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
    4008           solution_[sequenceOut_]=lowerOut_;
    4009         }
    4010         solution_[sequenceIn_]=valueIn_+movement;
    4011       }
    4012     }   
    4013     double objectiveChange = dualIn_*movement;
    4014     // update duals
    4015     if (pivotRow_>=0) {
    4016       alpha_ = rowArray_[1]->denseVector()[pivotRow_];
    4017       assert (fabs(alpha_)>1.0e-8);
    4018       double multiplier = dualIn_/alpha_;
    4019       rowArray_[0]->insert(pivotRow_,multiplier);
    4020       factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
    4021       // put row of tableau in rowArray[0] and columnArray[0]
    4022       matrix_->transposeTimes(this,-1.0,
    4023                               rowArray_[0],columnArray_[1],columnArray_[0]);
    4024       // update column djs
    4025       int i;
    4026       int * index = columnArray_[0]->getIndices();
    4027       int number = columnArray_[0]->getNumElements();
    4028       double * element = columnArray_[0]->denseVector();
    4029       for (i=0;i<number;i++) {
    4030         int ii = index[i];
    4031         dj_[ii] += element[ii];
    4032         element[ii]=0.0;
    4033       }
    4034       columnArray_[0]->setNumElements(0);
    4035       // and row djs
    4036       index = rowArray_[0]->getIndices();
    4037       number = rowArray_[0]->getNumElements();
    4038       element = rowArray_[0]->denseVector();
    4039       for (i=0;i<number;i++) {
    4040         int ii = index[i];
    4041         dj_[ii+numberColumns_] += element[ii];
    4042         dual_[ii] = dj_[ii+numberColumns_];
    4043         element[ii]=0.0;
    4044       }
    4045       rowArray_[0]->setNumElements(0);
    4046       // check incoming
    4047       assert (fabs(dj_[sequenceIn_])<1.0e-6);
    4048     }
    4049    
    4050     // if stable replace in basis
    4051     int updateStatus = factorization_->replaceColumn(rowArray_[2],
    4052                                                    pivotRow_,
    4053                                                      alpha_);
    4054     bool takePivot=true;
    4055     // if no pivots, bad update but reasonable alpha - take and invert
    4056     if (updateStatus==2&&
    4057         lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
    4058       updateStatus=4;
    4059     if (updateStatus==1||updateStatus==4) {
    4060       // slight error
    4061       if (factorization_->pivots()>5||updateStatus==4) {
    4062         returnCode=-1;
    4063       }
    4064     } else if (updateStatus==2) {
    4065       // major error
    4066       rowArray_[1]->clear();
    4067       takePivot=false;
    4068       if (factorization_->pivots()) {
    4069         // refactorize here
    4070         statusOfProblem();
    4071         roundAgain=true;
    4072       } else {
    4073         returnCode=1;
    4074       }
    4075     } else if (updateStatus==3) {
    4076       // out of memory
    4077       // increase space if not many iterations
    4078       if (factorization_->pivots()<
    4079           0.5*factorization_->maximumPivots()&&
    4080           factorization_->pivots()<200)
    4081         factorization_->areaFactor(
    4082                                    factorization_->areaFactor() * 1.1);
    4083       returnCode =-1; // factorize now
    4084     }
    4085     if (takePivot) {
    4086       int save = algorithm_;
    4087       // make simple so always primal
    4088       algorithm_=1;
    4089       housekeeping(objectiveChange);
    4090       algorithm_=save;
    4091     }
    4092   }
    4093   if (returnCode == -1) {
    4094     // refactorize here
    4095     statusOfProblem();
    4096   } else {
    4097     // just for now - recompute anyway
    4098     gutsOfSolution(NULL,NULL);
    4099   }
    4100   return returnCode;
    4101 }
    4102 
    4103 /* Pivot in a variable and choose an outgoing one.  Assumes primal
    4104    feasible - will not go through a bound.  Returns step length in theta
    4105    Returns ray in ray_ (or NULL if no pivot)
    4106    Return codes as before but -1 means no acceptable pivot
    4107 */
    4108 int ClpInterior::primalPivotResult()
    4109 {
    4110   assert (sequenceIn_>=0);
    4111   valueIn_=solution_[sequenceIn_];
    4112   lowerIn_=lower_[sequenceIn_];
    4113   upperIn_=upper_[sequenceIn_];
    4114   dualIn_=dj_[sequenceIn_];
    4115 
    4116   int returnCode = ((ClpInteriorPrimal *) this)->pivotResult();
    4117   if (returnCode<0&&returnCode>-4) {
    4118     return 0;
    4119   } else {
    4120     printf("Return code of %d from ClpInteriorPrimal::pivotResult\n",
    4121            returnCode);
    4122     return -1;
    4123   }
    4124 }
    4125  
    4126 /* Pivot out a variable and choose an incoing one.  Assumes dual
    4127    feasible - will not go through a reduced cost. 
    4128    Returns step length in theta
    4129    Returns ray in ray_ (or NULL if no pivot)
    4130    Return codes as before but -1 means no acceptable pivot
    4131 */
    4132 int
    4133 ClpInterior::dualPivotResult()
    4134 {
    4135   return ((ClpInteriorDual *) this)->pivotResult();
    4136 }
    4137 // Factorization frequency
    4138 int
    4139 ClpInterior::factorizationFrequency() const
    4140 {
    4141   if (factorization_)
    4142     return factorization_->maximumPivots();
    4143   else
    4144     return -1;
    4145 }
    4146 void
    4147 ClpInterior::setFactorizationFrequency(int value)
    4148 {
    4149   if (factorization_)
    4150     factorization_->maximumPivots(value);
    4151 }
    4152 // Common bits of coding for dual and primal
    4153 int
    4154 ClpInterior::startup(int ifValuesPass)
    4155 {
    4156   // sanity check
    4157   assert (numberRows_==matrix_->getNumRows());
    4158   assert (numberColumns_==matrix_->getNumCols());
    4159   // for moment all arrays must exist
    4160   assert(columnLower_);
    4161   assert(columnUpper_);
    4162   assert(rowLower_);
    4163   assert(rowUpper_);
    4164   pivotRow_=-1;
    4165   sequenceIn_=-1;
    4166   sequenceOut_=-1;
    4167   secondaryStatus_=0;
    4168 
    4169   primalTolerance_=dblParam_[ClpPrimalTolerance];
    4170   dualTolerance_=dblParam_[ClpDualTolerance];
    4171   if (problemStatus_!=10)
    4172     numberIterations_=0;
    4173 
    4174   // put in standard form (and make row copy)
    4175   // create modifiable copies of model rim and do optional scaling
    4176   bool goodMatrix=createRim(7+8+16,true);
    4177 
    4178   if (goodMatrix) {
    4179     // Model looks okay
    4180     // Do initial factorization
    4181     // and set certain stuff
    4182     // We can either set increasing rows so ...IsBasic gives pivot row
    4183     // or we can just increment iBasic one by one
    4184     // for now let ...iBasic give pivot row
    4185     factorization_->increasingRows(2);
    4186     // row activities have negative sign
    4187     factorization_->slackValue(-1.0);
    4188     factorization_->zeroTolerance(1.0e-13);
    4189     // Switch off dense (unless special option set)
    4190     int saveThreshold = factorization_->denseThreshold();
    4191     factorization_->setDenseThreshold(0);
    4192     // If values pass then perturb (otherwise may be optimal so leave a bit)
    4193     if (ifValuesPass) {
    4194       // do perturbation if asked for
    4195      
    4196       if (perturbation_<100) {
    4197         if (algorithm_>0) {
    4198           ((ClpInteriorPrimal *) this)->perturb(0);
    4199         } else if (algorithm_<0) {
    4200         ((ClpInteriorDual *) this)->perturb();
    4201         }
    4202       }
    4203     }
    4204     // for primal we will change bounds using infeasibilityCost_
    4205     if (nonLinearCost_==NULL&&algorithm_>0) {
    4206       // get a valid nonlinear cost function
    4207       delete nonLinearCost_;
    4208       nonLinearCost_= new ClpNonLinearCost(this);
    4209     }
    4210    
    4211     // loop round to clean up solution if values pass
    4212     int numberThrownOut = -1;
    4213     int totalNumberThrownOut=0;
    4214     while(numberThrownOut) {
    4215       int status = internalFactorize(0+10*ifValuesPass);
    4216       if (status<0)
    4217         return 1; // some error
    4218       else
    4219         numberThrownOut = status;
    4220      
    4221       // for this we need clean basis so it is after factorize
    4222       if (!numberThrownOut)
    4223         numberThrownOut = gutsOfSolution(  NULL,NULL,
    4224                                          ifValuesPass!=0);
    4225       totalNumberThrownOut+= numberThrownOut;
    4226      
    4227     }
    4228    
    4229     if (totalNumberThrownOut)
    4230       handler_->message(CLP_SINGULARITIES,messages_)
    4231         <<totalNumberThrownOut
    4232         <<CoinMessageEol;
    4233     // Switch back dense
    4234     factorization_->setDenseThreshold(saveThreshold);
    4235    
    4236     problemStatus_ = -1;
    4237    
    4238     // number of times we have declared optimality
    4239     numberTimesOptimal_=0;
    4240 
    4241     return 0;
    4242   } else {
    4243     // bad matrix
    4244     return 2;
    4245   }
    4246    
    4247 }
    4248 
    4249 
    4250 void
    4251 ClpInterior::finish()
    4252 {
    4253   // Get rid of some arrays and empty factorization
    4254   deleteRim();
    4255   // Skip message if changing algorithms
    4256   if (problemStatus_!=10) {
    4257     assert(problemStatus_>=0&&problemStatus_<5);
    4258     handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
    4259       <<objectiveValue()
    4260       <<CoinMessageEol;
    4261   }
    4262   factorization_->relaxAccuracyCheck(1.0);
    4263   // get rid of any network stuff - could do more
    4264   factorization_->cleanUp();
    4265 }
    4266 // Save data
    4267 ClpDataSave
    4268 ClpInterior::saveData()
    4269 {
    4270   ClpDataSave saved;
    4271   saved.dualBound_ = dualBound_;
    4272   saved.infeasibilityCost_ = infeasibilityCost_;
    4273   saved.sparseThreshold_ = factorization_->sparseThreshold();
    4274   saved.perturbation_ = perturbation_;
    4275   // Progress indicator
    4276   delete progress_;
    4277   progress_ = new ClpInteriorProgress (this);
    4278   return saved;
    4279 }
    4280 // Restore data
    4281 void
    4282 ClpInterior::restoreData(ClpDataSave saved)
    4283 {
    4284   factorization_->sparseThreshold(saved.sparseThreshold_);
    4285   perturbation_ = saved.perturbation_;
    4286   infeasibilityCost_ = saved.infeasibilityCost_;
    4287   dualBound_ = saved.dualBound_;
    4288   delete progress_;
    4289   progress_=NULL;
    4290 }
    4291 /* Factorizes and returns true if optimal.  Used by user */
    4292 bool
    4293 ClpInterior::statusOfProblem()
    4294 {
    4295   // is factorization okay?
    4296   assert (internalFactorize(1)==0);
    4297   // put back original costs and then check
    4298   // also move to work arrays
    4299   createRim(4+32);
    4300   //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));
    4301   //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));
    4302   gutsOfSolution(NULL,NULL);
    4303   //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));
    4304   //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));
    4305   //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));
    4306   deleteRim(-1);
    4307   return (primalFeasible()&&dualFeasible());
    4308 }
    4309 /* Return model - updates any scalars */
    4310 void
    4311 ClpInterior::returnModel(ClpInterior & otherModel)
    4312 {
    4313   ClpModel::returnModel(otherModel);
    4314   otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;
    4315   otherModel.columnPrimalSequence_ = columnPrimalSequence_;
    4316   otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;
    4317   otherModel.rowPrimalSequence_ = rowPrimalSequence_;
    4318   otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
    4319   otherModel.columnDualSequence_ = columnDualSequence_;
    4320   otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
    4321   otherModel.rowDualSequence_ = rowDualSequence_;
    4322   otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
    4323   otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
    4324   otherModel.largestPrimalError_ = largestPrimalError_;
    4325   otherModel.largestDualError_ = largestDualError_;
    4326   otherModel.largestSolutionError_ = largestSolutionError_;
    4327   otherModel.alpha_ = alpha_;
    4328   otherModel.theta_ = theta_;
    4329   otherModel.lowerIn_ = lowerIn_;
    4330   otherModel.valueIn_ = valueIn_;
    4331   otherModel.upperIn_ = upperIn_;
    4332   otherModel.dualIn_ = dualIn_;
    4333   otherModel.sequenceIn_ = sequenceIn_;
    4334   otherModel.directionIn_ = directionIn_;
    4335   otherModel.lowerOut_ = lowerOut_;
    4336   otherModel.valueOut_ = valueOut_;
    4337   otherModel.upperOut_ = upperOut_;
    4338   otherModel.dualOut_ = dualOut_;
    4339   otherModel.sequenceOut_ = sequenceOut_;
    4340   otherModel.directionOut_ = directionOut_;
    4341   otherModel.pivotRow_ = pivotRow_;
    4342   otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;
    4343   otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;
    4344   otherModel.numberDualInfeasibilitiesWithoutFree_ =
    4345     numberDualInfeasibilitiesWithoutFree_;
    4346   otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;
    4347   otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
    4348   otherModel.numberTimesOptimal_ = numberTimesOptimal_;
    4349   otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
    4350   otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
    4351 }
    4352 /* Constructs a non linear cost from list of non-linearities (columns only)
    4353    First lower of each column is taken as real lower
    4354    Last lower is taken as real upper and cost ignored
    4355    
    4356    Returns nonzero if bad data e.g. lowers not monotonic
    4357 */
    4358 int
    4359 ClpInterior::createPiecewiseLinearCosts(const int * starts,
    4360                                        const double * lower, const double * gradient)
    4361 {
    4362   delete nonLinearCost_;
    4363   // Set up feasible bounds and check monotonicity
    4364   int iColumn;
    4365   int returnCode=0;
    4366 
    4367   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    4368     int iIndex = starts[iColumn];
    4369     int end = starts[iColumn+1]-1;
    4370     columnLower_[iColumn] = lower[iIndex];
    4371     columnUpper_[iColumn] = lower[end];
    4372     double value = columnLower_[iColumn];
    4373     iIndex++;
    4374     for (;iIndex<end;iIndex++) {
    4375       if (lower[iIndex]<value)
    4376         returnCode++; // not monotonic
    4377       value=lower[iIndex];
    4378     }
    4379   }
    4380   nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
    4381   specialOptions_ |= 2; // say keep
    4382   return returnCode;
    4383 }
    4384 /* For advanced use.  When doing iterative solves things can get
    4385    nasty so on values pass if incoming solution has largest
    4386    infeasibility < incomingInfeasibility throw out variables
    4387    from basis until largest infeasibility < allowedInfeasibility
    4388    or incoming largest infeasibility.
    4389    If allowedInfeasibility>= incomingInfeasibility this is
    4390    always possible altough you may end up with an all slack basis.
    4391    
    4392    Defaults are 1.0,10.0
    4393 */
    4394 void
    4395 ClpInterior::setValuesPassAction(float incomingInfeasibility,
    4396                                 float allowedInfeasibility)
    4397 {
    4398   incomingInfeasibility_=incomingInfeasibility;
    4399   allowedInfeasibility_=allowedInfeasibility;
    4400   assert(incomingInfeasibility_>=0.0);
    4401   assert(allowedInfeasibility_>=incomingInfeasibility_);
    4402 }
    4403 //#############################################################################
    4404 
    4405 ClpInteriorProgress::ClpInteriorProgress ()
    4406 {
    4407   int i;
    4408   for (i=0;i<CLP_PROGRESS;i++) {
    4409     objective_[i] = COIN_DBL_MAX;
    4410     infeasibility_[i] = -1.0; // set to an impossible value
    4411     numberInfeasibilities_[i]=-1;
    4412     iterationNumber_[i]=-1;
    4413   }
    4414   for (i=0;i<CLP_CYCLE;i++) {
    4415     in_[i]=-1;
    4416     out_[i]=-1;
    4417     way_[i]=0;
    4418   }
    4419   numberTimes_ = 0;
    4420   numberBadTimes_ = 0;
    4421   model_ = NULL;
    4422 }
    4423 
    4424 
    4425 //-----------------------------------------------------------------------------
    4426 
    4427 ClpInteriorProgress::~ClpInteriorProgress ()
    4428 {
    4429 }
    4430 // Copy constructor.
    4431 ClpInteriorProgress::ClpInteriorProgress(const ClpInteriorProgress &rhs)
    4432 {
    4433   int i;
    4434   for (i=0;i<CLP_PROGRESS;i++) {
    4435     objective_[i] = rhs.objective_[i];
    4436     infeasibility_[i] = rhs.infeasibility_[i];
    4437     numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
    4438     iterationNumber_[i]=rhs.iterationNumber_[i];
    4439   }
    4440   for (i=0;i<CLP_CYCLE;i++) {
    4441     in_[i]=rhs.in_[i];
    4442     out_[i]=rhs.out_[i];
    4443     way_[i]=rhs.way_[i];
    4444   }
    4445   numberTimes_ = rhs.numberTimes_;
    4446   numberBadTimes_ = rhs.numberBadTimes_;
    4447   model_ = rhs.model_;
    4448 }
    4449 // Copy constructor.from model
    4450 ClpInteriorProgress::ClpInteriorProgress(ClpInterior * model)
    4451 {
    4452   model_ = model;
    4453   int i;
    4454   for (i=0;i<CLP_PROGRESS;i++) {
    4455     if (model_->algorithm()>=0)
    4456       objective_[i] = COIN_DBL_MAX;
    4457     else
    4458       objective_[i] = -COIN_DBL_MAX;
    4459     infeasibility_[i] = -1.0; // set to an impossible value
    4460     numberInfeasibilities_[i]=-1;
    4461     iterationNumber_[i]=-1;
    4462   }
    4463   for (i=0;i<CLP_CYCLE;i++) {
    4464     in_[i]=-1;
    4465     out_[i]=-1;
    4466     way_[i]=0;
    4467   }
    4468   numberTimes_ = 0;
    4469   numberBadTimes_ = 0;
    4470 }
    4471 // Assignment operator. This copies the data
    4472 ClpInteriorProgress &
    4473 ClpInteriorProgress::operator=(const ClpInteriorProgress & rhs)
    4474 {
    4475   if (this != &rhs) {
    4476     int i;
    4477     for (i=0;i<CLP_PROGRESS;i++) {
    4478       objective_[i] = rhs.objective_[i];
    4479       infeasibility_[i] = rhs.infeasibility_[i];
    4480       numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i];
    4481       iterationNumber_[i]=rhs.iterationNumber_[i];
    4482     }
    4483     for (i=0;i<CLP_CYCLE;i++) {
    4484       in_[i]=rhs.in_[i];
    4485       out_[i]=rhs.out_[i];
    4486       way_[i]=rhs.way_[i];
    4487     }
    4488     numberTimes_ = rhs.numberTimes_;
    4489     numberBadTimes_ = rhs.numberBadTimes_;
    4490     model_ = rhs.model_;
    4491   }
    4492   return *this;
    4493 }
    4494 // Seems to be something odd about exact comparison of doubles on linux
    4495 static bool equalDouble(double value1, double value2)
    4496 {
    4497   unsigned int *i1 = (unsigned int *) &value1;
    4498   unsigned int *i2 = (unsigned int *) &value2;
    4499   if (sizeof(unsigned int)*2==sizeof(double))
    4500     return (i1[0]==i2[0]&&i1[1]==i2[1]);
    4501   else
    4502     return (i1[0]==i2[0]);
    4503 }
    4504 int
    4505 ClpInteriorProgress::looping()
    4506 {
    4507   if (!model_)
    4508     return -1;
    4509   double objective = model_->rawObjectiveValue();
    4510   double infeasibility;
    4511   int numberInfeasibilities;
    4512   int iterationNumber = model_->numberIterations();
    4513   if (model_->algorithm()<0) {
    4514     // dual
    4515     infeasibility = model_->sumPrimalInfeasibilities();
    4516     numberInfeasibilities = model_->numberPrimalInfeasibilities();
    4517   } else {
    4518     //primal
    4519     infeasibility = model_->sumDualInfeasibilities();
    4520     numberInfeasibilities = model_->numberDualInfeasibilities();
    4521   }
    4522   int i;
    4523   int numberMatched=0;
    4524   int matched=0;
    4525   int nsame=0;
    4526   for (i=0;i<CLP_PROGRESS;i++) {
    4527     bool matchedOnObjective = equalDouble(objective,objective_[i]);
    4528     bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
    4529     bool matchedOnInfeasibilities =
    4530       (numberInfeasibilities==numberInfeasibilities_[i]);
    4531    
    4532     if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
    4533       matched |= (1<<i);
    4534       // Check not same iteration
    4535       if (iterationNumber!=iterationNumber_[i]) {
    4536         numberMatched++;
    4537         // here mainly to get over compiler bug?
    4538         if (model_->messageHandler()->logLevel()>10)
    4539           printf("%d %d %d %d %d loop check\n",i,numberMatched,
    4540                  matchedOnObjective, matchedOnInfeasibility,
    4541                  matchedOnInfeasibilities);
    4542       } else {
    4543         // stuck but code should notice
    4544         nsame++;
    4545       }
    4546     }
    4547     if (i) {
    4548       objective_[i-1] = objective_[i];
    4549       infeasibility_[i-1] = infeasibility_[i];
    4550       numberInfeasibilities_[i-1]=numberInfeasibilities_[i];
    4551       iterationNumber_[i-1]=iterationNumber_[i];
    4552     }
    4553   }
    4554   objective_[CLP_PROGRESS-1] = objective;
    4555   infeasibility_[CLP_PROGRESS-1] = infeasibility;
    4556   numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
    4557   iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
    4558   if (nsame==CLP_PROGRESS)
    4559     numberMatched=CLP_PROGRESS; // really stuck
    4560   if (model_->progressFlag())
    4561     numberMatched=0;
    4562   numberTimes_++;
    4563   if (numberTimes_<10)
    4564     numberMatched=0;
    4565   // skip if just last time as may be checking something
    4566   if (matched==(1<<(CLP_PROGRESS-1)))
    4567     numberMatched=0;
    4568   if (numberMatched) {
    4569     model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
    4570       <<numberMatched
    4571       <<matched
    4572       <<numberTimes_
    4573       <<CoinMessageEol;
    4574     numberBadTimes_++;
    4575     if (numberBadTimes_<10) {
    4576       // make factorize every iteration
    4577       model_->forceFactorization(1);
    4578       if (model_->algorithm()<0) {
    4579         // dual - change tolerance
    4580         model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
    4581         // if infeasible increase dual bound
    4582         if (model_->dualBound()<1.0e19) {
    4583           model_->setDualBound(model_->dualBound()*1.1);
    4584         }
    4585       } else {
    4586         // primal - change tolerance    model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
    4587         // if infeasible increase infeasibility cost
    4588         if (model_->nonLinearCost()->numberInfeasibilities()&&
    4589             model_->infeasibilityCost()<1.0e19) {
    4590           model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
    4591         }
    4592       }
    4593       return -2;
    4594     } else {
    4595       // look at solution and maybe declare victory
    4596       if (infeasibility<1.0e-4) {
    4597         return 0;
    4598       } else {
    4599         model_->messageHandler()->message(CLP_LOOP,model_->messages())
    4600           <<CoinMessageEol;
    4601 #ifndef NDEBUG
    4602         abort();
    4603 #endif
    4604         return 3;
    4605       }
    4606     }
    4607   }
    4608   return -1;
    4609 }
    4610 // Returns previous objective (if -1) - current if (0)
    4611 double
    4612 ClpInteriorProgress::lastObjective(int back) const
    4613 {
    4614   return objective_[CLP_PROGRESS-1-back];
    4615 }
    4616 // Modify objective e.g. if dual infeasible in dual
    4617 void
    4618 ClpInteriorProgress::modifyObjective(double value)
    4619 {
    4620   objective_[CLP_PROGRESS-1]=value;
    4621 }
    4622 // Returns previous iteration number (if -1) - current if (0)
    4623 int
    4624 ClpInteriorProgress::lastIterationNumber(int back) const
    4625 {
    4626   return iterationNumber_[CLP_PROGRESS-1-back];
    4627 }
    4628 // Start check at beginning of whileIterating
    4629 void
    4630 ClpInteriorProgress::startCheck()
    4631 {
    4632   int i;
    4633   for (i=0;i<CLP_CYCLE;i++) {
    4634     in_[i]=-1;
    4635     out_[i]=-1;
    4636     way_[i]=0;
    4637   }
    4638 }
    4639 // Returns cycle length in whileIterating
    4640 int
    4641 ClpInteriorProgress::cycle(int in, int out,int wayIn,int wayOut)
    4642 {
    4643   int i;
    4644   int matched=0;
    4645   return 0;
    4646   // first see if in matches any out
    4647   for (i=1;i<CLP_CYCLE;i++) {
    4648     if (in==out_[i]) {
    4649       // even if flip then suspicious
    4650       matched=-1;
    4651       break;
    4652     }
    4653   }
    4654   if (!matched) {
    4655     // can't be cycle
    4656     for (i=0;i<CLP_CYCLE-1;i++) {
    4657       in_[i]=in_[i+1];
    4658       out_[i]=out_[i+1];
    4659       way_[i]=way_[i+1];
    4660     }
    4661   } else {
    4662     // possible cycle
    4663     matched=0;
    4664     for (i=0;i<CLP_CYCLE-1;i++) {
    4665       int k;
    4666       char wayThis = way_[i];
    4667       int inThis = in_[i];
    4668       int outThis = out_[i];
    4669       for(k=i+1;k<CLP_CYCLE;k++) {
    4670         if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
    4671           matched=k-i;
    4672         }
    4673       }
    4674       in_[i]=in_[i+1];
    4675       out_[i]=out_[i+1];
    4676       way_[i]=way_[i+1];
    4677     }
    4678   }
    4679   char way = 1-wayIn+4*(1-wayOut);
    4680   in_[CLP_CYCLE-1]=in;
    4681   out_[CLP_CYCLE-1]=out;
    4682   way_[CLP_CYCLE-1]=way;
    4683   return matched;
    4684 }
    4685 // Allow initial dense factorization
    4686 void
    4687 ClpInterior::setInitialDenseFactorization(bool onOff)
    4688 {
    4689   if (onOff)
    4690     specialOptions_ |= 8;
    4691   else
    4692     specialOptions_ &= ~8;
    4693 }
    4694 bool
    4695 ClpInterior::initialDenseFactorization() const
    4696 {
    4697   return (specialOptions_&8)!=0;
    4698 }
  • branches/pre/ClpPackedMatrix.cpp

    r212 r213  
    15471547*/
    15481548bool
    1549 ClpPackedMatrix::allElementsInRange(ClpSimplex * model,
     1549ClpPackedMatrix::allElementsInRange(ClpModel * model,
    15501550                                    double smallest, double largest)
    15511551{
  • branches/pre/Makefile.Clp

    r212 r213  
    3434LIBSRC += ClpSimplexPrimalQuadratic.cpp
    3535LIBSRC += ClpSolve.cpp
     36LIBSRC += ClpInterior.cpp
    3637# and Presolve stuff
    3738LIBSRC += ClpPresolve.cpp
  • branches/pre/include/ClpInterior.hpp

    r212 r213  
    2323class ClpInteriorProgress;
    2424
    25 /** This solves LPs using the simplex method
     25/** This solves LPs using interior point methods
    2626
    2727    It inherits from ClpModel and all its arrays are created at
    28     algorithm time. Originally I tried to work with model arrays
    29     but for simplicity of coding I changed to single arrays with
    30     structural variables then row variables.  Some coding is still
    31     based on old style and needs cleaning up.
    32 
    33     For a description of algorithms:
    34 
    35     for dual see ClpInteriorDual.hpp and at top of ClpInteriorDual.cpp
    36     for primal see ClpInteriorPrimal.hpp and at top of ClpInteriorPrimal.cpp
    37 
    38     There is an algorithm data member.  + for primal variations
    39     and - for dual variations
    40 
    41     This file also includes (at end) a very simple class ClpInteriorProgress
    42     which is where anti-looping stuff should migrate to
     28    algorithm time.
    4329
    4430*/
     
    4935
    5036public:
    51   /** enums for status of various sorts.
    52       First 4 match CoinWarmStartBasis,
    53       isFixed means fixed at lower bound and out of basis
    54   */
    55   enum Status {
    56     isFree = 0x00,
    57     basic = 0x01,
    58     atUpperBound = 0x02,
    59     atLowerBound = 0x03,
    60     superBasic = 0x04,
    61     isFixed = 0x05
    62   };
    63 
    64   enum FakeBound {
    65     noFake = 0x00,
    66     bothFake = 0x01,
    67     upperFake = 0x02,
    68     lowerFake = 0x03
    69   };
    7037
    7138  /**@name Constructors and destructor and copy */
     
    134101              bool keepNames=false,
    135102              bool ignoreErrors = false);
    136   /** Borrow model.  This is so we dont have to copy large amounts
    137       of data around.  It assumes a derived class wants to overwrite
    138       an empty model with a real one - while it does an algorithm.
    139       This is same as ClpModel one, but sets scaling on etc. */
    140   void borrowModel(ClpModel & otherModel);
    141103  //@}
    142104
    143105  /**@name Functions most useful to user */
    144106  //@{
    145   /** General solve algorithm which can do presolve.
    146       See  ClpSolve.hpp for options
    147    */
    148   int initialSolve(ClpSolve & options);
    149   /** Dual algorithm - see ClpInteriorDual.hpp for method */
    150   int dual(int ifValuesPass=0);
    151   /** Primal algorithm - see ClpInteriorPrimal.hpp for method */
    152   int primal(int ifValuesPass=0);
    153   /** Solves quadratic problem using SLP - may be used as crash
    154       for other algorithms when number of iterations small.
    155       Also exits if all problematical variables are changing
    156       less than deltaTolerance
    157   */
    158   int quadraticSLP(int numberPasses,double deltaTolerance);
    159   /// Solves quadratic using Dantzig's algorithm - primal
    160   int quadraticPrimal(int phase=2);
    161   /// Passes in factorization
    162   void setFactorization( ClpFactorization & factorization);
    163   /// Sets or unsets scaling, 0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
    164   void scaling(int mode=1);
    165   /// Gets scalingFlag
    166   inline int scalingFlag() const {return scalingFlag_;};
    167   /** Tightens primal bounds to make dual faster.  Unless
    168       fixed, bounds are slightly looser than they could be.
    169       This is to make dual go faster and is probably not needed
    170       with a presolve.  Returns non-zero if problem infeasible.
    171 
    172       Fudge for branch and bound - put bounds on columns of factor *
    173       largest value (at continuous) - should improve stability
    174       in branch and bound on infeasible branches (0.0 is off)
    175   */
    176   int tightenPrimalBounds(double factor=0.0);
    177   /** Crash - at present just aimed at dual, returns
    178       -2 if dual preferred and crash basis created
    179       -1 if dual preferred and all slack basis preferred
    180        0 if basis going in was not all slack
    181        1 if primal preferred and all slack basis preferred
    182        2 if primal preferred and crash basis created.
    183        
    184        if gap between bounds <="gap" variables can be flipped
    185 
    186        If "pivot" is
    187        0 No pivoting (so will just be choice of algorithm)
    188        1 Simple pivoting e.g. gub
    189        2 Mini iterations
    190   */
    191   int crash(double gap,int pivot);
    192   /// Sets row pivot choice algorithm in dual
    193   void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
    194   /// Sets column pivot choice algorithm in primal
    195   void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
    196   /** For strong branching.  On input lower and upper are new bounds
    197       while on output they are change in objective function values
    198       (>1.0e50 infeasible).
    199       Return code is 0 if nothing interesting, -1 if infeasible both
    200       ways and +1 if infeasible one way (check values to see which one(s))
    201       Solutions are filled in as well - even down, odd up - also
    202       status and number of iterations
    203   */
    204   int strongBranching(int numberVariables,const int * variables,
    205                       double * newLower, double * newUpper,
    206                       double ** outputSolution,
    207                       int * outputStatus, int * outputIterations,
    208                       bool stopOnFirstInfeasible=true,
    209                       bool alwaysFinish=false);
     107  /** Pdco algorithm - see ClpPdco.hpp for method */
     108  int pdco();
    210109  //@}
    211110
    212111  /**@name Needed for functionality of OsiSimplexInterface */
    213112  //@{
    214   /** Pivot in a variable and out a variable.  Returns 0 if okay,
    215       1 if inaccuracy forced re-factorization, -1 if would be singular.
    216       Also updates primal/dual infeasibilities.
    217       Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
    218   */
    219   int pivot();
    220 
    221   /** Pivot in a variable and choose an outgoing one.  Assumes primal
    222       feasible - will not go through a bound.  Returns step length in theta
    223       Returns ray in ray_ (or NULL if no pivot)
    224       Return codes as before but -1 means no acceptable pivot
    225   */
    226   int primalPivotResult();
    227  
    228   /** Pivot out a variable and choose an incoing one.  Assumes dual
    229       feasible - will not go through a reduced cost. 
    230       Returns step length in theta
    231       Returns ray in ray_ (or NULL if no pivot)
    232       Return codes as before but -1 means no acceptable pivot
    233   */
    234   int dualPivotResult();
    235 
    236   /** Common bits of coding for dual and primal.  Return s0 if okay,
    237       1 if bad matrix, 2 if very bad factorization
    238   */
    239   int startup(int ifValuesPass);
    240   void finish();
    241  
    242   /** Factorizes and returns true if optimal.  Used by user */
    243   bool statusOfProblem();
     113  /// LSQR
     114  void lsqr();
    244115  //@}
    245116
     
    248119  /// If problem is primal feasible
    249120  inline bool primalFeasible() const
    250          { return (numberPrimalInfeasibilities_==0);};
     121         { return (sumPrimalInfeasibilities_<=1.0e-5);};
    251122  /// If problem is dual feasible
    252123  inline bool dualFeasible() const
    253          { return (numberDualInfeasibilities_==0);};
    254   /// factorization
    255   inline ClpFactorization * factorization() const
    256           { return factorization_;};
    257   /// Sparsity on or off
    258   bool sparseFactorization() const;
    259   void setSparseFactorization(bool value);
    260   /// Factorization frequency
    261   int factorizationFrequency() const;
    262   void setFactorizationFrequency(int value);
    263   /// Dual bound
    264   inline double dualBound() const
    265           { return dualBound_;};
    266   void setDualBound(double value);
    267   /// Infeasibility cost
    268   inline double infeasibilityCost() const
    269           { return infeasibilityCost_;};
    270   void setInfeasibilityCost(double value);
    271   /** Amount of print out:
    272       0 - none
    273       1 - just final
    274       2 - just factorizations
    275       3 - as 2 plus a bit more
    276       4 - verbose
    277       above that 8,16,32 etc just for selective debug
    278   */
    279   /** Perturbation:
    280       -50 to +50 - perturb by this power of ten (-6 sounds good)
    281       100 - auto perturb if takes too long (1.0e-6 largest nonzero)
    282       101 - we are perturbed
    283       102 - don't try perturbing again
    284       default is 100
    285   */
    286   inline int perturbation() const
    287     { return perturbation_;};
    288   void setPerturbation(int value);
     124         { return (sumDualInfeasibilities_<=1.0e-5);};
    289125  /// Current (or last) algorithm
    290126  inline int algorithm() const
     
    296132  inline double sumDualInfeasibilities() const
    297133          { return sumDualInfeasibilities_;} ;
    298   /// Number of dual infeasibilities
    299   inline int numberDualInfeasibilities() const
    300           { return numberDualInfeasibilities_;} ;
    301134  /// Sum of primal infeasibilities
    302135  inline double sumPrimalInfeasibilities() const
    303136          { return sumPrimalInfeasibilities_;} ;
    304   /// Number of primal infeasibilities
    305   inline int numberPrimalInfeasibilities() const
    306           { return numberPrimalInfeasibilities_;} ;
    307   /** Save model to file, returns 0 if success.  This is designed for
    308       use outside algorithms so does not save iterating arrays etc.
    309   It does not save any messaging information.
    310   Does not save scaling values.
    311   It does not know about all types of virtual functions.
    312   */
    313   int saveModel(const char * fileName);
    314   /** Restore model from file, returns 0 if success,
    315       deletes current model */
    316   int restoreModel(const char * fileName);
    317  
    318   /** Just check solution (for external use) - sets sum of
    319       infeasibilities etc */
    320   void checkSolution();
    321   /// Useful row length arrays (0,1,2,3,4,5)
    322   inline CoinIndexedVector * rowArray(int index) const
    323   { return rowArray_[index];};
    324   /// Useful column length arrays (0,1,2,3,4,5)
    325   inline CoinIndexedVector * columnArray(int index) const
    326   { return columnArray_[index];};
    327137  //@}
    328138
     
    330140  /**@name Functions less likely to be useful to casual user */
    331141  //@{
    332   /** Given an existing factorization computes and checks
    333       primal and dual solutions.  Uses input arrays for variables at
    334       bounds.  Returns feasibility states */
    335   int getSolution (  const double * rowActivities,
    336                      const double * columnActivities);
    337   /** Given an existing factorization computes and checks
    338       primal and dual solutions.  Uses current problem arrays for
    339       bounds.  Returns feasibility states */
    340   int getSolution ();
    341   /** Constructs a non linear cost from list of non-linearities (columns only)
    342       First lower of each column is taken as real lower
    343       Last lower is taken as real upper and cost ignored
    344 
    345       Returns nonzero if bad data e.g. lowers not monotonic
    346   */
    347   int createPiecewiseLinearCosts(const int * starts,
    348                    const double * lower, const double * gradient);
    349   /** Return model - updates any scalars */
    350   void returnModel(ClpInterior & otherModel);
    351   /** Factorizes using current basis. 
    352       solveType - 1 iterating, 0 initial, -1 external
    353       If 10 added then in primal values pass
    354       Return codes are as from ClpFactorization unless initial factorization
    355       when total number of singularities is returned
    356   */
    357     /// Save data
    358     ClpDataSave saveData() ;
    359     /// Restore data
    360     void restoreData(ClpDataSave saved);
    361   int internalFactorize(int solveType);
    362   /// Clean up status
    363   void cleanStatus();
    364   /// Factorizes using current basis. For external use
    365   int factorize();
    366   /** Computes duals from scratch. If givenDjs then
    367       allows for nonzero basic djs */
    368   void computeDuals(double * givenDjs);
    369   /// Computes primals from scratch
    370   void computePrimals (  const double * rowActivities,
    371                      const double * columnActivities);
    372   /**
    373      Unpacks one column of the matrix into indexed array
    374      Uses sequenceIn_
    375      Also applies scaling if needed
    376   */
    377   void unpack(CoinIndexedVector * rowArray) const ;
    378   /**
    379      Unpacks one column of the matrix into indexed array
    380      Slack if sequence>= numberColumns
    381      Also applies scaling if needed
    382   */
    383   void unpack(CoinIndexedVector * rowArray,int sequence) const;
    384   /**
    385      Unpacks one column of the matrix into indexed array
    386      ** as packed vector
    387      Uses sequenceIn_
    388      Also applies scaling if needed
    389   */
    390   void unpackPacked(CoinIndexedVector * rowArray) ;
    391   /**
    392      Unpacks one column of the matrix into indexed array
    393      ** as packed vector
    394      Slack if sequence>= numberColumns
    395      Also applies scaling if needed
    396   */
    397   void unpackPacked(CoinIndexedVector * rowArray,int sequence);
    398  
    399   /**
    400       This does basis housekeeping and does values for in/out variables.
    401       Can also decide to re-factorize
    402   */
    403   int housekeeping(double objectiveChange);
    404   /** This sets largest infeasibility and most infeasible and sum
    405       and number of infeasibilities (Primal) */
    406   void checkPrimalSolution(const double * rowActivities=NULL,
    407                            const double * columnActivies=NULL);
    408   /** This sets largest infeasibility and most infeasible and sum
    409       and number of infeasibilities (Dual) */
    410   void checkDualSolution();
    411   /** For advanced use.  When doing iterative solves things can get
    412       nasty so on values pass if incoming solution has largest
    413       infeasibility < incomingInfeasibility throw out variables
    414       from basis until largest infeasibility < allowedInfeasibility
    415       or incoming largest infeasibility.
    416       If allowedInfeasibility>= incomingInfeasibility this is
    417       always possible altough you may end up with an all slack basis.
    418 
    419       Defaults are 1.0,10.0
    420   */
    421   void setValuesPassAction(float incomingInfeasibility,
    422                            float allowedInfeasibility);
    423142  //@}
    424143  /**@name Matrix times vector methods
     
    441160  /**@name most useful gets and sets */
    442161  //@{
    443   /// Worst column primal infeasibility
    444   inline double columnPrimalInfeasibility() const
    445           { return columnPrimalInfeasibility_;} ;
    446   /// Sequence of worst (-1 if feasible)
    447   inline int columnPrimalSequence() const
    448           { return columnPrimalSequence_;} ;
    449   /// Worst row primal infeasibility
    450   inline double rowPrimalInfeasibility() const
    451           { return rowPrimalInfeasibility_;} ;
    452   /// Sequence of worst (-1 if feasible)
    453   inline int rowPrimalSequence() const
    454           { return rowPrimalSequence_;} ;
    455   /** Worst column dual infeasibility (note - these may not be as meaningful
    456       if the problem is primal infeasible */
    457   inline double columnDualInfeasibility() const
    458           { return columnDualInfeasibility_;} ;
    459   /// Sequence of worst (-1 if feasible)
    460   inline int columnDualSequence() const
    461           { return columnDualSequence_;} ;
    462   /// Worst row dual infeasibility
    463   inline double rowDualInfeasibility() const
    464           { return rowDualInfeasibility_;} ;
    465   /// Sequence of worst (-1 if feasible)
    466   inline int rowDualSequence() const
    467           { return rowDualSequence_;} ;
    468   /// Primal tolerance needed to make dual feasible (<largeTolerance)
    469   inline double primalToleranceToGetOptimal() const
    470           { return primalToleranceToGetOptimal_;} ;
    471   /// Remaining largest dual infeasibility
    472   inline double remainingDualInfeasibility() const
    473           { return remainingDualInfeasibility_;} ;
    474   /// Large bound value (for complementarity etc)
    475   inline double largeValue() const
    476           { return largeValue_;} ;
    477   void setLargeValue( double value) ;
    478162  /// Largest error on Ax-b
    479163  inline double largestPrimalError() const
     
    482166  inline double largestDualError() const
    483167          { return largestDualError_;} ;
    484   /// Largest difference between input primal solution and computed
    485   inline double largestSolutionError() const
    486           { return largestSolutionError_;} ;
    487   /// Basic variables pivoting on which rows
    488   inline int * pivotVariable() const
    489           { return pivotVariable_;};
    490   /// Scaling of objective 12345.0 (auto), 0.0 (off), other user
    491   inline double objectiveScale() const
    492           { return objectiveScale_;} ;
    493   inline void setObjectiveScale(double value)
    494           { objectiveScale_ = value;} ;
    495   /// Current dual tolerance
    496   inline double currentDualTolerance() const
    497           { return dualTolerance_;} ;
    498   inline void setCurrentDualTolerance(double value)
    499           { dualTolerance_ = value;} ;
    500   /// Current primal tolerance
    501   inline double currentPrimalTolerance() const
    502           { return primalTolerance_;} ;
    503   inline void setCurrentPrimalTolerance(double value)
    504           { primalTolerance_ = value;} ;
    505   /// How many iterative refinements to do
    506   inline int numberRefinements() const
    507           { return numberRefinements_;} ;
    508   void setNumberRefinements( int value) ;
    509   /// Alpha (pivot element) for use by classes e.g. steepestedge
    510   inline double alpha() const { return alpha_;};
    511   /// Reduced cost of last incoming for use by classes e.g. steepestedge
    512   inline double dualIn() const { return dualIn_;};
    513   /// Pivot Row for use by classes e.g. steepestedge
    514   inline int pivotRow() const{ return pivotRow_;};
    515   /// value of incoming variable (in Dual)
    516   double valueIncomingDual() const;
    517168  //@}
    518169
     
    520171  /**@name protected methods */
    521172  //@{
    522   /** May change basis and then returns number changed.
    523       Computation of solutions may be overriden by given pi and solution
    524   */
    525   int gutsOfSolution ( double * givenDuals,
    526                        const double * givenPrimals,
    527                        bool valuesPass=false);
    528   /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
    529   void gutsOfDelete(int type);
     173  /// Does most of deletion
     174  void gutsOfDelete();
    530175  /// Does most of copying
    531176  void gutsOfCopy(const ClpInterior & rhs);
    532   /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
    533       1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
    534       8 bit does solution scaling in
    535       16 bit does rowArray and columnArray indexed vectors
    536       and makes row copy if wanted, also sets columnStart_ etc
    537       Also creates scaling arrays if needed.  It does scaling if needed.
    538       16 also moves solutions etc in to work arrays
    539       On 16 returns false if problem "bad" i.e. matrix or bounds bad
    540   */
    541   bool createRim(int what,bool makeRowCopy=false);
    542   /** releases above arrays and does solution scaling out.  May also
    543       get rid of factorization data -
    544       0 get rid of nothing, 1 get rid of arrays, 2 also factorization
    545   */
    546   void deleteRim(int getRidOfFactorizationData=2);
    547   /// Sanity check on input rim data (after scaling) - returns true if okay
     177  /// Returns true if data looks okay, false if not
     178  bool createWorkingData();
     179  void deleteWorkingData();
     180  /// Sanity check on input rim data
    548181  bool sanityCheck();
     182  ///  This does housekeeping
     183  int housekeeping();
    549184  //@}
    550185  public:
    551186  /**@name public methods */
    552187  //@{
    553   /** Return row or column sections - not as much needed as it
    554       once was.  These just map into single arrays */
    555   inline double * solutionRegion(int section)
    556   { if (!section) return rowActivityWork_; else return columnActivityWork_;};
    557   inline double * djRegion(int section)
    558   { if (!section) return rowReducedCost_; else return reducedCostWork_;};
    559   inline double * lowerRegion(int section)
    560   { if (!section) return rowLowerWork_; else return columnLowerWork_;};
    561   inline double * upperRegion(int section)
    562   { if (!section) return rowUpperWork_; else return columnUpperWork_;};
    563   inline double * costRegion(int section)
    564   { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
    565   /// Return region as single array
    566   inline double * solutionRegion()
    567   { return solution_;};
    568   inline double * djRegion()
    569   { return dj_;};
    570   inline double * lowerRegion()
    571   { return lower_;};
    572   inline double * upperRegion()
    573   { return upper_;};
    574   inline double * costRegion()
    575   { return cost_;};
    576   inline Status getStatus(int sequence) const
    577   {return static_cast<Status> (status_[sequence]&7);};
    578   inline void setStatus(int sequence, Status status)
    579   {
    580     unsigned char & st_byte = status_[sequence];
    581     st_byte &= ~7;
    582     st_byte |= status;
    583   };
    584   /** Normally the first factorization does sparse coding because
    585       the factorization could be singular.  This allows initial dense
    586       factorization when it is known to be safe
    587   */
    588   void setInitialDenseFactorization(bool onOff);
    589   bool  initialDenseFactorization() const;
    590   /** Return sequence In or Out */
    591   inline int sequenceIn() const
    592   {return sequenceIn_;};
    593   inline int sequenceOut() const
    594   {return sequenceOut_;};
    595   /** Set sequenceIn or Out */
    596   inline void  setSequenceIn(int sequence)
    597   { sequenceIn_=sequence;};
    598   inline void  setSequenceOut(int sequence)
    599   { sequenceOut_=sequence;};
    600   /** Return direction In or Out */
    601   inline int directionIn() const
    602   {return directionIn_;};
    603   inline int directionOut() const
    604   {return directionOut_;};
    605   /** Set directionIn or Out */
    606   inline void  setDirectionIn(int direction)
    607   { directionIn_=direction;};
    608   inline void  setDirectionOut(int direction)
    609   { directionOut_=direction;};
     188  /// Raw objective value (so always minimize)
     189  inline double rawObjectiveValue() const
     190  { return objectiveValue_;};
    610191  /// Returns 1 if sequence indicates column
    611192  inline int isColumn(int sequence) const
     
    614195  inline int sequenceWithin(int sequence) const
    615196  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
    616   /// Return row or column values
    617   inline double solution(int sequence)
    618   { return solution_[sequence];};
    619   /// Return address of row or column values
    620   inline double & solutionAddress(int sequence)
    621   { return solution_[sequence];};
    622   inline double reducedCost(int sequence)
    623    { return dj_[sequence];};
    624   inline double & reducedCostAddress(int sequence)
    625    { return dj_[sequence];};
    626   inline double lower(int sequence)
    627   { return lower_[sequence];};
    628   /// Return address of row or column lower bound
    629   inline double & lowerAddress(int sequence)
    630   { return lower_[sequence];};
    631   inline double upper(int sequence)
    632   { return upper_[sequence];};
    633   /// Return address of row or column upper bound
    634   inline double & upperAddress(int sequence)
    635   { return upper_[sequence];};
    636   inline double cost(int sequence)
    637   { return cost_[sequence];};
    638   /// Return address of row or column cost
    639   inline double & costAddress(int sequence)
    640   { return cost_[sequence];};
    641   /// Return original lower bound
    642   inline double originalLower(int iSequence) const
    643   { if (iSequence<numberColumns_) return columnLower_[iSequence]; else
    644     return rowLower_[iSequence-numberColumns_];};
    645   /// Return original lower bound
    646   inline double originalUpper(int iSequence) const
    647   { if (iSequence<numberColumns_) return columnUpper_[iSequence]; else
    648     return rowUpper_[iSequence-numberColumns_];};
    649   /// Theta (pivot change)
    650   inline double theta() const
    651   { return theta_;};
    652   /// Scaling
    653   const double * rowScale() const {return rowScale_;};
    654   const double * columnScale() const {return columnScale_;};
    655   void setRowScale(double * scale) { rowScale_ = scale;};
    656   void setColumnScale(double * scale) { columnScale_ = scale;};
    657   /// Return pointer to details of costs
    658   inline ClpNonLinearCost * nonLinearCost() const
    659   { return nonLinearCost_;};
    660   //@}
    661   /**@name status methods */
    662   //@{
    663   inline void setFakeBound(int sequence, FakeBound fakeBound)
    664   {
    665     unsigned char & st_byte = status_[sequence];
    666     st_byte &= ~24;
    667     st_byte |= fakeBound<<3;
    668   };
    669   inline FakeBound getFakeBound(int sequence) const
    670   {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
    671   inline void setRowStatus(int sequence, Status status)
    672   {
    673     unsigned char & st_byte = status_[sequence+numberColumns_];
    674     st_byte &= ~7;
    675     st_byte |= status;
    676   };
    677   inline Status getRowStatus(int sequence) const
    678   {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
    679   inline void setColumnStatus(int sequence, Status status)
    680   {
    681     unsigned char & st_byte = status_[sequence];
    682     st_byte &= ~7;
    683     st_byte |= status;
    684   };
    685   inline Status getColumnStatus(int sequence) const
    686   {return static_cast<Status> (status_[sequence]&7);};
    687   inline void setPivoted( int sequence)
    688   { status_[sequence] |= 32;};
    689   inline void clearPivoted( int sequence)
    690   { status_[sequence] &= ~32; };
    691   inline bool pivoted(int sequence) const
    692   {return (((status_[sequence]>>5)&1)!=0);};
    693   inline void setFlagged( int sequence)
    694   {
    695     status_[sequence] |= 64;
    696   };
    697   inline void clearFlagged( int sequence)
    698   {
    699     status_[sequence] &= ~64;
    700   };
    701   inline bool flagged(int sequence) const
    702   {return (((status_[sequence]>>6)&1)!=0);};
    703   /** Set up status array (can be used by OsiClp).
    704       Also can be used to set up all slack basis */
    705   void createStatus() ;
    706   inline void allSlackBasis()
    707   { createStatus();};
    708    
    709   /// So we know when to be cautious
    710   inline int lastBadIteration() const
    711   {return lastBadIteration_;};
    712   /// Progress flag - at present 0 bit says artificials out
    713   inline int progressFlag() const
    714   {return progressFlag_;};
    715   /// Force re-factorization early
    716   inline void forceFactorization(int value)
    717   { forceFactorization_ = value;};
    718   /// Raw objective value (so always minimize in primal)
    719   inline double rawObjectiveValue() const
    720   { return objectiveValue_;};
    721197  //@}
    722198
     
    730206  */
    731207  //@{
    732   /// Worst column primal infeasibility
    733   double columnPrimalInfeasibility_;
    734   /// Worst row primal infeasibility
    735   double rowPrimalInfeasibility_;
    736   /// Sequence of worst (-1 if feasible)
    737   int columnPrimalSequence_;
    738   /// Sequence of worst (-1 if feasible)
    739   int rowPrimalSequence_;
    740   /// Worst column dual infeasibility
    741   double columnDualInfeasibility_;
    742   /// Worst row dual infeasibility
    743   double rowDualInfeasibility_;
    744   /// Sequence of worst (-1 if feasible)
    745   int columnDualSequence_;
    746   /// Sequence of worst (-1 if feasible)
    747   int rowDualSequence_;
    748   /// Primal tolerance needed to make dual feasible (<largeTolerance)
    749   double primalToleranceToGetOptimal_;
    750   /// Remaining largest dual infeasibility
    751   double remainingDualInfeasibility_;
    752   /// Large bound value (for complementarity etc)
    753   double largeValue_;
    754   /// Scaling of objective
    755   double objectiveScale_;
    756208  /// Largest error on Ax-b
    757209  double largestPrimalError_;
    758210  /// Largest error on basic duals
    759211  double largestDualError_;
    760   /// Largest difference between input primal solution and computed
    761   double largestSolutionError_;
    762   /// Dual bound
    763   double dualBound_;
    764   /// Alpha (pivot element)
    765   double alpha_;
    766   /// Theta (pivot change)
    767   double theta_;
    768   /// Lower Bound on In variable
    769   double lowerIn_;
    770   /// Value of In variable
    771   double valueIn_;
    772   /// Upper Bound on In variable
    773   double upperIn_;
    774   /// Reduced cost of In variable
    775   double dualIn_;
    776   /// Lower Bound on Out variable
    777   double lowerOut_;
    778   /// Value of Out variable
    779   double valueOut_;
    780   /// Upper Bound on Out variable
    781   double upperOut_;
    782   /// Infeasibility (dual) or ? (primal) of Out variable
    783   double dualOut_;
    784   /// Current dual tolerance for algorithm
    785   double dualTolerance_;
    786   /// Current primal tolerance for algorithm
    787   double primalTolerance_;
    788212  /// Sum of dual infeasibilities
    789213  double sumDualInfeasibilities_;
    790214  /// Sum of primal infeasibilities
    791215  double sumPrimalInfeasibilities_;
    792   /// Weight assigned to being infeasible in primal
    793   double infeasibilityCost_;
    794   /// Sum of Dual infeasibilities using tolerance based on error in duals
    795   double sumOfRelaxedDualInfeasibilities_;
    796   /// Sum of Primal infeasibilities using tolerance based on error in primals
    797   double sumOfRelaxedPrimalInfeasibilities_;
    798216  /// Working copy of lower bounds (Owner of arrays below)
    799217  double * lower_;
     
    808226  /// Column upper bounds - working copy
    809227  double * columnUpperWork_;
    810   /// Working copy of objective (Owner of arrays below)
     228  /// Working copy of objective
    811229  double * cost_;
    812   /// Row objective - working copy
    813   double * rowObjectiveWork_;
    814   /// Column objective - working copy
    815   double * objectiveWork_;
    816   /// Useful row length arrays
    817   CoinIndexedVector * rowArray_[6];
    818   /// Useful column length arrays
    819   CoinIndexedVector * columnArray_[6];
    820   /// Sequence of In variable
    821   int sequenceIn_;
    822   /// Direction of In, 1 going up, -1 going down, 0 not a clude
    823   int directionIn_;
    824   /// Sequence of Out variable
    825   int sequenceOut_;
    826   /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
    827   int directionOut_;
    828   /// Pivot Row
    829   int pivotRow_;
    830   /// Last good iteration (immediately after a re-factorization)
    831   int lastGoodIteration_;
    832   /// Working copy of reduced costs (Owner of arrays below)
    833   double * dj_;
    834   /// Reduced costs of slacks not same as duals (or - duals)
    835   double * rowReducedCost_;
    836   /// Possible scaled reduced costs
    837   double * reducedCostWork_;
    838   /// Working copy of primal solution (Owner of arrays below)
    839   double * solution_;
    840   /// Row activities - working copy
    841   double * rowActivityWork_;
    842   /// Column activities - working copy
    843   double * columnActivityWork_;
    844   /// Number of dual infeasibilities
    845   int numberDualInfeasibilities_;
    846   /// Number of dual infeasibilities (without free)
    847   int numberDualInfeasibilitiesWithoutFree_;
    848   /// Number of primal infeasibilities
    849   int numberPrimalInfeasibilities_;
    850   /// How many iterative refinements to do
    851   int numberRefinements_;
    852   /// dual row pivot choice
    853   ClpDualRowPivot * dualRowPivot_;
    854   /// primal column pivot choice
    855   ClpPrimalColumnPivot * primalColumnPivot_;
    856   /// Basic variables pivoting on which rows
    857   int * pivotVariable_;
    858   /// factorization
    859   ClpFactorization * factorization_;
    860   /// Row scale factors for matrix
    861   // ****** get working simply then make coding more efficient
    862   // on full matrix operations
    863   double * rowScale_;
    864   /// Saved version of solution
    865   double * savedSolution_;
    866   /// Column scale factors
    867   double * columnScale_;
    868   /// Scale flag, 0 none, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic
    869   int scalingFlag_;
    870   /// Number of times code has tentatively thought optimal
    871   int numberTimesOptimal_;
    872   /// If change has been made (first attempt at stopping looping)
    873   int changeMade_;
    874   /// Algorithm >0 == Primal, <0 == Dual
     230  /// Which algorithm being used
    875231  int algorithm_;
    876   /** Now for some reliability aids
    877       This forces re-factorization early */
    878   int forceFactorization_;
    879   /** Perturbation:
    880       -50 to +50 - perturb by this power of ten (-6 sounds good)
    881       100 - auto perturb if takes too long (1.0e-6 largest nonzero)
    882       101 - we are perturbed
    883       102 - don't try perturbing again
    884       default is 100
    885   */
    886   int perturbation_;
    887   /// Saved status regions
    888   unsigned char * saveStatus_;
    889   /** Very wasteful way of dealing with infeasibilities in primal.
    890       However it will allow non-linearities and use of dual
    891       analysis.  If it doesn't work it can easily be replaced.
    892   */
    893   ClpNonLinearCost * nonLinearCost_;
    894   /** For advanced options
    895       1 - Don't keep changing infeasibility weight
    896       2 - Keep nonLinearCost round solves
    897       4 - Force outgoing variables to exact bound (primal)
    898       8 - Safe to use dense initial factorization
    899   */
    900   unsigned int specialOptions_;
    901   /// So we know when to be cautious
    902   int lastBadIteration_;
    903   /// Can be used for count of fake bounds (dual) or fake costs (primal)
    904   int numberFake_;
    905   /// Progress flag - at present 0 bit says artificials out, 1 free in
    906   int progressFlag_;
    907   /// First free/super-basic variable (-1 if none)
    908   int firstFree_;
    909   /** For advanced use.  When doing iterative solves things can get
    910       nasty so on values pass if incoming solution has largest
    911       infeasibility < incomingInfeasibility throw out variables
    912       from basis until largest infeasibility < allowedInfeasibility.
    913       if allowedInfeasibility>= incomingInfeasibility this is
    914       always possible altough you may end up with an all slack basis.
    915 
    916       Defaults are 1.0,10.0
    917   */
    918   float incomingInfeasibility_;
    919   float allowedInfeasibility_;
    920   /// For dealing with all issues of cycling etc
    921   ClpInteriorProgress * progress_;
    922232  //@}
    923233};
     
    936246
    937247
    938 /// For saving extra information to see if looping.
    939 class ClpInteriorProgress {
    940 
    941 public:
    942 
    943 
    944   /**@name Constructors and destructor and copy */
    945   //@{
    946   /// Default constructor
    947     ClpInteriorProgress (  );
    948 
    949   /// Constructor from model
    950     ClpInteriorProgress ( ClpInterior * model );
    951 
    952   /// Copy constructor.
    953   ClpInteriorProgress(const ClpInteriorProgress &);
    954 
    955   /// Assignment operator. This copies the data
    956     ClpInteriorProgress & operator=(const ClpInteriorProgress & rhs);
    957   /// Destructor
    958    ~ClpInteriorProgress (  );
    959   //@}
    960 
    961   /**@name Check progress */
    962   //@{
    963   /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
    964       >=0 if give up and use as problem status
    965   */
    966     int looping (  );
    967   /// Start check at beginning of whileIterating
    968   void startCheck();
    969   /// Returns cycle length in whileIterating
    970   int cycle(int in, int out,int wayIn,int wayOut);
    971 
    972   /// Returns previous objective (if -1) - current if (0)
    973   double lastObjective(int back=1) const;
    974   /// Modify objective e.g. if dual infeasible in dual
    975   void modifyObjective(double value);
    976   /// Returns previous iteration number (if -1) - current if (0)
    977   int lastIterationNumber(int back=1) const;
    978 
    979   //@}
    980   /**@name Data  */
    981 #define CLP_PROGRESS 5
    982   //@{
    983   /// Objective values
    984   double objective_[CLP_PROGRESS];
    985   /// Sum of infeasibilities for algorithm
    986   double infeasibility_[CLP_PROGRESS];
    987   /// Pointer back to model so we can get information
    988   ClpInterior * model_;
    989   /// Number of infeasibilities
    990   int numberInfeasibilities_[CLP_PROGRESS];
    991   /// Iteration number at which occurred
    992   int iterationNumber_[CLP_PROGRESS];
    993   /// Number of times checked (so won't stop too early)
    994   int numberTimes_;
    995   /// Number of times it looked like loop
    996   int numberBadTimes_;
    997 #define CLP_CYCLE 12
    998   /// For cycle checking
    999   int in_[CLP_CYCLE];
    1000   int out_[CLP_CYCLE];
    1001   char way_[CLP_CYCLE];
    1002   //@}
    1003 };
    1004248#endif
  • branches/pre/include/ClpMatrixBase.hpp

    r210 r213  
    99class CoinIndexedVector;
    1010class ClpSimplex;
     11class ClpModel;
    1112
    1213/** Abstract base class for Clp Matrices
     
    9596      small elements.
    9697  */
    97   virtual bool allElementsInRange(ClpSimplex * model,
     98  virtual bool allElementsInRange(ClpModel * model,
    9899                                  double smallest, double largest)
    99100  { return true;};
  • branches/pre/include/ClpPackedMatrix.hpp

    r210 r213  
    9797      small elements.
    9898  */
    99   virtual bool allElementsInRange(ClpSimplex * model,
     99  virtual bool allElementsInRange(ClpModel * model,
    100100                                  double smallest, double largest);
    101101
Note: See TracChangeset for help on using the changeset viewer.