Changeset 673


Ignore:
Timestamp:
Oct 5, 2005 1:38:14 PM (14 years ago)
Author:
forrest
Message:

mods

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Samples/dualCuts.cpp

    r359 r673  
    1111
    1212#include "ClpSimplex.hpp"
     13#include "ClpPresolve.hpp"
     14#include "ClpFactorization.hpp"
    1315#include "CoinSort.hpp"
    1416#include "CoinHelperFunctions.hpp"
     
    2527    status=model.readMps("small.mps",true);
    2628  } else {
    27     status=model.readMps(argv[1],true);
     29    status=model.readMps(argv[1],false);
    2830  }
    2931  if (status)
     
    3537
    3638  double time1 = CoinCpuTime();
    37  
    38   int numberColumns = model.numberColumns();
    39   int originalNumberRows = model.numberRows();
     39  ClpSimplex * model2;
     40  ClpPresolve pinfo;
     41  int numberPasses=5; // can change this
     42  /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
     43     not being integer, do "numberpasses" passes and throw away names
     44     in presolved model */
     45  model2 = pinfo.presolvedModel(model,1.0e-8,false,numberPasses,false);
     46  if (!model2) {
     47    fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
     48            argv[1],1.0e-8);
     49    fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
     50            argv[1],1.0e-8);
     51    // model was infeasible - maybe try again with looser tolerances
     52    model2 = pinfo.presolvedModel(model,1.0e-7,false,numberPasses,false);
     53    if (!model2) {
     54      fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
     55              argv[1],1.0e-7);
     56      fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
     57              argv[1],1.0e-7);
     58      exit(2);
     59    }
     60  }
     61  // change factorization frequency from 200
     62  model2->setFactorizationFrequency(100+model2->numberRows()/50);
     63 
     64  int numberColumns = model2->numberColumns();
     65  int originalNumberRows = model2->numberRows();
    4066 
    4167  // We will need arrays to choose rows to add
     
    4571  char * take = new char [originalNumberRows];
    4672 
    47   const double * rowLower = model.rowLower();
    48   const double * rowUpper = model.rowUpper();
     73  const double * rowLower = model2->rowLower();
     74  const double * rowUpper = model2->rowUpper();
    4975  int iRow,iColumn;
    5076  // Set up initial list
     
    5783    }
    5884  }
     85  numberSort /= 2;
    5986  // Just add this number of rows each time in small problem
    6087  int smallNumberRows = 2*numberColumns;
     
    73100  */
    74101#if 0
    75   model.tightenPrimalBounds();
    76102  // However for some we need to do anyway
    77   double * columnLower = model.columnLower();
    78   double * columnUpper = model.columnUpper();
     103  double * columnLower = model2->columnLower();
     104  double * columnUpper = model2->columnUpper();
    79105  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    80106    columnLower[iColumn]=max(-1.0e6,columnLower[iColumn]);
     
    82108  }
    83109#endif
     110  model2->tightenPrimalBounds(-1.0e4,true);
    84111  printf("%d rows in initial problem\n",numberSort);
    85   double * fullSolution = model.primalRowSolution();
     112  double * fullSolution = model2->primalRowSolution();
    86113 
    87114  // Just do this number of passes
     
    91118  int iPass;
    92119 
    93   const int * start = model.clpMatrix()->getVectorStarts();
    94   const int * length = model.clpMatrix()->getVectorLengths();
    95   const int * row = model.clpMatrix()->getIndices();
     120  const int * start = model2->clpMatrix()->getVectorStarts();
     121  const int * length = model2->clpMatrix()->getVectorLengths();
     122  const int * row = model2->clpMatrix()->getIndices();
    96123  int * whichColumns = new int [numberColumns];
    97124  for (iColumn=0;iColumn<numberColumns;iColumn++)
     
    103130    std::sort(sort,sort+numberSort);
    104131    // Create small problem
    105     ClpSimplex small(&model,numberSort,sort,numberSmallColumns,whichColumns);
     132    ClpSimplex small(model2,numberSort,sort,numberSmallColumns,whichColumns);
    106133    small.setFactorizationFrequency(100+numberSort/200);
    107134    //small.setPerturbation(50);
     
    110137    //if (iPass)
    111138    //small.setMaximumIterations(100);
    112     // Solve
    113     small.dual();
     139    // Solve
     140    small.factorization()->messageLevel(8);
     141    if (iPass) {
     142      small.dual();
     143    } else {
     144      small.writeMps("continf.mps");
     145      ClpSolve solveOptions;
     146      solveOptions.setSolveType(ClpSolve::useDual);
     147      //solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
     148      //solveOptions.setSpecialOption(1,2,200); // idiot
     149      small.initialSolve(solveOptions);
     150    }
    114151    // move solution back
    115     double * solution = model.primalColumnSolution();
     152    double * solution = model2->primalColumnSolution();
    116153    const double * smallSolution = small.primalColumnSolution();
    117154    for (int j=0;j<numberSmallColumns;j++) {
    118155      iColumn = whichColumns[j];
    119156      solution[iColumn]=smallSolution[j];
    120       model.setColumnStatus(iColumn,small.getColumnStatus(j));
     157      model2->setColumnStatus(iColumn,small.getColumnStatus(j));
    121158    }
    122159    for (iRow=0;iRow<numberSort;iRow++) {
    123160      int kRow = sort[iRow];
    124       model.setRowStatus(kRow,small.getRowStatus(iRow));
     161      model2->setRowStatus(kRow,small.getRowStatus(iRow));
    125162    }
    126163    // compute full solution
    127164    memset(fullSolution,0,originalNumberRows*sizeof(double));
    128     model.times(1.0,model.primalColumnSolution(),fullSolution);
     165    model2->times(1.0,model2->primalColumnSolution(),fullSolution);
    129166    if (iPass!=maxPass-1) {
    130167      // Mark row as not looked at
     
    141178        iRow=sort[iSort];
    142179        //printf("%d %g %g\n",iRow,fullSolution[iRow],small.primalRowSolution()[iSort]);
    143         if (model.getRowStatus(iRow)==ClpSimplex::basic) {
     180        if (model2->getRowStatus(iRow)==ClpSimplex::basic) {
    144181          // Basic - we can get rid of if early on
    145182          if (iPass<takeOutPass) {
     
    223260  delete [] take;
    224261  // If problem is big you may wish to skip this
    225   model.dual();
     262  model2->dual();
    226263  int numberBinding=0;
    227264  for (iRow=0;iRow<originalNumberRows;iRow++) {
    228     if (model.getRowStatus(iRow)!=ClpSimplex::basic)
     265    if (model2->getRowStatus(iRow)!=ClpSimplex::basic)
    229266      numberBinding++;
    230267  }
    231268  printf("%d binding rows at end\n",numberBinding);
     269  pinfo.postsolve(true);
     270
     271  int numberIterations=model2->numberIterations();;
     272  delete model2;
     273  /* After this postsolve model should be optimal.
     274     We can use checkSolution and test feasibility */
     275  model.checkSolution();
     276  if (model.numberDualInfeasibilities()||
     277      model.numberPrimalInfeasibilities())
     278    printf("%g dual %g(%d) Primal %g(%d)\n",
     279           model.objectiveValue(),
     280           model.sumDualInfeasibilities(),
     281           model.numberDualInfeasibilities(),
     282           model.sumPrimalInfeasibilities(),
     283           model.numberPrimalInfeasibilities());
     284  // But resolve for safety
     285  model.primal(1);
     286
     287  numberIterations += model.numberIterations();;
    232288  printf("Solve took %g seconds\n",CoinCpuTime()-time1);
    233289  return 0;
Note: See TracChangeset for help on using the changeset viewer.