source: trunk/Samples/decompose.cpp @ 423

Last change on this file since 423 was 423, checked in by forrest, 15 years ago

various Cholesky

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