source: trunk/Clp/examples/testGub2.cpp @ 1552

Last change on this file since 1552 was 1552, checked in by mjs, 11 years ago

Format examples with 'astyle -A4 -p'.

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