Changeset 17 for branches


Ignore:
Timestamp:
Sep 3, 2002 1:43:51 PM (17 years ago)
Author:
forrest
Message:

Synchronizing

Location:
branches/devel-1
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • branches/devel-1/ClpDualRowSteepest.cpp

    r2 r17  
    3030    savedWeights_(NULL)
    3131{
    32   type_=2;
     32  type_=2+64*mode;
    3333}
    3434
     
    245245      }
    246246    }
     247#ifdef CLP_DEBUG
    247248    assert(work3[pivotRow]&&work[pivotRow]);
     249#endif
    248250    alternateWeights_->setNumElements(nSave);
    249251    if (norm < TRY_NORM)
     
    275277  const int * pivotVariable = model_->pivotVariable();
    276278  double * infeas = infeasible_->denseVector();
     279  int pivotRow = model_->pivotRow();
     280  double * solution = model_->solutionRegion();
    277281  for (i=0;i<number;i++) {
    278282    int iRow=which[i];
    279283    int iPivot=pivotVariable[iRow];
    280     double & value = model_->solutionAddress(iPivot);
     284    double value = solution[iPivot];
    281285    double cost = model_->cost(iPivot);
    282286    double change = primalRatio*work[iRow];
    283287    value -= change;
     288    changeObj -= change*cost;
     289    solution[iPivot] = value;
    284290    double lower = model_->lower(iPivot);
    285291    double upper = model_->upper(iPivot);
     292    // But if pivot row then use value of incoming
     293    if (iRow==pivotRow) {
     294      iPivot = model_->sequenceIn();
     295      // make last resort choice
     296      lower = 1.0e-6*model_->lower(iPivot);
     297      upper = 1.0e-6*model_->upper(iPivot);
     298      value = 1.0e-6*model_->valueIncomingDual();
     299    }
    286300    if (value>upper+tolerance) {
    287301      // store square in list
     
    301315        infeas[iRow] = 1.0e-100;
    302316    }
    303     changeObj -= change*cost;
    304317    work[iRow]=0.0;
    305318  }
     
    332345    }
    333346    state_=1;
    334   } else if (mode==2||mode==4) {
     347  } else if (mode==2||mode==4||mode==5) {
    335348    // restore
    336     if (!weights_||state_==-1) {
     349    if (!weights_||state_==-1||mode==5) {
    337350      // initialize weights
    338351      delete [] weights_;
     
    343356      alternateWeights_->reserve(numberRows+
    344357                                 model_->factorization()->maximumPivots());
    345       if (!mode_) {
     358      if (!mode_||mode==5) {
    346359        // initialize to 1.0 (can we do better?)
    347360        for (i=0;i<numberRows;i++) {
  • branches/devel-1/ClpModel.cpp

    r9 r17  
    9393  lengthNames_(0),
    9494  rowNames_(),
    95   columnNames_()
     95  columnNames_(),
     96  integerType_(NULL)
    9697{
    9798  intParam_[OsiMaxNumIteration] = 9999999;
     
    148149  delete [] ray_;
    149150  ray_ = NULL;
     151  delete [] integerType_;
     152  integerType_ = NULL;
    150153}
    151154//#############################################################################
     
    328331    rowNames_ = rhs.rowNames_;
    329332    columnNames_ = rhs.columnNames_;
     333    if (rhs.integerType_) {
     334      integerType_ = new char[numberColumns_];
     335      memcpy(integerType_,rhs.integerType_,numberColumns_*sizeof(char));
     336    } else {
     337      integerType_ = NULL;
     338    }
    330339    if (rhs.rowActivity_) {
    331340      rowActivity_=new double[numberRows_];
     
    384393    rowNames_ = std::vector<std::string> ();
    385394    columnNames_ = std::vector<std::string> ();
     395    integerType_ = NULL;
    386396  }
    387397}
     
    400410  numberRows_ = rhs.numberRows_;
    401411  numberColumns_ = rhs.numberColumns_;
     412  delete [] rhs.ray_;
     413  rhs.ray_=NULL;
    402414  gutsOfCopy(rhs,false);
    403415}
     
    421433  matrix_ = NULL;
    422434  rowCopy_ = NULL;
     435  delete [] otherModel.ray_;
     436  otherModel.ray_ = ray_;
    423437  ray_ = NULL;
    424438}
     
    528542    newSize = size-numberDeleted;
    529543    double * newArray = new double[newSize];
     544    int put=0;
     545    for (i=0;i<size;i++) {
     546      if (!deleted[i]) {
     547        newArray[put++]=array[i];
     548      }
     549    }
     550    delete [] array;
     551    array = newArray;
     552    delete [] deleted;
     553  }
     554  return array;
     555}
     556char * deleteChar(char * array , int size,
     557                      int number, const int * which,int & newSize)
     558{
     559  if (array) {
     560    int i ;
     561    char * deleted = new char[size];
     562    int numberDeleted=0;
     563    memset(deleted,0,size*sizeof(char));
     564    for (i=0;i<number;i++) {
     565      int j = which[i];
     566      if (j>=0&&j<size&&!deleted[j]) {
     567        numberDeleted++;
     568        deleted[j]=1;
     569      }
     570    }
     571    newSize = size-numberDeleted;
     572    char * newArray = new char[newSize];
    530573    int put=0;
    531574    for (i=0;i<size;i++) {
     
    587630    delete [] which;
    588631  }
     632  if (integerType_) {
     633    char * temp = new char [newNumberColumns];
     634    memset(temp,0,newNumberColumns*sizeof(char));
     635    memcpy(temp,integerType_,
     636           min(newNumberColumns,numberColumns_)*sizeof(char));
     637    delete [] integerType_;
     638    integerType_ = temp;
     639  }
    589640  numberColumns_ = newNumberColumns;
    590641  // for now gets rid of names
     
    644695  rowNames_ = std::vector<std::string> ();
    645696  columnNames_ = std::vector<std::string> ();
     697  integerType_ = deleteChar(integerType_,numberColumns_,
     698                            number, which, newSize);
    646699}
    647700// Infeasibility/unbounded ray (NULL returned if none/wrong)
     
    712765                m.getObjCoefficients(),
    713766                m.getRowLower(),m.getRowUpper());
     767    if (m.integerColumns()) {
     768      integerType_ = new char[numberColumns_];
     769      memcpy(integerType_,m.integerColumns(),numberColumns_*sizeof(char));
     770    } else {
     771      integerType_ = NULL;
     772    }
    714773    // set problem name
    715774    setStrParam(OsiProbName,m.getProblemName());
     
    749808  return status;
    750809}
     810bool ClpModel::isPrimalObjectiveLimitReached() const
     811{
     812  double limit = 0.0;
     813  getDblParam(OsiPrimalObjectiveLimit, limit);
     814  if (limit > 1e30) {
     815    // was not ever set
     816    return false;
     817  }
     818   
     819  const double obj = objectiveValue();
     820  const int maxmin = optimizationDirection();
     821
     822  if (problemStatus_ == 0) // optimal
     823    return maxmin > 0 ? (obj < limit) /*minim*/ : (obj > limit) /*maxim*/;
     824  else if (problemStatus_==2)
     825    return true;
     826  else
     827    return false;
     828}
     829
     830bool ClpModel::isDualObjectiveLimitReached() const
     831{
     832
     833  double limit = 0.0;
     834  getDblParam(OsiDualObjectiveLimit, limit);
     835  if (limit > 1e30) {
     836    // was not ever set
     837    return false;
     838  }
     839   
     840  const double obj = objectiveValue();
     841  const int maxmin = optimizationDirection();
     842
     843  if (problemStatus_ == 0) // optimal
     844    return maxmin > 0 ? (obj > limit) /*minim*/ : (obj < limit) /*maxim*/;
     845  else if (problemStatus_==1)
     846    return true;
     847  else
     848    return false;
     849
     850}
  • branches/devel-1/ClpPackedMatrix.cpp

    r8 r17  
    3636
    3737  matrix_ = new OsiPackedMatrix(*(rhs.matrix_));
     38 
     39}
     40
     41//-------------------------------------------------------------------
     42// assign matrix (for space reasons)
     43//-------------------------------------------------------------------
     44ClpPackedMatrix::ClpPackedMatrix (OsiPackedMatrix * rhs)
     45: ClpMatrixBase()
     46
     47  matrix_ = rhs;
     48  setType(1);
    3849 
    3950}
  • branches/devel-1/ClpPrimalColumnSteepest.cpp

    r15 r17  
    3535    devex_(0.0)
    3636{
    37   type_=2;
     37  type_=2+64*mode;
    3838}
    3939
  • branches/devel-1/ClpSimplex.cpp

    r15 r17  
    1414#include "ClpSimplex.hpp"
    1515#include "ClpFactorization.hpp"
    16 #include "OsiPackedMatrix.hpp"
     16#include "ClpPackedMatrix.hpp"
    1717#include "OsiIndexedVector.hpp"
    1818#include "OsiWarmStartBasis.hpp"
     19#include "ClpDualRowDantzig.hpp"
    1920#include "ClpDualRowSteepest.hpp"
    2021#include "ClpPrimalColumnDantzig.hpp"
     
    126127  nonLinearCost_(NULL),
    127128  specialOptions_(0),
    128   lastBadIteration_(-999999)
     129  lastBadIteration_(-999999),
     130  numberFake_(0)
    129131
    130132{
     
    845847        numberSlacks++;
    846848    }
    847     status= numberSlacks-totalSlacks;
     849    status= max(numberSlacks-totalSlacks,0);
    848850  }
    849851
     
    10721074  nonLinearCost_(NULL),
    10731075  specialOptions_(0),
    1074   lastBadIteration_(-999999)
     1076  lastBadIteration_(-999999),
     1077  numberFake_(0)
    10751078{
    10761079  int i;
     
    11571160  nonLinearCost_(NULL),
    11581161  specialOptions_(0),
    1159   lastBadIteration_(-999999)
     1162  lastBadIteration_(-999999),
     1163  numberFake_(0)
    11601164{
    11611165  int i;
     
    12851289  specialOptions_ = rhs.specialOptions_;
    12861290  lastBadIteration_ = rhs.lastBadIteration_;
     1291  numberFake_ = rhs.numberFake_;
    12871292  if (rhs.nonLinearCost_!=NULL)
    12881293    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
     
    21612166  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
    21622167}
     2168/* For strong branching.  On input lower and upper are new bounds
     2169   while on output they are objective function values (>1.0e50 infeasible).
     2170   Return code is 0 if nothing interesting, -1 if infeasible both
     2171   ways and +1 if infeasible one way (check values to see which one(s))
     2172*/
     2173int ClpSimplex::strongBranching(int numberVariables,const int * variables,
     2174                                double * newLower, double * newUpper,
     2175                                bool stopOnFirstInfeasible,
     2176                                bool alwaysFinish)
     2177{
     2178  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
     2179                                                    newLower,  newUpper,
     2180                                                    stopOnFirstInfeasible,
     2181                                                    alwaysFinish);
     2182}
     2183/* Borrow model.  This is so we dont have to copy large amounts
     2184   of data around.  It assumes a derived class wants to overwrite
     2185   an empty model with a real one - while it does an algorithm.
     2186   This is same as ClpModel one, but sets scaling on etc. */
     2187void
     2188ClpSimplex::borrowModel(ClpModel & otherModel)
     2189{
     2190  ClpModel::borrowModel(otherModel);
     2191  scaling();
     2192  ClpDualRowSteepest steep1;
     2193  setDualRowPivotAlgorithm(steep1);
     2194  ClpPrimalColumnSteepest steep2;
     2195  setPrimalColumnPivotAlgorithm(steep2);
     2196}
     2197typedef struct {
     2198  double optimizationDirection;
     2199  double dblParam[OsiLastDblParam];
     2200  double objectiveValue;
     2201  double dualBound;
     2202  double dualTolerance;
     2203  double primalTolerance;
     2204  double sumDualInfeasibilities;
     2205  double sumPrimalInfeasibilities;
     2206  double infeasibilityCost;
     2207  int numberRows;
     2208  int numberColumns;
     2209  int intParam[OsiLastIntParam];
     2210  int numberIterations;
     2211  int problemStatus;
     2212  int maximumIterations;
     2213  int lengthNames;
     2214  int numberDualInfeasibilities;
     2215  int numberDualInfeasibilitiesWithoutFree;
     2216  int numberPrimalInfeasibilities;
     2217  int numberRefinements;
     2218  int scalingFlag;
     2219  int algorithm;
     2220  unsigned int specialOptions;
     2221  int dualPivotChoice;
     2222  int primalPivotChoice;
     2223  int matrixStorageChoice;
     2224} Clp_scalars;
     2225
     2226int outDoubleArray(double * array, int length, FILE * fp)
     2227{
     2228  int numberWritten;
     2229  if (array&&length) {
     2230    numberWritten = fwrite(&length,sizeof(int),1,fp);
     2231    if (numberWritten!=1)
     2232      return 1;
     2233    numberWritten = fwrite(array,sizeof(double),length,fp);
     2234    if (numberWritten!=length)
     2235      return 1;
     2236  } else {
     2237    length = 0;
     2238    numberWritten = fwrite(&length,sizeof(int),1,fp);
     2239    if (numberWritten!=1)
     2240      return 1;
     2241  }
     2242  return 0;
     2243}
     2244// Save model to file, returns 0 if success
     2245int
     2246ClpSimplex::saveModel(const char * fileName)
     2247{
     2248  FILE * fp = fopen(fileName,"wb");
     2249  if (fp) {
     2250    Clp_scalars scalars;
     2251    int i, numberWritten;
     2252    // Fill in scalars
     2253    scalars.optimizationDirection = optimizationDirection_;
     2254    memcpy(scalars.dblParam, dblParam_,OsiLastDblParam * sizeof(double));
     2255    scalars.objectiveValue = objectiveValue_;
     2256    scalars.dualBound = dualBound_;
     2257    scalars.dualTolerance = dualTolerance_;
     2258    scalars.primalTolerance = primalTolerance_;
     2259    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
     2260    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
     2261    scalars.infeasibilityCost = infeasibilityCost_;
     2262    scalars.numberRows = numberRows_;
     2263    scalars.numberColumns = numberColumns_;
     2264    memcpy(scalars.intParam, intParam_,OsiLastIntParam * sizeof(double));
     2265    scalars.numberIterations = numberIterations_;
     2266    scalars.problemStatus = problemStatus_;
     2267    scalars.maximumIterations = maximumIterations_;
     2268    scalars.lengthNames = lengthNames_;
     2269    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
     2270    scalars.numberDualInfeasibilitiesWithoutFree
     2271      = numberDualInfeasibilitiesWithoutFree_;
     2272    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
     2273    scalars.numberRefinements = numberRefinements_;
     2274    scalars.scalingFlag = scalingFlag_;
     2275    scalars.algorithm = algorithm_;
     2276    scalars.specialOptions = specialOptions_;
     2277    scalars.dualPivotChoice = dualRowPivot_->type();
     2278    scalars.primalPivotChoice = primalColumnPivot_->type();
     2279    scalars.matrixStorageChoice = matrix_->type();
     2280
     2281    // put out scalars
     2282    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
     2283    if (numberWritten!=1)
     2284      return 1;
     2285    // strings
     2286    int length;
     2287    for (i=0;i<OsiLastStrParam;i++) {
     2288      length = strParam_[i].size();
     2289      numberWritten = fwrite(&length,sizeof(int),1,fp);
     2290      if (numberWritten!=1)
     2291        return 1;
     2292      numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
     2293      if (numberWritten!=1)
     2294        return 1;
     2295    }
     2296    // arrays - in no particular order
     2297    if (outDoubleArray(rowActivity_,numberRows_,fp))
     2298        return 1;
     2299    if (outDoubleArray(columnActivity_,numberColumns_,fp))
     2300        return 1;
     2301    if (outDoubleArray(dual_,numberRows_,fp))
     2302        return 1;
     2303    if (outDoubleArray(reducedCost_,numberColumns_,fp))
     2304        return 1;
     2305    if (outDoubleArray(rowLower_,numberRows_,fp))
     2306        return 1;
     2307    if (outDoubleArray(rowUpper_,numberRows_,fp))
     2308        return 1;
     2309    if (outDoubleArray(objective_,numberColumns_,fp))
     2310        return 1;
     2311    if (outDoubleArray(rowObjective_,numberRows_,fp))
     2312        return 1;
     2313    if (outDoubleArray(columnLower_,numberColumns_,fp))
     2314        return 1;
     2315    if (outDoubleArray(columnUpper_,numberColumns_,fp))
     2316        return 1;
     2317    if (ray_) {
     2318      if (problemStatus_==1)
     2319        if (outDoubleArray(ray_,numberRows_,fp))
     2320          return 1;
     2321      else if (problemStatus_==2)
     2322        if (outDoubleArray(ray_,numberColumns_,fp))
     2323          return 1;
     2324      else
     2325        if (outDoubleArray(NULL,0,fp))
     2326          return 1;
     2327    } else {
     2328      if (outDoubleArray(NULL,0,fp))
     2329        return 1;
     2330    }
     2331    if (status_&&(numberRows_+numberColumns_)>0) {
     2332      length = numberRows_+numberColumns_;
     2333      numberWritten = fwrite(&length,sizeof(int),1,fp);
     2334      if (numberWritten!=1)
     2335        return 1;
     2336      numberWritten = fwrite(status_,sizeof(char),length, fp);
     2337      if (numberWritten!=length)
     2338        return 1;
     2339    } else {
     2340      length = 0;
     2341      numberWritten = fwrite(&length,sizeof(int),1,fp);
     2342      if (numberWritten!=1)
     2343        return 1;
     2344    }
     2345    if (lengthNames_) {
     2346      assert (numberRows_ == (int) rowNames_.size());
     2347      for (i=0;i<numberRows_;i++) {
     2348        length = rowNames_[i].size();
     2349        numberWritten = fwrite(&length,sizeof(int),1,fp);
     2350        if (numberWritten!=1)
     2351          return 1;
     2352        numberWritten = fwrite(rowNames_[i].c_str(),length,1,fp);
     2353        if (numberWritten!=1)
     2354          return 1;
     2355      }
     2356      assert (numberColumns_ == (int) columnNames_.size());
     2357      for (i=0;i<numberColumns_;i++) {
     2358        length = columnNames_[i].size();
     2359        numberWritten = fwrite(&length,sizeof(int),1,fp);
     2360        if (numberWritten!=1)
     2361          return 1;
     2362        numberWritten = fwrite(columnNames_[i].c_str(),length,1,fp);
     2363        if (numberWritten!=1)
     2364          return 1;
     2365      }
     2366    }
     2367    // just standard type at present
     2368    assert (matrix_->type()==1);
     2369    assert (matrix_->getNumCols() == numberColumns_);
     2370    assert (matrix_->getNumRows() == numberRows_);
     2371    // we are going to save with gaps
     2372    length = matrix_->getVectorStarts()[numberColumns_-1]
     2373      + matrix_->getVectorLengths()[numberColumns_-1];
     2374    numberWritten = fwrite(&length,sizeof(int),1,fp);
     2375    if (numberWritten!=1)
     2376      return 1;
     2377    numberWritten = fwrite(matrix_->getElements(),
     2378                           sizeof(double),length,fp);
     2379    if (numberWritten!=length)
     2380      return 1;
     2381    numberWritten = fwrite(matrix_->getIndices(),
     2382                           sizeof(int),length,fp);
     2383    if (numberWritten!=length)
     2384      return 1;
     2385    numberWritten = fwrite(matrix_->getVectorStarts(),
     2386                           sizeof(int),numberColumns_,fp);
     2387    if (numberWritten!=numberColumns_)
     2388      return 1;
     2389    numberWritten = fwrite(matrix_->getVectorLengths(),
     2390                           sizeof(int),numberColumns_,fp);
     2391    if (numberWritten!=numberColumns_)
     2392      return 1;
     2393    // finished
     2394    fclose(fp);
     2395    return 0;
     2396  } else {
     2397    return -1;
     2398  }
     2399}
     2400
     2401int inDoubleArray(double * &array, int length, FILE * fp)
     2402{
     2403  int numberRead;
     2404  int length2;
     2405  numberRead = fread(&length2,sizeof(int),1,fp);
     2406  if (numberRead!=1)
     2407    return 1;
     2408  if (length2) {
     2409    // lengths must match
     2410    if (length!=length2)
     2411      return 2;
     2412    array = new double[length];
     2413    numberRead = fread(array,sizeof(double),length,fp);
     2414    if (numberRead!=length)
     2415      return 1;
     2416  }
     2417  return 0;
     2418}
     2419/* Restore model from file, returns 0 if success,
     2420   deletes current model */
     2421int
     2422ClpSimplex::restoreModel(const char * fileName)
     2423{
     2424  FILE * fp = fopen(fileName,"rb");
     2425  if (fp) {
     2426    // Get rid of current model
     2427    ClpModel::gutsOfDelete();
     2428    gutsOfDelete(0);
     2429    int i;
     2430    for (i=0;i<6;i++) {
     2431      rowArray_[i]=NULL;
     2432      columnArray_[i]=NULL;
     2433    }
     2434    // get an empty factorization so we can set tolerances etc
     2435    factorization_ = new ClpFactorization();
     2436    // Say sparse
     2437    factorization_->sparseThreshold(1);
     2438    Clp_scalars scalars;
     2439    int numberRead;
     2440
     2441    // get scalars
     2442    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
     2443    if (numberRead!=1)
     2444      return 1;
     2445    // Fill in scalars
     2446    optimizationDirection_ = scalars.optimizationDirection;
     2447    memcpy(dblParam_, scalars.dblParam, OsiLastDblParam * sizeof(double));
     2448    objectiveValue_ = scalars.objectiveValue;
     2449    dualBound_ = scalars.dualBound;
     2450    dualTolerance_ = scalars.dualTolerance;
     2451    primalTolerance_ = scalars.primalTolerance;
     2452    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
     2453    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
     2454    infeasibilityCost_ = scalars.infeasibilityCost;
     2455    numberRows_ = scalars.numberRows;
     2456    numberColumns_ = scalars.numberColumns;
     2457    memcpy(intParam_, scalars.intParam,OsiLastIntParam * sizeof(double));
     2458    numberIterations_ = scalars.numberIterations;
     2459    problemStatus_ = scalars.problemStatus;
     2460    maximumIterations_ = scalars.maximumIterations;
     2461    lengthNames_ = scalars.lengthNames;
     2462    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
     2463    numberDualInfeasibilitiesWithoutFree_
     2464      = scalars.numberDualInfeasibilitiesWithoutFree;
     2465    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
     2466    numberRefinements_ = scalars.numberRefinements;
     2467    scalingFlag_ = scalars.scalingFlag;
     2468    algorithm_ = scalars.algorithm;
     2469    specialOptions_ = scalars.specialOptions;
     2470    // strings
     2471    int length;
     2472    for (i=0;i<OsiLastStrParam;i++) {
     2473      numberRead = fread(&length,sizeof(int),1,fp);
     2474      if (numberRead!=1)
     2475        return 1;
     2476      char * array = new char[length+1];
     2477      numberRead = fread(array,length,1,fp);
     2478      if (numberRead!=1)
     2479        return 1;
     2480      array[length]='\0';
     2481      strParam_[i]=array;
     2482      delete [] array;
     2483    }
     2484    // arrays - in no particular order
     2485    if (inDoubleArray(rowActivity_,numberRows_,fp))
     2486        return 1;
     2487    if (inDoubleArray(columnActivity_,numberColumns_,fp))
     2488        return 1;
     2489    if (inDoubleArray(dual_,numberRows_,fp))
     2490        return 1;
     2491    if (inDoubleArray(reducedCost_,numberColumns_,fp))
     2492        return 1;
     2493    if (inDoubleArray(rowLower_,numberRows_,fp))
     2494        return 1;
     2495    if (inDoubleArray(rowUpper_,numberRows_,fp))
     2496        return 1;
     2497    if (inDoubleArray(objective_,numberColumns_,fp))
     2498        return 1;
     2499    if (inDoubleArray(rowObjective_,numberRows_,fp))
     2500        return 1;
     2501    if (inDoubleArray(columnLower_,numberColumns_,fp))
     2502        return 1;
     2503    if (inDoubleArray(columnUpper_,numberColumns_,fp))
     2504        return 1;
     2505    if (problemStatus_==1) {
     2506      if (inDoubleArray(ray_,numberRows_,fp))
     2507        return 1;
     2508    } else if (problemStatus_==2) {
     2509      if (inDoubleArray(ray_,numberColumns_,fp))
     2510        return 1;
     2511    } else {
     2512      // ray should be null
     2513      numberRead = fread(&length,sizeof(int),1,fp);
     2514      if (numberRead!=1)
     2515        return 1;
     2516      if (length)
     2517        return 2;
     2518    }
     2519    delete [] status_;
     2520    status_=NULL;
     2521    // status region
     2522    numberRead = fread(&length,sizeof(int),1,fp);
     2523    if (numberRead!=1)
     2524        return 1;
     2525    if (length) {
     2526      if (length!=numberRows_+numberColumns_)
     2527        return 1;
     2528      status_ = new char unsigned[length];
     2529      numberRead = fread(status_,sizeof(char),length, fp);
     2530      if (numberRead!=length)
     2531        return 1;
     2532    }
     2533    if (lengthNames_) {
     2534      char * array = new char[lengthNames_+1];
     2535      rowNames_.resize(0);
     2536      for (i=0;i<numberRows_;i++) {
     2537        numberRead = fread(&length,sizeof(int),1,fp);
     2538        if (numberRead!=1)
     2539          return 1;
     2540        numberRead = fread(array,length,1,fp);
     2541        if (numberRead!=1)
     2542          return 1;
     2543        rowNames_[i]=array;
     2544      }
     2545      columnNames_.resize(0);
     2546      for (i=0;i<numberColumns_;i++) {
     2547        numberRead = fread(&length,sizeof(int),1,fp);
     2548        if (numberRead!=1)
     2549          return 1;
     2550        numberRead = fread(array,length,1,fp);
     2551        if (numberRead!=1)
     2552          return 1;
     2553        columnNames_[i]=array;
     2554      }
     2555      delete [] array;
     2556    }
     2557    // Pivot choices
     2558    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
     2559    delete dualRowPivot_;
     2560    switch ((scalars.dualPivotChoice&63)) {
     2561    case 1:
     2562      // Dantzig
     2563      dualRowPivot_ = new ClpDualRowDantzig();
     2564      break;
     2565    case 2:
     2566      // Steepest - use mode
     2567      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
     2568      break;
     2569    default:
     2570      abort();
     2571    }
     2572    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
     2573    delete primalColumnPivot_;
     2574    switch ((scalars.primalPivotChoice&63)) {
     2575    case 1:
     2576      // Dantzig
     2577      primalColumnPivot_ = new ClpPrimalColumnDantzig();
     2578      break;
     2579    case 2:
     2580      // Steepest - use mode
     2581      primalColumnPivot_
     2582        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
     2583      break;
     2584    default:
     2585      abort();
     2586    }
     2587    assert(scalars.matrixStorageChoice==1);
     2588    delete matrix_;
     2589    // get arrays
     2590    numberRead = fread(&length,sizeof(int),1,fp);
     2591    if (numberRead!=1)
     2592      return 1;
     2593    double * elements = new double[length];
     2594    int * indices = new int[length];
     2595    int * starts = new int[numberColumns_];
     2596    int * lengths = new int[numberColumns_];
     2597    numberRead = fread(elements, sizeof(double),length,fp);
     2598    if (numberRead!=length)
     2599      return 1;
     2600    numberRead = fread(indices, sizeof(int),length,fp);
     2601    if (numberRead!=length)
     2602      return 1;
     2603    numberRead = fread(starts, sizeof(int),numberColumns_,fp);
     2604    if (numberRead!=numberColumns_)
     2605      return 1;
     2606    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
     2607    if (numberRead!=numberColumns_)
     2608      return 1;
     2609    // assign matrix
     2610    OsiPackedMatrix * matrix = new OsiPackedMatrix();
     2611    matrix->assignMatrix(true, numberRows_, numberColumns_,
     2612                         length, elements, indices, starts, lengths);
     2613    // and transfer to Clp
     2614    matrix_ = new ClpPackedMatrix(matrix);
     2615    // finished
     2616    fclose(fp);
     2617    return 0;
     2618  } else {
     2619    return -1;
     2620  }
     2621  return 0;
     2622}
     2623// value of incoming variable (in Dual)
     2624double
     2625ClpSimplex::valueIncomingDual() const
     2626{
     2627  // Need value of incoming for list of infeasibilities as may be infeasible
     2628  double valueIncoming = (dualOut_/alpha_)*directionOut_;
     2629  if (directionIn_==-1)
     2630    valueIncoming = upperIn_-valueIncoming;
     2631  else
     2632    valueIncoming = lowerIn_-valueIncoming;
     2633  return valueIncoming;
     2634}
  • branches/devel-1/ClpSimplexDual.cpp

    r15 r17  
    278278  // for this we need clean basis so it is after factorize
    279279  gutsOfSolution(rowActivityWork_,columnActivityWork_);
     280
     281  numberFake_ =0; // Number of variables at fake bounds
    280282  changeBounds(true,NULL,objectiveChange);
    281283
     
    334336  }
    335337
     338  assert(!numberFake_||problemStatus_); // all bounds should be okay
    336339  // at present we are leaving factorization around
    337340  // maybe we should empty it
     
    346349  return problemStatus_;
    347350}
    348 void
     351/* Reasons to come out:
     352   -1 iterations etc
     353   -2 inaccuracy
     354   -3 slight inaccuracy (and done iterations)
     355   +0 looks optimal (might be unbounded - but we will investigate)
     356   +1 looks infeasible
     357   +3 max iterations
     358 */
     359int
    349360ClpSimplexDual::whileIterating()
    350361{
    351362  // status stays at -1 while iterating, >=0 finished, -2 to invert
    352363  // status -3 to go to top without an invert
     364  int returnCode = -1;
    353365  while (problemStatus_==-1) {
    354366#ifdef CLP_DEBUG
     
    470482            rowArray_[1]->clear();
    471483            columnArray_[0]->clear();
     484            returnCode=-2;
    472485            break;
    473486          } else {
     
    525538            dualRowPivot_->unrollWeights();
    526539            problemStatus_=-2; // factorize now
     540            returnCode=-2;
    527541            break;
    528542          }
     
    543557            dualRowPivot_->unrollWeights();
    544558            problemStatus_=-2; // factorize now
     559            returnCode=-2;
    545560            break;
    546561          }
     
    552567        if (updateStatus==1) {
    553568          // slight error
    554           if (factorization_->pivots()>5)
     569          if (factorization_->pivots()>5) {
    555570            problemStatus_=-2; // factorize now
     571            returnCode=-3;
     572          }
    556573        } else if (updateStatus==2) {
    557574          // major error
     
    560577          if (factorization_->pivots()) {
    561578            problemStatus_=-2; // factorize now
     579            returnCode=-2;
    562580            break;
    563581          } else {
     
    613631          // maximum iterations or equivalent
    614632          problemStatus_= 3;
     633          returnCode=3;
    615634          break;
    616635        }
     
    630649        rowArray_[0]->clear();
    631650        columnArray_[0]->clear();
     651        returnCode=1;
    632652        break;
    633653      }
     
    660680        }
    661681      }
     682      returnCode=0;
    662683      break;
    663684    }
    664685  }
     686  return returnCode;
    665687}
    666688/* The duals are updated by the given arrays.
     
    838860    // do actual flips
    839861    flipBounds(rowArray,columnArray,theta);
     862    numberFake_ = numberAtFake;
    840863  }
    841864  objectiveChange += changeObj;
     
    17061729    // save dual weights
    17071730    dualRowPivot_->saveWeights(this,1);
    1708     // is factorization okay?
    1709     if (internalFactorize(1)) {
    1710       // no - restore previous basis
    1711       assert (type==1);
    1712       changeMade_++; // say something changed
    1713       memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    1714       memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
    1715              numberRows_*sizeof(double));
    1716       memcpy(columnActivityWork_,savedSolution_ ,
    1717              numberColumns_*sizeof(double));
    1718       // get correct bounds on all variables
    1719       double dummyChangeCost=0.0;
    1720       changeBounds(true,rowArray_[2],dummyChangeCost);
    1721       // throw away change
    1722       rowArray_[2]->clear();
    1723       forceFactorization_=1; // a bit drastic but ..
    1724       type = 2;
    1725       assert (internalFactorize(1)==0);
     1731    if (type) {
     1732      // is factorization okay?
     1733      if (internalFactorize(1)) {
     1734        // no - restore previous basis
     1735        assert (type==1);
     1736        changeMade_++; // say something changed
     1737        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
     1738        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
     1739               numberRows_*sizeof(double));
     1740        memcpy(columnActivityWork_,savedSolution_ ,
     1741               numberColumns_*sizeof(double));
     1742        // get correct bounds on all variables
     1743        double dummyChangeCost=0.0;
     1744        changeBounds(true,rowArray_[2],dummyChangeCost);
     1745        // throw away change
     1746        rowArray_[2]->clear();
     1747        forceFactorization_=1; // a bit drastic but ..
     1748        type = 2;
     1749        assert (internalFactorize(1)==0);
     1750      }
    17261751    }
    17271752    problemStatus_=-3;
     
    19331958}
    19341959/* While updateDualsInDual sees what effect is of flip
    1935    this does actuall flipping.
     1960   this does actual flipping.
    19361961   If change >0.0 then value in array >0.0 => from lower to upper
    19371962*/
     
    21292154
    21302155}
    2131 
     2156/* For strong branching.  On input lower and upper are new bounds
     2157   while on output they are change in objective function values
     2158   (>1.0e50 infeasible).
     2159   Return code is 0 if nothing interesting, -1 if infeasible both
     2160   ways and +1 if infeasible one way (check values to see which one(s))
     2161*/
     2162int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
     2163                                    double * newLower, double * newUpper,
     2164                                    bool stopOnFirstInfeasible,
     2165                                    bool alwaysFinish)
     2166{
     2167  int i;
     2168  int returnCode=0;
     2169  double saveObjective = objectiveValue_;
     2170#if 1
     2171  algorithm_ = -1;
     2172
     2173  //scaling(false);
     2174
     2175  // put in standard form (and make row copy)
     2176  // create modifiable copies of model rim and do optional scaling
     2177  createRim(7+8+16,true);
     2178
     2179  // change newLower and newUpper if scaled
     2180
     2181  // Do initial factorization
     2182  // and set certain stuff
     2183  // We can either set increasing rows so ...IsBasic gives pivot row
     2184  // or we can just increment iBasic one by one
     2185  // for now let ...iBasic give pivot row
     2186  factorization_->increasingRows(2);
     2187  // row activities have negative sign
     2188  factorization_->slackValue(-1.0);
     2189  factorization_->zeroTolerance(1.0e-13);
     2190  // save if sparse factorization wanted
     2191  int saveSparse = factorization_->sparseThreshold();
     2192
     2193  int factorizationStatus = internalFactorize(0);
     2194  if (factorizationStatus<0)
     2195    return 1; // some error
     2196  else if (factorizationStatus)
     2197    handler_->message(CLP_SINGULARITIES,messages_)
     2198      <<factorizationStatus
     2199      <<OsiMessageEol;
     2200  if (saveSparse) {
     2201    // use default at present
     2202    factorization_->sparseThreshold(0);
     2203    factorization_->goSparse();
     2204  }
     2205 
     2206  // save stuff
     2207  ClpFactorization saveFactorization(*factorization_);
     2208  // save basis and solution
     2209  double * saveSolution = new double[numberRows_+numberColumns_];
     2210  memcpy(saveSolution,solution_,
     2211         (numberRows_+numberColumns_)*sizeof(double));
     2212  unsigned char * saveStatus =
     2213    new unsigned char [numberRows_+numberColumns_];
     2214  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
     2215  // save bounds as createRim makes clean copies
     2216  double * saveLower = new double[numberRows_+numberColumns_];
     2217  memcpy(saveLower,lower_,
     2218         (numberRows_+numberColumns_)*sizeof(double));
     2219  double * saveUpper = new double[numberRows_+numberColumns_];
     2220  memcpy(saveUpper,upper_,
     2221         (numberRows_+numberColumns_)*sizeof(double));
     2222  int * savePivot = new int [numberRows_];
     2223  memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int));
     2224  // need to save/restore weights.
     2225
     2226  for (i=0;i<numberVariables;i++) {
     2227    int iColumn = variables[i];
     2228    double objectiveChange;
     2229    double saveBound;
     2230   
     2231    // try down
     2232   
     2233    saveBound = columnUpper_[iColumn];
     2234    // external view - in case really getting optimal
     2235    columnUpper_[iColumn] = newUpper[i];
     2236    if (scalingFlag_<=0)
     2237      upper_[iColumn] = newUpper[i];
     2238    else
     2239      upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale
     2240    // Start of fast iterations
     2241    int status = fastDual(alwaysFinish);
     2242
     2243    // restore
     2244    memcpy(solution_,saveSolution,
     2245           (numberRows_+numberColumns_)*sizeof(double));
     2246    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
     2247    memcpy(lower_,saveLower,
     2248           (numberRows_+numberColumns_)*sizeof(double));
     2249    memcpy(upper_,saveUpper,
     2250           (numberRows_+numberColumns_)*sizeof(double));
     2251    columnUpper_[iColumn] = saveBound;
     2252    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
     2253    delete factorization_;
     2254    factorization_ = new ClpFactorization(saveFactorization);
     2255
     2256    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
     2257      objectiveChange = objectiveValue_-saveObjective;
     2258    } else {
     2259      objectiveChange = 1.0e100;
     2260    }
     2261    newUpper[i]=objectiveChange;
     2262#ifdef CLP_DEBUG
     2263    printf("down on %d costs %g\n",iColumn,objectiveChange);
     2264#endif
     2265         
     2266    // try up
     2267   
     2268    saveBound = columnLower_[iColumn];
     2269    // external view - in case really getting optimal
     2270    columnLower_[iColumn] = newLower[i];
     2271    if (scalingFlag_<=0)
     2272      lower_[iColumn] = newLower[i];
     2273    else
     2274      lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale
     2275    // Start of fast iterations
     2276    status = fastDual(alwaysFinish);
     2277
     2278    // restore
     2279    memcpy(solution_,saveSolution,
     2280           (numberRows_+numberColumns_)*sizeof(double));
     2281    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
     2282    memcpy(lower_,saveLower,
     2283           (numberRows_+numberColumns_)*sizeof(double));
     2284    memcpy(upper_,saveUpper,
     2285           (numberRows_+numberColumns_)*sizeof(double));
     2286    columnLower_[iColumn] = saveBound;
     2287    memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int));
     2288    delete factorization_;
     2289    factorization_ = new ClpFactorization(saveFactorization);
     2290
     2291    if (status||(problemStatus_==0&&!isDualObjectiveLimitReached())) {
     2292      objectiveChange = objectiveValue_-saveObjective;
     2293    } else {
     2294      objectiveChange = 1.0e100;
     2295    }
     2296    newLower[i]=objectiveChange;
     2297#ifdef CLP_DEBUG
     2298    printf("up on %d costs %g\n",iColumn,objectiveChange);
     2299#endif
     2300         
     2301    /* Possibilities are:
     2302       Both sides feasible - store
     2303       Neither side feasible - set objective high and exit
     2304       One side feasible - change bounds and resolve
     2305    */
     2306    if (newUpper[i]<1.0e100) {
     2307      if(newLower[i]<1.0e100) {
     2308        // feasible - no action
     2309      } else {
     2310        // up feasible, down infeasible
     2311        returnCode=1;
     2312        if (stopOnFirstInfeasible)
     2313          break;
     2314      }
     2315    } else {
     2316      if(newLower[i]<1.0e100) {
     2317        // down feasible, up infeasible
     2318        returnCode=1;
     2319        if (stopOnFirstInfeasible)
     2320          break;
     2321      } else {
     2322        // neither side feasible
     2323        returnCode=-1;
     2324        break;
     2325      }
     2326    }
     2327  }
     2328  delete [] saveSolution;
     2329  delete [] saveLower;
     2330  delete [] saveUpper;
     2331  delete [] saveStatus;
     2332  delete [] savePivot;
     2333
     2334  // at present we are leaving factorization around
     2335  // maybe we should empty it
     2336  deleteRim();
     2337  factorization_->sparseThreshold(saveSparse);
     2338#else
     2339  // save basis and solution
     2340  double * rowSolution = new double[numberRows_];
     2341  memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double));
     2342  double * columnSolution = new double[numberColumns_];
     2343  memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double));
     2344  unsigned char * saveStatus =
     2345    new unsigned char [numberRows_+numberColumns_];
     2346  if (!status_) {
     2347    // odd but anyway setup all slack basis
     2348    status_ = new unsigned char [numberColumns_+numberRows_];
     2349    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
     2350  }
     2351  memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
     2352  for (i=0;i<numberVariables;i++) {
     2353    int iColumn = variables[i];
     2354    double objectiveChange;
     2355   
     2356    // try down
     2357   
     2358    double saveUpper = columnUpper_[iColumn];
     2359    columnUpper_[iColumn] = newUpper[i];
     2360    dual();
     2361
     2362    // restore
     2363    columnUpper_[iColumn] = saveUpper;
     2364    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
     2365    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
     2366    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
     2367
     2368    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
     2369      objectiveChange = objectiveValue_-saveObjective;
     2370    } else {
     2371      objectiveChange = 1.0e100;
     2372    }
     2373    newUpper[i]=objectiveChange;
     2374#ifdef CLP_DEBUG
     2375    printf("down on %d costs %g\n",iColumn,objectiveChange);
     2376#endif
     2377         
     2378    // try up
     2379   
     2380    double saveLower = columnLower_[iColumn];
     2381    columnLower_[iColumn] = newLower[i];
     2382    dual();
     2383
     2384    // restore
     2385    columnLower_[iColumn] = saveLower;
     2386    memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double));
     2387    memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double));
     2388    memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char));
     2389
     2390    if (problemStatus_==0&&!isDualObjectiveLimitReached()) {
     2391      objectiveChange = objectiveValue_-saveObjective;
     2392    } else {
     2393      objectiveChange = 1.0e100;
     2394    }
     2395    newLower[i]=objectiveChange;
     2396#ifdef CLP_DEBUG
     2397    printf("up on %d costs %g\n",iColumn,objectiveChange);
     2398#endif
     2399         
     2400    /* Possibilities are:
     2401       Both sides feasible - store
     2402       Neither side feasible - set objective high and exit
     2403       One side feasible - change bounds and resolve
     2404    */
     2405    if (newUpper[i]<1.0e100) {
     2406      if(newLower[i]<1.0e100) {
     2407        // feasible - no action
     2408      } else {
     2409        // up feasible, down infeasible
     2410        returnCode=1;
     2411        if (stopOnFirstInfeasible)
     2412          break;
     2413      }
     2414    } else {
     2415      if(newLower[i]<1.0e100) {
     2416        // down feasible, up infeasible
     2417        returnCode=1;
     2418        if (stopOnFirstInfeasible)
     2419          break;
     2420      } else {
     2421        // neither side feasible
     2422        returnCode=-1;
     2423        break;
     2424      }
     2425    }
     2426  }
     2427  delete [] rowSolution;
     2428  delete [] columnSolution;
     2429  delete [] saveStatus;
     2430#endif
     2431  objectiveValue_ = saveObjective;
     2432  return returnCode;
     2433}
     2434// treat no pivot as finished (unless interesting)
     2435int ClpSimplexDual::fastDual(bool alwaysFinish)
     2436{
     2437  algorithm_ = -1;
     2438  dualTolerance_=dblParam_[OsiDualTolerance];
     2439  primalTolerance_=dblParam_[OsiPrimalTolerance];
     2440
     2441  // save dual bound
     2442  double saveDualBound = dualBound_;
     2443
     2444  double objectiveChange;
     2445  // for dual we will change bounds using dualBound_
     2446  // for this we need clean basis so it is after factorize
     2447  gutsOfSolution(rowActivityWork_,columnActivityWork_);
     2448  numberFake_ =0; // Number of variables at fake bounds
     2449  changeBounds(true,NULL,objectiveChange);
     2450
     2451  problemStatus_ = -1;
     2452  numberIterations_=0;
     2453
     2454  int lastCleaned=0; // last time objective or bounds cleaned up
     2455
     2456  // number of times we have declared optimality
     2457  numberTimesOptimal_=0;
     2458
     2459  // This says whether to restore things etc
     2460  int factorType=0;
     2461  /*
     2462    Status of problem:
     2463    0 - optimal
     2464    1 - infeasible
     2465    2 - unbounded
     2466    -1 - iterating
     2467    -2 - factorization wanted
     2468    -3 - redo checking without factorization
     2469    -4 - looks infeasible
     2470
     2471    BUT also from whileIterating return code is:
     2472
     2473   -1 iterations etc
     2474   -2 inaccuracy
     2475   -3 slight inaccuracy (and done iterations)
     2476   +0 looks optimal (might be unbounded - but we will investigate)
     2477   +1 looks infeasible
     2478   +3 max iterations
     2479
     2480  */
     2481
     2482  int returnCode = 0;
     2483
     2484  while (problemStatus_<0) {
     2485    int iRow,iColumn;
     2486    // clear
     2487    for (iRow=0;iRow<4;iRow++) {
     2488      rowArray_[iRow]->clear();
     2489    }   
     2490   
     2491    for (iColumn=0;iColumn<2;iColumn++) {
     2492      columnArray_[iColumn]->clear();
     2493    }   
     2494
     2495    // give matrix (and model costs and bounds a chance to be
     2496    // refreshed (normally null)
     2497    matrix_->refresh(this);
     2498    // may factorize, checks if problem finished
     2499    // should be able to speed this up on first time
     2500    statusOfProblemInDual(lastCleaned,factorType);
     2501
     2502    // Say good factorization
     2503    factorType=1;
     2504
     2505    // Do iterations
     2506    if (problemStatus_<0) {
     2507#if 1
     2508      returnCode = whileIterating();
     2509      if (!alwaysFinish&&returnCode<1) {
     2510        double limit = 0.0;
     2511        getDblParam(OsiDualObjectiveLimit, limit);
     2512        if(objectiveValue()*optimizationDirection_<limit||
     2513           numberAtFakeBound()) {
     2514          // can't say anything interesting - might as well return
     2515#ifdef CLP_DEBUG
     2516          printf("returning from fastDual after %d iterations with code %d\n",
     2517                 numberIterations_,returnCode);
     2518#endif
     2519          returnCode=1;
     2520          break;
     2521        }
     2522      }
     2523      returnCode=0;
     2524#else
     2525      whileIterating();
     2526#endif
     2527    }
     2528  }
     2529
     2530  assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
     2531  dualBound_ = saveDualBound;
     2532  return returnCode;
     2533}
     2534/* Checks number of variables at fake bounds.  This is used by fastDual
     2535   so can exit gracefully before end */
     2536int
     2537ClpSimplexDual::numberAtFakeBound()
     2538{
     2539  int iSequence;
     2540  int numberFake=0;
     2541 
     2542  for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     2543    FakeBound bound = getFakeBound(iSequence);
     2544    switch(getStatus(iSequence)) {
     2545
     2546    case ClpSimplex::basic:
     2547      break;
     2548    case ClpSimplex::isFree:
     2549    case ClpSimplex::superBasic:
     2550      break;
     2551    case ClpSimplex::atUpperBound:
     2552      if (bound==ClpSimplex::upperFake||bound==ClpSimplex::bothFake)
     2553        numberFake++;
     2554      break;
     2555    case ClpSimplex::atLowerBound:
     2556      if (bound==ClpSimplex::lowerFake||bound==ClpSimplex::bothFake)
     2557        numberFake++;
     2558      break;
     2559    }
     2560  }
     2561  return numberFake;
     2562}
  • branches/devel-1/ClpSimplexPrimal.cpp

    r15 r17  
    376376  return problemStatus_;
    377377}
    378 void
     378/*
     379  Reasons to come out:
     380  -1 iterations etc
     381  -2 inaccuracy
     382  -3 slight inaccuracy (and done iterations)
     383  -4 end of values pass and done iterations
     384  +0 looks optimal (might be infeasible - but we will investigate)
     385  +2 looks unbounded
     386  +3 max iterations
     387*/
     388int
    379389ClpSimplexPrimal::whileIterating(int & firstSuperBasic)
    380390{
     
    382392  // Say if values pass
    383393  int ifValuesPass=0;
     394  int returnCode=-1;
    384395  if (firstSuperBasic<numberRows_+numberColumns_)
    385396    ifValuesPass=1;
     
    451462          problemStatus_=-2; // factorize now
    452463          pivotRow_=-1; // say no weights update
     464          returnCode=-4;
    453465          break;
    454466        }
     
    516528          rowArray_[1]->clear();
    517529          pivotRow_=-1; // say no weights update
     530          returnCode=-2;
    518531          break;
    519532        } else {
     
    541554        if (updateStatus==1) {
    542555          // slight error
    543           if (factorization_->pivots()>5)
     556          if (factorization_->pivots()>5) {
    544557            problemStatus_=-2; // factorize now
     558            returnCode=-3;
     559          }
    545560        } else if (updateStatus==2) {
    546561          // major error
     
    548563          if(saveNumber != numberIterations_) {
    549564            problemStatus_=-2; // factorize now
     565            returnCode=-2;
    550566            break;
    551567          } else {
     
    603619          }
    604620          rowArray_[0]->clear();
     621          returnCode=2;
    605622          break;
    606623        } else {
     
    628645        // maximum iterations or equivalent
    629646        problemStatus_= 3;
     647        returnCode=3;
    630648        break;
    631649      }
     
    638656      if (nonLinearCost_->numberInfeasibilities())
    639657        problemStatus_=-4; // might be infeasible
     658      returnCode=0;
    640659      break;
    641660    }
    642661  }
     662  return returnCode;
    643663}
    644664/* Checks if finished.  Updates status */
     
    658678  }
    659679  int tentativeStatus = problemStatus_;
    660 
    661680  if (problemStatus_>-3||problemStatus_==-4) {
    662681    // factorize
     
    666685    // This may save pivotRow_ for use
    667686    primalColumnPivot_->saveWeights(this,1);
    668     // is factorization okay?
    669     if (internalFactorize(1)) {
    670       // no - restore previous basis
    671       assert (type==1);
    672       memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
    673       memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
    674              numberRows_*sizeof(double));
    675       memcpy(columnActivityWork_,savedSolution_ ,
    676              numberColumns_*sizeof(double));
    677       forceFactorization_=1; // a bit drastic but ..
    678       type = 2;
    679       assert (internalFactorize(1)==0);
    680       changeMade_++; // say change made
     687
     688    if (type) {
     689      // is factorization okay?
     690      if (internalFactorize(1)) {
     691        // no - restore previous basis
     692        assert (type==1);
     693        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
     694        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
     695               numberRows_*sizeof(double));
     696        memcpy(columnActivityWork_,savedSolution_ ,
     697               numberColumns_*sizeof(double));
     698        forceFactorization_=1; // a bit drastic but ..
     699        type = 2;
     700        assert (internalFactorize(1)==0);
     701        changeMade_++; // say change made
     702      }
    681703    }
    682704    if (problemStatus_!=-4)
     
    9841006      index[numberRemaining++]=iRow;
    9851007      double value=oldValue-upperTheta*fabs(alpha);
    986       if (value<-primalTolerance_)
     1008      if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
    9871009        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
    9881010    }
     
    11191141              // back to previous one
    11201142              goBackOne = true;
     1143              //break;
    11211144            } else {
    1122               // can only get here if all pivots too small
     1145              // can only get here if all pivots so far too small
    11231146            }
    11241147            break;
  • branches/devel-1/Makefile

    r16 r17  
    4242include ${MakefileDir}/Makefile.location
    4343# This modification seems to be needed
    44 # export ExtraIncDir := ../Osi/include
     44export ExtraIncDir := ../Osi/include
    4545export ExtraLibDir :=
    4646export ExtraLibName :=
  • branches/devel-1/include/ClpDualRowPivot.hpp

    r2 r17  
    4040  virtual void updatePrimalSolution(OsiIndexedVector * input,
    4141                                    double theta,
    42                                     double & changeInObjective)=0;
     42                                    double & changeInObjective) =0;
    4343  /** Saves any weights round factorization as pivot rows may change
    4444      Will be empty unless steepest edge (will save model)
     
    4949         (e.g. check for infeasible)
    5050      4) as 2 but restore weights from previous snapshot
     51      5) for strong branching - initialize  , infeasibilities
    5152  */
    5253  virtual void saveWeights(ClpSimplex * model,int mode);
     
    8384  { return model_;};
    8485 
    85   /// Returns type
     86  /// Returns type (above 63 is extra information)
    8687  inline int type()
    8788  { return type_;};
  • branches/devel-1/include/ClpDualRowSteepest.hpp

    r2 r17  
    4747         (e.g. check for infeasible)
    4848      4) as 2 but restore weights from previous snapshot
     49      5) for strong branching - initialize (uninitialized) , infeasibilities
    4950  */
    5051  virtual void saveWeights(ClpSimplex * model, int mode);
     
    7374  virtual ClpDualRowPivot * clone(bool copyData = true) const;
    7475 
    75   //@}
     76   //@}
     77  /**@name gets and sets */
     78  //@{
     79  /// Mode
     80  inline int mode() const
     81    { return mode_;};
     82 //@}
    7683
    7784  //---------------------------------------------------------------------------
  • branches/devel-1/include/ClpModel.hpp

    r14 r17  
    149149  { return problemStatus_==2;};
    150150  /// Is the given primal objective limit reached?
    151   bool isPrimalObjectiveLimitReached() const
    152   { return false;}; // not implemented yet
     151  bool isPrimalObjectiveLimitReached() const ;
    153152  /// Is the given dual objective limit reached?
    154   bool isDualObjectiveLimitReached() const
    155   { return false;}; // not implememented yet
     153  bool isDualObjectiveLimitReached() const ;
    156154  /// Iteration limit reached?
    157155  bool isIterationLimitReached() const
     
    227225  inline double getObjValue() const
    228226  { return objectiveValue_*optimizationDirection_ - dblParam_[OsiObjOffset];};
     227  /// Integer information
     228  inline char * integerInformation() const
     229  {return integerType_;};
    229230  /** Infeasibility/unbounded ray (NULL returned if none/wrong)
    230231      Up to user to use delete [] on these arrays.  */
     
    321322  /**@name private or protected methods */
    322323  //@{
    323 private:
     324protected:
    324325  /// Does most of deletion
    325326  void gutsOfDelete();
     
    327328      If trueCopy false then just points to arrays */
    328329  void gutsOfCopy(const ClpModel & rhs, bool trueCopy=true);
    329 protected:
    330330  /// gets lower and upper bounds on rows
    331331  void getRowBound(int iRow, double& lower, double& upper) const;
     
    381381  double dblParam_[OsiLastDblParam];
    382382  /// Array of string parameters
    383   std::string strParam_[OsiLastDblParam];
     383  std::string strParam_[OsiLastStrParam];
    384384  /// Objective value
    385385  double objectiveValue_;
     
    396396  /// Messages
    397397  OsiMessages messages_;
    398   /// length of names (0 means no names0
     398  /// length of names (0 means no names)
    399399  int lengthNames_;
    400400  /// Row names
     
    402402  /// Column names
    403403  std::vector<std::string> columnNames_;
     404  /// Integer information
     405  char * integerType_;
    404406  //@}
    405407};
  • branches/devel-1/include/ClpPackedMatrix.hpp

    r2 r17  
    152152   ClpPackedMatrix(const OsiPackedMatrix&);
    153153
     154  /** This takes over ownership (for space reasons) */
     155   ClpPackedMatrix(OsiPackedMatrix * matrix);
     156
    154157   ClpPackedMatrix& operator=(const ClpPackedMatrix&);
    155158  /// Clone
  • branches/devel-1/include/ClpPrimalColumnPivot.hpp

    r2 r17  
    8484  { return model_;};
    8585 
    86   /// Returns type
     86  /// Returns type (above 63 is extra information)
    8787  inline int type()
    8888  { return type_;};
  • branches/devel-1/include/ClpPrimalColumnSteepest.hpp

    r2 r17  
    5050  //@}
    5151 
     52  /**@name gets and sets */
     53  //@{
     54  /// Mode
     55  inline int mode() const
     56    { return mode_;};
     57 //@}
     58
    5259 
    5360  ///@name Constructors and destructors
  • branches/devel-1/include/ClpSimplex.hpp

    r15 r17  
    7575  /// Destructor
    7676   ~ClpSimplex (  );
     77  /** Borrow model.  This is so we dont have to copy large amounts
     78      of data around.  It assumes a derived class wants to overwrite
     79      an empty model with a real one - while it does an algorithm.
     80      This is same as ClpModel one, but sets scaling on etc. */
     81  void borrowModel(ClpModel & otherModel);
    7782  //@}
    7883
     
    101106  /// Sets column pivot choice algorithm in primal
    102107  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
     108  /** For strong branching.  On input lower and upper are new bounds
     109      while on output they are change in objective function values
     110      (>1.0e50 infeasible).
     111      Return code is 0 if nothing interesting, -1 if infeasible both
     112      ways and +1 if infeasible one way (check values to see which one(s))
     113  */
     114  int strongBranching(int numberVariables,const int * variables,
     115                      double * newLower, double * newUpper,
     116                      bool stopOnFirstInfeasible=true,
     117                      bool alwaysFinish=false);
    103118  //@}
    104119
     
    160175  /// Warm start
    161176  OsiWarmStartBasis getBasis() const;
     177  /** Save model to file, returns 0 if success.  This is designed for
     178      use outside algorithms so does not save iterating arrays etc.
     179  It does not save any messaging information.
     180  Does not save scaling values.
     181  It does not know about all types of virtual functions.
     182  */
     183  int saveModel(const char * fileName);
     184  /** Restore model from file, returns 0 if success,
     185      deletes current model */
     186  int restoreModel(const char * fileName);
     187 
    162188  //@}
    163189
     
    296322  /// Pivot Row for use by classes e.g. steepestedge
    297323  inline int pivotRow() const{ return pivotRow_;};
     324  /// value of incoming variable (in Dual)
     325  double valueIncomingDual() const;
    298326  //@}
    299327
     
    618646  /// So we know when to be cautious
    619647  int lastBadIteration_;
     648  /// Can be used for count of fake bounds (dual) or fake costs (primal)
     649  int numberFake_;
    620650  //@}
    621651};
  • branches/devel-1/include/ClpSimplexDual.hpp

    r14 r17  
    115115
    116116  int dual();
     117  /** For strong branching.  On input lower and upper are new bounds
     118      while on output they are change in objective function values
     119      (>1.0e50 infeasible).
     120      Return code is 0 if nothing interesting, -1 if infeasible both
     121      ways and +1 if infeasible one way (check values to see which one(s))
     122  */
     123  int strongBranching(int numberVariables,const int * variables,
     124                      double * newLower, double * newUpper,
     125                      bool stopOnFirstInfeasible=true,
     126                      bool alwaysFinish=false);
    117127  //@}
    118128
     
    121131  /** This has the flow between re-factorizations
    122132      Broken out for clarity and will be used by strong branching
     133
     134      Reasons to come out:
     135      -1 iterations etc
     136      -2 inaccuracy
     137      -3 slight inaccuracy (and done iterations)
     138      +0 looks optimal (might be unbounded - but we will investigate)
     139      +1 looks infeasible
     140      +3 max iterations
    123141   */
    124   void whileIterating();
     142  int whileIterating();
    125143  /** The duals are updated by the given arrays.
    126144      Returns number of infeasibilities.
     
    194212  /// Perturbs problem (method depends on perturbation())
    195213  void perturb();
     214  /** Fast iterations.  Misses out a lot of initialization.
     215      Normally stops on maximum iterations, first re-factorization
     216      or tentative optimum.  If looks interesting then continues as
     217      normal.  Returns 0 if finished properly, 1 otherwise.
     218  */
     219  int fastDual(bool alwaysFinish=false);
     220  /** Checks number of variables at fake bounds.  This is used by fastDual
     221      so can exit gracefully before end */
     222  int numberAtFakeBound();
    196223  //@}
    197224};
  • branches/devel-1/include/ClpSimplexPrimal.hpp

    r15 r17  
    125125      firstSuperBasic == number rows + columns normally,
    126126      otherwise first super basic variable
     127
     128      Reasons to come out:
     129      -1 iterations etc
     130      -2 inaccuracy
     131      -3 slight inaccuracy (and done iterations)
     132      -4 end of values pass and done iterations
     133      +0 looks optimal (might be infeasible - but we will investigate)
     134      +2 looks unbounded
     135      +3 max iterations
    127136   */
    128   void whileIterating(int & firstSuperBasic);
     137  int whileIterating(int & firstSuperBasic);
    129138  /** The primals are updated by the given array.
    130139      Returns number of infeasibilities.
Note: See TracChangeset for help on using the changeset viewer.