Changeset 383


Ignore:
Timestamp:
Jun 1, 2004 2:38:34 PM (15 years ago)
Author:
forrest
Message:

add barrier

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpSimplex.cpp

    r372 r383  
    34873487{
    34883488  return ((ClpSimplexPrimalQuadratic *) this)->primalQuadratic(phase);
     3489}
     3490#include "ClpPredictorCorrector.hpp"
     3491#ifdef REAL_BARRIER
     3492#include "ClpCholeskyWssmp.hpp"
     3493#include "ClpCholeskyWssmpKKT.hpp"
     3494//#include "ClpCholeskyTaucs.hpp"
     3495#endif
     3496#include "ClpPresolve.hpp"
     3497/* Solves using barrier (assumes you have good cholesky factor code).
     3498   Does crossover to simplex if asked*/
     3499int
     3500ClpSimplex::barrier(bool crossover)
     3501{
     3502  ClpSimplex * model2 = this;
     3503  int savePerturbation=perturbation_;
     3504  ClpInterior barrier;
     3505  barrier.borrowModel(*model2);
     3506#ifdef REAL_BARRIER
     3507  // uncomment this if you have Anshul Gupta's wsmp package
     3508  ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(max(100,model2->numberRows()/10));
     3509  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
     3510  //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(max(100,model2->numberRows()/10));
     3511  // uncomment this if you have Sivan Toledo's Taucs package
     3512  //ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
     3513  barrier.setCholesky(cholesky);
     3514#endif
     3515  int numberRows = model2->numberRows();
     3516  int numberColumns = model2->numberColumns();
     3517  int saveMaxIts = model2->maximumIterations();
     3518  if (saveMaxIts<1000) {
     3519    barrier.setMaximumBarrierIterations(saveMaxIts);
     3520    model2->setMaximumIterations(1000000);
     3521  }
     3522  barrier.primalDual();
     3523  int barrierStatus=barrier.status();
     3524  double gap = barrier.complementarityGap();
     3525  // get which variables are fixed
     3526  double * saveLower=NULL;
     3527  double * saveUpper=NULL;
     3528  ClpPresolve pinfo2;
     3529  ClpSimplex * saveModel2=NULL;
     3530  int numberFixed = barrier.numberFixed();
     3531  if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover) {
     3532    // may as well do presolve
     3533    int numberRows = barrier.numberRows();
     3534    int numberColumns = barrier.numberColumns();
     3535    int numberTotal = numberRows+numberColumns;
     3536    saveLower = new double [numberTotal];
     3537    saveUpper = new double [numberTotal];
     3538    memcpy(saveLower,barrier.columnLower(),numberColumns*sizeof(double));
     3539    memcpy(saveLower+numberColumns,barrier.rowLower(),numberRows*sizeof(double));
     3540    memcpy(saveUpper,barrier.columnUpper(),numberColumns*sizeof(double));
     3541    memcpy(saveUpper+numberColumns,barrier.rowUpper(),numberRows*sizeof(double));
     3542    barrier.fixFixed();
     3543    saveModel2=model2;
     3544  }
     3545  barrier.returnModel(*model2);
     3546  double * rowPrimal = new double [numberRows];
     3547  double * columnPrimal = new double [numberColumns];
     3548  double * rowDual = new double [numberRows];
     3549  double * columnDual = new double [numberColumns];
     3550  // move solutions other way
     3551  CoinMemcpyN(model2->primalRowSolution(),
     3552              numberRows,rowPrimal);
     3553  CoinMemcpyN(model2->dualRowSolution(),
     3554              numberRows,rowDual);
     3555  CoinMemcpyN(model2->primalColumnSolution(),
     3556              numberColumns,columnPrimal);
     3557  CoinMemcpyN(model2->dualColumnSolution(),
     3558              numberColumns,columnDual);
     3559  if (saveModel2) {
     3560    // do presolve
     3561    model2 = pinfo2.presolvedModel(*model2,1.0e-8,
     3562                                   false,5,true);
     3563  }
     3564  if (barrierStatus<4&&crossover) {
     3565    // make sure no status left
     3566    model2->createStatus();
     3567    // solve
     3568    model2->setPerturbation(100);
     3569    // throw some into basis
     3570    {
     3571      int numberRows = model2->numberRows();
     3572      int numberColumns = model2->numberColumns();
     3573      double * dsort = new double[numberColumns];
     3574      int * sort = new int[numberColumns];
     3575      int n=0;
     3576      const double * columnLower = model2->columnLower();
     3577      const double * columnUpper = model2->columnUpper();
     3578      const double * primalSolution = model2->primalColumnSolution();
     3579      double tolerance = 10.0*primalTolerance_;
     3580      int i;
     3581      for ( i=0;i<numberRows;i++)
     3582        model2->setRowStatus(i,superBasic);
     3583      for ( i=0;i<numberColumns;i++) {
     3584        double distance = min(columnUpper[i]-primalSolution[i],
     3585                              primalSolution[i]-columnLower[i]);
     3586        if (distance>tolerance) {
     3587          dsort[n]=-distance;
     3588          sort[n++]=i;
     3589          model2->setStatus(i,superBasic);
     3590        } else if (distance>primalTolerance_) {
     3591          model2->setStatus(i,superBasic);
     3592        } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
     3593          model2->setStatus(i,atLowerBound);
     3594        } else {
     3595          model2->setStatus(i,atUpperBound);
     3596        }
     3597      }
     3598      CoinSort_2(dsort,dsort+n,sort);
     3599      n = min(numberRows,n);
     3600      for ( i=0;i<n;i++) {
     3601        int iColumn = sort[i];
     3602        model2->setStatus(iColumn,basic);
     3603      }
     3604      delete [] sort;
     3605      delete [] dsort;
     3606    }
     3607    if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
     3608      int numberRows = model2->numberRows();
     3609      int numberColumns = model2->numberColumns();
     3610      // just primal values pass
     3611      model2->primal(2);
     3612      // save primal solution and copy back dual
     3613      CoinMemcpyN(model2->primalRowSolution(),
     3614                  numberRows,rowPrimal);
     3615      CoinMemcpyN(rowDual,
     3616                  numberRows,model2->dualRowSolution());
     3617      CoinMemcpyN(model2->primalColumnSolution(),
     3618                  numberColumns,columnPrimal);
     3619      CoinMemcpyN(columnDual,
     3620                  numberColumns,model2->dualColumnSolution());
     3621      //model2->primal(1);
     3622      // clean up reduced costs and flag variables
     3623      {
     3624        double * dj = model2->dualColumnSolution();
     3625        double * cost = model2->objective();
     3626        double * saveCost = new double[numberColumns];
     3627        memcpy(saveCost,cost,numberColumns*sizeof(double));
     3628        double * saveLower = new double[numberColumns];
     3629        double * lower = model2->columnLower();
     3630        memcpy(saveLower,lower,numberColumns*sizeof(double));
     3631        double * saveUpper = new double[numberColumns];
     3632        double * upper = model2->columnUpper();
     3633        memcpy(saveUpper,upper,numberColumns*sizeof(double));
     3634        int i;
     3635        double tolerance = 10.0*dualTolerance_;
     3636        for ( i=0;i<numberColumns;i++) {
     3637          if (model2->getStatus(i)==basic) {
     3638            dj[i]=0.0;
     3639          } else if (model2->getStatus(i)==atLowerBound) {
     3640            if (optimizationDirection_*dj[i]<tolerance) {
     3641              if (optimizationDirection_*dj[i]<0.0) {
     3642                //if (dj[i]<-1.0e-3)
     3643                //printf("bad dj at lb %d %g\n",i,dj[i]);
     3644                cost[i] -= dj[i];
     3645                dj[i]=0.0;
     3646              }
     3647            } else {
     3648              upper[i]=lower[i];
     3649            }
     3650          } else if (model2->getStatus(i)==atUpperBound) {
     3651            if (optimizationDirection_*dj[i]>tolerance) {
     3652              if (optimizationDirection_*dj[i]>0.0) {
     3653                //if (dj[i]>1.0e-3)
     3654                //printf("bad dj at ub %d %g\n",i,dj[i]);
     3655                cost[i] -= dj[i];
     3656                dj[i]=0.0;
     3657              }
     3658            } else {
     3659              lower[i]=upper[i];
     3660            }
     3661          }
     3662        }
     3663        // just dual values pass
     3664        //model2->setLogLevel(63);
     3665        //model2->setFactorizationFrequency(1);
     3666        model2->dual(2);
     3667        memcpy(cost,saveCost,numberColumns*sizeof(double));
     3668        delete [] saveCost;
     3669        memcpy(lower,saveLower,numberColumns*sizeof(double));
     3670        delete [] saveLower;
     3671        memcpy(upper,saveUpper,numberColumns*sizeof(double));
     3672        delete [] saveUpper;
     3673      }
     3674      // and finish
     3675      // move solutions
     3676      CoinMemcpyN(rowPrimal,
     3677                  numberRows,model2->primalRowSolution());
     3678      CoinMemcpyN(columnPrimal,
     3679                  numberColumns,model2->primalColumnSolution());
     3680    }
     3681    model2->primal(1);
     3682  } else if (barrierStatus==4&&crossover) {
     3683    // memory problems
     3684    model2->setPerturbation(savePerturbation);
     3685    model2->createStatus();
     3686    model2->dual();
     3687  }
     3688  model2->setMaximumIterations(saveMaxIts);
     3689  delete [] rowPrimal;
     3690  delete [] columnPrimal;
     3691  delete [] rowDual;
     3692  delete [] columnDual;
     3693  if (saveLower) {
     3694    pinfo2.postsolve(true);
     3695    delete model2;
     3696    model2=saveModel2;
     3697    int numberRows = model2->numberRows();
     3698    int numberColumns = model2->numberColumns();
     3699    memcpy(model2->columnLower(),saveLower,numberColumns*sizeof(double));
     3700    memcpy(model2->rowLower(),saveLower+numberColumns,numberRows*sizeof(double));
     3701    delete [] saveLower;
     3702    memcpy(model2->columnUpper(),saveUpper,numberColumns*sizeof(double));
     3703    memcpy(model2->rowUpper(),saveUpper+numberColumns,numberRows*sizeof(double));
     3704    delete [] saveUpper;
     3705    model2->primal(1);
     3706  }
     3707  model2->setPerturbation(savePerturbation);
     3708  return model2->status();
    34893709}
    34903710/* For strong branching.  On input lower and upper are new bounds
  • trunk/Test/ClpMain.cpp

    r374 r383  
    1414#include "CoinPragma.hpp"
    1515#include "CoinHelperFunctions.hpp"
    16 #define CLPVERSION "0.99.6"
     16#define CLPVERSION "0.99.7"
    1717
    1818#include "CoinMpsIO.hpp"
     
    5050#endif
    5151
    52 int mainTest (int argc, const char *argv[],bool doDual,
     52int mainTest (int argc, const char *argv[],int algorithm,
    5353              ClpSimplex empty, bool doPresolve,int doIdiot);
    5454enum ClpParameterType {
    55   GENERALQUERY=-100,
     55  GENERALQUERY=-100,FULLGENERALQUERY,
    5656 
    5757  DUALTOLERANCE=1,PRIMALTOLERANCE,DUALBOUND,PRIMALWEIGHT,MAXTIME,OBJSCALE,
     
    6666  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,
    6767  MAXIMIZE,MINIMIZE,EXIT,STDIN,UNITTEST,NETLIB_DUAL,NETLIB_PRIMAL,SOLUTION,
    68   TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER,
     68  TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,REVERSE,BARRIER,NETLIB_BARRIER,
    6969
    7070  INVALID=1000
     
    897897      ClpItem("?","For help",GENERALQUERY,-1,false);
    898898    parameters[numberParameters++]=
     899      ClpItem("???","For help",FULLGENERALQUERY,-1,false);
     900    parameters[numberParameters++]=
    899901      ClpItem("-","From stdin",
    900902              STDIN,299,false);
     
    922924      (
    923925       "This command solves the current model using the  primal dual predictor \
    924 corrector algorithm.  This is not a sophisticated version just something JJF\
    925 knocked up\
    926 ** another slight drawback (early December) is that it does not work"
     926corrector algorithm.  This is not a sophisticated version just something JJF \
     927knocked up"
    927928
    928929       );
     
    11501151The user can set options before e.g. clp -presolve off -netlib"
    11511152       );
     1153#ifdef REAL_BARRIER
     1154    parameters[numberParameters++]=
     1155      ClpItem("netlibB!arrier","Solve entire netlib test set with barrier",
     1156              NETLIB_BARRIER);
     1157    parameters[numberParameters-1].setLonghelp
     1158      (
     1159       "This exercises the unit test for clp and then solves the netlib test set using barrier.\
     1160The user can set options before e.g. clp -presolve off -netlib"
     1161       );
     1162#endif
    11521163    parameters[numberParameters++]=
    11531164      ClpItem("netlibP!rimal","Solve entire netlib test set (primal)",
     
    14441455      // see if ? at end
    14451456      int numberQuery=0;
    1446       if (field!="?") {
     1457      if (field!="?"&&field!="???") {
    14471458        int length = field.length();
    14481459        int i;
     
    15001511              int type = parameters[iParam].indexNumber();
    15011512              if (parameters[iParam].displayThis()&&type>=limits[iType]
     1513                  &&type<limits[iType+1]) {
     1514                if (!across)
     1515                  std::cout<<"  ";
     1516                std::cout<<parameters[iParam].matchName()<<"  ";
     1517                across++;
     1518                if (across==maxAcross) {
     1519                  std::cout<<std::endl;
     1520                  across=0;
     1521                }
     1522              }
     1523            }
     1524            if (across)
     1525              std::cout<<std::endl;
     1526          }
     1527        } else if (type==FULLGENERALQUERY) {
     1528          std::cout<<"Full list of ommands is:"<<std::endl;
     1529          int maxAcross=5;
     1530          int limits[]={1,101,201,301,401};
     1531          std::vector<std::string> types;
     1532          types.push_back("Double parameters:");
     1533          types.push_back("Int parameters:");
     1534          types.push_back("Keyword parameters and others:");
     1535          types.push_back("Actions:");
     1536          int iType;
     1537          for (iType=0;iType<4;iType++) {
     1538            int across=0;
     1539            std::cout<<types[iType]<<std::endl;
     1540            for ( iParam=0; iParam<numberParameters; iParam++ ) {
     1541              int type = parameters[iParam].indexNumber();
     1542              if (type>=limits[iType]
    15021543                  &&type<limits[iType+1]) {
    15031544                if (!across)
     
    21302171            break;
    21312172          case NETLIB_DUAL:
     2173          case NETLIB_BARRIER:
    21322174          case NETLIB_PRIMAL:
    21332175            {
     
    21422184                nFields=4;
    21432185              }
    2144               if (type==NETLIB_DUAL)
     2186              int algorithm;
     2187              if (type==NETLIB_DUAL) {
    21452188                std::cerr<<"Doing netlib with dual agorithm"<<std::endl;
    2146               else
     2189                algorithm =0;
     2190              } else if (type==NETLIB_BARRIER) {
     2191                std::cerr<<"Doing netlib with barrier agorithm"<<std::endl;
     2192                algorithm =2;
     2193              } else {
    21472194                std::cerr<<"Doing netlib with primal agorithm"<<std::endl;
    2148               mainTest(nFields,fields,(type==NETLIB_DUAL),models[iModel],
     2195                algorithm=1;
     2196              }
     2197              mainTest(nFields,fields,algorithm,models[iModel],
    21492198                       (preSolve!=0),doIdiot);
    21502199            }
  • trunk/Test/unitTest.cpp

    r345 r383  
    5656// All parameters are optional.
    5757//----------------------------------------------------------------
    58 int mainTest (int argc, const char *argv[],bool doDual,
     58int mainTest (int argc, const char *argv[],int algorithm,
    5959              ClpSimplex empty, bool doPresolve,int doIdiot)
    6060{
     
    275275            model2->numberRows()>-5000)
    276276          model2->factorization()->maximumPivots(100+model2->numberRows()/100);
    277         if (doDual) {
     277        if (algorithm==0) {
    278278          // faster if bounds tightened
    279279          int numberInfeasibilities = model2->tightenPrimalBounds();
     
    284284          //model2->crash(1000,1);
    285285          model2->dual();
     286        } else if (algorithm==2) {
     287          model2->barrier();
    286288        } else {
    287289          if (doIdiot>0) {
     
    322324        }
    323325      } else {
    324         if (doDual) {
     326        if (algorithm==0) {
    325327          if (doIdiot<0)
    326328            solution.crash(1000,1);
    327329          solution.dual();
     330        } else if (algorithm==2) {
     331          solution.barrier();
    328332        } else {
    329333          if (doIdiot>0) {
  • trunk/include/ClpSimplex.hpp

    r372 r383  
    178178  /// Solves quadratic using Dantzig's algorithm - primal
    179179  int quadraticPrimal(int phase=2);
     180  /** Solves using barrier (assumes you have good cholesky factor code).
     181      Does crossover to simplex if asked*/
     182  int barrier(bool crossover=true);
    180183  /** Dual ranging.
    181184      This computes increase/decrease in cost for each given variable and corresponding
Note: See TracChangeset for help on using the changeset viewer.