Changeset 109


Ignore:
Timestamp:
Jan 31, 2003 3:32:13 PM (17 years ago)
Author:
forrest
Message:

To help OsiSimplexInterface?

Location:
trunk
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpModel.cpp

    r82 r109  
    6464  objectiveValue_(0.0),
    6565  numberIterations_(0),
     66  solveType_(0),
    6667  problemStatus_(-1),
    6768  defaultHandler_(true),
     
    324325  objectiveValue_=rhs.objectiveValue_;
    325326  numberIterations_ = rhs.numberIterations_;
     327  solveType_ = rhs.solveType_;
    326328  problemStatus_ = rhs.problemStatus_;
    327329  if (trueCopy) {
     
    952954  }
    953955}
     956//#############################################################################
     957// Constructors / Destructor / Assignment
     958//#############################################################################
     959
     960//-------------------------------------------------------------------
     961// Default Constructor
     962//-------------------------------------------------------------------
     963ClpDataSave::ClpDataSave ()
     964{
     965  dualBound_ = 0.0;
     966  infeasibilityCost_ = 0.0;
     967  sparseThreshold_ = 0;
     968  perturbation_ = 0;
     969}
     970
     971//-------------------------------------------------------------------
     972// Copy constructor
     973//-------------------------------------------------------------------
     974ClpDataSave::ClpDataSave (const ClpDataSave & rhs)
     975
     976  dualBound_ = rhs.dualBound_;
     977  infeasibilityCost_ = rhs.infeasibilityCost_;
     978  sparseThreshold_ = rhs.sparseThreshold_;
     979  perturbation_ = rhs.perturbation_;
     980}
     981
     982//-------------------------------------------------------------------
     983// Destructor
     984//-------------------------------------------------------------------
     985ClpDataSave::~ClpDataSave ()
     986{
     987
     988}
     989
     990//----------------------------------------------------------------
     991// Assignment operator
     992//-------------------------------------------------------------------
     993ClpDataSave &
     994ClpDataSave::operator=(const ClpDataSave& rhs)
     995{
     996  if (this != &rhs) {
     997    dualBound_ = rhs.dualBound_;
     998    infeasibilityCost_ = rhs.infeasibilityCost_;
     999    sparseThreshold_ = rhs.sparseThreshold_;
     1000    perturbation_ = rhs.perturbation_;
     1001  }
     1002  return *this;
     1003}
  • trunk/ClpPackedMatrix.cpp

    r93 r109  
    843843      double value = fabs(elementByColumn[j]);
    844844      int iRow = row[j];
     845      if (iRow<0||iRow>=numberRows) {
     846        printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
     847        return false;
     848      }
    845849      if (mark[iRow]==-1) {
    846850        mark[iRow]=j;
  • trunk/ClpSimplex.cpp

    r108 r109  
    946946  }
    947947  // change of incoming
    948   if (algorithm_<0) {
    949     dualOut_ /= alpha_;
    950     dualOut_ *= -directionOut_;
    951   }
    952948  char rowcol[]={'R','C'};
    953   double cost = cost_[sequenceIn_];
    954   double value=valueIn_;
    955949  if (pivotRow_>=0)
    956950    pivotVariable_[pivotRow_]=sequenceIn();
    957   if (directionIn_==-1) {
    958     // as if from upper bound
    959     if (sequenceIn_!=sequenceOut_) {
    960       // variable becoming basic
    961       setStatus(sequenceIn_,basic);
    962       if (algorithm_<0) {
    963         value = upperIn_+dualOut_;
    964         dj_[sequenceIn_]=0.0;
    965       } else {
    966         value = valueIn_-fabs(theta_);
    967       }
    968     } else {
    969       value=lowerIn_;
    970       setStatus(sequenceIn_, atLowerBound);
    971     }
    972   } else {
    973     // as if from lower bound
    974     if (sequenceIn_!=sequenceOut_) {
    975       // variable becoming basic
    976       setStatus(sequenceIn_,basic);
    977       if (algorithm_<0) {
    978         value = lowerIn_+dualOut_;
    979         dj_[sequenceIn_]=0.0;
    980       } else {
    981         value = valueIn_+fabs(theta_);
    982       }
    983     } else {
    984       value=upperIn_;
    985       setStatus(sequenceIn_, atUpperBound);
    986     }
    987   }
    988951  if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<-1.0e20)
    989952    progressFlag_ |= 2; // making real progress
    990   if (algorithm_<0)
    991     objectiveChange += cost*(value-valueIn_);
    992   else
    993     objectiveChange += dualIn_*(value-valueIn_);
    994   solution_[sequenceIn_]=value;
    995 
    996   // outgoing
     953  solution_[sequenceIn_]=valueIn_;
     954  if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
     955    progressFlag_ |= 1; // making real progress
    997956  if (sequenceIn_!=sequenceOut_) {
    998     assert( getStatus(sequenceOut_)== basic);
    999     if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
    1000       progressFlag_ |= 1; // making real progress
    1001     if (algorithm_<0) {
    1002       if (directionOut_>0) {
    1003         value = lowerOut_;
    1004         setStatus(sequenceOut_,atLowerBound);
    1005         dj_[sequenceOut_]=theta_;
    1006       } else {
    1007         value = upperOut_;
    1008         setStatus(sequenceOut_,atUpperBound);
    1009         dj_[sequenceOut_]=-theta_;
    1010       }
    1011       solution_[sequenceOut_]=value;
     957    //assert( getStatus(sequenceOut_)== basic);
     958    setStatus(sequenceIn_,basic);
     959    if (directionOut_>0) {
     960      // going to lower
     961      setStatus(sequenceOut_,atLowerBound);
    1012962    } else {
    1013       if (directionOut_>0) {
    1014         value = lowerOut_;
    1015       } else {
    1016         value = upperOut_;
    1017       }
    1018       double lowerValue = lower_[sequenceOut_];
    1019       double upperValue = upper_[sequenceOut_];
    1020       assert(value>=lowerValue-primalTolerance_&&
    1021              value<=upperValue+primalTolerance_);
    1022       // may not be exactly at bound and bounds may have changed
    1023       if (value<=lowerValue+primalTolerance_) {
    1024         setStatus(sequenceOut_,atLowerBound);
    1025       } else if (value>=upperValue-primalTolerance_) {
    1026         setStatus(sequenceOut_,atUpperBound);
    1027       } else {
    1028         printf("*** variable wandered off bound %g %g %g!\n",
    1029                lowerValue,value,upperValue);
    1030         if (value-lowerValue<=upperValue-value) {
    1031           setStatus(sequenceOut_,atLowerBound);
    1032           value=lowerValue;
    1033         } else {
    1034           setStatus(sequenceOut_,atUpperBound);
    1035           value=upperValue;
    1036         }
    1037       }
    1038       solution_[sequenceOut_]=value;
    1039       nonLinearCost_->setOne(sequenceOut_,value);
    1040     }
    1041   }
    1042   // change cost and bounds on incoming if primal
    1043   if (algorithm_>0)
    1044     nonLinearCost_->setOne(sequenceIn_,solution_[sequenceIn_]);
     963      // going to upper
     964      setStatus(sequenceOut_,atUpperBound);
     965    }
     966    solution_[sequenceOut_]=valueOut_;
     967  } else {
     968    // flip from bound to bound
     969    if (directionIn_==-1) {
     970      // as if from upper bound
     971      setStatus(sequenceIn_, atLowerBound);
     972    } else {
     973      // as if from lower bound
     974      setStatus(sequenceIn_, atUpperBound);
     975    }
     976  }
    1045977  objectiveValue_ += objectiveChange;
    1046978  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
     
    11591091  gutsOfDelete(0);
    11601092  gutsOfCopy(rhs);
     1093  solveType_=1; // say simplex based life form
    11611094}
    11621095// Copy constructor from model
     
    23362269  ClpModel::borrowModel(otherModel);
    23372270  createStatus();
    2338   scaling();
    23392271  ClpDualRowSteepest steep1;
    23402272  setDualRowPivotAlgorithm(steep1);
     
    36723604          } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
    36733605            setStatus(iColumn,isFree);
     3606            break;
    36743607          } else {
    36753608            break;
     
    36933626  }
    36943627}
     3628/* Pivot in a variable and out a variable.  Returns 0 if okay,
     3629   1 if inaccuracy forced re-factorization, -1 if would be singular.
     3630   Also updates primal/dual infeasibilities.
     3631   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
     3632*/
     3633int ClpSimplex::pivot()
     3634{
     3635  // scaling not allowed
     3636  assert (!scalingFlag_);
     3637  // assume In_ and Out_ are correct and directionOut_ set
     3638  // (or In_ if flip
     3639  lowerIn_ = lower_[sequenceIn_];
     3640  valueIn_ = solution_[sequenceIn_];
     3641  upperIn_ = upper_[sequenceIn_];
     3642  dualIn_ = dj_[sequenceIn_];
     3643  if (sequenceOut_>=0) {
     3644    assert (pivotRow_>=0&&pivotRow_<numberRows_);
     3645    assert (pivotVariable_[pivotRow_]==sequenceOut_);
     3646    lowerOut_ = lower_[sequenceOut_];
     3647    valueOut_ = solution_[sequenceOut_];
     3648    upperOut_ = upper_[sequenceOut_];
     3649    // for now assume primal is feasible (or in dual)
     3650    dualOut_ = dj_[sequenceOut_];
     3651    assert(fabs(dualOut_)<1.0e-6);
     3652  } else {
     3653    assert (pivotRow_<0);
     3654  }
     3655  unpack(rowArray_[1]);
     3656  factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
     3657  // we are going to subtract movement from current basic
     3658  double movement;
     3659  // see where incoming will go to
     3660  if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
     3661    // flip so go to bound
     3662    movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
     3663  } else {
     3664    // get where outgoing needs to get to
     3665    double outValue = (directionOut_>0) ? upperOut_ : lowerOut_;
     3666    // solutionOut_ - movement*alpha_ == outValue
     3667    movement = (outValue-valueOut_)/alpha_;
     3668    // set directionIn_ correctly
     3669    directionIn_ = (movement>0) ? 1 :-1;
     3670  }
     3671  // update primal solution
     3672  {
     3673    int i;
     3674    int * index = rowArray_[1]->getIndices();
     3675    int number = rowArray_[1]->getNumElements();
     3676    double * element = rowArray_[1]->denseVector();
     3677    for (i=0;i<number;i++) {
     3678      int ii = index[i];
     3679      // get column
     3680      ii = pivotVariable_[ii];
     3681      solution_[ii] -= movement*element[i];
     3682      element[i]=0.0;
     3683    }
     3684    // see where something went to
     3685    if (sequenceOut_<0) {
     3686      if (directionIn_<0) {
     3687        assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
     3688        solution_[sequenceIn_]=upperIn_;
     3689      } else {
     3690        assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
     3691        solution_[sequenceIn_]=lowerIn_;
     3692      }
     3693    } else {
     3694      if (directionOut_<0) {
     3695        assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
     3696        solution_[sequenceOut_]=upperOut_;
     3697      } else {
     3698        assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
     3699        solution_[sequenceOut_]=lowerOut_;
     3700      }
     3701      solution_[sequenceIn_]=valueIn_+movement;
     3702    }
     3703  }   
     3704  // update duals
     3705  if (pivotRow_>=0) {
     3706    alpha_ = rowArray_[1]->denseVector()[pivotRow_];
     3707    assert (fabs(alpha_)>1.0e-8);
     3708    double multiplier = dualIn_/alpha_;
     3709    rowArray_[0]->insert(pivotRow_,multiplier);
     3710    factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
     3711    // put row of tableau in rowArray[0] and columnArray[0]
     3712    matrix_->transposeTimes(this,-1.0,
     3713                            rowArray_[0],columnArray_[1],columnArray_[0]);
     3714    // update column djs
     3715    int i;
     3716    int * index = columnArray_[0]->getIndices();
     3717    int number = columnArray_[0]->getNumElements();
     3718    double * element = columnArray_[0]->denseVector();
     3719    for (i=0;i<number;i++) {
     3720      int ii = index[i];
     3721      dj_[ii] -= element[i];
     3722      element[i]=0.0;
     3723    }
     3724    columnArray_[0]->setNumElements(0);
     3725    // and row djs
     3726    index = rowArray_[0]->getIndices();
     3727    number = rowArray_[0]->getNumElements();
     3728    element = rowArray_[0]->denseVector();
     3729    for (i=0;i<number;i++) {
     3730      int ii = index[i];
     3731      dj_[ii+numberRows_] -= element[i];
     3732      element[i]=0.0;
     3733    }
     3734    rowArray_[0]->setNumElements(0);
     3735    // check incoming
     3736    assert (fabs(dj_[sequenceIn_])<1.0e-6);
     3737  }
     3738#if 0
     3739  sol update
     3740    eta update
     3741  int save = algorithm_;
     3742  // make simple so always primal
     3743  algorithm_=1;
     3744  housekeeping();
     3745  algorithm_=save;
     3746  checksol
     3747#endif
     3748  abort();
     3749  return 0;
     3750}
     3751
     3752/* Pivot in a variable and choose an outgoing one.  Assumes primal
     3753   feasible - will not go through a bound.  Returns step length in theta
     3754   Returns ray in ray_ (or NULL if no pivot)
     3755   Return codes as before but -1 means no acceptable pivot
     3756*/
     3757int ClpSimplex::primalPivotResult()
     3758{
     3759  return ((ClpSimplexPrimal *) this)->pivotResult();
     3760}
     3761 
     3762/* Pivot out a variable and choose an incoing one.  Assumes dual
     3763   feasible - will not go through a reduced cost. 
     3764   Returns step length in theta
     3765   Returns ray in ray_ (or NULL if no pivot)
     3766   Return codes as before but -1 means no acceptable pivot
     3767*/
     3768int
     3769ClpSimplex::dualPivotResult()
     3770{
     3771  return ((ClpSimplexDual *) this)->pivotResult();
     3772}
     3773// Common bits of coding for dual and primal
     3774int
     3775ClpSimplex::startup(int ifValuesPass)
     3776{
     3777  // sanity check
     3778  assert (numberRows_==matrix_->getNumRows());
     3779  assert (numberColumns_==matrix_->getNumCols());
     3780  // for moment all arrays must exist
     3781  assert(columnLower_);
     3782  assert(columnUpper_);
     3783  assert(rowLower_);
     3784  assert(rowUpper_);
     3785
     3786
     3787  primalTolerance_=dblParam_[ClpPrimalTolerance];
     3788  dualTolerance_=dblParam_[ClpDualTolerance];
     3789
     3790  // put in standard form (and make row copy)
     3791  // create modifiable copies of model rim and do optional scaling
     3792  bool goodMatrix=createRim(7+8+16,true);
     3793
     3794  if (goodMatrix) {
     3795    // Model looks okay
     3796    // Do initial factorization
     3797    // and set certain stuff
     3798    // We can either set increasing rows so ...IsBasic gives pivot row
     3799    // or we can just increment iBasic one by one
     3800    // for now let ...iBasic give pivot row
     3801    factorization_->increasingRows(2);
     3802    // row activities have negative sign
     3803    factorization_->slackValue(-1.0);
     3804    factorization_->zeroTolerance(1.0e-13);
     3805
     3806    // do perturbation if asked for
     3807
     3808    if (perturbation_<100) {
     3809      if (algorithm_>0) {
     3810        ((ClpSimplexPrimal *) this)->perturb();
     3811      } else if (algorithm_<0) {
     3812        ((ClpSimplexDual *) this)->perturb();
     3813      }
     3814    }
     3815    // for primal we will change bounds using infeasibilityCost_
     3816    if (nonLinearCost_==NULL&&algorithm_>0) {
     3817      // get a valid nonlinear cost function
     3818      delete nonLinearCost_;
     3819      nonLinearCost_= new ClpNonLinearCost(this);
     3820    }
     3821   
     3822    // loop round to clean up solution if values pass
     3823    int numberThrownOut = -1;
     3824    int totalNumberThrownOut=0;
     3825    while(numberThrownOut) {
     3826      int status = internalFactorize(0+10*ifValuesPass);
     3827      if (status<0)
     3828        return 1; // some error
     3829      else
     3830        totalNumberThrownOut+= status;
     3831     
     3832      // for this we need clean basis so it is after factorize
     3833      numberThrownOut=gutsOfSolution(rowActivityWork_,columnActivityWork_,
     3834                                     ifValuesPass!=0);
     3835      totalNumberThrownOut+= numberThrownOut;
     3836     
     3837    }
     3838   
     3839    if (totalNumberThrownOut)
     3840      handler_->message(CLP_SINGULARITIES,messages_)
     3841        <<totalNumberThrownOut
     3842        <<CoinMessageEol;
     3843   
     3844    problemStatus_ = -1;
     3845    numberIterations_=0;
     3846   
     3847    // number of times we have declared optimality
     3848    numberTimesOptimal_=0;
     3849
     3850    return 0;
     3851  } else {
     3852    // bad matrix
     3853    return 2;
     3854  }
     3855   
     3856}
     3857
     3858
     3859void
     3860ClpSimplex::finish()
     3861{
     3862  // Get rid of some arrays and empty factorization
     3863  deleteRim();
     3864  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
     3865    <<objectiveValue()
     3866    <<CoinMessageEol;
     3867  factorization_->relaxAccuracyCheck(1.0);
     3868}
     3869// Save data
     3870ClpDataSave
     3871ClpSimplex::saveData() const
     3872{
     3873  ClpDataSave saved;
     3874  saved.dualBound_ = dualBound_;
     3875  saved.infeasibilityCost_ = infeasibilityCost_;
     3876  saved.sparseThreshold_ = factorization_->sparseThreshold();
     3877  saved.perturbation_ = perturbation_;
     3878  return saved;
     3879}
     3880// Restore data
     3881void
     3882ClpSimplex::restoreData(ClpDataSave saved)
     3883{
     3884  factorization_->sparseThreshold(saved.sparseThreshold_);
     3885  perturbation_ = saved.perturbation_;
     3886  infeasibilityCost_ = saved.infeasibilityCost_;
     3887  dualBound_ = saved.dualBound_;
     3888}
  • trunk/ClpSimplexDual.cpp

    r105 r109  
    197197  */
    198198
     199  algorithm_ = -1;
     200
     201  // save data
     202  ClpDataSave data = saveData();
    199203
    200204  // sanity check
    201   assert (numberRows_==matrix_->getNumRows());
    202   assert (numberColumns_==matrix_->getNumCols());
    203   // for moment all arrays must exist
    204   assert(columnLower_);
    205   assert(columnUpper_);
    206   assert(rowLower_);
    207   assert(rowUpper_);
    208 
    209 
    210   algorithm_ = -1;
    211   dualTolerance_=dblParam_[ClpDualTolerance];
    212   primalTolerance_=dblParam_[ClpPrimalTolerance];
    213 
     205  // initialize - no values pass and algorithm_ is -1
    214206  // put in standard form (and make row copy)
    215207  // create modifiable copies of model rim and do optional scaling
    216   bool goodMatrix = createRim(7+8+16,true);
    217 
    218   // save dual bound
    219   double saveDualBound_ = dualBound_;
    220   // save perturbation
    221   int savePerturbation = perturbation_;
    222   // save if sparse factorization wanted
    223   int saveSparse = factorization_->sparseThreshold();
    224 
    225   if (goodMatrix) {
    226     // Problem looks okay
    227     int iRow,iColumn;
    228     // Do initial factorization
    229     // and set certain stuff
    230     // We can either set increasing rows so ...IsBasic gives pivot row
    231     // or we can just increment iBasic one by one
    232     // for now let ...iBasic give pivot row
    233     factorization_->increasingRows(2);
    234     // row activities have negative sign
    235     factorization_->slackValue(-1.0);
    236     factorization_->zeroTolerance(1.0e-13);
    237     // might want to do first time - factorization_->pivotTolerance(0.99);
    238 
    239     // Do Initial factorization
    240     int factorizationStatus = internalFactorize(0);
    241 
    242     if (factorizationStatus<0)
    243       return 1; // some error
    244     else if (factorizationStatus)
    245       handler_->message(CLP_SINGULARITIES,messages_)
    246         <<factorizationStatus
    247         <<CoinMessageEol;
    248    
    249     // If user asked for perturbation - do it
    250    
    251     if (perturbation_<100)
    252       perturb();
    253    
    254     // for dual we will change bounds using dualBound_
    255     // for this we need clean basis so it is after factorize
    256     gutsOfSolution(rowActivityWork_,columnActivityWork_);
     208  // If problem looks okay
     209  // Do initial factorization
     210  // If user asked for perturbation - do it
     211  if (!startup(0)) {
     212    // looks okay
    257213   
    258214    double objectiveChange;
     
    260216    changeBounds(true,NULL,objectiveChange);
    261217   
    262     problemStatus_ = -1;
    263     numberIterations_=0;
    264    
    265218    int lastCleaned=0; // last time objective or bounds cleaned up
    266219   
    267     // number of times we have declared optimality
    268     numberTimesOptimal_=0;
    269220
    270221    // Progress indicator
     
    284235    */
    285236    while (problemStatus_<0) {
     237      int iRow, iColumn;
    286238      // clear
    287239      for (iRow=0;iRow<4;iRow++) {
     
    307259      // Say good factorization
    308260      factorType=1;
    309       if (saveSparse) {
     261      if (data.sparseThreshold_) {
    310262        // use default at present
    311263        factorization_->sparseThreshold(0);
     
    319271
    320272  assert (problemStatus_||!sumPrimalInfeasibilities_);
    321   //assert(!numberFake_||problemStatus_); // all bounds should be okay
    322   factorization_->sparseThreshold(saveSparse);
    323   // Get rid of some arrays and empty factorization
    324   deleteRim();
    325   handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
    326     <<objectiveValue()
    327     <<CoinMessageEol;
     273
     274  // clean up
     275  finish();
     276
    328277  // Restore any saved stuff
    329   perturbation_ = savePerturbation;
    330   dualBound_ = saveDualBound_;
    331   factorization_->relaxAccuracyCheck(1.0);
     278  restoreData(data);
    332279  return problemStatus_;
    333280}
     
    630577        double oldobj=objectiveValue_;
    631578#endif
     579        // modify dualout
     580        dualOut_ /= alpha_;
     581        dualOut_ *= -directionOut_;
     582        //setStatus(sequenceIn_,basic);
     583        dj_[sequenceIn_]=0.0;
     584        double oldValue=valueIn_;
     585        if (directionIn_==-1) {
     586          // as if from upper bound
     587          valueIn_ = upperIn_+dualOut_;
     588        } else {
     589          // as if from lower bound
     590          valueIn_ = lowerIn_+dualOut_;
     591        }
     592        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
     593        // outgoing
     594        if (directionOut_>0) {
     595          valueOut_ = lowerOut_;
     596          //setStatus(sequenceOut_,atLowerBound);
     597          dj_[sequenceOut_]=theta_;
     598        } else {
     599          valueOut_ = upperOut_;
     600          //setStatus(sequenceOut_,atUpperBound);
     601          dj_[sequenceOut_]=-theta_;
     602        }
     603        solution_[sequenceOut_]=valueOut_;
    632604        int whatNext=housekeeping(objectiveChange);
    633605        // and set bounds correctly
     
    896868{
    897869  // get pivot row using whichever method it is
    898 #if 1
    899   // Doesnt seem to help - think harder about free variables
    900870  // first see if any free variables and put them in basis
    901871  int nextFree = nextSuperBasic();
     
    904874  if (nextFree>=0) {
    905875    // unpack vector and find a good pivot
    906     int saveIn=sequenceIn_;
    907     sequenceIn_=nextFree;
    908     unpack(rowArray_[1]);
     876    unpack(rowArray_[1],nextFree);
    909877    factorization_->updateColumn(rowArray_[2],rowArray_[1],false);
    910878
     
    949917    if (chosenRow>=0)
    950918      pivotRow_=chosenRow;
    951     sequenceIn_=saveIn;
    952919    rowArray_[1]->clear();
    953920  }
    954921  if (chosenRow<0)
    955 #endif
    956922    pivotRow_=dualRowPivot_->pivotRow();
    957923
     
    18191785  int tentativeStatus = problemStatus_;
    18201786  double changeCost;
     1787  bool unflagVariables = true;
    18211788
    18221789  if (problemStatus_>-3||factorization_->pivots()) {
     
    18301797      if (internalFactorize(1)) {
    18311798        // no - restore previous basis
     1799        unflagVariables = false;
    18321800        assert (type==1);
    18331801        changeMade_++; // say something changed
     
    21432111    dualRowPivot_->saveWeights(this,3);
    21442112  // unflag all variables (we may want to wait a bit?)
    2145   if (tentativeStatus!=-2) {
     2113  if (tentativeStatus!=-2&&unflagVariables) {
    21462114    int iRow;
    21472115    for (iRow=0;iRow<numberRows_;iRow++) {
     
    27752743  return numberFake;
    27762744}
     2745/* Pivot out a variable and choose an incoing one.  Assumes dual
     2746   feasible - will not go through a reduced cost. 
     2747   Returns step length in theta
     2748   Returns ray in ray_ (or NULL if no pivot)
     2749   Return codes as before but -1 means no acceptable pivot
     2750*/
     2751int
     2752ClpSimplexDual::pivotResult()
     2753{
     2754  abort();
     2755  return 0;
     2756}
  • trunk/ClpSimplexPrimal.cpp

    r105 r109  
    180180  */
    181181
    182   // sanity check
    183   assert (numberRows_==matrix_->getNumRows());
    184   assert (numberColumns_==matrix_->getNumCols());
    185   // for moment all arrays must exist
    186   assert(columnLower_);
    187   assert(columnUpper_);
    188   assert(rowLower_);
    189   assert(rowUpper_);
    190 
    191 
    192182  algorithm_ = +1;
    193   primalTolerance_=dblParam_[ClpPrimalTolerance];
    194   dualTolerance_=dblParam_[ClpDualTolerance];
    195 
    196   // put in standard form (and make row copy)
    197   // create modifiable copies of model rim and do optional scaling
    198   bool goodMatrix=createRim(7+8+16,true);
    199 
    200   // save infeasibility cost
    201   double saveInfeasibilityCost = infeasibilityCost_;
    202   // Save perturbation
    203   int savePerturbation = perturbation_;
    204   // save if sparse factorization wanted
    205   int saveSparse = factorization_->sparseThreshold();
    206 
    207   if (goodMatrix) {
    208     // Model looks okay
    209     int iRow,iColumn;
    210     // Do initial factorization
    211     // and set certain stuff
    212     // We can either set increasing rows so ...IsBasic gives pivot row
    213     // or we can just increment iBasic one by one
    214     // for now let ...iBasic give pivot row
    215     factorization_->increasingRows(2);
    216     // row activities have negative sign
    217     factorization_->slackValue(-1.0);
    218     factorization_->zeroTolerance(1.0e-13);
    219    
    220    
    221    
    222     // If user asked for perturbation - do it
    223    
    224     if (perturbation_<100)
    225       perturb();
    226 
    227     // for primal we will change bounds using infeasibilityCost_
    228     if (nonLinearCost_==NULL) {
    229       // get a valid nonlinear cost function
    230       delete nonLinearCost_;
    231       nonLinearCost_= new ClpNonLinearCost(this);
    232     }
    233    
    234     // loop round to clean up solution if values pass
    235     int numberThrownOut = -1;
    236     int firstSuperBasic=numberRows_+numberColumns_;
    237     int totalNumberThrownOut=0;
    238     while(numberThrownOut) {
    239       int status = internalFactorize(0+10*ifValuesPass);
    240       if (status<0)
    241         return 1; // some error
    242       else
    243         totalNumberThrownOut+= status;
    244      
    245       // for this we need clean basis so it is after factorize
    246       numberThrownOut=gutsOfSolution(rowActivityWork_,columnActivityWork_,
    247                                      ifValuesPass!=0);
    248       totalNumberThrownOut+= numberThrownOut;
    249      
    250       // find first superbasic - columns, then rows
    251       if (ifValuesPass) {
    252         nextSuperBasic(firstSuperBasic);
    253         if (firstSuperBasic==numberRows_+numberColumns_)
    254           ifValuesPass=0; // signal no values pass
    255       }
    256     }
    257    
    258     if (totalNumberThrownOut)
    259       handler_->message(CLP_SINGULARITIES,messages_)
    260         <<totalNumberThrownOut
    261         <<CoinMessageEol;
    262    
    263     problemStatus_ = -1;
    264     numberIterations_=0;
    265    
     183
     184  // save data
     185  ClpDataSave data = saveData();
     186
     187
     188    // initialize - maybe values pass and algorithm_ is +1
     189    if (!startup(ifValuesPass)) {
     190
    266191    int lastCleaned=0; // last time objective or bounds cleaned up
    267    
    268     // number of times we have declared optimality
    269     numberTimesOptimal_=0;
    270192   
    271193    // Progress indicator
     
    291213    */
    292214    while (problemStatus_<0) {
     215      int iRow,iColumn;
    293216      // clear
    294217      for (iRow=0;iRow<4;iRow++) {
     
    317240      // Say good factorization
    318241      factorType=1;
    319       if (saveSparse) {
     242      if (data.sparseThreshold_) {
    320243        // use default at present
    321244        factorization_->sparseThreshold(0);
     
    330253     
    331254      // Iterate
    332       whileIterating(firstSuperBasic);
     255      whileIterating();
    333256    }
    334257  }
     
    343266    computeDuals();
    344267  }
    345   factorization_->sparseThreshold(saveSparse);
    346   // Get rid of some arrays and empty factorization
    347   deleteRim();
    348   handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
    349     <<objectiveValue()
    350     <<CoinMessageEol;
    351   // Restore any saved stuff
    352   perturbation_ = savePerturbation;
    353   infeasibilityCost_ = saveInfeasibilityCost;
     268  // clean up
     269  finish();
     270  restoreData(data);
    354271  return problemStatus_;
    355272}
     
    365282*/
    366283int
    367 ClpSimplexPrimal::whileIterating(int & firstSuperBasic)
     284ClpSimplexPrimal::whileIterating()
    368285{
    369286
    370287  // Say if values pass
    371   int ifValuesPass=0;
     288  int ifValuesPass=(firstFree_>=0) ? 1 : 0;
    372289  int returnCode=-1;
    373290  int startIteration = numberIterations_;
    374   if (firstSuperBasic<numberRows_+numberColumns_)
    375     ifValuesPass=1;
    376291  int saveNumber = numberIterations_;
    377292  // status stays at -1 while iterating, >=0 finished, -2 to invert
     
    431346      // in values pass
    432347      if (ifValuesPass>0) {
    433         nextSuperBasic(firstSuperBasic);
    434         if (firstSuperBasic==numberRows_+numberColumns_)
     348        int sequenceIn=nextSuperBasic();
     349        if (sequenceIn<0) {
    435350          ifValuesPass=-1; // signal end of values pass after this
     351        } else {
     352          // normal
     353          sequenceIn_ = sequenceIn;
     354          valueIn_=solution_[sequenceIn_];
     355          lowerIn_=lower_[sequenceIn_];
     356          upperIn_=upper_[sequenceIn_];
     357          dualIn_=dj_[sequenceIn_];
     358        }
    436359      } else {
    437360        // end of values pass - initialize weights etc
     
    625548      if (pivotRow_>=0)
    626549        dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
    627      
     550      double oldValue = valueIn_;
     551      if (directionIn_==-1) {
     552        // as if from upper bound
     553        if (sequenceIn_!=sequenceOut_) {
     554          // variable becoming basic
     555          //setStatus(sequenceIn_,basic);
     556          valueIn_ -= fabs(theta_);
     557        } else {
     558          valueIn_=lowerIn_;
     559          //setStatus(sequenceIn_, atLowerBound);
     560        }
     561      } else {
     562        // as if from lower bound
     563        if (sequenceIn_!=sequenceOut_) {
     564          // variable becoming basic
     565          //setStatus(sequenceIn_,basic);
     566          valueIn_ += fabs(theta_);
     567        } else {
     568          valueIn_=upperIn_;
     569          //setStatus(sequenceIn_, atUpperBound);
     570        }
     571      }
     572      objectiveChange += dualIn_*(valueIn_-oldValue);
     573      // outgoing
     574      if (sequenceIn_!=sequenceOut_) {
     575        if (directionOut_>0) {
     576          valueOut_ = lowerOut_;
     577        } else {
     578          valueOut_ = upperOut_;
     579        }
     580        double lowerValue = lower_[sequenceOut_];
     581        double upperValue = upper_[sequenceOut_];
     582        assert(valueOut_>=lowerValue-primalTolerance_&&
     583               valueOut_<=upperValue+primalTolerance_);
     584        // may not be exactly at bound and bounds may have changed
     585        if (valueOut_<=lowerValue+primalTolerance_) {
     586          //setStatus(sequenceOut_,atLowerBound);
     587          directionOut_=1;
     588        } else if (valueOut_>=upperValue-primalTolerance_) {
     589          //setStatus(sequenceOut_,atUpperBound);
     590          directionOut_=-1;
     591        } else {
     592          printf("*** variable wandered off bound %g %g %g!\n",
     593                 lowerValue,valueOut_,upperValue);
     594          if (valueOut_-lowerValue<=upperValue-valueOut_) {
     595            //setStatus(sequenceOut_,atLowerBound);
     596            valueOut_=lowerValue;
     597            directionOut_=1;
     598          } else {
     599            //setStatus(sequenceOut_,atUpperBound);
     600            valueOut_=upperValue;
     601            directionOut_=-1;
     602          }
     603        }
     604        solution_[sequenceOut_]=valueOut_;
     605        nonLinearCost_->setOne(sequenceOut_,valueOut_);
     606      }
     607      // change cost and bounds on incoming if primal
     608      nonLinearCost_->setOne(sequenceIn_,valueIn_);
    628609      int whatNext=housekeeping(objectiveChange);
    629610     
     
    14021383  return 0;
    14031384}
    1404 void
    1405 ClpSimplexPrimal::nextSuperBasic(int & firstSuperBasic)
    1406 {
    1407   int iColumn;
    1408   if (firstSuperBasic==numberRows_+numberColumns_) {
    1409     // initialization
    1410     iColumn=0;
    1411   } else {
    1412     // normal
    1413     sequenceIn_=firstSuperBasic;
    1414     valueIn_=solution_[sequenceIn_];
    1415     lowerIn_=lower_[sequenceIn_];
    1416     upperIn_=upper_[sequenceIn_];
    1417     dualIn_=dj_[sequenceIn_];
    1418     iColumn=firstSuperBasic+1;
    1419   }
    1420   for (;iColumn<numberRows_+numberColumns_;iColumn++) {
    1421     if (getStatus(iColumn)==superBasic) {
    1422       // is it really super basic
    1423       if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
    1424         solution_[iColumn]=lower_[iColumn];
    1425         setStatus(iColumn,atLowerBound);
    1426       } else if (fabs(solution_[iColumn]-upper_[iColumn])
    1427                  <=primalTolerance_) {
    1428         solution_[iColumn]=upper_[iColumn];
    1429         setStatus(iColumn,atUpperBound);
    1430       } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
    1431         setStatus(iColumn,isFree);
    1432       } else {
    1433         break;
    1434       }
    1435     }
    1436   }
    1437   firstSuperBasic = iColumn;
    1438 }
    14391385// Perturbs problem
    14401386void
     
    15301476    specialOptions_ &= ~1;
    15311477}
    1532 
     1478/* Pivot in a variable and choose an outgoing one.  Assumes primal
     1479   feasible - will not go through a bound.  Returns step length in theta
     1480   Returns ray in ray_ (or NULL if no pivot)
     1481   Return codes as before but -1 means no acceptable pivot
     1482*/
     1483int
     1484ClpSimplexPrimal::pivotResult()
     1485{
     1486  abort();
     1487  return 0;
     1488}
     1489 
     1490
  • trunk/Makefile.Clp

    r95 r109  
    6363ifeq ($(OptLevel),-O2)
    6464#     CXXFLAGS += -DNDEBUG
     65#    CXXFLAGS += -DPRESOLVE_SUMMARY=1 -DDEBUG_PRESOLVE -DCHECK_CONSISTENCY=1
    6566endif
    6667
  • trunk/Presolve.cpp

    r102 r109  
    427427  originalModel_=model;
    428428}
    429 
     429#if 0
    430430// A lazy way to restrict which transformations are applied
    431431// during debugging.
     
    448448#endif
    449449}
     450#endif
    450451//#define DEBUG_PRESOLVE 1
    451452#if DEBUG_PRESOLVE
     
    521522
    522523  paction_ = make_fixed(prob, paction_);
     524  // if integers then switch off dual stuff
     525  // later just do individually
     526  bool doDualStuff = (presolvedModel_->integerInformation()!=NULL);
    523527
    524528#if     CHECK_CONSISTENCY
     
    527531
    528532  if (!prob->status_) {
     533#if 0
    529534    const bool slackd = ATOI("SLACKD")!=0;
    530535    //const bool forcing = ATOI("FORCING")!=0;
     
    536541    const bool duprow = ATOI("off")!=0;
    537542    const bool dual = ATOI("off")!=0;
     543#else
     544    const bool slackd = doDualStuff;
     545    const bool doubleton = true;
     546    const bool forcing = true;
     547    const bool ifree = true;
     548    const bool zerocost = true;
     549    const bool dupcol = true;
     550    const bool duprow = true;
     551    const bool dual = doDualStuff;
     552#endif
    538553   
    539554    // some things are expensive so just do once (normally)
  • trunk/PresolveMatrix.cpp

    r60 r109  
    22// Corporation and others.  All Rights Reserved.
    33
    4 #define CHECK_CONSISTENCY       1
     4//#define       CHECK_CONSISTENCY       1
    55
    66#include <stdio.h>
     
    545545                                     CoinBigIndex nelems0)
    546546{
     547#if 0
     548  // out as it was hardware error !!
     549  {
     550    // temporary sanity check (for funny nw04 bug)
     551    int icol;
     552    for (icol=0;icol<ncols_;icol++) {
     553      int j;
     554      for (j=mcstrt_[icol];j<mcstrt_[icol]+hincol_[icol];j++) {
     555        assert(hrow_[j]>=0&&hrow_[j]<nrows_);
     556      }
     557    }
     558  }
     559#endif
    547560  si.loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
    548561                 clo_, cup_, cost_, rlo_, rup_);
     
    569582}
    570583
     584#if     CHECK_CONSISTENCY
    571585
    572586// The matrix is represented redundantly in both row and column format,
     
    591605                              const char *ROW, const char *COL)
    592606{
    593 #if     CHECK_CONSISTENCY
    594607  for (int irow=0; irow<nrows; irow++) {
    595608    if (hinrow[irow] > 0) {
     
    618631    }
    619632  }
    620 #endif
    621 }
     633}
     634#endif
    622635
    623636
  • trunk/PresolveSubst.cpp

    r95 r109  
    10241024#if     PRESOLVE_SUMMARY
    10251025    printf("NSUBSTS:  %d\n", nactions);
    1026     printf("NT: %d  NGOOD:  %d FILL_LEVEL:  %d\n", nt, ngood, fill_level);
     1026    // printf("NT: %d  NGOOD:  %d FILL_LEVEL:  %d\n", nt, ngood, fill_level);
    10271027#endif
    10281028    next = new subst_constraint_action(nactions, copyOfArray(actions,nactions), next);
  • trunk/include/ClpModel.hpp

    r82 r109  
    145145  inline void setNumberIterations(int numberIterations)
    146146  { numberIterations_ = numberIterations;};
     147  /// Solve type - 1 simplex, 2 simplex interface
     148  inline int solveType() const
     149  { return solveType_;};
     150  inline void setSolveType(int type)
     151  { solveType_=type;};
    147152   /// Maximum number of iterations
    148153   inline int maximumIterations() const { return intParam_[ClpMaxNumIteration]; }
     
    405410  /// Number of iterations
    406411  int numberIterations_;
     412  /// Solve type - 1 simplex, 2 simplex interface
     413  int solveType_;
    407414  /// Status of problem
    408415  int problemStatus_;
     
    431438  //@}
    432439};
     440/** This is a tiny class where data can be saved round calls */
     441class ClpDataSave {
     442
     443public:
     444  /**@name Constructors and destructor
     445   */
     446  //@{
     447    /// Default constructor
     448    ClpDataSave (  );
     449
     450    /// Copy constructor.
     451    ClpDataSave(const ClpDataSave &);
     452    /// Assignment operator. This copies the data
     453    ClpDataSave & operator=(const ClpDataSave & rhs);
     454    /// Destructor
     455    ~ClpDataSave (  );
     456
     457  //@}
     458
     459////////////////// data //////////////////
     460public:
     461
     462  /**@name data - with same names as in other classes*/
     463  //@{
     464  double dualBound_;
     465  double infeasibilityCost_;
     466  int sparseThreshold_;
     467  int perturbation_;
     468
     469  //@}
     470};
     471
    433472#endif
  • trunk/include/ClpSimplex.hpp

    r104 r109  
    177177  //@}
    178178
     179  /**@name Needed for functionality of OsiSimplexInterface */
     180  //@{
     181  /** Pivot in a variable and out a variable.  Returns 0 if okay,
     182      1 if inaccuracy forced re-factorization, -1 if would be singular.
     183      Also updates primal/dual infeasibilities.
     184      Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
     185  */
     186  int pivot();
     187
     188  /** Pivot in a variable and choose an outgoing one.  Assumes primal
     189      feasible - will not go through a bound.  Returns step length in theta
     190      Returns ray in ray_ (or NULL if no pivot)
     191      Return codes as before but -1 means no acceptable pivot
     192  */
     193  int primalPivotResult();
     194 
     195  /** Pivot out a variable and choose an incoing one.  Assumes dual
     196      feasible - will not go through a reduced cost. 
     197      Returns step length in theta
     198      Returns ray in ray_ (or NULL if no pivot)
     199      Return codes as before but -1 means no acceptable pivot
     200  */
     201  int dualPivotResult();
     202
     203  /** Common bits of coding for dual and primal.  Return s0 if okay,
     204      1 if bad matrix, 2 if very bad factorization
     205  */
     206  int startup(int ifValuesPass);
     207  void finish();
     208 
     209  //@}
     210
    179211  /**@name most useful gets and sets */
    180212  //@{
     
    246278      infeasibilities etc */
    247279  void checkSolution();
     280  /// Useful row length arrays (0,1,2,3,4,5)
     281  inline CoinIndexedVector * rowArray(int index) const
     282  { return rowArray_[index];};
     283  /// Useful column length arrays (0,1,2,3,4,5)
     284  inline CoinIndexedVector * columnArray(int index) const
     285  { return columnArray_[index];};
    248286  //@}
    249287
     
    266304      when total number of singularities is returned
    267305  */
     306    /// Save data
     307    ClpDataSave saveData() const;
     308    /// Restore data
     309    void restoreData(ClpDataSave saved);
    268310  int internalFactorize(int solveType);
    269311  /// Factorizes using current basis. For external use
  • trunk/include/ClpSimplexDual.hpp

    r50 r109  
    222222      so can exit gracefully before end */
    223223  int numberAtFakeBound();
     224
     225  /** Pivot in a variable and choose an outgoing one.  Assumes dual
     226      feasible - will not go through a reduced cost.  Returns step length in theta
     227      Returns ray in ray_ (or NULL if no pivot)
     228      Return codes as before but -1 means no acceptable pivot
     229  */
     230  int pivotResult();
     231 
    224232  //@}
    225233};
  • trunk/include/ClpSimplexPrimal.hpp

    r54 r109  
    123123  //@{
    124124  /** This has the flow between re-factorizations
    125       firstSuperBasic == number rows + columns normally,
    126       otherwise first super basic variable
    127125
    128126      Reasons to come out:
     
    135133      +3 max iterations
    136134   */
    137   int whileIterating(int & firstSuperBasic);
     135  int whileIterating();
    138136  /** The primals are updated by the given array.
    139137      Returns number of infeasibilities.
     
    191189  /// Unflag all variables and return number unflagged
    192190  int unflag();
    193   /// Sets sequenceIn_ to next superBasic (input by first..) and updates
    194   void nextSuperBasic(int & firstSuperBasic);
     191
     192  /** Pivot in a variable and choose an outgoing one.  Assumes primal
     193      feasible - will not go through a bound.  Returns step length in theta
     194      Returns ray in ray_ (or NULL if no pivot)
     195      Return codes as before but -1 means no acceptable pivot
     196  */
     197  int pivotResult();
     198 
    195199  //@}
    196200};
Note: See TracChangeset for help on using the changeset viewer.