source: trunk/Clp/examples/testGub2.cpp

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

COIN_BIG_INDEX 2 changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 16.7 KB
Line 
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"
13int 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}
Note: See TracBrowser for help on using the repository browser.