source: stable/1.13/Clp/examples/testGub2.cpp @ 1625

Last change on this file since 1625 was 1561, checked in by stefan, 10 years ago

make check for sample dir work without including ClpConfig?.h

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