Changeset 383 for trunk/ClpSimplex.cpp


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

add barrier

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.