source: trunk/Clp/examples/testGub.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: 19.9 KB
Line 
1/* $Id: testGub.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 "ClpDynamicExampleMatrix.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 maxFactor = 100;
17     if (argc < 2) {
18          status = model.readMps("../../Data/Netlib/czprob.mps");
19          if (status) {
20               printf("Unable to read matrix - trying gzipped version\n");
21               status = model.readMps("../../Data/Netlib/czprob.mps.gz");
22          }
23     } else {
24          status = model.readMps(argv[1]);
25     }
26     if (status) {
27          printf("errors on input\n");
28          exit(77);
29     }
30     if (argc > 2) {
31          maxFactor = atoi(argv[2]);
32          printf("max factor %d\n", maxFactor);
33     }
34     if (argc > 3) {
35          printf("Using ClpDynamicMatrix\n");
36     } else {
37          printf("Using ClpDynamicExampleMatrix\n");
38     }
39     // find gub
40     int numberRows = model.numberRows();
41     int * gubStart = new int[numberRows+1];
42     int * gubEnd = new int[numberRows];
43     int * which = new int[numberRows];
44     int * whichGub = new int[numberRows];
45     int numberColumns = model.numberColumns();
46     int * mark = new int[numberColumns];
47     int iRow, iColumn;
48     // delete variables fixed to zero
49     const double * columnLower = model.columnLower();
50     const double * columnUpper = model.columnUpper();
51     int numberDelete = 0;
52     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
53          if (columnUpper[iColumn] == 0.0 && columnLower[iColumn] == 0.0)
54               mark[numberDelete++] = iColumn;
55     }
56     if (numberDelete) {
57          model.deleteColumns(numberDelete, mark);
58          numberColumns -= numberDelete;
59          columnLower = model.columnLower();
60          columnUpper = model.columnUpper();
61     }
62     double * lower = new double[numberRows];
63     double * upper = new double[numberRows];
64     const double * rowLower = model.rowLower();
65     const double * rowUpper = model.rowUpper();
66     for (iColumn = 0; iColumn < numberColumns; iColumn++)
67          mark[iColumn] = -1;
68     CoinPackedMatrix * matrix = model.matrix();
69     // get row copy
70     CoinPackedMatrix rowCopy = *matrix;
71     rowCopy.reverseOrdering();
72     const int * column = rowCopy.getIndices();
73     const int * rowLength = rowCopy.getVectorLengths();
74     const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
75     const double * element = rowCopy.getElements();
76     int putGub = numberRows;
77     int putNonGub = numberRows;
78     int * rowIsGub = new int [numberRows];
79     for (iRow = numberRows - 1; iRow >= 0; iRow--) {
80          bool gubRow = true;
81          int first = numberColumns + 1;
82          int last = -1;
83          for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
84               if (element[j] != 1.0) {
85                    gubRow = false;
86                    break;
87               } else {
88                    int iColumn = column[j];
89                    if (mark[iColumn] >= 0) {
90                         gubRow = false;
91                         break;
92                    } else {
93                         last = CoinMax(last, iColumn);
94                         first = CoinMin(first, iColumn);
95                    }
96               }
97          }
98          if (last - first + 1 != rowLength[iRow] || !gubRow) {
99               which[--putNonGub] = iRow;
100               rowIsGub[iRow] = 0;
101          } else {
102               for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
103                    int iColumn = column[j];
104                    mark[iColumn] = iRow;
105               }
106               rowIsGub[iRow] = -1;
107               putGub--;
108               gubStart[putGub] = first;
109               gubEnd[putGub] = last + 1;
110               lower[putGub] = rowLower[iRow];
111               upper[putGub] = rowUpper[iRow];
112               whichGub[putGub] = iRow;
113          }
114     }
115     int numberNonGub = numberRows - putNonGub;
116     int numberGub = numberRows - putGub;
117     if (numberGub > 0) {
118          printf("** %d gub rows\n", numberGub);
119          int numberNormal = 0;
120          const int * row = matrix->getIndices();
121          const int * columnLength = matrix->getVectorLengths();
122          const CoinBigIndex * columnStart = matrix->getVectorStarts();
123          const double * elementByColumn = matrix->getElements();
124          int numberElements = 0;
125          bool doLower = false;
126          bool doUpper = false;
127          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
128               if (mark[iColumn] < 0) {
129                    mark[numberNormal++] = iColumn;
130               } else {
131                    numberElements += columnLength[iColumn];
132                    if (columnLower[iColumn] != 0.0)
133                         doLower = true;
134                    if (columnUpper[iColumn] < 1.0e20)
135                         doUpper = true;
136               }
137          }
138          if (!numberNormal) {
139               printf("Putting back one gub row to make non-empty\n");
140               for (iColumn = gubStart[putGub]; iColumn < gubEnd[putGub]; iColumn++)
141                    mark[numberNormal++] = iColumn;
142               putGub++;
143               numberGub--;
144          }
145          ClpSimplex model2(&model, numberNonGub, which + putNonGub, numberNormal, mark);
146          // and copy for restart test
147          ClpSimplex model3 = model2;
148          int numberGubColumns = numberColumns - numberNormal;
149          // sort gubs so monotonic
150          int * which = new int[numberGub];
151          int i;
152          for (i = 0; i < numberGub; i++)
153               which[i] = i;
154          CoinSort_2(gubStart + putGub, gubStart + putGub + numberGub, which);
155          // move to bottom if we want to use later
156          memmove(gubStart, gubStart + putGub, numberGub * sizeof(int));
157          int * temp1 = new int [numberGub];
158          for (i = 0; i < numberGub; i++) {
159               int k = which[i];
160               temp1[i] = gubEnd[putGub+k];
161          }
162          memcpy(gubEnd, temp1, numberGub * sizeof(int));
163          delete [] temp1;
164          double * temp2 = new double [numberGub];
165          for (i = 0; i < numberGub; i++) {
166               int k = which[i];
167               temp2[i] = lower[putGub+k];
168          }
169          memcpy(lower, temp2, numberGub * sizeof(double));
170          for (i = 0; i < numberGub; i++) {
171               int k = which[i];
172               temp2[i] = upper[putGub+k];
173          }
174          memcpy(upper, temp2, numberGub * sizeof(double));
175          delete [] temp2;
176          delete [] which;
177          numberElements -= numberGubColumns;
178          int * start2 = new int[numberGubColumns+1];
179          int * row2 = new int[numberElements];
180          double * element2 = new double[numberElements];
181          double * cost2 = new double [numberGubColumns];
182          double * lowerColumn2 = NULL;
183          if (doLower) {
184               lowerColumn2 = new double [numberGubColumns];
185               CoinFillN(lowerColumn2, numberGubColumns, 0.0);
186          }
187          double * upperColumn2 = NULL;
188          if (doUpper) {
189               upperColumn2 = new double [numberGubColumns];
190               CoinFillN(upperColumn2, numberGubColumns, COIN_DBL_MAX);
191          }
192          numberElements = 0;
193          int numberNonGubRows = 0;
194          for (iRow = 0; iRow < numberRows; iRow++) {
195               if (!rowIsGub[iRow])
196                    rowIsGub[iRow] = numberNonGubRows++;
197          }
198          numberColumns = 0;
199          int iStart = gubStart[0];
200          gubStart[0] = 0;
201          start2[0] = 0;
202          const double * cost = model.objective();
203          for (int iSet = 0; iSet < numberGub; iSet++) {
204               int iEnd = gubEnd[iSet];
205               for (int k = iStart; k < iEnd; k++) {
206                    cost2[numberColumns] = cost[k];
207                    if (columnLower[k])
208                         lowerColumn2[numberColumns] = columnLower[k];
209                    if (columnUpper[k] < 1.0e20)
210                         upperColumn2[numberColumns] = columnUpper[k];
211                    for (int j = columnStart[k]; j < columnStart[k] + columnLength[k]; j++) {
212                         int iRow = rowIsGub[row[j]];
213                         if (iRow >= 0) {
214                              row2[numberElements] = iRow;
215                              element2[numberElements++] = elementByColumn[j];
216                         }
217                    }
218                    start2[++numberColumns] = numberElements;
219               }
220               if (iSet < numberGub - 1)
221                    iStart = gubStart[iSet+1];
222               gubStart[iSet+1] = numberColumns;
223          }
224          printf("** Before adding matrix there are %d rows and %d columns\n",
225                 model2.numberRows(), model2.numberColumns());
226          if (argc > 3) {
227               ClpDynamicMatrix * newMatrix = new ClpDynamicMatrix(&model2, numberGub,
228                         numberColumns, gubStart,
229                         lower, upper,
230                         start2, row2, element2, cost2,
231                         lowerColumn2, upperColumn2);
232               model2.replaceMatrix(newMatrix);
233               newMatrix->switchOffCheck();
234               newMatrix->setRefreshFrequency(1000);
235          } else {
236               ClpDynamicExampleMatrix * newMatrix = new ClpDynamicExampleMatrix(&model2, numberGub,
237                         numberColumns, gubStart,
238                         lower, upper,
239                         start2, row2, element2, cost2,
240                         lowerColumn2, upperColumn2);
241               model2.replaceMatrix(newMatrix);
242               newMatrix->switchOffCheck();
243               newMatrix->setRefreshFrequency(1000);
244          }
245          printf("** While after adding matrix there are %d rows and %d columns\n",
246                 model2.numberRows(), model2.numberColumns());
247          model2.setSpecialOptions(4); // exactly to bound
248          // For now scaling off
249          model2.scaling(0);
250          ClpPrimalColumnSteepest steepest(5);
251          model2.setPrimalColumnPivotAlgorithm(steepest);
252          //model2.messageHandler()->setLogLevel(63);
253          model2.setFactorizationFrequency(maxFactor);
254          model2.setMaximumIterations(4000000);
255          double time1 = CoinCpuTime();
256          model2.primal();
257          // can't use values pass
258          model2.primal(0);
259          // test proper restart
260          if (argc > 3) {
261               ClpDynamicMatrix * oldMatrix =
262                    dynamic_cast< ClpDynamicMatrix*>(model2.clpMatrix());
263               assert (oldMatrix);
264               ClpDynamicMatrix * newMatrix = new
265               ClpDynamicMatrix(&model3, numberGub,
266                                numberColumns, gubStart,
267                                lower, upper,
268                                start2, row2, element2, cost2,
269                                lowerColumn2, upperColumn2,
270                                oldMatrix->gubRowStatus(), oldMatrix->dynamicStatus());
271               model3.replaceMatrix(newMatrix);
272               // and ordinary status (but only NON gub rows)
273               memcpy(model3.statusArray(), model2.statusArray(),
274                      (newMatrix->numberStaticRows() + model3.numberColumns())*sizeof(unsigned char));
275               newMatrix->switchOffCheck();
276               newMatrix->setRefreshFrequency(1000);
277          } else {
278               ClpDynamicExampleMatrix * oldMatrix =
279                    dynamic_cast< ClpDynamicExampleMatrix*>(model2.clpMatrix());
280               assert (oldMatrix);
281               ClpDynamicExampleMatrix * newMatrix = new
282               ClpDynamicExampleMatrix(&model3, numberGub,
283                                       numberColumns, gubStart,
284                                       lower, upper,
285                                       start2, row2, element2, cost2,
286                                       lowerColumn2, upperColumn2,
287                                       oldMatrix->gubRowStatus(), oldMatrix->dynamicStatus(),
288                                       oldMatrix->numberGubColumns(), oldMatrix->idGen());
289               model3.replaceMatrix(newMatrix);
290               // and ordinary status (but only NON gub rows)
291               memcpy(model3.statusArray(), model2.statusArray(),
292                      (newMatrix->numberStaticRows() + model3.numberColumns())*sizeof(unsigned char));
293               newMatrix->switchOffCheck();
294               newMatrix->setRefreshFrequency(1000);
295          }
296          model3.setSpecialOptions(4); // exactly to bound
297          // For now scaling off
298          model3.scaling(0);
299          model3.setPrimalColumnPivotAlgorithm(steepest);
300          model3.messageHandler()->setLogLevel(63);
301          model3.setFactorizationFrequency(maxFactor);
302          model3.setMaximumIterations(4000000);
303          delete [] rowIsGub;
304          delete [] start2;
305          delete [] row2;
306          delete [] element2;
307          delete [] cost2;
308          delete [] lowerColumn2;
309          delete [] upperColumn2;
310          model3.primal();
311          // this code expects non gub first in original matrix
312          // and only works at present for ClpDynamicMatrix
313          ClpDynamicMatrix * gubMatrix =
314               dynamic_cast< ClpDynamicMatrix*>(model2.clpMatrix());
315          assert (gubMatrix);
316          ClpDynamicExampleMatrix * gubMatrix2 =
317               dynamic_cast< ClpDynamicExampleMatrix*>(model2.clpMatrix());
318          if (!gubMatrix2) {
319               const double * solution = model2.primalColumnSolution();
320               const double * cost = model.objective();
321               int numberGubColumns = gubMatrix->numberGubColumns();
322               int firstOdd = gubMatrix->firstDynamic();
323               int lastOdd = gubMatrix->firstAvailable();
324               int numberTotalColumns = firstOdd + numberGubColumns;
325               int originalNumberRows = model.numberRows();
326               int numberStaticRows = gubMatrix->numberStaticRows();
327               char * status = new char [numberTotalColumns];
328               double * gubSolution = new double [numberTotalColumns];
329               int numberSets = gubMatrix->numberSets();
330               const int * id = gubMatrix->id();
331               int i;
332               const float * columnLower = gubMatrix->columnLower();
333               const float * columnUpper = gubMatrix->columnUpper();
334               for (i = 0; i < numberGubColumns; i++) {
335                    if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) {
336                         gubSolution[i+firstOdd] = columnUpper[i];
337                         status[i+firstOdd] = 2;
338                    } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound && columnLower) {
339                         gubSolution[i+firstOdd] = columnLower[i];
340                         status[i+firstOdd] = 1;
341                    } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) {
342                         int iSet = gubMatrix->whichSet(i);
343                         gubSolution[i+firstOdd] = gubMatrix->keyValue(iSet);
344                         status[i+firstOdd] = 0;
345                    } else {
346                         gubSolution[i+firstOdd] = 0.0;
347                         status[i+firstOdd] = 1;
348                    }
349               }
350               for (i = 0; i < firstOdd; i++) {
351                    ClpSimplex::Status thisStatus = model2.getStatus(i);
352                    if (thisStatus == ClpSimplex::basic)
353                         status[i] = 0;
354                    else if (thisStatus == ClpSimplex::atLowerBound)
355                         status[i] = 1;
356                    else if (thisStatus == ClpSimplex::atUpperBound)
357                         status[i] = 2;
358                    else if (thisStatus == ClpSimplex::isFixed)
359                         status[i] = 3;
360                    else
361                         abort();
362                    gubSolution[i] = solution[i];
363               }
364               for (i = firstOdd; i < lastOdd; i++) {
365                    int iBig = id[i-firstOdd] + firstOdd;
366                    ClpSimplex::Status thisStatus = model2.getStatus(i);
367                    if (thisStatus == ClpSimplex::basic)
368                         status[iBig] = 0;
369                    else if (thisStatus == ClpSimplex::atLowerBound)
370                         status[iBig] = 1;
371                    else if (thisStatus == ClpSimplex::atUpperBound)
372                         status[iBig] = 2;
373                    else if (thisStatus == ClpSimplex::isFixed)
374                         status[iBig] = 3;
375                    else
376                         abort();
377                    gubSolution[iBig] = solution[i];
378               }
379               char * rowStatus = new char[originalNumberRows];
380               for (i = 0; i < numberStaticRows; i++) {
381                    ClpSimplex::Status thisStatus = model2.getRowStatus(i);
382                    if (thisStatus == ClpSimplex::basic)
383                         rowStatus[i] = 0;
384                    else if (thisStatus == ClpSimplex::atLowerBound)
385                         rowStatus[i] = 1;
386                    else if (thisStatus == ClpSimplex::atUpperBound)
387                         rowStatus[i] = 2;
388                    else if (thisStatus == ClpSimplex::isFixed)
389                         rowStatus[i] = 3;
390                    else
391                         abort();
392               }
393               double objValue = 0.0;
394               for (i = 0; i < numberTotalColumns; i++)
395                    objValue += cost[i] * gubSolution[i];
396               printf("objective value is %g\n", objValue);
397               for (i = 0; i < numberSets; i++) {
398                    ClpSimplex::Status thisStatus = gubMatrix->getStatus(i);
399                    if (thisStatus == ClpSimplex::basic)
400                         rowStatus[numberStaticRows+i] = 0;
401                    else if (thisStatus == ClpSimplex::atLowerBound)
402                         rowStatus[numberStaticRows+i] = 1;
403                    else if (thisStatus == ClpSimplex::atUpperBound)
404                         rowStatus[numberStaticRows+i] = 2;
405                    else if (thisStatus == ClpSimplex::isFixed)
406                         rowStatus[numberStaticRows+i] = 3;
407                    else
408                         abort();
409               }
410               // Coding below may not work if gub rows not at end
411               FILE * fp = fopen ("xx.sol", "w");
412               fwrite(gubSolution, sizeof(double), numberTotalColumns, fp);
413               fwrite(status, sizeof(char), numberTotalColumns, fp);
414               const double * rowsol = model2.primalRowSolution();
415               double * rowsol2 = new double[originalNumberRows];
416               memset(rowsol2, 0, originalNumberRows * sizeof(double));
417               model.times(1.0, gubSolution, rowsol2);
418               for (i = 0; i < numberStaticRows; i++)
419                    assert (fabs(rowsol[i] - rowsol2[i]) < 1.0e-3);
420               for (; i < originalNumberRows; i++)
421                    assert (rowsol2[i] > lower[i-numberStaticRows] - 1.0e-3 &&
422                            rowsol2[i] < upper[i-numberStaticRows] + 1.0e-3);
423               //for (;i<originalNumberRows;i++)
424               //printf("%d %g\n",i,rowsol2[i]);
425               fwrite(rowsol2, sizeof(double), originalNumberRows, fp);
426               delete [] rowsol2;
427               fwrite(rowStatus, sizeof(char), originalNumberRows, fp);
428               fclose(fp);
429               delete [] status;
430               delete [] gubSolution;
431               // ** if going to rstart as dynamic need id_
432               // also copy coding in useEf.. from ClpGubMatrix (i.e. test for basis)
433          }
434          printf("obj offset is %g\n", model2.objectiveOffset());
435          printf("Primal took %g seconds\n", CoinCpuTime() - time1);
436     }
437     delete [] mark;
438     delete [] gubStart;
439     delete [] gubEnd;
440     delete [] which;
441     delete [] whichGub;
442     delete [] lower;
443     delete [] upper;
444     return 0;
445}
Note: See TracBrowser for help on using the repository browser.