source: trunk/Clp/examples/sprint.cpp

Last change on this file was 2278, checked in by forrest, 12 months ago

COIN_BIG_INDEX 2 changes

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