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

Last change on this file since 1552 was 1552, checked in by mjs, 11 years ago

Format examples with 'astyle -A4 -p'.

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