source: trunk/Clp/examples/sprint.cpp @ 800

Last change on this file since 800 was 230, checked in by forrest, 16 years ago

3000 min

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.6 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpSimplex.hpp"
5#include "CoinSort.hpp"
6#include <iomanip>
7
8int main (int argc, const char *argv[])
9{
10  ClpSimplex  model;
11  int status;
12  // Keep names
13  if (argc<2) {
14    status=model.readMps("small.mps",true);
15  } else {
16    status=model.readMps(argv[1],true);
17  }
18  if (status)
19    exit(10);
20  /*
21    This driver implements what I called Sprint.  Cplex calls it
22    "sifting" which is just as silly.  When I thought of this trivial idea
23    it reminded me of an LP code of the 60's called sprint which after
24    every factorization took a subset of the matrix into memory (all
25    64K words!) and then iterated very fast on that subset.  On the
26    problems of those days it did not work very well, but it worked very
27    well on aircrew scheduling problems where there were very large numbers
28    of columns all with the same flavor.
29  */
30
31  /* The idea works best if you can get feasible easily.  To make it
32     more general we can add in costed slacks */
33
34  int originalNumberColumns = model.numberColumns();
35  int numberRows = model.numberRows();
36
37  // We will need arrays to choose variables.  These are too big but ..
38  double * weight = new double [numberRows+originalNumberColumns];
39  int * sort = new int [numberRows+originalNumberColumns];
40  int numberSort=0;
41  // Say we are going to add slacks - if you can get a feasible
42  // solution then do that at the comment - Add in your own coding here
43  bool addSlacks = true;
44
45  if (addSlacks) {
46    // initial list will just be artificials
47    // first we will set all variables as close to zero as possible
48    int iColumn;
49    const double * columnLower = model.columnLower();
50    const double * columnUpper = model.columnUpper();
51    double * columnSolution = model.primalColumnSolution();
52
53    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
54      double value =0.0;
55      if (columnLower[iColumn]>0.0)
56        value = columnLower[iColumn];
57      else if (columnUpper[iColumn]<0.0)
58        value = columnUpper[iColumn];
59      columnSolution[iColumn]=value;
60    }
61    // now see what that does to row solution
62    double * rowSolution = model.primalRowSolution();
63    memset (rowSolution,0,numberRows*sizeof(double));
64    model.times(1.0,columnSolution,rowSolution);
65
66    int * addStarts = new int [numberRows+1];
67    int * addRow = new int[numberRows];
68    double * addElement = new double[numberRows];
69    const double * lower = model.rowLower();
70    const double * upper = model.rowUpper();
71    addStarts[0]=0;
72    int numberArtificials=0;
73    double * addCost = new double [numberRows];
74    const double penalty=1.0e8;
75    int iRow;
76    for (iRow=0;iRow<numberRows;iRow++) {
77      if (lower[iRow]>rowSolution[iRow]) {
78        addRow[numberArtificials]=iRow;
79        addElement[numberArtificials]=1.0;
80        addCost[numberArtificials]=penalty;
81        numberArtificials++;
82        addStarts[numberArtificials]=numberArtificials;
83      } else if (upper[iRow]<rowSolution[iRow]) {
84        addRow[numberArtificials]=iRow;
85        addElement[numberArtificials]=-1.0;
86        addCost[numberArtificials]=penalty;
87        numberArtificials++;
88        addStarts[numberArtificials]=numberArtificials;
89      }
90    }
91    model.addColumns(numberArtificials,NULL,NULL,addCost,
92                      addStarts,addRow,addElement);
93    delete [] addStarts;
94    delete [] addRow;
95    delete [] addElement;
96    delete [] addCost;
97    // Set up initial list
98    numberSort=numberArtificials;
99    int i;
100    for (i=0;i<numberSort;i++)
101      sort[i] = i+originalNumberColumns;
102  } else {
103    // Get initial list in some magical way
104    // Add in your own coding here
105    abort();
106  }
107
108  int numberColumns = model.numberColumns();
109  const double * columnLower = model.columnLower();
110  const double * columnUpper = model.columnUpper();
111  double * fullSolution = model.primalColumnSolution();
112
113  // Just do this number of passes
114  int maxPass=100;
115  int iPass;
116  double lastObjective=1.0e31;
117
118  // Just take this number of columns in small problem
119  int smallNumberColumns = min(3*numberRows,numberColumns);
120  smallNumberColumns = max(smallNumberColumns,3000);
121  // We will be using all rows
122  int * whichRows = new int [numberRows];
123  for (int iRow=0;iRow<numberRows;iRow++)
124    whichRows[iRow]=iRow;
125  double originalOffset;
126  model.getDblParam(ClpObjOffset,originalOffset);
127
128  for (iPass=0;iPass<maxPass;iPass++) {
129    printf("Start of pass %d\n",iPass);
130    //printf("Bug until submodel new version\n");
131    CoinSort_2(sort,sort+numberSort,weight);
132    // Create small problem
133    ClpSimplex small(&model,numberRows,whichRows,numberSort,sort);
134    // now see what variables left out do to row solution
135    double * rowSolution = model.primalRowSolution();
136    memset (rowSolution,0,numberRows*sizeof(double));
137    int iRow,iColumn;
138    // zero out ones in small problem
139    for (iColumn=0;iColumn<numberSort;iColumn++) {
140      int kColumn = sort[iColumn];
141      fullSolution[kColumn]=0.0;
142    }
143    // Get objective offset
144    double offset=0.0;
145    const double * objective = model.objective();
146    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) 
147      offset += fullSolution[iColumn]*objective[iColumn];
148    small.setDblParam(ClpObjOffset,originalOffset-offset);
149    model.times(1.0,fullSolution,rowSolution);
150
151    double * lower = small.rowLower();
152    double * upper = small.rowUpper();
153    for (iRow=0;iRow<numberRows;iRow++) {
154      if (lower[iRow]>-1.0e50) 
155        lower[iRow] -= rowSolution[iRow];
156      if (upper[iRow]<1.0e50)
157        upper[iRow] -= rowSolution[iRow];
158    }
159    /* For some problems a useful variant is to presolve problem.
160       In this case you need to adjust smallNumberColumns to get
161       right size problem.  Also you can dispense with creating
162       small problem and fix variables in large problem and do presolve
163       on that. */
164    // Solve
165    small.primal();
166    // move solution back
167    const double * solution = small.primalColumnSolution();
168    for (iColumn=0;iColumn<numberSort;iColumn++) {
169      int kColumn = sort[iColumn];
170      model.setColumnStatus(kColumn,small.getColumnStatus(iColumn));
171      fullSolution[kColumn]=solution[iColumn];
172    }
173    for (iRow=0;iRow<numberRows;iRow++) 
174      model.setRowStatus(iRow,small.getRowStatus(iRow));
175    memcpy(model.primalRowSolution(),small.primalRowSolution(),
176           numberRows*sizeof(double));
177    if ((small.objectiveValue()>lastObjective-1.0e-7&&iPass>5)||
178        !small.numberIterations()||
179        iPass==maxPass-1) {
180
181      break; // finished
182    } else {
183      lastObjective = small.objectiveValue();
184      // get reduced cost for large problem
185      // this assumes minimization
186      memcpy(weight,model.objective(),numberColumns*sizeof(double));
187      model.transposeTimes(-1.0,small.dualRowSolution(),weight);
188      // now massage weight so all basic in plus good djs
189      for (iColumn=0;iColumn<numberColumns;iColumn++) {
190        double dj = weight[iColumn];
191        double value = fullSolution[iColumn];
192        if (model.getColumnStatus(iColumn)==ClpSimplex::basic)
193          dj = -1.0e50;
194        else if (dj<0.0&&value<columnUpper[iColumn])
195          dj = dj;
196        else if (dj>0.0&&value>columnLower[iColumn])
197          dj = -dj;
198        else if (columnUpper[iColumn]>columnLower[iColumn])
199          dj = fabs(dj);
200        else
201          dj = 1.0e50;
202        weight[iColumn] = dj;
203        sort[iColumn] = iColumn;
204      }
205      // sort
206      CoinSort_2(weight,weight+numberColumns,sort);
207      numberSort = smallNumberColumns;
208    }
209  }
210  if (addSlacks) {
211    int i;
212    int numberArtificials = numberColumns-originalNumberColumns;
213    for (i=0;i<numberArtificials;i++)
214      sort[i] = i + originalNumberColumns;
215    model.deleteColumns(numberArtificials,sort);
216  }
217  delete [] weight;
218  delete [] sort;
219  delete [] whichRows;
220  model.primal(1);
221  return 0;
222}   
Note: See TracBrowser for help on using the repository browser.