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

Last change on this file since 1370 was 1370, checked in by forrest, 11 years ago

add ids

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 16.3 KB
Line 
1/* $Id: decompose.cpp 1370 2009-06-04 09:37:13Z forrest $ */
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.