Ignore:
Timestamp:
May 24, 2010 9:03:59 PM (10 years ago)
Author:
mjs
Message:

Format examples with 'astyle -A4 -p'.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/examples/decompose.cpp

    r1370 r1552  
    99int main (int argc, const char *argv[])
    1010{
    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);
     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);
    171171#if 0
    172     // temp
    173     double * upper = sub[iBlock].columnUpper();
    174     for (iColumn=0;iColumn<numberColumn2;iColumn++)
    175       upper[iColumn]=100.0;
     172          // temp
     173          double * upper = sub[iBlock].columnUpper();
     174          for (iColumn = 0; iColumn < numberColumn2; iColumn++)
     175               upper[iColumn] = 100.0;
    176176#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;
     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;
    201201
    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);
     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);
    242242
    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     }
     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          }
    250250
    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 }   
     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 TracChangeset for help on using the changeset viewer.