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
main(int argc,const char * argv[])10 int 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 }
227