source: trunk/Clp/examples/testGub.cpp

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

COIN_BIG_INDEX 2 changes

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