Changeset 111


Ignore:
Timestamp:
Feb 5, 2003 11:49:19 AM (17 years ago)
Author:
forrest
Message:

Changes to make OsiSimplex? stuff easier

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpSimplex.cpp

    r109 r111  
    6969  directionOut_(-1),
    7070  pivotRow_(-1),
     71  lastGoodIteration_(-100),
    7172  dj_(NULL),
    7273  rowReducedCost_(NULL),
     
    722723            // but put to bound if close
    723724            if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow])
    724                 <=primalTolerance_)
     725                <=primalTolerance_) {
    725726              rowActivityWork_[iRow]=rowLowerWork_[iRow];
    726             else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
    727                 <=primalTolerance_)
     727              setRowStatus(iRow,atLowerBound);
     728            } else if (fabs(rowActivityWork_[iRow]-rowUpperWork_[iRow])
     729                       <=primalTolerance_) {
    728730              rowActivityWork_[iRow]=rowUpperWork_[iRow];
     731              setRowStatus(iRow,atUpperBound);
     732            }
    729733          }
    730734         
     
    738742              // but put to bound if close
    739743              if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    740                   <=primalTolerance_)
     744                  <=primalTolerance_) {
    741745                columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    742               else if (fabs(columnActivityWork_[iColumn]
     746                setColumnStatus(iColumn,atLowerBound);
     747              } else if (fabs(columnActivityWork_[iColumn]
    743748                            -columnUpperWork_[iColumn])
    744                        <=primalTolerance_)
     749                         <=primalTolerance_) {
    745750                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     751                setColumnStatus(iColumn,atUpperBound);
     752              }
    746753            } else
    747754              numberBasic++;
     
    750757            // but put to bound if close
    751758            if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    752                 <=primalTolerance_)
     759                <=primalTolerance_) {
    753760              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    754             else if (fabs(columnActivityWork_[iColumn]
     761              setColumnStatus(iColumn,atLowerBound);
     762            } else if (fabs(columnActivityWork_[iColumn]
    755763                          -columnUpperWork_[iColumn])
    756                      <=primalTolerance_)
     764                       <=primalTolerance_) {
    757765              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     766              setColumnStatus(iColumn,atUpperBound);
     767            }
    758768          }
    759769        }
     
    772782          // but put to bound if close
    773783          if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn])
    774               <=primalTolerance_)
     784              <=primalTolerance_) {
    775785            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    776           else if (fabs(columnActivityWork_[iColumn]
     786            setColumnStatus(iColumn,atLowerBound);
     787          } else if (fabs(columnActivityWork_[iColumn]
    777788                        -columnUpperWork_[iColumn])
    778                    <=primalTolerance_)
     789                     <=primalTolerance_) {
    779790            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
     791            setColumnStatus(iColumn,atUpperBound);
     792          }
    780793        }
    781794      }
     
    10451058  directionOut_(-1),
    10461059  pivotRow_(-1),
     1060  lastGoodIteration_(-100),
    10471061  dj_(NULL),
    10481062  rowReducedCost_(NULL),
     
    11351149  directionOut_(-1),
    11361150  pivotRow_(-1),
     1151  lastGoodIteration_(-100),
    11371152  dj_(NULL),
    11381153  rowReducedCost_(NULL),
     
    12731288  directionOut_ = rhs.directionOut_;
    12741289  pivotRow_ = rhs.pivotRow_;
     1290  lastGoodIteration_ = rhs.lastGoodIteration_;
    12751291  numberRefinements_ = rhs.numberRefinements_;
    12761292  dualTolerance_ = rhs.dualTolerance_;
     
    19902006   This is to make dual go faster and is probably not needed
    19912007   with a presolve.  Returns non-zero if problem infeasible
     2008
     2009   Fudge for branch and bound - put bounds on columns of factor *
     2010   largest value (at continuous) - should improve stability
     2011   in branch and bound on infeasible branches (0.0 is off)
    19922012*/
    19932013int
    1994 ClpSimplex::tightenPrimalBounds()
     2014ClpSimplex::tightenPrimalBounds(double factor)
    19952015{
    19962016 
     
    20162036  double * saveUpper = new double [numberColumns_];
    20172037  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
     2038
     2039  int iRow, iColumn;
     2040
     2041  // If wanted - tighten column bounds using solution
     2042  if (factor) {
     2043    assert (factor>1.0);
     2044    double largest=0.0;
     2045    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     2046      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
     2047        largest = max(largest,fabs(columnActivity_[iColumn]));
     2048      }
     2049    }
     2050    largest *= factor;
     2051    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     2052      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
     2053        columnUpper_[iColumn] = min(columnUpper_[iColumn],largest);
     2054        columnLower_[iColumn] = max(columnLower_[iColumn],-largest);
     2055      }
     2056    }
     2057  }
    20182058#define MAXPASS 10 
    20192059
     
    20212061  // Would be faster to have a stack of possible rows
    20222062  // and we put altered rows back on stack
    2023 
    2024   int iRow, iColumn;
    20252063
    20262064  while(numberChanged) {
     
    36533691    assert (pivotRow_<0);
    36543692  }
    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_;
     3693  bool roundAgain = true;
     3694  int returnCode=0;
     3695  while (roundAgain) {
     3696    roundAgain=false;
     3697    unpack(rowArray_[1]);
     3698    factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
     3699    // we are going to subtract movement from current basic
     3700    double movement;
     3701    // see where incoming will go to
     3702    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
     3703      // flip so go to bound
     3704      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
     3705    } else {
     3706      // get where outgoing needs to get to
     3707      double outValue = (directionOut_>0) ? upperOut_ : lowerOut_;
     3708      // solutionOut_ - movement*alpha_ == outValue
     3709      movement = (outValue-valueOut_)/alpha_;
     3710      // set directionIn_ correctly
     3711      directionIn_ = (movement>0) ? 1 :-1;
     3712    }
     3713    // update primal solution
     3714    {
     3715      int i;
     3716      int * index = rowArray_[1]->getIndices();
     3717      int number = rowArray_[1]->getNumElements();
     3718      double * element = rowArray_[1]->denseVector();
     3719      for (i=0;i<number;i++) {
     3720        int ii = index[i];
     3721        // get column
     3722        ii = pivotVariable_[ii];
     3723        solution_[ii] -= movement*element[i];
     3724        element[i]=0.0;
     3725      }
     3726      // see where something went to
     3727      if (sequenceOut_<0) {
     3728        if (directionIn_<0) {
     3729          assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
     3730          solution_[sequenceIn_]=upperIn_;
     3731        } else {
     3732          assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
     3733          solution_[sequenceIn_]=lowerIn_;
     3734        }
     3735      } else {
     3736        if (directionOut_<0) {
     3737          assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
     3738          solution_[sequenceOut_]=upperOut_;
     3739        } else {
     3740          assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
     3741          solution_[sequenceOut_]=lowerOut_;
     3742        }
     3743        solution_[sequenceIn_]=valueIn_+movement;
     3744      }
     3745    }   
     3746    double objectiveChange = dualIn_*movement;
     3747    // update duals
     3748    if (pivotRow_>=0) {
     3749      alpha_ = rowArray_[1]->denseVector()[pivotRow_];
     3750      assert (fabs(alpha_)>1.0e-8);
     3751      double multiplier = dualIn_/alpha_;
     3752      rowArray_[0]->insert(pivotRow_,multiplier);
     3753      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
     3754      // put row of tableau in rowArray[0] and columnArray[0]
     3755      matrix_->transposeTimes(this,-1.0,
     3756                              rowArray_[0],columnArray_[1],columnArray_[0]);
     3757      // update column djs
     3758      int i;
     3759      int * index = columnArray_[0]->getIndices();
     3760      int number = columnArray_[0]->getNumElements();
     3761      double * element = columnArray_[0]->denseVector();
     3762      for (i=0;i<number;i++) {
     3763        int ii = index[i];
     3764        dj_[ii] -= element[i];
     3765        element[i]=0.0;
     3766      }
     3767      columnArray_[0]->setNumElements(0);
     3768      // and row djs
     3769      index = rowArray_[0]->getIndices();
     3770      number = rowArray_[0]->getNumElements();
     3771      element = rowArray_[0]->denseVector();
     3772      for (i=0;i<number;i++) {
     3773        int ii = index[i];
     3774        dj_[ii+numberRows_] -= element[i];
     3775        element[i]=0.0;
     3776      }
     3777      rowArray_[0]->setNumElements(0);
     3778      // check incoming
     3779      assert (fabs(dj_[sequenceIn_])<1.0e-6);
     3780    }
     3781   
     3782    // if stable replace in basis
     3783    int updateStatus = factorization_->replaceColumn(rowArray_[2],
     3784                                                   pivotRow_,
     3785                                                     alpha_);
     3786    bool takePivot=true;
     3787    // if no pivots, bad update but reasonable alpha - take and invert
     3788    if (updateStatus==2&&
     3789        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
     3790      updateStatus=4;
     3791    if (updateStatus==1||updateStatus==4) {
     3792      // slight error
     3793      if (factorization_->pivots()>5||updateStatus==4) {
     3794        returnCode=-1;
     3795      }
     3796    } else if (updateStatus==2) {
     3797      // major error
     3798      rowArray_[1]->clear();
     3799      takePivot=false;
     3800      if (factorization_->pivots()) {
     3801        // refactorize here
     3802        statusOfProblem();
     3803        roundAgain=true;
     3804      } else {
     3805        returnCode=1;
     3806      }
     3807    } else if (updateStatus==3) {
     3808      // out of memory
     3809      // increase space if not many iterations
     3810      if (factorization_->pivots()<
     3811          0.5*factorization_->maximumPivots()&&
     3812          factorization_->pivots()<200)
     3813        factorization_->areaFactor(
     3814                                   factorization_->areaFactor() * 1.1);
     3815      returnCode =-1; // factorize now
     3816    }
     3817    if (takePivot) {
     3818      int save = algorithm_;
     3819      // make simple so always primal
     3820      algorithm_=1;
     3821      housekeeping(objectiveChange);
     3822      algorithm_=save;
     3823    }
     3824  }
     3825  if (returnCode == -1) {
     3826    // refactorize here
     3827    statusOfProblem();
    36633828  } 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;
     3829    // just for now - recompute anyway
     3830    gutsOfSolution(rowActivityWork_,columnActivityWork_);
     3831  }
     3832  return returnCode;
    37503833}
    37513834
     
    38873970  dualBound_ = saved.dualBound_;
    38883971}
     3972/* Factorizes and returns true if optimal.  Used by user */
     3973bool
     3974ClpSimplex::statusOfProblem()
     3975{
     3976  // is factorization okay?
     3977  assert (internalFactorize(1)==0);
     3978  // put back original costs and then check
     3979  createRim(4);
     3980  gutsOfSolution(rowActivityWork_, columnActivityWork_);
     3981  return (primalFeasible()&&dualFeasible());
     3982}
  • trunk/ClpSimplexDual.cpp

    r109 r111  
    21292129    problemStatus_=4; // unknown
    21302130  }
     2131  lastGoodIteration_ = numberIterations_;
    21312132
    21322133}
  • trunk/ClpSimplexPrimal.cpp

    r109 r111  
    199199    // This says whether to restore things etc
    200200    int factorType=0;
    201     // Save iteration number
    202     int saveNumber = -1;
    203201    /*
    204202      Status of problem:
     
    233231#endif
    234232      // If we have done no iterations - special
    235       if (saveNumber==numberIterations_)
     233      if (lastGoodIteration_==numberIterations_)
    236234        factorType=3;
    237235      // may factorize, checks if problem finished
     
    240238      // Say good factorization
    241239      factorType=1;
    242       if (data.sparseThreshold_) {
    243         // use default at present
    244         factorization_->sparseThreshold(0);
    245         factorization_->goSparse();
    246       }
    247240     
    248241      // Say no pivot has occurred (for steepest edge and updates)
    249242      pivotRow_=-2;
    250      
    251       // Save iteration number
    252       saveNumber = numberIterations_;
    253243     
    254244      // Iterate
     
    288278  int ifValuesPass=(firstFree_>=0) ? 1 : 0;
    289279  int returnCode=-1;
    290   int startIteration = numberIterations_;
    291   int saveNumber = numberIterations_;
    292280  // status stays at -1 while iterating, >=0 finished, -2 to invert
    293281  // status -3 to go to top without an invert
     
    361349        primalColumnPivot_->saveWeights(this,5);
    362350        ifValuesPass=0;
    363         if(saveNumber != numberIterations_) {
     351        if(lastGoodIteration_ != numberIterations_) {
    364352          problemStatus_=-2; // factorize now
    365353          pivotRow_=-1; // say no weights update
     
    385373      }
    386374#endif
    387       // update the incoming column
    388       unpack(rowArray_[1]);
    389       // save reduced cost
    390       double saveDj = dualIn_;
    391       factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
    392       // do ratio test and re-compute dj
    393       primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
    394                 ifValuesPass);
    395       if (ifValuesPass) {
    396         saveDj=dualIn_;
    397         if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
    398           if(fabs(dualIn_)<1.0e2*dualTolerance_) {
    399             // try other way
    400             directionIn_=-directionIn_;
    401             primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
    402                       0);
    403           }
    404           if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
    405             // reject it
    406             char x = isColumn(sequenceIn_) ? 'C' :'R';
    407             handler_->message(CLP_SIMPLEX_FLAG,messages_)
    408               <<x<<sequenceWithin(sequenceIn_)
    409               <<CoinMessageEol;
    410             setFlagged(sequenceIn_);
    411             lastBadIteration_ = numberIterations_; // say be more cautious
    412             rowArray_[1]->clear();
    413             pivotRow_=-1;
    414             continue;
    415           }
    416         }
    417       }
    418      
    419 #ifdef CLP_DEBUG
    420       if ((handler_->logLevel()&32))
    421         printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_);
    422 #endif
    423       if ((saveDj*dualIn_<1.0e-20&&!ifValuesPass)||
    424           fabs(saveDj-dualIn_)>1.0e-3*(1.0+fabs(saveDj))) {
    425         handler_->message(CLP_PRIMAL_DJ,messages_)
    426           <<saveDj<<dualIn_
    427           <<CoinMessageEol;
    428         if(saveNumber != numberIterations_) {
    429           problemStatus_=-2; // factorize now
    430           rowArray_[1]->clear();
    431           pivotRow_=-1; // say no weights update
    432           returnCode=-2;
    433           if(saveNumber+1 == numberIterations_) {
    434             // not looking wonderful - try cleaning bounds
    435             // put non-basics to bounds in case tolerance moved
    436             nonLinearCost_->checkInfeasibilities(true);
    437           }
    438           break;
    439         } else {
    440           // take on more relaxed criterion
    441           if (saveDj*dualIn_<1.0e-20||
    442               fabs(saveDj-dualIn_)>1.0e-1*(1.0+fabs(dualIn_))) {
    443             // need to reject something
    444             char x = isColumn(sequenceIn_) ? 'C' :'R';
    445             handler_->message(CLP_SIMPLEX_FLAG,messages_)
    446               <<x<<sequenceWithin(sequenceIn_)
    447               <<CoinMessageEol;
    448             setFlagged(sequenceIn_);
    449             lastBadIteration_ = numberIterations_; // say be more cautious
    450             rowArray_[1]->clear();
    451             pivotRow_=-1;
    452             continue;
    453           }
    454         }
    455       }
    456       if (pivotRow_>=0) {
    457         // if stable replace in basis
    458         int updateStatus = factorization_->replaceColumn(rowArray_[2],
    459                                                          pivotRow_,
    460                                                          alpha_);
    461         // if no pivots, bad update but reasonable alpha - take and invert
    462         if (updateStatus==2&&
    463                    saveNumber==numberIterations_&&fabs(alpha_)>1.0e-5)
    464           updateStatus=4;
    465         if (updateStatus==1||updateStatus==4) {
    466           // slight error
    467           if (factorization_->pivots()>5||updateStatus==4) {
    468             problemStatus_=-2; // factorize now
    469             returnCode=-3;
    470           }
    471         } else if (updateStatus==2) {
    472           // major error
    473           // later we may need to unwind more e.g. fake bounds
    474           if(saveNumber != numberIterations_) {
    475             problemStatus_=-2; // factorize now
    476             returnCode=-2;
    477             break;
    478           } else {
    479             // need to reject something
    480             char x = isColumn(sequenceIn_) ? 'C' :'R';
    481             handler_->message(CLP_SIMPLEX_FLAG,messages_)
    482               <<x<<sequenceWithin(sequenceIn_)
    483               <<CoinMessageEol;
    484             setFlagged(sequenceIn_);
    485             lastBadIteration_ = numberIterations_; // say be more cautious
    486             rowArray_[1]->clear();
    487             pivotRow_=-1;
    488             continue;
    489           }
    490         } else if (updateStatus==3) {
    491           // out of memory
    492           // increase space if not many iterations
    493           if (factorization_->pivots()<
    494               0.5*factorization_->maximumPivots()&&
    495               factorization_->pivots()<200)
    496             factorization_->areaFactor(
    497                                        factorization_->areaFactor() * 1.1);
    498           problemStatus_=-2; // factorize now
    499         }
    500         // here do part of steepest - ready for next iteration
    501         primalColumnPivot_->updateWeights(rowArray_[1]);
    502       } else {
    503         if (pivotRow_==-1) {
    504           // no outgoing row is valid
    505 #ifdef CLP_DEBUG
    506           if (handler_->logLevel()&32)
    507             printf("** no row pivot\n");
    508 #endif
    509           if (!factorization_->pivots()) {
    510             problemStatus_=-5; //say looks unbounded
    511             // do ray
    512             delete [] ray_;
    513             ray_ = new double [numberColumns_];
    514             ClpFillN(ray_,numberColumns_,0.0);
    515             int number=rowArray_[1]->getNumElements();
    516             int * index = rowArray_[1]->getIndices();
    517             double * array = rowArray_[1]->denseVector();
    518             double way=-directionIn_;
    519             int i;
    520             double zeroTolerance=1.0e-12;
    521             if (sequenceIn_<numberColumns_)
    522               ray_[sequenceIn_]=directionIn_;
    523             for (i=0;i<number;i++) {
    524               int iRow=index[i];
    525               int iPivot=pivotVariable_[iRow];
    526               double arrayValue = array[iRow];
    527               if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
    528                 ray_[iPivot] = way* array[iRow];
    529             }
    530           }
    531           rowArray_[0]->clear();
    532           returnCode=2;
    533           break;
    534         } else {
    535           // flipping from bound to bound
    536         }
    537       }
    538      
    539       // update primal solution
    540      
    541       double objectiveChange=0.0;
    542       // Cost on pivot row may change - may need to change dualIn
    543       double oldCost=0.0;
    544       if (pivotRow_>=0)
    545         oldCost = cost(pivotVariable_[pivotRow_]);
    546       // rowArray_[1] is not empty - used to update djs
    547       updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
    548       if (pivotRow_>=0)
    549         dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
    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_);
    609       int whatNext=housekeeping(objectiveChange);
    610      
    611       if (whatNext==1) {
    612         problemStatus_ =-2; // refactorize
    613       } else if (whatNext==2) {
    614         // maximum iterations or equivalent
    615         problemStatus_= 3;
    616         returnCode=3;
    617         break;
    618       } else if(numberIterations_ == startIteration
    619                 + 2 * factorization_->maximumPivots()) {
    620         // done a lot of flips - be safe
    621         problemStatus_ =-2; // refactorize
     375      // do second half of iteration
     376      returnCode = pivotResult();
     377      if (returnCode<-1&&returnCode>-5) {
     378        problemStatus_=-2; //
     379      } else if (returnCode==-5) {
     380        // something flagged - continue;
     381      } else if (returnCode==2) {
     382        problemStatus_=-5; // looks unbounded
     383      } else if (returnCode==4) {
     384        problemStatus_=-2; // looks unbounded but has iterated
     385      } else if (returnCode!=-1) {
     386        assert(returnCode==3);
     387        problemStatus_=3;
    622388      }
    623389    } else {
     
    651417    changeMade_++; // say change made
    652418  }
     419  int saveThreshold = factorization_->sparseThreshold();
    653420  int tentativeStatus = problemStatus_;
    654421  if (problemStatus_>-3||problemStatus_==-4) {
     
    883650    problemStatus_=4; // unknown
    884651  }
     652  if (saveThreshold) {
     653    // use default at present
     654    factorization_->sparseThreshold(0);
     655    factorization_->goSparse();
     656  }
     657  lastGoodIteration_ = numberIterations_;
    885658}
    886659/*
     
    14761249    specialOptions_ &= ~1;
    14771250}
    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
     1251/*
     1252  Reasons to come out (normal mode/user mode):
     1253  -1 normal
     1254  -2 factorize now - good iteration/ NA
     1255  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
     1256  -4 inaccuracy - refactorize - no iteration/ NA
     1257  -5 something flagged - go round again/ pivot not possible
     1258  +2 looks unbounded
     1259  +3 max iterations (iteration done)
    14821260*/
    1483 int 
     1261int
    14841262ClpSimplexPrimal::pivotResult()
    14851263{
    1486   abort();
    1487   return 0;
     1264
     1265  // Say if values pass
     1266  int ifValuesPass=(firstFree_>=0) ? 1 : 0;
     1267  bool roundAgain=true;
     1268  int returnCode=-1;
     1269
     1270  // loop round if user setting and doing refactorization
     1271  while (roundAgain) {
     1272    roundAgain=false;
     1273    returnCode=-1;
     1274    pivotRow_=-1;
     1275    sequenceOut_=-1;
     1276    rowArray_[1]->clear();
     1277    // we found a pivot column
     1278    // update the incoming column
     1279    unpack(rowArray_[1]);
     1280    // save reduced cost
     1281    double saveDj = dualIn_;
     1282    factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
     1283    // do ratio test and re-compute dj
     1284    primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
     1285              ifValuesPass);
     1286    if (ifValuesPass) {
     1287      saveDj=dualIn_;
     1288      if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
     1289        if(fabs(dualIn_)<1.0e2*dualTolerance_) {
     1290          // try other way
     1291          directionIn_=-directionIn_;
     1292          primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
     1293                    0);
     1294        }
     1295        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
     1296          if (solveType_==1) {
     1297            // reject it
     1298            char x = isColumn(sequenceIn_) ? 'C' :'R';
     1299            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     1300              <<x<<sequenceWithin(sequenceIn_)
     1301              <<CoinMessageEol;
     1302            setFlagged(sequenceIn_);
     1303            lastBadIteration_ = numberIterations_; // say be more cautious
     1304            rowArray_[1]->clear();
     1305            pivotRow_=-1;
     1306          }
     1307          returnCode=-5;
     1308          break;
     1309        }
     1310      }
     1311    }
     1312    if ((saveDj*dualIn_<1.0e-20&&!ifValuesPass&&solveType_==1)||
     1313        fabs(saveDj-dualIn_)>1.0e-3*(1.0+fabs(saveDj))) {
     1314      handler_->message(CLP_PRIMAL_DJ,messages_)
     1315        <<saveDj<<dualIn_
     1316        <<CoinMessageEol;
     1317      if(lastGoodIteration_ != numberIterations_) {
     1318        rowArray_[1]->clear();
     1319        pivotRow_=-1; // say no weights update
     1320        returnCode=-4;
     1321        if(lastGoodIteration_+1 == numberIterations_) {
     1322          // not looking wonderful - try cleaning bounds
     1323          // put non-basics to bounds in case tolerance moved
     1324          nonLinearCost_->checkInfeasibilities(true);
     1325        }
     1326        break;
     1327      } else {
     1328        // take on more relaxed criterion
     1329        if (saveDj*dualIn_<1.0e-20||
     1330            fabs(saveDj-dualIn_)>1.0e-1*(1.0+fabs(dualIn_))) {
     1331          // need to reject something
     1332          char x = isColumn(sequenceIn_) ? 'C' :'R';
     1333          handler_->message(CLP_SIMPLEX_FLAG,messages_)
     1334            <<x<<sequenceWithin(sequenceIn_)
     1335            <<CoinMessageEol;
     1336          setFlagged(sequenceIn_);
     1337          lastBadIteration_ = numberIterations_; // say be more cautious
     1338          rowArray_[1]->clear();
     1339          pivotRow_=-1;
     1340          returnCode=-5;
     1341          break;
     1342        }
     1343      }
     1344    }
     1345    if (pivotRow_>=0) {
     1346      if (solveType_==2) {
     1347        // do ray
     1348        primalRay(rowArray_[1]);
     1349      }
     1350      // if stable replace in basis
     1351      int updateStatus = factorization_->replaceColumn(rowArray_[2],
     1352                                                       pivotRow_,
     1353                                                       alpha_);
     1354      // if no pivots, bad update but reasonable alpha - take and invert
     1355      if (updateStatus==2&&
     1356          lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
     1357        updateStatus=4;
     1358      if (updateStatus==1||updateStatus==4) {
     1359        // slight error
     1360        if (factorization_->pivots()>5||updateStatus==4) {
     1361          returnCode=-3;
     1362        }
     1363      } else if (updateStatus==2) {
     1364        // major error
     1365        // later we may need to unwind more e.g. fake bounds
     1366        if(lastGoodIteration_ != numberIterations_) {
     1367          rowArray_[1]->clear();
     1368          pivotRow_=-1;
     1369          if (solveType_==1) {
     1370            returnCode=-4;
     1371            break;
     1372          } else {
     1373            // user in charge - re-factorize
     1374            int lastCleaned;
     1375            ClpSimplexProgress dummyProgress;
     1376            statusOfProblemInPrimal(lastCleaned,1,dummyProgress);
     1377            roundAgain=true;
     1378            continue;
     1379          }
     1380        } else {
     1381          // need to reject something
     1382          if (solveType_==1) {
     1383            char x = isColumn(sequenceIn_) ? 'C' :'R';
     1384            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     1385              <<x<<sequenceWithin(sequenceIn_)
     1386              <<CoinMessageEol;
     1387            setFlagged(sequenceIn_);
     1388          }
     1389          lastBadIteration_ = numberIterations_; // say be more cautious
     1390          rowArray_[1]->clear();
     1391          pivotRow_=-1;
     1392          returnCode = -5;
     1393          break;
     1394
     1395        }
     1396      } else if (updateStatus==3) {
     1397        // out of memory
     1398        // increase space if not many iterations
     1399        if (factorization_->pivots()<
     1400            0.5*factorization_->maximumPivots()&&
     1401            factorization_->pivots()<200)
     1402          factorization_->areaFactor(
     1403                                     factorization_->areaFactor() * 1.1);
     1404        returnCode =-2; // factorize now
     1405      }
     1406      // here do part of steepest - ready for next iteration
     1407      primalColumnPivot_->updateWeights(rowArray_[1]);
     1408    } else {
     1409      if (pivotRow_==-1) {
     1410        // no outgoing row is valid
     1411        rowArray_[0]->clear();
     1412        if (!factorization_->pivots()) {
     1413          returnCode = 2; //say looks unbounded
     1414          // do ray
     1415          primalRay(rowArray_[1]);
     1416        } else if (solveType_==2) {
     1417          // refactorize
     1418          int lastCleaned;
     1419          ClpSimplexProgress dummyProgress;
     1420          statusOfProblemInPrimal(lastCleaned,1,dummyProgress);
     1421          roundAgain=true;
     1422          continue;
     1423        } else {
     1424          returnCode = 4; //say looks unbounded but has iterated
     1425        }
     1426        break;
     1427      } else {
     1428        // flipping from bound to bound
     1429      }
     1430    }
     1431   
     1432    // update primal solution
     1433   
     1434    double objectiveChange=0.0;
     1435    // Cost on pivot row may change - may need to change dualIn
     1436    double oldCost=0.0;
     1437    if (pivotRow_>=0)
     1438      oldCost = cost(pivotVariable_[pivotRow_]);
     1439    // rowArray_[1] is not empty - used to update djs
     1440    updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
     1441    if (pivotRow_>=0)
     1442      dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
     1443    double oldValue = valueIn_;
     1444    if (directionIn_==-1) {
     1445      // as if from upper bound
     1446      if (sequenceIn_!=sequenceOut_) {
     1447        // variable becoming basic
     1448        valueIn_ -= fabs(theta_);
     1449      } else {
     1450        valueIn_=lowerIn_;
     1451      }
     1452    } else {
     1453      // as if from lower bound
     1454      if (sequenceIn_!=sequenceOut_) {
     1455        // variable becoming basic
     1456        valueIn_ += fabs(theta_);
     1457      } else {
     1458        valueIn_=upperIn_;
     1459      }
     1460    }
     1461    objectiveChange += dualIn_*(valueIn_-oldValue);
     1462    // outgoing
     1463    if (sequenceIn_!=sequenceOut_) {
     1464      if (directionOut_>0) {
     1465        valueOut_ = lowerOut_;
     1466      } else {
     1467        valueOut_ = upperOut_;
     1468      }
     1469      double lowerValue = lower_[sequenceOut_];
     1470      double upperValue = upper_[sequenceOut_];
     1471      assert(valueOut_>=lowerValue-primalTolerance_&&
     1472             valueOut_<=upperValue+primalTolerance_);
     1473      // may not be exactly at bound and bounds may have changed
     1474      if (valueOut_<=lowerValue+primalTolerance_) {
     1475        directionOut_=1;
     1476      } else if (valueOut_>=upperValue-primalTolerance_) {
     1477        directionOut_=-1;
     1478      } else {
     1479        printf("*** variable wandered off bound %g %g %g!\n",
     1480               lowerValue,valueOut_,upperValue);
     1481        if (valueOut_-lowerValue<=upperValue-valueOut_) {
     1482          valueOut_=lowerValue;
     1483          directionOut_=1;
     1484        } else {
     1485          valueOut_=upperValue;
     1486          directionOut_=-1;
     1487        }
     1488      }
     1489      solution_[sequenceOut_]=valueOut_;
     1490      nonLinearCost_->setOne(sequenceOut_,valueOut_);
     1491    }
     1492    // change cost and bounds on incoming if primal
     1493    nonLinearCost_->setOne(sequenceIn_,valueIn_);
     1494    int whatNext=housekeeping(objectiveChange);
     1495   
     1496    if (whatNext==1) {
     1497        returnCode =-2; // refactorize
     1498    } else if (whatNext==2) {
     1499      // maximum iterations or equivalent
     1500      returnCode=3;
     1501    } else if(numberIterations_ == lastGoodIteration_
     1502              + 2 * factorization_->maximumPivots()) {
     1503      // done a lot of flips - be safe
     1504      returnCode =-2; // refactorize
     1505    }
     1506  }
     1507  if (solveType_==2&&(returnCode == -2||returnCode==-3)) {
     1508    // refactorize here
     1509    int lastCleaned;
     1510    ClpSimplexProgress dummyProgress;
     1511    statusOfProblemInPrimal(lastCleaned,1,dummyProgress);
     1512  }
     1513#ifdef CLP_DEBUG
     1514  {
     1515    int i;
     1516    // not [1] as may have information
     1517    for (i=0;i<4;i++) {
     1518      if (i!=1)
     1519        rowArray_[i]->checkClear();
     1520    }   
     1521    for (i=0;i<2;i++) {
     1522      columnArray_[i]->checkClear();
     1523    }   
     1524  }     
     1525#endif
     1526  return returnCode;
     1527}
     1528// Create primal ray
     1529void
     1530ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
     1531{
     1532  delete [] ray_;
     1533  ray_ = new double [numberColumns_];
     1534  ClpFillN(ray_,numberColumns_,0.0);
     1535  int number=rowArray->getNumElements();
     1536  int * index = rowArray->getIndices();
     1537  double * array = rowArray->denseVector();
     1538  double way=-directionIn_;
     1539  int i;
     1540  double zeroTolerance=1.0e-12;
     1541  if (sequenceIn_<numberColumns_)
     1542    ray_[sequenceIn_]=directionIn_;
     1543  for (i=0;i<number;i++) {
     1544    int iRow=index[i];
     1545    int iPivot=pivotVariable_[iRow];
     1546    double arrayValue = array[iRow];
     1547    if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
     1548      ray_[iPivot] = way* array[iRow];
     1549  }
    14881550}
    14891551 
  • trunk/Presolve.cpp

    r110 r111  
    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 = true;
     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/PresolveSubst.cpp

    r109 r111  
    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/ClpSimplex.hpp

    r109 r111  
    143143      fixed, bounds are slightly looser than they could be.
    144144      This is to make dual go faster and is probably not needed
    145       with a presolve.  Returns non-zero if problem infeasible
    146   */
    147   int tightenPrimalBounds();
     145      with a presolve.  Returns non-zero if problem infeasible.
     146
     147      Fudge for branch and bound - put bounds on columns of factor *
     148      largest value (at continuous) - should improve stability
     149      in branch and bound on infeasible branches (0.0 is off)
     150  */
     151  int tightenPrimalBounds(double factor=0.0);
    148152  /** Crash - at present just aimed at dual, returns
    149153      -2 if dual preferred and crash basis created
     
    207211  void finish();
    208212 
     213  /** Factorizes and returns true if optimal.  Used by user */
     214  bool statusOfProblem();
    209215  //@}
    210216
     
    698704  /// Pivot Row
    699705  int pivotRow_;
     706  /// Last good iteration (immediately after a re-factorization)
     707  int lastGoodIteration_;
    700708  /// Working copy of reduced costs (Owner of arrays below)
    701709  double * dj_;
  • trunk/include/ClpSimplexPrimal.hpp

    r109 r111  
    124124  /** This has the flow between re-factorizations
    125125
    126       Reasons to come out:
    127       -1 iterations etc
    128       -2 inaccuracy
    129       -3 slight inaccuracy (and done iterations)
    130       -4 end of values pass and done iterations
    131       +0 looks optimal (might be infeasible - but we will investigate)
    132       +2 looks unbounded
     126      Returns a code to say where decision to exit was made
     127      Problem status set to:
     128
     129      -2 re-factorize
     130      -4 Looks optimal/infeasible
     131      -5 Looks unbounded
    133132      +3 max iterations
    134133   */
    135134  int whileIterating();
     135
     136  /** Do last half of an iteration.  This is split out so people can
     137      force incoming variable.  If solveType_ is 2 then this may
     138      re-factorize while normally it would exit to re-factorize.
     139      Return codes
     140      Reasons to come out (normal mode/user mode):
     141      -1 normal
     142      -2 factorize now - good iteration/ NA
     143      -3 slight inaccuracy - refactorize - iteration done/ same but factor done
     144      -4 inaccuracy - refactorize - no iteration/ NA
     145      -5 something flagged - go round again/ pivot not possible
     146      +2 looks unbounded
     147      +3 max iterations (iteration done)
     148
     149      With solveType_ ==2 this should
     150      Pivot in a variable and choose an outgoing one.  Assumes primal
     151      feasible - will not go through a bound.  Returns step length in theta
     152      Returns ray in ray_
     153  */
     154  int pivotResult();
     155
     156
    136157  /** The primals are updated by the given array.
    137158      Returns number of infeasibilities.
     
    190211  int unflag();
    191212
    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();
     213  /// Create primal ray
     214  void primalRay(CoinIndexedVector * rowArray);
    198215 
    199216  //@}
Note: See TracChangeset for help on using the changeset viewer.