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

Breaking out whileIterating

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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 */
Note: See TracChangeset for help on using the changeset viewer.