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

Last change on this file since 1559 was 1559, checked in by stefan, 10 years ago

merge split branch into trunk

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