Changeset 14 for branches


Ignore:
Timestamp:
Aug 7, 2002 11:55:35 AM (17 years ago)
Author:
forrest
Message:

Breaking out whileIterating

Location:
branches/devel-1
Files:
6 edited

Legend:

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

    r8 r14  
    534534  index = infeasible_->getIndices();
    535535  number = infeasible_->getNumElements();
    536   // we can't really trust infeasibilities if there is dual error
    537   double checkTolerance = 1.0e-8;
    538   if (!model_->factorization()->pivots()&&
    539       model_->numberIterations()<model_->lastBadIteration()+200)
    540     checkTolerance = 1.0e-6;
    541   if (model_->largestDualError()>checkTolerance)
    542     tolerance *= model_->largestDualError()/checkTolerance;
    543 
     536  if(model_->numberIterations()<model_->lastBadIteration()+200) {
     537    // we can't really trust infeasibilities if there is dual error
     538    double checkTolerance = 1.0e-8;
     539    if (!model_->factorization()->pivots())
     540      checkTolerance = 1.0e-6;
     541    if (model_->largestDualError()>checkTolerance)
     542      tolerance *= model_->largestDualError()/checkTolerance;
     543  }
    544544  tolerance *= tolerance; // as we are using squares
    545545  for (i=0;i<number;i++) {
  • branches/devel-1/ClpSimplexDual.cpp

    r12 r14  
    325325    }
    326326
    327     // status stays at -1 while iterating, >=0 finished, -2 to invert
    328     // status -3 to go to top without an invert
    329     while (problemStatus_==-1) {
     327    // Do iterations
     328    whileIterating();
     329  }
     330
     331  // at present we are leaving factorization around
     332  // maybe we should empty it
     333  deleteRim();
     334  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
     335    <<objectiveValue()
     336    <<OsiMessageEol;
     337  // Restore any saved stuff
     338  perturbation_ = savePerturbation;
     339  factorization_->sparseThreshold(saveSparse);
     340  dualBound_ = saveDualBound_;
     341  return problemStatus_;
     342}
     343void
     344ClpSimplexDual::whileIterating()
     345{
     346  // status stays at -1 while iterating, >=0 finished, -2 to invert
     347  // status -3 to go to top without an invert
     348  while (problemStatus_==-1) {
    330349#ifdef CLP_DEBUG
    331       {
    332         int i;
    333         for (i=0;i<4;i++) {
    334           rowArray_[i]->checkClear();
    335         }   
    336         for (i=0;i<2;i++) {
    337           columnArray_[i]->checkClear();
    338         }   
    339       }     
     350    {
     351      int i;
     352      for (i=0;i<4;i++) {
     353        rowArray_[i]->checkClear();
     354      }   
     355      for (i=0;i<2;i++) {
     356        columnArray_[i]->checkClear();
     357      }   
     358    }     
    340359#endif
    341360#if CLP_DEBUG>2
    342       // very expensive
    343       if (numberIterations_>0&&numberIterations_<-801) {
    344         handler_->setLogLevel(63);
    345         double saveValue = objectiveValue_;
    346         double * saveRow1 = new double[numberRows_];
    347         double * saveRow2 = new double[numberRows_];
    348         memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
    349         memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
    350         double * saveColumn1 = new double[numberColumns_];
    351         double * saveColumn2 = new double[numberColumns_];
    352         memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
    353         memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
    354         gutsOfSolution(rowActivityWork_,columnActivityWork_);
    355         printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
    356                numberIterations_,
    357                saveValue,objectiveValue_,sumDualInfeasibilities_);
    358         if (saveValue>objectiveValue_+1.0e-2)
    359           printf("**bad**\n");
    360         memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
    361         memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
    362         memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
    363         memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
    364         delete [] saveRow1;
    365         delete [] saveRow2;
    366         delete [] saveColumn1;
    367         delete [] saveColumn2;
    368         objectiveValue_=saveValue;
     361    // very expensive
     362    if (numberIterations_>0&&numberIterations_<-801) {
     363      handler_->setLogLevel(63);
     364      double saveValue = objectiveValue_;
     365      double * saveRow1 = new double[numberRows_];
     366      double * saveRow2 = new double[numberRows_];
     367      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
     368      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
     369      double * saveColumn1 = new double[numberColumns_];
     370      double * saveColumn2 = new double[numberColumns_];
     371      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
     372      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
     373      gutsOfSolution(rowActivityWork_,columnActivityWork_);
     374      printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
     375             numberIterations_,
     376             saveValue,objectiveValue_,sumDualInfeasibilities_);
     377      if (saveValue>objectiveValue_+1.0e-2)
     378        printf("**bad**\n");
     379      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
     380      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
     381      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
     382      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
     383      delete [] saveRow1;
     384      delete [] saveRow2;
     385      delete [] saveColumn1;
     386      delete [] saveColumn2;
     387      objectiveValue_=saveValue;
     388    }
     389#endif
     390#ifdef CLP_DEBUG
     391    {
     392      int iSequence, number=numberRows_+numberColumns_;
     393      for (iSequence=0;iSequence<number;iSequence++) {
     394        double lowerValue=lower_[iSequence];
     395        double upperValue=upper_[iSequence];
     396        double value=solution_[iSequence];
     397        if(getStatus(iSequence)!=ClpSimplex::basic) {
     398          assert(lowerValue>-1.0e20);
     399          assert(upperValue<1.0e20);
     400        }
     401        switch(getStatus(iSequence)) {
     402         
     403        case ClpSimplex::basic:
     404          break;
     405        case ClpSimplex::isFree:
     406        case ClpSimplex::superBasic:
     407          break;
     408        case ClpSimplex::atUpperBound:
     409          assert (fabs(value-upperValue)<=primalTolerance_) ;
     410          break;
     411        case ClpSimplex::atLowerBound:
     412          assert (fabs(value-lowerValue)<=primalTolerance_) ;
     413          break;
     414        }
    369415      }
    370 #endif
     416    }
     417    if(numberIterations_==debugIteration) {
     418      printf("dodgy iteration coming up\n");
     419    }
     420#endif
     421    // choose row to go out
     422    dualRow();
     423    if (pivotRow_>=0) {
     424      // we found a pivot row
     425      handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
     426        <<pivotRow_
     427        <<OsiMessageEol;
     428      // check accuracy of weights
     429      dualRowPivot_->checkAccuracy();
     430      // get sign for finding row of tableau
     431      rowArray_[0]->insert(pivotRow_,directionOut_);
     432      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
     433      // put row of tableau in rowArray[0] and columnArray[0]
     434      matrix_->transposeTimes(this,-1.0,
     435                              rowArray_[0],columnArray_[1],columnArray_[0]);
     436      // rowArray has pi equivalent
     437      // do ratio test
     438      dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
     439                 rowArray_[3]);
     440      if (sequenceIn_>=0) {
     441        // normal iteration
     442        // update the incoming column
     443        unpack(rowArray_[1]);
     444        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
     445        // and update dual weights (can do in parallel - with extra array)
     446        dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
     447                                     rowArray_[1]);
     448        // see if update stable
     449        double btranAlpha = -alpha_*directionOut_; // for check
     450        alpha_=(*rowArray_[1])[pivotRow_];
    371451#ifdef CLP_DEBUG
    372       {
    373         int iSequence, number=numberRows_+numberColumns_;
    374         for (iSequence=0;iSequence<number;iSequence++) {
    375           double lowerValue=lower_[iSequence];
    376           double upperValue=upper_[iSequence];
    377           double value=solution_[iSequence];
    378           if(getStatus(iSequence)!=ClpSimplex::basic) {
    379             assert(lowerValue>-1.0e20);
    380             assert(upperValue<1.0e20);
    381           }
    382           switch(getStatus(iSequence)) {
    383            
    384           case ClpSimplex::basic:
     452        if ((handler_->logLevel()&32))
     453          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
     454#endif
     455        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
     456            fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
     457          handler_->message(CLP_DUAL_CHECK,messages_)
     458            <<btranAlpha
     459            <<alpha_
     460            <<OsiMessageEol;
     461          dualRowPivot_->unrollWeights();
     462          if (factorization_->pivots()) {
     463            problemStatus_=-2; // factorize now
     464            rowArray_[0]->clear();
     465            rowArray_[1]->clear();
     466            columnArray_[0]->clear();
    385467            break;
    386           case ClpSimplex::isFree:
    387           case ClpSimplex::superBasic:
    388             break;
    389           case ClpSimplex::atUpperBound:
    390             assert (fabs(value-upperValue)<=primalTolerance_) ;
    391             break;
    392           case ClpSimplex::atLowerBound:
    393             assert (fabs(value-lowerValue)<=primalTolerance_) ;
    394               break;
    395           }
    396         }
    397       }
    398       if(numberIterations_==debugIteration) {
    399         printf("dodgy iteration coming up\n");
    400       }
    401 #endif
    402       // choose row to go out
    403       dualRow();
    404       if (pivotRow_>=0) {
    405         // we found a pivot row
    406         handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
    407           <<pivotRow_
    408           <<OsiMessageEol;
    409         // check accuracy of weights
    410         dualRowPivot_->checkAccuracy();
    411         // get sign for finding row of tableau
    412         rowArray_[0]->insert(pivotRow_,directionOut_);
    413         factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
    414         // put row of tableau in rowArray[0] and columnArray[0]
    415         matrix_->transposeTimes(this,-1.0,
    416                                 rowArray_[0],columnArray_[1],columnArray_[0]);
    417         // rowArray has pi equivalent
    418         // do ratio test
    419         dualColumn(rowArray_[0],columnArray_[0],columnArray_[1],
    420                    rowArray_[3]);
    421         if (sequenceIn_>=0) {
    422           // normal iteration
    423           // update the incoming column
    424           unpack(rowArray_[1]);
    425           factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
    426           // and update dual weights (can do in parallel - with extra array)
    427           dualRowPivot_->updateWeights(rowArray_[0],rowArray_[2],
    428                                        rowArray_[1]);
    429           // see if update stable
    430           double btranAlpha = -alpha_*directionOut_; // for check
    431           alpha_=(*rowArray_[1])[pivotRow_];
    432 #ifdef CLP_DEBUG
    433           if ((handler_->logLevel()&32))
    434             printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
    435 #endif
    436           if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    437               fabs(btranAlpha-alpha_)>1.0e-7*(1.0+fabs(alpha_))) {
    438             handler_->message(CLP_DUAL_CHECK,messages_)
    439               <<btranAlpha
    440               <<alpha_
    441               <<OsiMessageEol;
    442             dualRowPivot_->unrollWeights();
    443             if (factorization_->pivots()) {
    444               problemStatus_=-2; // factorize now
    445               rowArray_[0]->clear();
    446               rowArray_[1]->clear();
    447               columnArray_[0]->clear();
    448               break;
    449             } else {
    450               // take on more relaxed criterion
    451               if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
    452                   fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
    453                 // need to reject something
    454                 char x = isColumn(sequenceOut_) ? 'C' :'R';
    455                 handler_->message(CLP_SIMPLEX_FLAG,messages_)
    456                   <<x<<sequenceWithin(sequenceOut_)
    457                   <<OsiMessageEol;
    458                 setFlagged(sequenceOut_);
    459                 lastBadIteration_ = numberIterations_; // say be more cautious
    460                 rowArray_[0]->clear();
    461                 rowArray_[1]->clear();
    462                 columnArray_[0]->clear();
    463                 continue;
    464               }
    465             }
    466           }
    467           // update duals BEFORE replaceColumn so can do updateColumn
    468           objectiveChange=0.0;
    469           // do duals first as variables may flip bounds
    470           // rowArray_[0] and columnArray_[0] may have flips
    471           // so use rowArray_[3] for work array from here on
    472           int nswapped =
    473             updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
    474                         objectiveChange);
    475           // which will change basic solution
    476           if (nswapped) {
    477 #ifdef CLP_DEBUG
    478           if ((handler_->logLevel()&16))
    479             printf("old dualOut_ %g, v %g, l %g, u %g - new ",
    480                    dualOut_,valueOut_,lowerOut_,upperOut_);
    481             double oldOut=dualOut_;
    482 #endif
    483             factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
    484             dualRowPivot_->updatePrimalSolution(rowArray_[2],
    485                                                 1.0,objectiveChange);
    486 
    487             // recompute dualOut_
    488             valueOut_ = solution_[sequenceOut_];
    489             if (directionOut_<0) {
    490               dualOut_ = valueOut_ - upperOut_;
    491             } else {
    492               dualOut_ = lowerOut_ - valueOut_;
    493             }
    494 #ifdef CLP_DEBUG
    495             if ((handler_->logLevel()&16))
    496               printf("%g\n",dualOut_);
    497             assert(dualOut_<=oldOut);
    498 #endif
    499             if(dualOut_<0.0&&factorization_->pivots()) {
    500               // going backwards - factorize
    501               dualRowPivot_->unrollWeights();
    502               problemStatus_=-2; // factorize now
    503               break;
    504             }
    505           }
    506           // amount primal will move
    507           double movement = -dualOut_*directionOut_/alpha_;
    508           // so objective should increase by fabs(dj)*movement
    509           // but we already have objective change - so check will be good
    510           if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
    511 #ifdef CLP_DEBUG
    512             if (handler_->logLevel()&32)
    513               printf("movement %g, swap change %g, rest %g  * %g\n",
    514                      objectiveChange+fabs(movement*dualIn_),
    515                      objectiveChange,movement,dualIn_);
    516 #endif
    517             if(factorization_->pivots()>5) {
    518               // going backwards - factorize
    519               dualRowPivot_->unrollWeights();
    520               problemStatus_=-2; // factorize now
    521               break;
    522             }
    523           }
    524           // if stable replace in basis
    525           int updateStatus = factorization_->replaceColumn(rowArray_[2],
    526                                                            pivotRow_,
    527                                                            alpha_);
    528           if (updateStatus==1) {
    529             // slight error
    530             if (factorization_->pivots()>5)
    531               problemStatus_=-2; // factorize now
    532           } else if (updateStatus==2) {
    533             // major error
    534             dualRowPivot_->unrollWeights();
    535             // later we may need to unwind more e.g. fake bounds
    536             if (factorization_->pivots()) {
    537               problemStatus_=-2; // factorize now
    538               break;
    539             } else {
     468          } else {
     469            // take on more relaxed criterion
     470            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
     471                fabs(btranAlpha-alpha_)>1.0e-4*(1.0+fabs(alpha_))) {
    540472              // need to reject something
    541473              char x = isColumn(sequenceOut_) ? 'C' :'R';
     
    550482              continue;
    551483            }
    552           } else if (updateStatus==3) {
    553             // out of memory
    554             // increase space if not many iterations
    555             if (factorization_->pivots()<
    556                 0.5*factorization_->maximumPivots()&&
    557                 factorization_->pivots()<200)
    558               factorization_->areaFactor(
    559                                          factorization_->areaFactor() * 1.1);
     484          }
     485        }
     486        // update duals BEFORE replaceColumn so can do updateColumn
     487        double objectiveChange=0.0;
     488        // do duals first as variables may flip bounds
     489        // rowArray_[0] and columnArray_[0] may have flips
     490        // so use rowArray_[3] for work array from here on
     491        int nswapped =
     492          updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[2],theta_,
     493                            objectiveChange);
     494        // which will change basic solution
     495        if (nswapped) {
     496#ifdef CLP_DEBUG
     497          if ((handler_->logLevel()&16))
     498            printf("old dualOut_ %g, v %g, l %g, u %g - new ",
     499                   dualOut_,valueOut_,lowerOut_,upperOut_);
     500          double oldOut=dualOut_;
     501#endif
     502          factorization_->updateColumn(rowArray_[3],rowArray_[2],false);
     503          dualRowPivot_->updatePrimalSolution(rowArray_[2],
     504                                              1.0,objectiveChange);
     505         
     506          // recompute dualOut_
     507          valueOut_ = solution_[sequenceOut_];
     508          if (directionOut_<0) {
     509            dualOut_ = valueOut_ - upperOut_;
     510          } else {
     511            dualOut_ = lowerOut_ - valueOut_;
     512          }
     513#ifdef CLP_DEBUG
     514          if ((handler_->logLevel()&16))
     515            printf("%g\n",dualOut_);
     516          assert(dualOut_<=oldOut);
     517#endif
     518          if(dualOut_<0.0&&factorization_->pivots()) {
     519            // going backwards - factorize
     520            dualRowPivot_->unrollWeights();
    560521            problemStatus_=-2; // factorize now
    561           }
    562           // update primal solution
    563           if (theta_<0.0) {
    564 #ifdef CLP_DEBUG
    565             if (handler_->logLevel()&32)
    566               printf("negative theta %g\n",theta_);
    567 #endif
    568             theta_=0.0;
    569           }
    570           // do actual flips
    571           flipBounds(rowArray_[0],columnArray_[0],theta_);
    572           dualRowPivot_->updatePrimalSolution(rowArray_[1],
    573                                               movement,
    574                                               objectiveChange);
    575 #ifdef CLP_DEBUG
    576           double oldobj=objectiveValue_;
    577 #endif
    578           int whatNext=housekeeping(objectiveChange);
    579           // and set bounds correctly
    580           originalBound(sequenceIn_);
    581           changeBound(sequenceOut_);
    582 #ifdef CLP_DEBUG
    583           if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
    584             printf("obj backwards %g %g\n",objectiveValue_,oldobj);
    585 #endif
    586           if (whatNext==1) {
    587             problemStatus_ =-2; // refactorize
    588           } else if (whatNext==2) {
    589             // maximum iterations or equivalent
    590             problemStatus_= 3;
    591522            break;
    592523          }
    593         } else {
    594           // no incoming column is valid
     524        }
     525        // amount primal will move
     526        double movement = -dualOut_*directionOut_/alpha_;
     527        // so objective should increase by fabs(dj)*movement
     528        // but we already have objective change - so check will be good
     529        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
    595530#ifdef CLP_DEBUG
    596531          if (handler_->logLevel()&32)
    597             printf("** no column pivot\n");
    598 #endif
    599           if (factorization_->pivots()<5) {
    600             problemStatus_=-4; //say looks infeasible
    601             // create ray anyway
    602             delete [] ray_;
    603             ray_ = new double [ numberRows_];
    604             CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
    605           }
    606           rowArray_[0]->clear();
    607           columnArray_[0]->clear();
     532            printf("movement %g, swap change %g, rest %g  * %g\n",
     533                   objectiveChange+fabs(movement*dualIn_),
     534                   objectiveChange,movement,dualIn_);
     535#endif
     536          if(factorization_->pivots()>5) {
     537            // going backwards - factorize
     538            dualRowPivot_->unrollWeights();
     539            problemStatus_=-2; // factorize now
     540            break;
     541          }
     542        }
     543        // if stable replace in basis
     544        int updateStatus = factorization_->replaceColumn(rowArray_[2],
     545                                                         pivotRow_,
     546                                                         alpha_);
     547        if (updateStatus==1) {
     548          // slight error
     549          if (factorization_->pivots()>5)
     550            problemStatus_=-2; // factorize now
     551        } else if (updateStatus==2) {
     552          // major error
     553          dualRowPivot_->unrollWeights();
     554          // later we may need to unwind more e.g. fake bounds
     555          if (factorization_->pivots()) {
     556            problemStatus_=-2; // factorize now
     557            break;
     558          } else {
     559            // need to reject something
     560            char x = isColumn(sequenceOut_) ? 'C' :'R';
     561            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     562              <<x<<sequenceWithin(sequenceOut_)
     563              <<OsiMessageEol;
     564            setFlagged(sequenceOut_);
     565            lastBadIteration_ = numberIterations_; // say be more cautious
     566            rowArray_[0]->clear();
     567            rowArray_[1]->clear();
     568            columnArray_[0]->clear();
     569            continue;
     570          }
     571        } else if (updateStatus==3) {
     572          // out of memory
     573          // increase space if not many iterations
     574          if (factorization_->pivots()<
     575              0.5*factorization_->maximumPivots()&&
     576              factorization_->pivots()<200)
     577            factorization_->areaFactor(
     578                                       factorization_->areaFactor() * 1.1);
     579          problemStatus_=-2; // factorize now
     580        }
     581        // update primal solution
     582        if (theta_<0.0) {
     583#ifdef CLP_DEBUG
     584          if (handler_->logLevel()&32)
     585            printf("negative theta %g\n",theta_);
     586#endif
     587          theta_=0.0;
     588        }
     589        // do actual flips
     590        flipBounds(rowArray_[0],columnArray_[0],theta_);
     591        dualRowPivot_->updatePrimalSolution(rowArray_[1],
     592                                            movement,
     593                                            objectiveChange);
     594#ifdef CLP_DEBUG
     595        double oldobj=objectiveValue_;
     596#endif
     597        int whatNext=housekeeping(objectiveChange);
     598        // and set bounds correctly
     599        originalBound(sequenceIn_);
     600        changeBound(sequenceOut_);
     601#ifdef CLP_DEBUG
     602        if (objectiveValue_<oldobj-1.0e-5&&(handler_->logLevel()&16))
     603          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
     604#endif
     605        if (whatNext==1) {
     606          problemStatus_ =-2; // refactorize
     607        } else if (whatNext==2) {
     608          // maximum iterations or equivalent
     609          problemStatus_= 3;
    608610          break;
    609611        }
    610612      } else {
    611         // no pivot row
     613        // no incoming column is valid
    612614#ifdef CLP_DEBUG
    613615        if (handler_->logLevel()&32)
    614           printf("** no row pivot\n");
    615 #endif
    616         if (!factorization_->pivots()) {
    617           // may have crept through - so may be optimal
    618           //problemStatus_=-5; //say looks unbounded
    619           problemStatus_=0;
    620           // check any flagged variables
    621           int iRow;
    622           for (iRow=0;iRow<numberRows_;iRow++) {
    623             int iPivot=pivotVariable_[iRow];
    624             if (flagged(iPivot))
    625               break;
    626           }
    627           if (iRow<numberRows_) {
    628 #ifdef CLP_DEBUG
    629             std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
    630 #endif
    631             problemStatus_=-4; //say looks infeasible
    632             // create ray anyway
    633             delete [] ray_;
    634             ray_ = new double [ numberRows_];
    635             CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
    636           }
    637         }
     616          printf("** no column pivot\n");
     617#endif
     618        if (factorization_->pivots()<5) {
     619          problemStatus_=-4; //say looks infeasible
     620          // create ray anyway
     621          delete [] ray_;
     622          ray_ = new double [ numberRows_];
     623          CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
     624        }
     625        rowArray_[0]->clear();
     626        columnArray_[0]->clear();
    638627        break;
    639628      }
    640     }
    641   }
    642 
    643   // at present we are leaving factorization around
    644   // maybe we should empty it
    645   deleteRim();
    646   handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
    647     <<objectiveValue()
    648     <<OsiMessageEol;
    649   // Restore any saved stuff
    650   perturbation_ = savePerturbation;
    651   factorization_->sparseThreshold(saveSparse);
    652   dualBound_ = saveDualBound_;
    653   return problemStatus_;
     629    } else {
     630      // no pivot row
     631#ifdef CLP_DEBUG
     632      if (handler_->logLevel()&32)
     633        printf("** no row pivot\n");
     634#endif
     635      if (!factorization_->pivots()) {
     636        // may have crept through - so may be optimal
     637        //problemStatus_=-5; //say looks unbounded
     638        problemStatus_=0;
     639        // check any flagged variables
     640        int iRow;
     641        for (iRow=0;iRow<numberRows_;iRow++) {
     642          int iPivot=pivotVariable_[iRow];
     643          if (flagged(iPivot))
     644            break;
     645        }
     646        if (iRow<numberRows_) {
     647#ifdef CLP_DEBUG
     648          std::cerr<<"Flagged variables at end - infeasible?"<<std::endl;
     649#endif
     650          problemStatus_=-4; //say looks infeasible
     651          // create ray anyway
     652          delete [] ray_;
     653          ray_ = new double [ numberRows_];
     654          CoinDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
     655        }
     656      }
     657      break;
     658    }
     659  }
    654660}
    655661/* The duals are updated by the given arrays.
    656     Returns number of infeasibilities.
    657     rowArray and columnarray will have flipped
    658     The output vector has movement (row length array) */
     662   Returns number of infeasibilities.
     663   rowArray and columnarray will have flipped
     664   The output vector has movement (row length array) */
    659665int
    660666ClpSimplexDual::updateDualsInDual(OsiIndexedVector * rowArray,
    661                                 OsiIndexedVector * columnArray,
    662                                 OsiIndexedVector * outputArray,
    663                                 double theta,
    664                                 double & objectiveChange)
     667                                  OsiIndexedVector * columnArray,
     668                                  OsiIndexedVector * outputArray,
     669                                  double theta,
     670                                  double & objectiveChange)
    665671{
    666 
     672 
    667673  outputArray->clear();
    668674 
     
    678684                       columnArray->getNumElements()==numberColumns_);
    679685  int numberAtFake=0;
    680 
     686 
    681687  // use a tighter tolerance except for all being okay
    682688  double tolerance = dualTolerance_;
    683 
     689 
    684690  double changeObj=0.0;
    685 
     691 
    686692  int iSection;
    687693 
  • branches/devel-1/ClpSimplexPrimal.cpp

    r13 r14  
    216216  assert(rowUpper_);
    217217
    218 #ifdef CLP_DEBUG
    219   int debugIteration=-1;
    220 #endif
    221218
    222219  algorithm_ = +1;
     
    249246    perturb();
    250247
    251   double objectiveChange;
    252248  // for primal we will change bounds using infeasibilityCost_
    253249  if (nonLinearCost_==NULL) {
     
    343339    saveNumber = numberIterations_;
    344340
    345     // status stays at -1 while iterating, >=0 finished, -2 to invert
    346     // status -3 to go to top without an invert
    347     while (problemStatus_==-1) {
     341    // Iterate
     342    whileIterating(firstSuperBasic);
     343  }
     344
     345  // if infeasible get real values
     346  if (problemStatus_) {
     347    infeasibilityCost_=0.0;
     348    createRim(7);
     349    nonLinearCost_->checkInfeasibilities(true);
     350    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
     351    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
     352    // and get good feasible duals
     353    computeDuals();
     354  }
     355  // at present we are leaving factorization around
     356  // maybe we should empty it
     357  deleteRim();
     358  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
     359    <<objectiveValue()
     360    <<OsiMessageEol;
     361  // Restore any saved stuff
     362  perturbation_ = savePerturbation;
     363  factorization_->sparseThreshold(saveSparse);
     364  infeasibilityCost_ = saveInfeasibilityCost;
     365  return problemStatus_;
     366}
     367void
     368ClpSimplexPrimal::whileIterating(int firstSuperBasic)
     369{
     370
     371  // Say if values pass
     372  int ifValuesPass=0;
     373  if (firstSuperBasic<numberRows_+numberColumns_)
     374    ifValuesPass=1;
     375  int saveNumber = numberIterations_;
     376  // status stays at -1 while iterating, >=0 finished, -2 to invert
     377  // status -3 to go to top without an invert
     378  while (problemStatus_==-1) {
    348379#ifdef CLP_DEBUG
    349       {
    350         int i;
    351         // not [1] as has information
    352         for (i=0;i<4;i++) {
    353           if (i!=1)
    354             rowArray_[i]->checkClear();
    355         }   
    356         for (i=0;i<2;i++) {
    357           columnArray_[i]->checkClear();
    358         }   
    359       }     
     380    {
     381      int i;
     382      // not [1] as has information
     383      for (i=0;i<4;i++) {
     384        if (i!=1)
     385          rowArray_[i]->checkClear();
     386      }   
     387      for (i=0;i<2;i++) {
     388        columnArray_[i]->checkClear();
     389      }   
     390    }     
    360391#endif
    361392#if CLP_DEBUG>2
    362       // very expensive
    363       if (numberIterations_>0&&numberIterations_<-2534) {
    364         handler_->setLogLevel(63);
    365         double saveValue = objectiveValue_;
    366         double * saveRow1 = new double[numberRows_];
    367         double * saveRow2 = new double[numberRows_];
    368         memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
    369         memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
    370         double * saveColumn1 = new double[numberColumns_];
    371         double * saveColumn2 = new double[numberColumns_];
    372         memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
    373         memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
    374         createRim(7);
    375         gutsOfSolution(rowActivityWork_,columnActivityWork_);
    376         printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
    377                numberIterations_,
    378                saveValue,objectiveValue_,sumPrimalInfeasibilities_);
    379         memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
    380         memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
    381         memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
    382         memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
    383         delete [] saveRow1;
    384         delete [] saveRow2;
    385         delete [] saveColumn1;
    386         delete [] saveColumn2;
    387         objectiveValue_=saveValue;
    388       }
    389 #endif
    390 #ifdef CLP_DEBUG
    391       if(numberIterations_==debugIteration) {
    392         printf("dodgy iteration coming up\n");
    393       }
    394 #endif
    395       if (!ifValuesPass) {
    396         // choose column to come in
    397         // can use pivotRow_ to update weights
    398         // pass in list of cost changes so can do row updates (rowArray_[1])
    399         // NOTE rowArray_[0] is used by computeDuals which is a
    400         // slow way of getting duals but might be used
    401         primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
    402                      columnArray_[0],columnArray_[1]);
    403       } else {
    404         // in values pass
    405         if (ifValuesPass>0) {
    406           nextSuperBasic(firstSuperBasic);
    407           if (firstSuperBasic==numberRows_+numberColumns_)
    408             ifValuesPass=-1; // signal end of values pass after this
    409         } else {
    410           // end of values pass - initialize weights etc
    411           primalColumnPivot_->saveWeights(this,5);
    412           ifValuesPass=0;
    413           if(saveNumber != numberIterations_) {
    414             problemStatus_=-2; // factorize now
    415             pivotRow_=-1; // say no weights update
    416             break;
    417           }
    418            
    419           // and get variable
    420           primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
    421                      columnArray_[0],columnArray_[1]);
    422         }
    423       }
    424       pivotRow_=-1;
    425       sequenceOut_=-1;
    426       rowArray_[1]->clear();
    427       if (sequenceIn_>=0) {
    428         // we found a pivot column
    429 #ifdef CLP_DEBUG
    430         if ((handler_->logLevel()&32)) {
    431           char x = isColumn(sequenceIn_) ? 'C' :'R';
    432           std::cout<<"pivot column "<<
    433                   x<<sequenceWithin(sequenceIn_)<<std::endl;
    434         }
    435 #endif
    436         // update the incoming column
    437         unpack(rowArray_[1]);
    438         // save reduced cost
    439         double saveDj = dualIn_;
    440         factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
    441         // do ratio test and re-compute dj
    442         primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
    443                   ifValuesPass);
    444         if (ifValuesPass) {
    445           saveDj=dualIn_;
    446           if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
    447             if(fabs(dualIn_)<1.0e2*dualTolerance_) {
    448               // try other way
    449               directionIn_=-directionIn_;
    450               primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
    451                         0);
    452             }
    453             if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
    454               // reject it
    455               char x = isColumn(sequenceIn_) ? 'C' :'R';
    456               handler_->message(CLP_SIMPLEX_FLAG,messages_)
    457                 <<x<<sequenceWithin(sequenceIn_)
    458                 <<OsiMessageEol;
    459               setFlagged(sequenceIn_);
    460               lastBadIteration_ = numberIterations_; // say be more cautious
    461               rowArray_[1]->clear();
    462               pivotRow_=-1;
    463               continue;
    464             }
    465           }
    466         }
    467              
    468 #ifdef CLP_DEBUG
    469         if ((handler_->logLevel()&32))
    470           printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_);
    471 #endif
    472         if (saveDj*dualIn_<1.0e-20||
    473             fabs(saveDj-dualIn_)>1.0e-5*(1.0+fabs(saveDj))) {
    474           handler_->message(CLP_PRIMAL_DJ,messages_)
    475             <<saveDj<<dualIn_
    476             <<OsiMessageEol;
    477           if(saveNumber != numberIterations_) {
    478             problemStatus_=-2; // factorize now
    479             rowArray_[1]->clear();
    480             pivotRow_=-1; // say no weights update
    481             break;
    482           } else {
    483             // take on more relaxed criterion
    484             if (saveDj*dualIn_<1.0e-20||
    485                 fabs(saveDj-dualIn_)>1.0e-4*(1.0+fabs(dualIn_))) {
    486               // need to reject something
    487               char x = isColumn(sequenceIn_) ? 'C' :'R';
    488               handler_->message(CLP_SIMPLEX_FLAG,messages_)
    489                 <<x<<sequenceWithin(sequenceIn_)
    490                 <<OsiMessageEol;
    491               setFlagged(sequenceIn_);
    492               lastBadIteration_ = numberIterations_; // say be more cautious
    493               rowArray_[1]->clear();
    494               pivotRow_=-1;
    495               continue;
    496             }
    497           }
    498         }
    499         if (pivotRow_>=0) {
    500           // if stable replace in basis
    501           int updateStatus = factorization_->replaceColumn(rowArray_[2],
    502                                                            pivotRow_,
    503                                                            alpha_);
    504           if (updateStatus==1) {
    505             // slight error
    506             if (factorization_->pivots()>5)
    507               problemStatus_=-2; // factorize now
    508           } else if (updateStatus==2) {
    509             // major error
    510             // later we may need to unwind more e.g. fake bounds
    511             if(saveNumber != numberIterations_) {
    512               problemStatus_=-2; // factorize now
    513               break;
    514             } else {
    515               // need to reject something
    516               char x = isColumn(sequenceIn_) ? 'C' :'R';
    517               handler_->message(CLP_SIMPLEX_FLAG,messages_)
    518                 <<x<<sequenceWithin(sequenceIn_)
    519                 <<OsiMessageEol;
    520               setFlagged(sequenceIn_);
    521               lastBadIteration_ = numberIterations_; // say be more cautious
    522               rowArray_[1]->clear();
    523               pivotRow_=-1;
    524               continue;
    525             }
    526           } else if (updateStatus==3) {
    527             // out of memory
    528             // increase space if not many iterations
    529             if (factorization_->pivots()<
    530                 0.5*factorization_->maximumPivots()&&
    531                 factorization_->pivots()<200)
    532               factorization_->areaFactor(
    533                                          factorization_->areaFactor() * 1.1);
    534             problemStatus_=-2; // factorize now
    535           }
    536           // here do part of steepest - ready for next iteration
    537           primalColumnPivot_->updateWeights(rowArray_[1]);
    538         } else {
    539           if (pivotRow_==-1) {
    540             // no outgoing row is valid
    541 #ifdef CLP_DEBUG
    542             if (handler_->logLevel()&32)
    543               printf("** no row pivot\n");
    544 #endif
    545             if (!factorization_->pivots()) {
    546               problemStatus_=-5; //say looks unbounded
    547               // do ray
    548               delete [] ray_;
    549               ray_ = new double [numberColumns_];
    550               CoinFillN(ray_,numberColumns_,0.0);
    551               int number=rowArray_[1]->getNumElements();
    552               int * index = rowArray_[1]->getIndices();
    553               double * array = rowArray_[1]->denseVector();
    554               double way=-directionIn_;
    555               int i;
    556               double zeroTolerance=1.0e-12;
    557               if (sequenceIn_<numberColumns_)
    558                 ray_[sequenceIn_]=directionIn_;
    559               for (i=0;i<number;i++) {
    560                 int iRow=index[i];
    561                 int iPivot=pivotVariable_[iRow];
    562                 double arrayValue = array[iRow];
    563                 if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
    564                   ray_[iPivot] = way* array[iRow];
    565               }
    566             }
    567             rowArray_[0]->clear();
    568             break;
    569           } else {
    570             // flipping from bound to bound
    571           }
    572         }
    573 
    574         // update primal solution
    575        
    576         objectiveChange=0.0;
    577         // Cost on pivot row may change - may need to change dualIn
    578         double oldCost=0.0;
    579         if (pivotRow_>=0)
    580           oldCost = cost(pivotVariable_[pivotRow_]);
    581         // rowArray_[1] is not empty - used to update djs
    582         updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
    583         if (pivotRow_>=0)
    584           dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
    585        
    586         int whatNext=housekeeping(objectiveChange);
    587        
    588         if (whatNext==1) {
    589           problemStatus_ =-2; // refactorize
    590         } else if (whatNext==2) {
    591           // maximum iterations or equivalent
    592           problemStatus_= 3;
    593           break;
    594         }
    595       } else {
    596         // no pivot column
    597 #ifdef CLP_DEBUG
    598         if (handler_->logLevel()&32)
    599           printf("** no column pivot\n");
    600 #endif
    601         if (nonLinearCost_->numberInfeasibilities())
    602           problemStatus_=-4; // might be infeasible
    603         break;
    604       }
    605     }
    606 #if CLP_DEBUG>2
    607     if (numberIterations_>620&&numberIterations_<-2534) {
     393    // very expensive
     394    if (numberIterations_>0&&numberIterations_<-2534) {
    608395      handler_->setLogLevel(63);
    609396      double saveValue = objectiveValue_;
     
    632419    }
    633420#endif
    634   }
    635 
    636   // if infeasible get real values
    637   if (problemStatus_) {
    638     infeasibilityCost_=0.0;
    639     createRim(7);
    640     nonLinearCost_->checkInfeasibilities(true);
    641     sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
    642     numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
    643   }
    644   // at present we are leaving factorization around
    645   // maybe we should empty it
    646   deleteRim();
    647   handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
    648     <<objectiveValue()
    649     <<OsiMessageEol;
    650   // Restore any saved stuff
    651   perturbation_ = savePerturbation;
    652   factorization_->sparseThreshold(saveSparse);
    653   infeasibilityCost_ = saveInfeasibilityCost;
    654   return problemStatus_;
     421    if (!ifValuesPass) {
     422      // choose column to come in
     423      // can use pivotRow_ to update weights
     424      // pass in list of cost changes so can do row updates (rowArray_[1])
     425      // NOTE rowArray_[0] is used by computeDuals which is a
     426      // slow way of getting duals but might be used
     427      primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
     428                   columnArray_[0],columnArray_[1]);
     429    } else {
     430      // in values pass
     431      if (ifValuesPass>0) {
     432        nextSuperBasic(firstSuperBasic);
     433        if (firstSuperBasic==numberRows_+numberColumns_)
     434          ifValuesPass=-1; // signal end of values pass after this
     435      } else {
     436        // end of values pass - initialize weights etc
     437        primalColumnPivot_->saveWeights(this,5);
     438        ifValuesPass=0;
     439        if(saveNumber != numberIterations_) {
     440          problemStatus_=-2; // factorize now
     441          pivotRow_=-1; // say no weights update
     442          break;
     443        }
     444       
     445        // and get variable
     446        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
     447                     columnArray_[0],columnArray_[1]);
     448      }
     449    }
     450    pivotRow_=-1;
     451    sequenceOut_=-1;
     452    rowArray_[1]->clear();
     453    if (sequenceIn_>=0) {
     454      // we found a pivot column
     455#ifdef CLP_DEBUG
     456      if ((handler_->logLevel()&32)) {
     457        char x = isColumn(sequenceIn_) ? 'C' :'R';
     458        std::cout<<"pivot column "<<
     459          x<<sequenceWithin(sequenceIn_)<<std::endl;
     460      }
     461#endif
     462      // update the incoming column
     463      unpack(rowArray_[1]);
     464      // save reduced cost
     465      double saveDj = dualIn_;
     466      factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
     467      // do ratio test and re-compute dj
     468      primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
     469                ifValuesPass);
     470      if (ifValuesPass) {
     471        saveDj=dualIn_;
     472        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
     473          if(fabs(dualIn_)<1.0e2*dualTolerance_) {
     474            // try other way
     475            directionIn_=-directionIn_;
     476            primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
     477                      0);
     478          }
     479          if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
     480            // reject it
     481            char x = isColumn(sequenceIn_) ? 'C' :'R';
     482            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     483              <<x<<sequenceWithin(sequenceIn_)
     484              <<OsiMessageEol;
     485            setFlagged(sequenceIn_);
     486            lastBadIteration_ = numberIterations_; // say be more cautious
     487            rowArray_[1]->clear();
     488            pivotRow_=-1;
     489            continue;
     490          }
     491        }
     492      }
     493     
     494#ifdef CLP_DEBUG
     495      if ((handler_->logLevel()&32))
     496        printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_);
     497#endif
     498      if (saveDj*dualIn_<1.0e-20||
     499          fabs(saveDj-dualIn_)>1.0e-5*(1.0+fabs(saveDj))) {
     500        handler_->message(CLP_PRIMAL_DJ,messages_)
     501          <<saveDj<<dualIn_
     502          <<OsiMessageEol;
     503        if(saveNumber != numberIterations_) {
     504          problemStatus_=-2; // factorize now
     505          rowArray_[1]->clear();
     506          pivotRow_=-1; // say no weights update
     507          break;
     508        } else {
     509          // take on more relaxed criterion
     510          if (saveDj*dualIn_<1.0e-20||
     511              fabs(saveDj-dualIn_)>1.0e-4*(1.0+fabs(dualIn_))) {
     512            // need to reject something
     513            char x = isColumn(sequenceIn_) ? 'C' :'R';
     514            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     515              <<x<<sequenceWithin(sequenceIn_)
     516              <<OsiMessageEol;
     517            setFlagged(sequenceIn_);
     518            lastBadIteration_ = numberIterations_; // say be more cautious
     519            rowArray_[1]->clear();
     520            pivotRow_=-1;
     521            continue;
     522          }
     523        }
     524      }
     525      if (pivotRow_>=0) {
     526        // if stable replace in basis
     527        int updateStatus = factorization_->replaceColumn(rowArray_[2],
     528                                                         pivotRow_,
     529                                                         alpha_);
     530        if (updateStatus==1) {
     531          // slight error
     532          if (factorization_->pivots()>5)
     533            problemStatus_=-2; // factorize now
     534        } else if (updateStatus==2) {
     535          // major error
     536          // later we may need to unwind more e.g. fake bounds
     537          if(saveNumber != numberIterations_) {
     538            problemStatus_=-2; // factorize now
     539            break;
     540          } else {
     541            // need to reject something
     542            char x = isColumn(sequenceIn_) ? 'C' :'R';
     543            handler_->message(CLP_SIMPLEX_FLAG,messages_)
     544              <<x<<sequenceWithin(sequenceIn_)
     545              <<OsiMessageEol;
     546            setFlagged(sequenceIn_);
     547            lastBadIteration_ = numberIterations_; // say be more cautious
     548            rowArray_[1]->clear();
     549            pivotRow_=-1;
     550            continue;
     551          }
     552        } else if (updateStatus==3) {
     553          // out of memory
     554          // increase space if not many iterations
     555          if (factorization_->pivots()<
     556              0.5*factorization_->maximumPivots()&&
     557              factorization_->pivots()<200)
     558            factorization_->areaFactor(
     559                                       factorization_->areaFactor() * 1.1);
     560          problemStatus_=-2; // factorize now
     561        }
     562        // here do part of steepest - ready for next iteration
     563        primalColumnPivot_->updateWeights(rowArray_[1]);
     564      } else {
     565        if (pivotRow_==-1) {
     566          // no outgoing row is valid
     567#ifdef CLP_DEBUG
     568          if (handler_->logLevel()&32)
     569            printf("** no row pivot\n");
     570#endif
     571          if (!factorization_->pivots()) {
     572            problemStatus_=-5; //say looks unbounded
     573            // do ray
     574            delete [] ray_;
     575            ray_ = new double [numberColumns_];
     576            CoinFillN(ray_,numberColumns_,0.0);
     577            int number=rowArray_[1]->getNumElements();
     578            int * index = rowArray_[1]->getIndices();
     579            double * array = rowArray_[1]->denseVector();
     580            double way=-directionIn_;
     581            int i;
     582            double zeroTolerance=1.0e-12;
     583            if (sequenceIn_<numberColumns_)
     584              ray_[sequenceIn_]=directionIn_;
     585            for (i=0;i<number;i++) {
     586              int iRow=index[i];
     587              int iPivot=pivotVariable_[iRow];
     588              double arrayValue = array[iRow];
     589              if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
     590                ray_[iPivot] = way* array[iRow];
     591            }
     592          }
     593          rowArray_[0]->clear();
     594          break;
     595        } else {
     596          // flipping from bound to bound
     597        }
     598      }
     599     
     600      // update primal solution
     601     
     602      double objectiveChange=0.0;
     603      // Cost on pivot row may change - may need to change dualIn
     604      double oldCost=0.0;
     605      if (pivotRow_>=0)
     606        oldCost = cost(pivotVariable_[pivotRow_]);
     607      // rowArray_[1] is not empty - used to update djs
     608      updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
     609      if (pivotRow_>=0)
     610        dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
     611     
     612      int whatNext=housekeeping(objectiveChange);
     613     
     614      if (whatNext==1) {
     615        problemStatus_ =-2; // refactorize
     616      } else if (whatNext==2) {
     617        // maximum iterations or equivalent
     618        problemStatus_= 3;
     619        break;
     620      }
     621    } else {
     622      // no pivot column
     623#ifdef CLP_DEBUG
     624      if (handler_->logLevel()&32)
     625        printf("** no column pivot\n");
     626#endif
     627      if (nonLinearCost_->numberInfeasibilities())
     628        problemStatus_=-4; // might be infeasible
     629      break;
     630    }
     631  }
    655632}
    656633/* Checks if finished.  Updates status */
  • branches/devel-1/include/ClpModel.hpp

    r8 r14  
    183183  inline const double * getReducedCost() const
    184184          { return reducedCost_;} ;
    185   /// Row lower or rhs
     185  /// Row lower
    186186  inline double* rowLower() const
    187187          { return rowLower_;};
     
    200200  /// Row Objective
    201201  inline double * rowObjective() const
     202          { return rowObjective_;};
     203  inline const double * getRowObjCoefficients() const
    202204          { return rowObjective_;};
    203205  /// Column Lower
     
    225227  inline double getObjValue() const
    226228  { return objectiveValue_*optimizationDirection_ - dblParam_[OsiObjOffset];};
    227   /// Infeasibility/unbounded ray (NULL returned if none/wrong)
     229  /** Infeasibility/unbounded ray (NULL returned if none/wrong)
     230      Up to user to use delete [] on these arrays.  */
    228231  double * infeasibilityRay() const;
    229232  double * unboundedRay() const;
  • branches/devel-1/include/ClpSimplexDual.hpp

    r2 r14  
    119119  /**@name Functions used in dual */
    120120  //@{
     121  /** This has the flow between re-factorizations
     122      Broken out for clarity and will be used by strong branching
     123   */
     124  void whileIterating();
    121125  /** The duals are updated by the given arrays.
    122126      Returns number of infeasibilities.
  • branches/devel-1/include/ClpSimplexPrimal.hpp

    r2 r14  
    122122  /**@name Functions used in primal */
    123123  //@{
     124  /** This has the flow between re-factorizations
     125      firstSuperBasic == number rows + columns normally,
     126      otherwise first super basic variable
     127   */
     128  void whileIterating(int firstSuperBasic);
    124129  /** The primals are updated by the given array.
    125130      Returns number of infeasibilities.
Note: See TracChangeset for help on using the changeset viewer.