source: stable/1.15/Clp/examples/testGub.cpp @ 1941

Last change on this file since 1941 was 1941, checked in by stefan, 6 years ago

sync with trunk rev 1940

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