source: trunk/Clp/examples/decompose.cpp @ 2278

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

COIN_BIG_INDEX 2 changes

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