source: stable/1.13/Clp/examples/decompose.cpp @ 1898

Last change on this file since 1898 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: 22.2 KB
Line 
1/* $Id: decompose.cpp 1561 2010-06-13 13:08:03Z stefan $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "ClpSimplex.hpp"
6#include "CoinMpsIO.hpp"
7#include <iomanip>
8
9int main(int argc, const char *argv[])
10{
11     ClpSimplex  model;
12     int status;
13     // Keep names
14     if (argc < 2) {
15#if defined(NETLIBDIR)
16          status = model.readMps(NETLIBDIR "/czprob.mps", true);
17#else
18          fprintf(stderr, "Do not know where to find netlib MPS files.\n");
19          return 1;
20#endif
21     } else {
22          status = model.readMps(argv[1], true);
23     }
24     if (status)
25          exit(10);
26     /*
27       This driver does a simple Dantzig Wolfe decomposition
28     */
29     // Get master rows in some mysterious way
30     int numberRows = model.numberRows();
31     int * rowBlock = new int[numberRows];
32     int iRow;
33     for (iRow = 0; iRow < numberRows; iRow++)
34          rowBlock[iRow] = -2;
35     // these are master rows
36     if (numberRows == 105127) {
37          // ken-18
38          for (iRow = 104976; iRow < numberRows; iRow++)
39               rowBlock[iRow] = -1;
40     } else if (numberRows == 2426) {
41          // ken-7
42          for (iRow = 2401; iRow < numberRows; iRow++)
43               rowBlock[iRow] = -1;
44     } else if (numberRows == 810) {
45          for (iRow = 81; iRow < 84; iRow++)
46               rowBlock[iRow] = -1;
47     } else if (numberRows == 5418) {
48          for (iRow = 564; iRow < 603; iRow++)
49               rowBlock[iRow] = -1;
50     } else if (numberRows == 10280) {
51          // osa-60
52          for (iRow = 10198; iRow < 10280; iRow++)
53               rowBlock[iRow] = -1;
54     } else if (numberRows == 1503) {
55          // degen3
56          for (iRow = 0; iRow < 561; iRow++)
57               rowBlock[iRow] = -1;
58     } else if (numberRows == 929) {
59          // czprob
60          for (iRow = 0; iRow < 39; iRow++)
61               rowBlock[iRow] = -1;
62     }
63     CoinPackedMatrix * matrix = model.matrix();
64     // get row copy
65     CoinPackedMatrix rowCopy = *matrix;
66     rowCopy.reverseOrdering();
67     const int * row = matrix->getIndices();
68     const int * columnLength = matrix->getVectorLengths();
69     const CoinBigIndex * columnStart = matrix->getVectorStarts();
70     //const double * elementByColumn = matrix->getElements();
71     const int * column = rowCopy.getIndices();
72     const int * rowLength = rowCopy.getVectorLengths();
73     const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
74     //const double * elementByRow = rowCopy.getElements();
75     int numberBlocks = 0;
76     int * stack = new int [numberRows];
77     // to say if column looked at
78     int numberColumns = model.numberColumns();
79     int * columnBlock = new int[numberColumns];
80     int iColumn;
81     for (iColumn = 0; iColumn < numberColumns; iColumn++)
82          columnBlock[iColumn] = -2;
83     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
84          int kstart = columnStart[iColumn];
85          int kend = columnStart[iColumn] + columnLength[iColumn];
86          if (columnBlock[iColumn] == -2) {
87               // column not allocated
88               int j;
89               int nstack = 0;
90               for (j = kstart; j < kend; j++) {
91                    int iRow = row[j];
92                    if (rowBlock[iRow] != -1) {
93                         assert(rowBlock[iRow] == -2);
94                         rowBlock[iRow] = numberBlocks; // mark
95                         stack[nstack++] = iRow;
96                    }
97               }
98               if (nstack) {
99                    // new block - put all connected in
100                    numberBlocks++;
101                    columnBlock[iColumn] = numberBlocks - 1;
102                    while (nstack) {
103                         int iRow = stack[--nstack];
104                         int k;
105                         for (k = rowStart[iRow]; k < rowStart[iRow] + rowLength[iRow]; k++) {
106                              int iColumn = column[k];
107                              int kkstart = columnStart[iColumn];
108                              int kkend = kkstart + columnLength[iColumn];
109                              if (columnBlock[iColumn] == -2) {
110                                   columnBlock[iColumn] = numberBlocks - 1; // mark
111                                   // column not allocated
112                                   int jj;
113                                   for (jj = kkstart; jj < kkend; jj++) {
114                                        int jRow = row[jj];
115                                        if (rowBlock[jRow] == -2) {
116                                             rowBlock[jRow] = numberBlocks - 1;
117                                             stack[nstack++] = jRow;
118                                        }
119                                   }
120                              } else {
121                                   assert(columnBlock[iColumn] == numberBlocks - 1);
122                              }
123                         }
124                    }
125               } else {
126                    // Only in master
127                    columnBlock[iColumn] = -1;
128               }
129          }
130     }
131     printf("%d blocks found\n", numberBlocks);
132     if (numberBlocks > 50) {
133          int iBlock;
134          for (iRow = 0; iRow < numberRows; iRow++) {
135               iBlock = rowBlock[iRow];
136               if (iBlock >= 0)
137                    rowBlock[iRow] = iBlock % 50;
138          }
139          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
140               iBlock = columnBlock[iColumn];
141               if (iBlock >= 0)
142                    columnBlock[iColumn] = iBlock % 50;
143          }
144          numberBlocks = 50;
145     }
146     delete [] stack;
147     // make up problems
148     CoinPackedMatrix * top = new CoinPackedMatrix [numberBlocks];
149     ClpSimplex * sub = new ClpSimplex [numberBlocks];
150     ClpSimplex master;
151     // Create all sub problems
152     // Could do much faster - but do that later
153     int * whichRow = new int [numberRows];
154     int * whichColumn = new int [numberColumns];
155     // get top matrix
156     CoinPackedMatrix topMatrix = *model.matrix();
157     int numberRow2, numberColumn2;
158     numberRow2 = 0;
159     for (iRow = 0; iRow < numberRows; iRow++)
160          if (rowBlock[iRow] >= 0)
161               whichRow[numberRow2++] = iRow;
162     topMatrix.deleteRows(numberRow2, whichRow);
163     int iBlock;
164     for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
165          numberRow2 = 0;
166          numberColumn2 = 0;
167          for (iRow = 0; iRow < numberRows; iRow++)
168               if (iBlock == rowBlock[iRow])
169                    whichRow[numberRow2++] = iRow;
170          for (iColumn = 0; iColumn < numberColumns; iColumn++)
171               if (iBlock == columnBlock[iColumn])
172                    whichColumn[numberColumn2++] = iColumn;
173          sub[iBlock] = ClpSimplex(&model, numberRow2, whichRow,
174                                   numberColumn2, whichColumn);
175#if 0
176          // temp
177          double * upper = sub[iBlock].columnUpper();
178          for (iColumn = 0; iColumn < numberColumn2; iColumn++)
179               upper[iColumn] = 100.0;
180#endif
181          // and top matrix
182          CoinPackedMatrix matrix = topMatrix;
183          // and delete bits
184          numberColumn2 = 0;
185          for (iColumn = 0; iColumn < numberColumns; iColumn++)
186               if (iBlock != columnBlock[iColumn])
187                    whichColumn[numberColumn2++] = iColumn;
188          matrix.deleteCols(numberColumn2, whichColumn);
189          top[iBlock] = matrix;
190     }
191     // and master
192     numberRow2 = 0;
193     numberColumn2 = 0;
194     for (iRow = 0; iRow < numberRows; iRow++)
195          if (rowBlock[iRow] < 0)
196               whichRow[numberRow2++] = iRow;
197     for (iColumn = 0; iColumn < numberColumns; iColumn++)
198          if (columnBlock[iColumn] == -1)
199               whichColumn[numberColumn2++] = iColumn;
200     ClpModel masterModel(&model, numberRow2, whichRow,
201                          numberColumn2, whichColumn);
202     master = ClpSimplex(masterModel);
203     delete [] whichRow;
204     delete [] whichColumn;
205
206     // Overkill in terms of space
207     int numberMasterRows = master.numberRows();
208     int * columnAdd = new int[numberBlocks+1];
209     int * rowAdd = new int[numberBlocks*(numberMasterRows+1)];
210     double * elementAdd = new double[numberBlocks*(numberMasterRows+1)];
211     double * objective = new double[numberBlocks];
212     int maxPass = 500;
213     int iPass;
214     double lastObjective = 1.0e31;
215     // Create convexity rows for proposals
216     int numberMasterColumns = master.numberColumns();
217     master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
218     // Arrays to say which block and when created
219     int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
220     int * whichBlock = new int[maximumColumns];
221     int * when = new int[maximumColumns];
222     int numberColumnsGenerated = numberBlocks;
223     // fill in rhs and add in artificials
224     {
225          double * rowLower = master.rowLower();
226          double * rowUpper = master.rowUpper();
227          int iBlock;
228          columnAdd[0] = 0;
229          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
230               int iRow = iBlock + numberMasterRows;;
231               rowLower[iRow] = 1.0;
232               rowUpper[iRow] = 1.0;
233               rowAdd[iBlock] = iRow;
234               elementAdd[iBlock] = 1.0;
235               objective[iBlock] = 1.0e9;
236               columnAdd[iBlock+1] = iBlock + 1;
237               when[iBlock] = -1;
238               whichBlock[iBlock] = iBlock;
239          }
240          master.addColumns(numberBlocks, NULL, NULL, objective,
241                            columnAdd, rowAdd, elementAdd);
242     }
243     // and resize matrix to double check clp will be happy
244     //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
245     //                  numberMasterColumns+numberBlocks);
246
247     for (iPass = 0; iPass < maxPass; iPass++) {
248          printf("Start of pass %d\n", iPass);
249          // Solve master - may be infeasible
250          master.scaling(false);
251          if (0) {
252               master.writeMps("yy.mps");
253          }
254
255          master.primal();
256          int problemStatus = master.status(); // do here as can change (delcols)
257          if (master.numberIterations() == 0 && iPass)
258               break; // finished
259          if (master.objectiveValue() > lastObjective - 1.0e-7 && iPass > 555)
260               break; // finished
261          lastObjective = master.objectiveValue();
262          // mark basic ones and delete if necessary
263          int iColumn;
264          numberColumnsGenerated = master.numberColumns() - numberMasterColumns;
265          for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
266               if (master.getStatus(iColumn + numberMasterColumns) == ClpSimplex::basic)
267                    when[iColumn] = iPass;
268          }
269          if (numberColumnsGenerated + numberBlocks > maximumColumns) {
270               // delete
271               int numberKeep = 0;
272               int numberDelete = 0;
273               int * whichDelete = new int[numberColumns];
274               for (iColumn = 0; iColumn < numberColumnsGenerated; iColumn++) {
275                    if (when[iColumn] > iPass - 7) {
276                         // keep
277                         when[numberKeep] = when[iColumn];
278                         whichBlock[numberKeep++] = whichBlock[iColumn];
279                    } else {
280                         // delete
281                         whichDelete[numberDelete++] = iColumn + numberMasterColumns;
282                    }
283               }
284               numberColumnsGenerated -= numberDelete;
285               master.deleteColumns(numberDelete, whichDelete);
286               delete [] whichDelete;
287          }
288          const double * dual = NULL;
289          bool deleteDual = false;
290          if (problemStatus == 0) {
291               dual = master.dualRowSolution();
292          } else if (problemStatus == 1) {
293               // could do composite objective
294               dual = master.infeasibilityRay();
295               deleteDual = true;
296               printf("The sum of infeasibilities is %g\n",
297                      master.sumPrimalInfeasibilities());
298          } else if (!master.numberColumns()) {
299               assert(!iPass);
300               dual = master.dualRowSolution();
301               memset(master.dualRowSolution(),
302                      0, (numberMasterRows + numberBlocks) *sizeof(double));
303          } else {
304               abort();
305          }
306          // Create objective for sub problems and solve
307          columnAdd[0] = 0;
308          int numberProposals = 0;
309          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
310               int numberColumns2 = sub[iBlock].numberColumns();
311               double * saveObj = new double [numberColumns2];
312               double * objective2 = sub[iBlock].objective();
313               memcpy(saveObj, objective2, numberColumns2 * sizeof(double));
314               // new objective
315               top[iBlock].transposeTimes(dual, objective2);
316               int i;
317               if (problemStatus == 0) {
318                    for (i = 0; i < numberColumns2; i++)
319                         objective2[i] = saveObj[i] - objective2[i];
320               } else {
321                    for (i = 0; i < numberColumns2; i++)
322                         objective2[i] = -objective2[i];
323               }
324               sub[iBlock].primal();
325               memcpy(objective2, saveObj, numberColumns2 * sizeof(double));
326               // get proposal
327               if (sub[iBlock].numberIterations() || !iPass) {
328                    double objValue = 0.0;
329                    int start = columnAdd[numberProposals];
330                    // proposal
331                    if (sub[iBlock].isProvenOptimal()) {
332                         const double * solution = sub[iBlock].primalColumnSolution();
333                         top[iBlock].times(solution, elementAdd + start);
334                         for (i = 0; i < numberColumns2; i++)
335                              objValue += solution[i] * saveObj[i];
336                         // See if good dj and pack down
337                         int number = start;
338                         double dj = objValue;
339                         if (problemStatus)
340                              dj = 0.0;
341                         double smallest = 1.0e100;
342                         double largest = 0.0;
343                         for (i = 0; i < numberMasterRows; i++) {
344                              double value = elementAdd[start+i];
345                              if (fabs(value) > 1.0e-15) {
346                                   dj -= dual[i] * value;
347                                   smallest = CoinMin(smallest, fabs(value));
348                                   largest = CoinMax(largest, fabs(value));
349                                   rowAdd[number] = i;
350                                   elementAdd[number++] = value;
351                              }
352                         }
353                         // and convexity
354                         dj -= dual[numberMasterRows+iBlock];
355                         rowAdd[number] = numberMasterRows + iBlock;
356                         elementAdd[number++] = 1.0;
357                         // if elements large then scale?
358                         //if (largest>1.0e8||smallest<1.0e-8)
359                         printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
360                                iBlock, smallest, largest, dj);
361                         if (dj < -1.0e-6 || !iPass) {
362                              // take
363                              objective[numberProposals] = objValue;
364                              columnAdd[++numberProposals] = number;
365                              when[numberColumnsGenerated] = iPass;
366                              whichBlock[numberColumnsGenerated++] = iBlock;
367                         }
368                    } else if (sub[iBlock].isProvenDualInfeasible()) {
369                         // use ray
370                         const double * solution = sub[iBlock].unboundedRay();
371                         top[iBlock].times(solution, elementAdd + start);
372                         for (i = 0; i < numberColumns2; i++)
373                              objValue += solution[i] * saveObj[i];
374                         // See if good dj and pack down
375                         int number = start;
376                         double dj = objValue;
377                         double smallest = 1.0e100;
378                         double largest = 0.0;
379                         for (i = 0; i < numberMasterRows; i++) {
380                              double value = elementAdd[start+i];
381                              if (fabs(value) > 1.0e-15) {
382                                   dj -= dual[i] * value;
383                                   smallest = CoinMin(smallest, fabs(value));
384                                   largest = CoinMax(largest, fabs(value));
385                                   rowAdd[number] = i;
386                                   elementAdd[number++] = value;
387                              }
388                         }
389                         // if elements large or small then scale?
390                         //if (largest>1.0e8||smallest<1.0e-8)
391                         printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
392                                iBlock, smallest, largest, dj);
393                         if (dj < -1.0e-6) {
394                              // take
395                              objective[numberProposals] = objValue;
396                              columnAdd[++numberProposals] = number;
397                              when[numberColumnsGenerated] = iPass;
398                              whichBlock[numberColumnsGenerated++] = iBlock;
399                         }
400                    } else {
401                         abort();
402                    }
403               }
404               delete [] saveObj;
405          }
406          if (deleteDual)
407               delete [] dual;
408          if (numberProposals)
409               master.addColumns(numberProposals, NULL, NULL, objective,
410                                 columnAdd, rowAdd, elementAdd);
411     }
412     // now put back a good solution
413     double * lower = new double[numberMasterRows+numberBlocks];
414     double * upper = new double[numberMasterRows+numberBlocks];
415     numberColumnsGenerated  += numberMasterColumns;
416     double * sol = new double[numberColumnsGenerated];
417     const double * solution = master.primalColumnSolution();
418     const double * masterLower = master.rowLower();
419     const double * masterUpper = master.rowUpper();
420     double * fullSolution = model.primalColumnSolution();
421     const double * fullLower = model.columnLower();
422     const double * fullUpper = model.columnUpper();
423     const double * rowSolution = master.primalRowSolution();
424     double * fullRowSolution = model.primalRowSolution();
425     int kRow = 0;
426     for (iRow = 0; iRow < numberRows; iRow++) {
427          if (rowBlock[iRow] == -1) {
428               model.setRowStatus(iRow, master.getRowStatus(kRow));
429               fullRowSolution[iRow] = rowSolution[kRow++];
430          }
431     }
432     int kColumn = 0;
433     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
434          if (columnBlock[iColumn] == -1) {
435               model.setStatus(iColumn, master.getStatus(kColumn));
436               fullSolution[iColumn] = solution[kColumn++];
437          }
438     }
439     for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
440          // convert top bit to by rows
441          top[iBlock].reverseOrdering();
442          // zero solution
443          memset(sol, 0, numberColumnsGenerated * sizeof(double));
444          int i;
445          for (i = numberMasterColumns; i < numberColumnsGenerated; i++)
446               if (whichBlock[i-numberMasterColumns] == iBlock)
447                    sol[i] = solution[i];
448          memset(lower, 0, (numberMasterRows + numberBlocks) *sizeof(double));
449          master.times(1.0, sol, lower);
450          for (iRow = 0; iRow < numberMasterRows; iRow++) {
451               double value = lower[iRow];
452               if (masterUpper[iRow] < 1.0e20)
453                    upper[iRow] = value;
454               else
455                    upper[iRow] = COIN_DBL_MAX;
456               if (masterLower[iRow] > -1.0e20)
457                    lower[iRow] = value;
458               else
459                    lower[iRow] = -COIN_DBL_MAX;
460          }
461          sub[iBlock].addRows(numberMasterRows, lower, upper,
462                              top[iBlock].getVectorStarts(),
463                              top[iBlock].getVectorLengths(),
464                              top[iBlock].getIndices(),
465                              top[iBlock].getElements());
466          sub[iBlock].primal();
467          const double * subSolution = sub[iBlock].primalColumnSolution();
468          const double * subRowSolution = sub[iBlock].primalRowSolution();
469          // move solution
470          kColumn = 0;
471          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
472               if (columnBlock[iColumn] == iBlock) {
473                    model.setStatus(iColumn, sub[iBlock].getStatus(kColumn));
474                    fullSolution[iColumn] = subSolution[kColumn++];
475               }
476          }
477          assert(kColumn == sub[iBlock].numberColumns());
478          kRow = 0;
479          for (iRow = 0; iRow < numberRows; iRow++) {
480               if (rowBlock[iRow] == iBlock) {
481                    model.setRowStatus(iRow, sub[iBlock].getRowStatus(kRow));
482                    fullRowSolution[iRow] = subRowSolution[kRow++];
483               }
484          }
485          assert(kRow == sub[iBlock].numberRows() - numberMasterRows);
486     }
487     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
488          if (fullSolution[iColumn] < fullUpper[iColumn] - 1.0e-8 &&
489                    fullSolution[iColumn] > fullLower[iColumn] + 1.0e-8) {
490               assert(model.getStatus(iColumn) == ClpSimplex::basic);
491          } else if (fullSolution[iColumn] >= fullUpper[iColumn] - 1.0e-8) {
492               // may help to make rest non basic
493               model.setStatus(iColumn, ClpSimplex::atUpperBound);
494          } else if (fullSolution[iColumn] <= fullLower[iColumn] + 1.0e-8) {
495               // may help to make rest non basic
496               model.setStatus(iColumn, ClpSimplex::atLowerBound);
497          }
498     }
499     for (iRow = 0; iRow < numberRows; iRow++)
500          model.setRowStatus(iRow, ClpSimplex::superBasic);
501     model.primal(1);
502     delete [] sol;
503     delete [] lower;
504     delete [] upper;
505     delete [] whichBlock;
506     delete [] when;
507     delete [] columnAdd;
508     delete [] rowAdd;
509     delete [] elementAdd;
510     delete [] objective;
511     delete [] top;
512     delete [] sub;
513     delete [] rowBlock;
514     delete [] columnBlock;
515     return 0;
516}
Note: See TracBrowser for help on using the repository browser.