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