1 /* $Id: testGub2.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 "ClpGubDynamicMatrix.hpp"
8 #include "ClpPrimalColumnSteepest.hpp"
9 #include "CoinSort.hpp"
10 #include "CoinHelperFunctions.hpp"
11 #include "CoinTime.hpp"
12 #include "CoinMpsIO.hpp"
main(int argc,const char * argv[])13 int main(int argc, const char *argv[])
14 {
15 #if COIN_BIG_INDEX<2
16 ClpSimplex model;
17 int status;
18 int maxIts = 0;
19 int maxFactor = 100;
20 if (argc < 2) {
21 #if defined(SAMPLEDIR)
22 status = model.readMps(SAMPLEDIR "/p0033.mps", true);
23 #else
24 fprintf(stderr, "Do not know where to find sample MPS files.\n");
25 exit(1);
26 #endif
27 } else
28 status = model.readMps(argv[1]);
29 if (status) {
30 printf("errors on input\n");
31 exit(77);
32 }
33 if (argc > 2) {
34 maxFactor = atoi(argv[2]);
35 printf("max factor %d\n", maxFactor);
36 }
37 if (argc > 3) {
38 maxIts = atoi(argv[3]);
39 printf("max its %d\n", maxIts);
40 }
41 // For now scaling off
42 model.scaling(0);
43 if (maxIts) {
44 // Do partial dantzig
45 ClpPrimalColumnSteepest dantzig(5);
46 model.setPrimalColumnPivotAlgorithm(dantzig);
47 //model.messageHandler()->setLogLevel(63);
48 model.setFactorizationFrequency(maxFactor);
49 model.setMaximumIterations(maxIts);
50 model.primal();
51 if (!model.status())
52 exit(1);
53 }
54 // find gub
55 int numberRows = model.numberRows();
56 int * gubStart = new int[numberRows+1];
57 int * gubEnd = new int[numberRows];
58 int * which = new int[numberRows];
59 int * whichGub = new int[numberRows];
60 int numberColumns = model.numberColumns();
61 int * mark = new int[numberColumns];
62 int iRow, iColumn;
63 // delete variables fixed to zero
64 const double * columnLower = model.columnLower();
65 const double * columnUpper = model.columnUpper();
66 int numberDelete = 0;
67 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
68 if (columnUpper[iColumn] == 0.0 && columnLower[iColumn] == 0.0)
69 mark[numberDelete++] = iColumn;
70 }
71 if (numberDelete) {
72 model.deleteColumns(numberDelete, mark);
73 numberColumns -= numberDelete;
74 columnLower = model.columnLower();
75 columnUpper = model.columnUpper();
76 #if 0
77 CoinMpsIO writer;
78 writer.setMpsData(*model.matrix(), COIN_DBL_MAX,
79 model.getColLower(), model.getColUpper(),
80 model.getObjCoefficients(),
81 (const char*) 0 /*integrality*/,
82 model.getRowLower(), model.getRowUpper(),
83 NULL, NULL);
84 writer.writeMps("cza.mps", 0, 0, 1);
85 #endif
86 }
87 double * lower = new double[numberRows];
88 double * upper = new double[numberRows];
89 const double * rowLower = model.rowLower();
90 const double * rowUpper = model.rowUpper();
91 for (iColumn = 0; iColumn < numberColumns; iColumn++)
92 mark[iColumn] = -1;
93 CoinPackedMatrix * matrix = model.matrix();
94 // get row copy
95 CoinPackedMatrix rowCopy = *matrix;
96 rowCopy.reverseOrdering();
97 const int * column = rowCopy.getIndices();
98 const int * rowLength = rowCopy.getVectorLengths();
99 const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
100 const double * element = rowCopy.getElements();
101 int putGub = numberRows;
102 int putNonGub = numberRows;
103 int * rowIsGub = new int [numberRows];
104 for (iRow = numberRows - 1; iRow >= 0; iRow--) {
105 bool gubRow = true;
106 int first = numberColumns + 1;
107 int last = -1;
108 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
109 if (element[j] != 1.0) {
110 gubRow = false;
111 break;
112 } else {
113 int iColumn = column[j];
114 if (mark[iColumn] >= 0) {
115 gubRow = false;
116 break;
117 } else {
118 last = CoinMax(last, iColumn);
119 first = CoinMin(first, iColumn);
120 }
121 }
122 }
123 if (last - first + 1 != rowLength[iRow] || !gubRow) {
124 which[--putNonGub] = iRow;
125 rowIsGub[iRow] = 0;
126 } else {
127 for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
128 int iColumn = column[j];
129 mark[iColumn] = iRow;
130 }
131 rowIsGub[iRow] = -1;
132 putGub--;
133 gubStart[putGub] = first;
134 gubEnd[putGub] = last + 1;
135 lower[putGub] = rowLower[iRow];
136 upper[putGub] = rowUpper[iRow];
137 whichGub[putGub] = iRow;
138 }
139 }
140 int numberNonGub = numberRows - putNonGub;
141 int numberGub = numberRows - putGub;
142 if (numberGub > 0) {
143 printf("** %d gub rows\n", numberGub);
144 int numberNormal = 0;
145 const int * row = matrix->getIndices();
146 const int * columnLength = matrix->getVectorLengths();
147 const CoinBigIndex * columnStart = matrix->getVectorStarts();
148 const double * elementByColumn = matrix->getElements();
149 int numberElements = 0;
150 bool doLower = false;
151 bool doUpper = false;
152 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
153 if (mark[iColumn] < 0) {
154 mark[numberNormal++] = iColumn;
155 } else {
156 numberElements += columnLength[iColumn];
157 if (columnLower[iColumn] != 0.0)
158 doLower = true;
159 if (columnUpper[iColumn] < 1.0e20)
160 doUpper = true;
161 }
162 }
163 if (!numberNormal) {
164 printf("Putting back one gub row to make non-empty\n");
165 for (iColumn = gubStart[putGub]; iColumn < gubEnd[putGub]; iColumn++)
166 mark[numberNormal++] = iColumn;
167 putGub++;
168 numberGub--;
169 }
170 ClpSimplex model2(&model, numberNonGub, which + putNonGub, numberNormal, mark);
171 int numberGubColumns = numberColumns - numberNormal;
172 // sort gubs so monotonic
173 int * which = new int[numberGub];
174 int i;
175 for (i = 0; i < numberGub; i++)
176 which[i] = i;
177 CoinSort_2(gubStart + putGub, gubStart + putGub + numberGub, which);
178 int * temp1 = new int [numberGub];
179 for (i = 0; i < numberGub; i++) {
180 int k = which[i];
181 temp1[i] = gubEnd[putGub+k];
182 }
183 memcpy(gubEnd + putGub, temp1, numberGub * sizeof(int));
184 delete [] temp1;
185 double * temp2 = new double [numberGub];
186 for (i = 0; i < numberGub; i++) {
187 int k = which[i];
188 temp2[i] = lower[putGub+k];
189 }
190 memcpy(lower + putGub, temp2, numberGub * sizeof(double));
191 for (i = 0; i < numberGub; i++) {
192 int k = which[i];
193 temp2[i] = upper[putGub+k];
194 }
195 memcpy(upper + putGub, temp2, numberGub * sizeof(double));
196 delete [] temp2;
197 delete [] which;
198 numberElements -= numberGubColumns;
199 int * start2 = new int[numberGubColumns+1];
200 int * row2 = new int[numberElements];
201 double * element2 = new double[numberElements];
202 double * cost2 = new double [numberGubColumns];
203 double * lowerColumn2 = NULL;
204 if (doLower) {
205 lowerColumn2 = new double [numberGubColumns];
206 CoinFillN(lowerColumn2, numberGubColumns, 0.0);
207 }
208 double * upperColumn2 = NULL;
209 if (doUpper) {
210 upperColumn2 = new double [numberGubColumns];
211 CoinFillN(upperColumn2, numberGubColumns, COIN_DBL_MAX);
212 }
213 numberElements = 0;
214 int numberNonGubRows = 0;
215 for (iRow = 0; iRow < numberRows; iRow++) {
216 if (!rowIsGub[iRow])
217 rowIsGub[iRow] = numberNonGubRows++;
218 }
219 numberColumns = 0;
220 gubStart[0] = 0;
221 start2[0] = 0;
222 const double * cost = model.objective();
223 for (int iSet = 0; iSet < numberGub; iSet++) {
224 int iStart = gubStart[iSet+putGub];
225 int iEnd = gubEnd[iSet+putGub];
226 for (int k = iStart; k < iEnd; k++) {
227 cost2[numberColumns] = cost[k];
228 if (columnLower[k])
229 lowerColumn2[numberColumns] = columnLower[k];
230 if (columnUpper[k] < 1.0e20)
231 upperColumn2[numberColumns] = columnUpper[k];
232 for (int j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
233 int iRow = rowIsGub[row[j]];
234 if (iRow >= 0) {
235 row2[numberElements] = iRow;
236 element2[numberElements++] = elementByColumn[j];
237 }
238 }
239 start2[++numberColumns] = numberElements;
240 }
241 gubStart[iSet+1] = numberColumns;
242 }
243 model2.replaceMatrix(new ClpGubDynamicMatrix(&model2, numberGub,
244 numberColumns, gubStart,
245 lower + putGub, upper + putGub,
246 start2, row2, element2, cost2,
247 lowerColumn2, upperColumn2));
248 delete [] rowIsGub;
249 delete [] start2;
250 delete [] row2;
251 delete [] element2;
252 delete [] cost2;
253 delete [] lowerColumn2;
254 delete [] upperColumn2;
255 // For now scaling off
256 model2.scaling(0);
257 // Do partial dantzig
258 ClpPrimalColumnSteepest dantzig(5);
259 model2.setPrimalColumnPivotAlgorithm(dantzig);
260 //model2.messageHandler()->setLogLevel(63);
261 model2.setFactorizationFrequency(maxFactor);
262 model2.setMaximumIterations(4000000);
263 double time1 = CoinCpuTime();
264 model2.primal();
265 {
266 ClpGubDynamicMatrix * gubMatrix =
267 dynamic_cast< ClpGubDynamicMatrix*>(model2.clpMatrix());
268 assert(gubMatrix);
269 const double * solution = model2.primalColumnSolution();
270 int numberGubColumns = gubMatrix->numberGubColumns();
271 int firstOdd = gubMatrix->firstDynamic();
272 int lastOdd = gubMatrix->firstAvailable();
273 int numberTotalColumns = firstOdd + numberGubColumns;
274 int numberRows = model2.numberRows();
275 char * status = new char [numberTotalColumns];
276 double * gubSolution = new double [numberTotalColumns];
277 int numberSets = gubMatrix->numberSets();
278 const int * id = gubMatrix->id();
279 int i;
280 const double * lowerColumn = gubMatrix->lowerColumn();
281 const double * upperColumn = gubMatrix->upperColumn();
282 for (i = 0; i < numberGubColumns; i++) {
283 if (gubMatrix->getDynamicStatus(i) == ClpGubDynamicMatrix::atUpperBound) {
284 gubSolution[i+firstOdd] = upperColumn[i];
285 status[i+firstOdd] = 2;
286 } else if (gubMatrix->getDynamicStatus(i) == ClpGubDynamicMatrix::atLowerBound && lowerColumn) {
287 gubSolution[i+firstOdd] = lowerColumn[i];
288 status[i+firstOdd] = 1;
289 } else {
290 gubSolution[i+firstOdd] = 0.0;
291 status[i+firstOdd] = 1;
292 }
293 }
294 for (i = 0; i < firstOdd; i++) {
295 ClpSimplex::Status thisStatus = model2.getStatus(i);
296 if (thisStatus == ClpSimplex::basic)
297 status[i] = 0;
298 else if (thisStatus == ClpSimplex::atLowerBound)
299 status[i] = 1;
300 else if (thisStatus == ClpSimplex::atUpperBound)
301 status[i] = 2;
302 else if (thisStatus == ClpSimplex::isFixed)
303 status[i] = 3;
304 else
305 abort();
306 gubSolution[i] = solution[i];
307 }
308 for (i = firstOdd; i < lastOdd; i++) {
309 int iBig = id[i-firstOdd] + firstOdd;
310 ClpSimplex::Status thisStatus = model2.getStatus(i);
311 if (thisStatus == ClpSimplex::basic)
312 status[iBig] = 0;
313 else if (thisStatus == ClpSimplex::atLowerBound)
314 status[iBig] = 1;
315 else if (thisStatus == ClpSimplex::atUpperBound)
316 status[iBig] = 2;
317 else if (thisStatus == ClpSimplex::isFixed)
318 status[iBig] = 3;
319 else
320 abort();
321 gubSolution[iBig] = solution[i];
322 }
323 char * rowStatus = new char[numberRows];
324 for (i = 0; i < numberRows; i++) {
325 ClpSimplex::Status thisStatus = model2.getRowStatus(i);
326 if (thisStatus == ClpSimplex::basic)
327 rowStatus[i] = 0;
328 else if (thisStatus == ClpSimplex::atLowerBound)
329 rowStatus[i] = 1;
330 else if (thisStatus == ClpSimplex::atUpperBound)
331 rowStatus[i] = 2;
332 else if (thisStatus == ClpSimplex::isFixed)
333 rowStatus[i] = 3;
334 else
335 abort();
336 }
337 char * setStatus = new char[numberSets];
338 int * keyVariable = new int[numberSets];
339 memcpy(keyVariable, gubMatrix->keyVariable(), numberSets * sizeof(int));
340 for (i = 0; i < numberSets; i++) {
341 int iKey = keyVariable[i];
342 if (iKey > lastOdd)
343 iKey = numberTotalColumns + i;
344 else
345 iKey = id[iKey-firstOdd] + firstOdd;
346 keyVariable[i] = iKey;
347 ClpSimplex::Status thisStatus = gubMatrix->getStatus(i);
348 if (thisStatus == ClpSimplex::basic)
349 setStatus[i] = 0;
350 else if (thisStatus == ClpSimplex::atLowerBound)
351 setStatus[i] = 1;
352 else if (thisStatus == ClpSimplex::atUpperBound)
353 setStatus[i] = 2;
354 else if (thisStatus == ClpSimplex::isFixed)
355 setStatus[i] = 3;
356 else
357 abort();
358 }
359 FILE * fp = fopen("xx.sol", "w");
360 fwrite(gubSolution, sizeof(double), numberTotalColumns, fp);
361 fwrite(status, sizeof(char), numberTotalColumns, fp);
362 const double * rowsol = model2.primalRowSolution();
363 int originalNumberRows = model.numberRows();
364 double * rowsol2 = new double[originalNumberRows];
365 memset(rowsol2, 0, originalNumberRows * sizeof(double));
366 model.times(1.0, gubSolution, rowsol2);
367 for (i = 0; i < numberRows; i++)
368 assert(fabs(rowsol[i] - rowsol2[i]) < 1.0e-3);
369 //for (;i<originalNumberRows;i++)
370 //printf("%d %g\n",i,rowsol2[i]);
371 delete [] rowsol2;
372 fwrite(rowsol, sizeof(double), numberRows, fp);
373 fwrite(rowStatus, sizeof(char), numberRows, fp);
374 fwrite(setStatus, sizeof(char), numberSets, fp);
375 fwrite(keyVariable, sizeof(int), numberSets, fp);
376 fclose(fp);
377 delete [] status;
378 delete [] gubSolution;
379 delete [] setStatus;
380 delete [] keyVariable;
381 // ** if going to rstart as dynamic need id_
382 // also copy coding in useEf.. from ClpGubMatrix (i.e. test for basis)
383 }
384 printf("obj offset is %g\n", model2.objectiveOffset());
385 printf("Primal took %g seconds\n", CoinCpuTime() - time1);
386 //model2.primal(1);
387 }
388 delete [] mark;
389 delete [] gubStart;
390 delete [] gubEnd;
391 delete [] which;
392 delete [] whichGub;
393 delete [] lower;
394 delete [] upper;
395 #else
396 printf("testGub2 not available with COIN_BIG_INDEX=2\n");
397 #endif
398 return 0;
399 }
400