Changeset 1552 for trunk


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

Format examples with 'astyle -A4 -p'.

Location:
trunk/Clp/examples
Files:
31 edited

Legend:

Unmodified
Added
Removed

  • trunk/Clp/examples/addBits.cpp

    r1370 r1552  
    3434int main (int argc, const char *argv[])
    3535{
    36   // Empty model
    37   ClpSimplex  model;
    38   std::string mpsFileName = "../../Data/Netlib/25fv47.mps";
    39   if (argc>=2) mpsFileName = argv[1];
    40   int status=model.readMps(mpsFileName.c_str(),true);
     36     // Empty model
     37     ClpSimplex  model;
     38     std::string mpsFileName = "../../Data/Netlib/25fv47.mps";
     39     if (argc >= 2) mpsFileName = argv[1];
     40     int status = model.readMps(mpsFileName.c_str(), true);
    4141
    42   if (status) {
    43     fprintf(stderr,"Bad readMps %s\n",mpsFileName.c_str());
    44     fprintf(stdout,"Bad readMps %s\n",mpsFileName.c_str());
    45     exit(1);
    46   }
    47   // Point to data
    48   int numberRows = model.numberRows();
    49   const double * rowLower = model.rowLower();
    50   const double * rowUpper = model.rowUpper();
    51   int numberColumns = model.numberColumns();
    52   const double * columnLower = model.columnLower();
    53   const double * columnUpper = model.columnUpper();
    54   const double * columnObjective = model.objective();
    55   CoinPackedMatrix * matrix = model.matrix();
    56   // get row copy
    57   CoinPackedMatrix rowCopy = *matrix;
    58   rowCopy.reverseOrdering();
    59   const int * column = rowCopy.getIndices();
    60   const int * rowLength = rowCopy.getVectorLengths();
    61   const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
    62   const double * element = rowCopy.getElements();
    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  
    68   // solve
    69   model.dual();
    70   // Now build new model
    71   CoinModel build;
    72   double time1 = CoinCpuTime();
    73   // Row bounds
    74   int iRow;
    75   for (iRow=0;iRow<numberRows;iRow++) {
    76     build.setRowBounds(iRow,rowLower[iRow],rowUpper[iRow]);
    77     // optional name
    78     build.setRowName(iRow,model.rowName(iRow).c_str());
    79   }
    80   // Column bounds and objective
    81   int iColumn;
    82   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    83     build.setColumnLower(iColumn,columnLower[iColumn]);
    84     build.setColumnUpper(iColumn,columnUpper[iColumn]);
    85     build.setObjective(iColumn,columnObjective[iColumn]);
    86     // optional name
    87     build.setColumnName(iColumn,model.columnName(iColumn).c_str());
    88   }
    89   // Adds elements one by one by row (backwards by row)
    90   for (iRow=numberRows-1;iRow>=0;iRow--) {
    91     int start = rowStart[iRow];
    92     for (int j=start;j<start+rowLength[iRow];j++)
    93       build(iRow,column[j],element[j]);
    94   }
    95   double time2 = CoinCpuTime();
    96   // Now create clpsimplex
    97   ClpSimplex model2;
    98   model2.loadProblem(build);
    99   double time3 = CoinCpuTime();
    100   printf("Time for build using CoinModel is %g (%g for loadproblem)\n",time3-time1,
    101          time3-time2);
    102   model2.dual();
    103   // Now do with strings attached
    104   // Save build to show how to go over rows
    105   CoinModel saveBuild = build;
    106   build = CoinModel();
    107   time1 = CoinCpuTime();
    108   // Column bounds
    109   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    110     build.setColumnLower(iColumn,columnLower[iColumn]);
    111     build.setColumnUpper(iColumn,columnUpper[iColumn]);
    112   }
    113   // Objective - half the columns as is and half with multiplier of "1.0+multiplier"
    114   // Pick up from saveBuild (for no reason at all)
    115   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    116     double value = saveBuild.objective(iColumn);
    117     if (iColumn*2<numberColumns) {
    118       build.setObjective(iColumn,columnObjective[iColumn]);
    119     } else {
    120       // create as string
    121       char temp[100];
    122       sprintf(temp,"%g + abs(%g*multiplier)",value,value);
    123       build.setObjective(iColumn,temp);
    124     }
    125   }
    126   // It then adds rows one by one but for half the rows sets their values
    127   //      with multiplier of "1.0+1.5*multiplier"
    128   for (iRow=0;iRow<numberRows;iRow++) {
    129     if (iRow*2<numberRows) {
    130       // add row in simple way
    131       int start = rowStart[iRow];
    132       build.addRow(rowLength[iRow],column+start,element+start,
    133                    rowLower[iRow],rowUpper[iRow]);
    134     } else {
    135       // As we have to add one by one let's get from saveBuild
    136       CoinModelLink triple=saveBuild.firstInRow(iRow);
    137       while (triple.column()>=0) {
    138         int iColumn = triple.column();
    139         if (iColumn*2<numberColumns) {
    140           // just value as normal
    141           build(iRow,triple.column(),triple.value());
    142         } else {
    143           // create as string
    144           char temp[100];
    145           sprintf(temp,"%g + (1.5*%g*multiplier)",triple.value(), triple.value());
    146           build(iRow,iColumn,temp);
    147         }
    148         triple=saveBuild.next(triple);
    149       }
    150       // but remember to do rhs
    151       build.setRowLower(iRow,rowLower[iRow]);
    152       build.setRowUpper(iRow,rowUpper[iRow]);
    153     }
    154   }
    155   time2 = CoinCpuTime();
    156   // Now create ClpSimplex
    157   // If small switch on error printing
    158   if (numberColumns<50)
    159     build.setLogLevel(1);
    160   int numberErrors=model2.loadProblem(build);
    161   // should fail as we never set multiplier
    162   assert (numberErrors);
    163   time3 = CoinCpuTime()-time2;
    164   // subtract out unsuccessful times
    165   time1 += time3;
    166   time2 += time3;
    167   build.associateElement("multiplier",0.0);
    168   numberErrors=model2.loadProblem(build);
    169   assert (!numberErrors);
    170   time3 = CoinCpuTime();
    171   printf("Time for build using CoinModel is %g (%g for successful loadproblem)\n",time3-time1,
    172          time3-time2);
    173   build.writeMps("zero.mps");
    174   // It then loops with multiplier going from 0.0 to 2.0 in increments of 0.1
    175   for (double multiplier=0.0;multiplier<2.0;multiplier+= 0.1) {
    176     build.associateElement("multiplier",multiplier);
    177     numberErrors=model2.loadProblem(build,true);
    178     assert (!numberErrors);
    179     model2.dual();
    180   }
     42     if (status) {
     43          fprintf(stderr, "Bad readMps %s\n", mpsFileName.c_str());
     44          fprintf(stdout, "Bad readMps %s\n", mpsFileName.c_str());
     45          exit(1);
     46     }
     47     // Point to data
     48     int numberRows = model.numberRows();
     49     const double * rowLower = model.rowLower();
     50     const double * rowUpper = model.rowUpper();
     51     int numberColumns = model.numberColumns();
     52     const double * columnLower = model.columnLower();
     53     const double * columnUpper = model.columnUpper();
     54     const double * columnObjective = model.objective();
     55     CoinPackedMatrix * matrix = model.matrix();
     56     // get row copy
     57     CoinPackedMatrix rowCopy = *matrix;
     58     rowCopy.reverseOrdering();
     59     const int * column = rowCopy.getIndices();
     60     const int * rowLength = rowCopy.getVectorLengths();
     61     const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
     62     const double * element = rowCopy.getElements();
     63     //const int * row = matrix->getIndices();
     64     //const int * columnLength = matrix->getVectorLengths();
     65     //const CoinBigIndex * columnStart = matrix->getVectorStarts();
     66     //const double * elementByColumn = matrix->getElements();
    18167
    182   return 0;
    183 }   
     68     // solve
     69     model.dual();
     70     // Now build new model
     71     CoinModel build;
     72     double time1 = CoinCpuTime();
     73     // Row bounds
     74     int iRow;
     75     for (iRow = 0; iRow < numberRows; iRow++) {
     76          build.setRowBounds(iRow, rowLower[iRow], rowUpper[iRow]);
     77          // optional name
     78          build.setRowName(iRow, model.rowName(iRow).c_str());
     79     }
     80     // Column bounds and objective
     81     int iColumn;
     82     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     83          build.setColumnLower(iColumn, columnLower[iColumn]);
     84          build.setColumnUpper(iColumn, columnUpper[iColumn]);
     85          build.setObjective(iColumn, columnObjective[iColumn]);
     86          // optional name
     87          build.setColumnName(iColumn, model.columnName(iColumn).c_str());
     88     }
     89     // Adds elements one by one by row (backwards by row)
     90     for (iRow = numberRows - 1; iRow >= 0; iRow--) {
     91          int start = rowStart[iRow];
     92          for (int j = start; j < start + rowLength[iRow]; j++)
     93               build(iRow, column[j], element[j]);
     94     }
     95     double time2 = CoinCpuTime();
     96     // Now create clpsimplex
     97     ClpSimplex model2;
     98     model2.loadProblem(build);
     99     double time3 = CoinCpuTime();
     100     printf("Time for build using CoinModel is %g (%g for loadproblem)\n", time3 - time1,
     101            time3 - time2);
     102     model2.dual();
     103     // Now do with strings attached
     104     // Save build to show how to go over rows
     105     CoinModel saveBuild = build;
     106     build = CoinModel();
     107     time1 = CoinCpuTime();
     108     // Column bounds
     109     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     110          build.setColumnLower(iColumn, columnLower[iColumn]);
     111          build.setColumnUpper(iColumn, columnUpper[iColumn]);
     112     }
     113     // Objective - half the columns as is and half with multiplier of "1.0+multiplier"
     114     // Pick up from saveBuild (for no reason at all)
     115     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     116          double value = saveBuild.objective(iColumn);
     117          if (iColumn * 2 < numberColumns) {
     118               build.setObjective(iColumn, columnObjective[iColumn]);
     119          } else {
     120               // create as string
     121               char temp[100];
     122               sprintf(temp, "%g + abs(%g*multiplier)", value, value);
     123               build.setObjective(iColumn, temp);
     124          }
     125     }
     126     // It then adds rows one by one but for half the rows sets their values
     127     //      with multiplier of "1.0+1.5*multiplier"
     128     for (iRow = 0; iRow < numberRows; iRow++) {
     129          if (iRow * 2 < numberRows) {
     130               // add row in simple way
     131               int start = rowStart[iRow];
     132               build.addRow(rowLength[iRow], column + start, element + start,
     133                            rowLower[iRow], rowUpper[iRow]);
     134          } else {
     135               // As we have to add one by one let's get from saveBuild
     136               CoinModelLink triple = saveBuild.firstInRow(iRow);
     137               while (triple.column() >= 0) {
     138                    int iColumn = triple.column();
     139                    if (iColumn * 2 < numberColumns) {
     140                         // just value as normal
     141                         build(iRow, triple.column(), triple.value());
     142                    } else {
     143                         // create as string
     144                         char temp[100];
     145                         sprintf(temp, "%g + (1.5*%g*multiplier)", triple.value(), triple.value());
     146                         build(iRow, iColumn, temp);
     147                    }
     148                    triple = saveBuild.next(triple);
     149               }
     150               // but remember to do rhs
     151               build.setRowLower(iRow, rowLower[iRow]);
     152               build.setRowUpper(iRow, rowUpper[iRow]);
     153          }
     154     }
     155     time2 = CoinCpuTime();
     156     // Now create ClpSimplex
     157     // If small switch on error printing
     158     if (numberColumns < 50)
     159          build.setLogLevel(1);
     160     int numberErrors = model2.loadProblem(build);
     161     // should fail as we never set multiplier
     162     assert (numberErrors);
     163     time3 = CoinCpuTime() - time2;
     164     // subtract out unsuccessful times
     165     time1 += time3;
     166     time2 += time3;
     167     build.associateElement("multiplier", 0.0);
     168     numberErrors = model2.loadProblem(build);
     169     assert (!numberErrors);
     170     time3 = CoinCpuTime();
     171     printf("Time for build using CoinModel is %g (%g for successful loadproblem)\n", time3 - time1,
     172            time3 - time2);
     173     build.writeMps("zero.mps");
     174     // It then loops with multiplier going from 0.0 to 2.0 in increments of 0.1
     175     for (double multiplier = 0.0; multiplier < 2.0; multiplier += 0.1) {
     176          build.associateElement("multiplier", multiplier);
     177          numberErrors = model2.loadProblem(build, true);
     178          assert (!numberErrors);
     179          model2.dual();
     180     }
     181
     182     return 0;
     183}

  • trunk/Clp/examples/addColumns.cpp

    r1370 r1552  
    1515int main (int argc, const char *argv[])
    1616{
    17   {
    18     // Empty model
    19     ClpSimplex  model;
    20    
    21     // Bounds on rows - as dense vector
    22     double lower[]={2.0,1.0};
    23     double upper[]={COIN_DBL_MAX,1.0};
    24    
    25     // Create space for 2 rows
    26     model.resize(2,0);
    27     // Fill in
    28     int i;
    29     // Now row bounds
    30     for (i=0;i<2;i++) {
    31       model.setRowLower(i,lower[i]);
    32       model.setRowUpper(i,upper[i]);
    33     }
    34     // Now add column 1
    35     int column1Index[] = {0,1};
    36     double column1Value[]={1.0,1.0};
    37     model.addColumn(2,column1Index,column1Value,
    38                     0.0,2,1.0);
    39     // Now add column 2
    40     int column2Index[] = {1};
    41     double column2Value[]={-5.0};
    42     model.addColumn(1,column2Index,column2Value,
    43                     0.0,COIN_DBL_MAX,0.0);
    44     // Now add column 3
    45     int column3Index[] = {0,1};
    46     double column3Value[]={1.0,1.0};
    47     model.addColumn(2,column3Index,column3Value,
    48                     0.0,4.0,4.0);
    49     // solve
    50     model.dual();
    51    
    52     /*
    53       Adding one column at a time has a significant overhead so let's
    54       try a more complicated but faster way
    55      
    56       First time adding in 10000 columns one by one
    57      
    58     */
    59     model.allSlackBasis();
    60     ClpSimplex modelSave=model;
    61     double time1 = CoinCpuTime();
    62     int k;
    63     for ( k=0;k<10000;k++) {
    64       int column2Index[] = {0,1};
    65       double column2Value[]={1.0,-5.0};
    66       model.addColumn(2,column2Index,column2Value,
    67                       0.0,1.0,10000.0);
    68     }
    69     printf("Time for 10000 addColumn is %g\n",CoinCpuTime()-time1);
    70     model.dual();
    71     model=modelSave;
    72     // Now use build
    73     CoinBuild buildObject;
    74     time1 = CoinCpuTime();
    75     for ( k=0;k<100000;k++) {
    76       int column2Index[] = {0,1};
    77       double column2Value[]={1.0,-5.0};
    78       buildObject.addColumn(2,column2Index,column2Value,
    79                             0.0,1.0,10000.0);
    80     }
    81     model.addColumns(buildObject);
    82     printf("Time for 100000 addColumn using CoinBuild is %g\n",CoinCpuTime()-time1);
    83     model.dual();
    84     model=modelSave;
    85     // Now use build +-1
    86     int del[]={0,1,2};
    87     model.deleteColumns(3,del);
    88     CoinBuild buildObject2;
    89     time1 = CoinCpuTime();
    90     for ( k=0;k<10000;k++) {
    91       int column2Index[] = {0,1};
    92       double column2Value[]={1.0,1.0,-1.0};
    93       int bias = k&1;
    94       buildObject2.addColumn(2,column2Index,column2Value+bias,
    95                              0.0,1.0,10000.0);
    96     }
    97     model.addColumns(buildObject2,true);
    98     printf("Time for 10000 addColumn using CoinBuild+-1 is %g\n",CoinCpuTime()-time1);
    99     model.dual();
    100     model=modelSave;
    101     // Now use build +-1
    102     model.deleteColumns(3,del);
    103     CoinModel modelObject2;
    104     time1 = CoinCpuTime();
    105     for ( k=0;k<10000;k++) {
    106       int column2Index[] = {0,1};
    107       double column2Value[]={1.0,1.0,-1.0};
    108       int bias = k&1;
    109       modelObject2.addColumn(2,column2Index,column2Value+bias,
    110                              0.0,1.0,10000.0);
    111     }
    112     model.addColumns(modelObject2,true);
    113     printf("Time for 10000 addColumn using CoinModel+-1 is %g\n",CoinCpuTime()-time1);
    114     //model.writeMps("xx.mps");
    115     model.dual();
    116     model=modelSave;
    117     // Now use model
    118     CoinModel modelObject;
    119     time1 = CoinCpuTime();
    120     for ( k=0;k<100000;k++) {
    121       int column2Index[] = {0,1};
    122       double column2Value[]={1.0,-5.0};
    123       modelObject.addColumn(2,column2Index,column2Value,
    124                             0.0,1.0,10000.0);
    125     }
    126     model.addColumns(modelObject);
    127     printf("Time for 100000 addColumn using CoinModel is %g\n",CoinCpuTime()-time1);
    128     model.dual();
    129     // Print column solution Just first 3 columns
    130     int numberColumns = model.numberColumns();
    131     numberColumns = CoinMin(3,numberColumns);
    132    
    133     // Alternatively getColSolution()
    134     double * columnPrimal = model.primalColumnSolution();
    135     // Alternatively getReducedCost()
    136     double * columnDual = model.dualColumnSolution();
    137     // Alternatively getColLower()
    138     double * columnLower = model.columnLower();  // Alternatively getColUpper()
    139     double * columnUpper = model.columnUpper();
    140     // Alternatively getObjCoefficients()
    141     double * columnObjective = model.objective();
    142    
    143     int iColumn;
    144    
    145     std::cout<<"               Primal          Dual         Lower         Upper          Cost"
    146              <<std::endl;
    147    
    148     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    149       double value;
    150       std::cout<<std::setw(6)<<iColumn<<" ";
    151       value = columnPrimal[iColumn];
    152       if (fabs(value)<1.0e5)
    153         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    154       else
    155         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    156       value = columnDual[iColumn];
    157       if (fabs(value)<1.0e5)
    158         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    159       else
    160         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    161       value = columnLower[iColumn];
    162       if (fabs(value)<1.0e5)
    163         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    164       else
    165         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    166       value = columnUpper[iColumn];
    167       if (fabs(value)<1.0e5)
    168         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    169       else
    170         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    171       value = columnObjective[iColumn];
    172       if (fabs(value)<1.0e5)
    173         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    174       else
    175         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    176      
    177       std::cout<<std::endl;
    178     }
    179     std::cout<<"--------------------------------------"<<std::endl;
    180   }
    181   {
    182     // Now copy a model
    183     ClpSimplex  model;
    184     int status;
    185     if (argc<2)
    186       status=model.readMps("../../Data/Sample/p0033.mps");
    187     else
    188       status=model.readMps(argv[1]);
    189     if (status) {
    190       printf("errors on input\n");
    191       exit(77);
    192     }
    193     model.initialSolve();
    194     int numberRows = model.numberRows();
    195     int numberColumns = model.numberColumns();
    196     const double * rowLower = model.rowLower();
    197     const double * rowUpper = model.rowUpper();
    198 
    199     // Start off model2
    200     ClpSimplex model2;
    201     model2.addRows(numberRows,rowLower,rowUpper,NULL);
    202 
    203     // Build object
    204     CoinBuild buildObject;
    205     // Add columns
    206     const double * columnLower = model.columnLower();
    207     const double * columnUpper = model.columnUpper();
    208     const double * objective = model.objective();
    209     CoinPackedMatrix * matrix = model.matrix();
    210     const int * row = matrix->getIndices();
    211     const int * columnLength = matrix->getVectorLengths();
    212     const CoinBigIndex * columnStart = matrix->getVectorStarts();
    213     const double * elementByColumn = matrix->getElements();
    214     for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    215       CoinBigIndex start = columnStart[iColumn];
    216       buildObject.addColumn(columnLength[iColumn],row+start,elementByColumn+start,
    217                             columnLower[iColumn],columnUpper[iColumn],
    218                             objective[iColumn]);
    219     }
    220 
    221     // add in
    222     model2.addColumns(buildObject);
    223     model2.initialSolve();
    224   }
    225   {
    226     // and again
    227     ClpSimplex  model;
    228     int status;
    229     if (argc<2)
    230       status=model.readMps("../../Data/Sample/p0033.mps");
    231     else
    232       status=model.readMps(argv[1]);
    233     if (status) {
    234       printf("errors on input\n");
    235       exit(77);
    236     }
    237     model.initialSolve();
    238     int numberRows = model.numberRows();
    239     int numberColumns = model.numberColumns();
    240     const double * rowLower = model.rowLower();
    241     const double * rowUpper = model.rowUpper();
    242 
    243     // Build object
    244     CoinModel buildObject;
    245     for (int iRow=0;iRow<numberRows;iRow++)
    246       buildObject.setRowBounds(iRow,rowLower[iRow],rowUpper[iRow]);
    247 
    248     // Add columns
    249     const double * columnLower = model.columnLower();
    250     const double * columnUpper = model.columnUpper();
    251     const double * objective = model.objective();
    252     CoinPackedMatrix * matrix = model.matrix();
    253     const int * row = matrix->getIndices();
    254     const int * columnLength = matrix->getVectorLengths();
    255     const CoinBigIndex * columnStart = matrix->getVectorStarts();
    256     const double * elementByColumn = matrix->getElements();
    257     for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    258       CoinBigIndex start = columnStart[iColumn];
    259       buildObject.addColumn(columnLength[iColumn],row+start,elementByColumn+start,
    260                             columnLower[iColumn],columnUpper[iColumn],
    261                             objective[iColumn]);
    262     }
    263 
    264     // add in
    265     ClpSimplex model2;
    266     model2.loadProblem(buildObject);
    267     model2.initialSolve();
    268   }
    269   return 0;
    270 }   
     17     {
     18          // Empty model
     19          ClpSimplex  model;
     20
     21          // Bounds on rows - as dense vector
     22          double lower[] = {2.0, 1.0};
     23          double upper[] = {COIN_DBL_MAX, 1.0};
     24
     25          // Create space for 2 rows
     26          model.resize(2, 0);
     27          // Fill in
     28          int i;
     29          // Now row bounds
     30          for (i = 0; i < 2; i++) {
     31               model.setRowLower(i, lower[i]);
     32               model.setRowUpper(i, upper[i]);
     33          }
     34          // Now add column 1
     35          int column1Index[] = {0, 1};
     36          double column1Value[] = {1.0, 1.0};
     37          model.addColumn(2, column1Index, column1Value,
     38                          0.0, 2, 1.0);
     39          // Now add column 2
     40          int column2Index[] = {1};
     41          double column2Value[] = { -5.0};
     42          model.addColumn(1, column2Index, column2Value,
     43                          0.0, COIN_DBL_MAX, 0.0);
     44          // Now add column 3
     45          int column3Index[] = {0, 1};
     46          double column3Value[] = {1.0, 1.0};
     47          model.addColumn(2, column3Index, column3Value,
     48                          0.0, 4.0, 4.0);
     49          // solve
     50          model.dual();
     51
     52          /*
     53            Adding one column at a time has a significant overhead so let's
     54            try a more complicated but faster way
     55
     56            First time adding in 10000 columns one by one
     57
     58          */
     59          model.allSlackBasis();
     60          ClpSimplex modelSave = model;
     61          double time1 = CoinCpuTime();
     62          int k;
     63          for ( k = 0; k < 10000; k++) {
     64               int column2Index[] = {0, 1};
     65               double column2Value[] = {1.0, -5.0};
     66               model.addColumn(2, column2Index, column2Value,
     67                               0.0, 1.0, 10000.0);
     68          }
     69          printf("Time for 10000 addColumn is %g\n", CoinCpuTime() - time1);
     70          model.dual();
     71          model = modelSave;
     72          // Now use build
     73          CoinBuild buildObject;
     74          time1 = CoinCpuTime();
     75          for ( k = 0; k < 100000; k++) {
     76               int column2Index[] = {0, 1};
     77               double column2Value[] = {1.0, -5.0};
     78               buildObject.addColumn(2, column2Index, column2Value,
     79                                     0.0, 1.0, 10000.0);
     80          }
     81          model.addColumns(buildObject);
     82          printf("Time for 100000 addColumn using CoinBuild is %g\n", CoinCpuTime() - time1);
     83          model.dual();
     84          model = modelSave;
     85          // Now use build +-1
     86          int del[] = {0, 1, 2};
     87          model.deleteColumns(3, del);
     88          CoinBuild buildObject2;
     89          time1 = CoinCpuTime();
     90          for ( k = 0; k < 10000; k++) {
     91               int column2Index[] = {0, 1};
     92               double column2Value[] = {1.0, 1.0, -1.0};
     93               int bias = k & 1;
     94               buildObject2.addColumn(2, column2Index, column2Value + bias,
     95                                      0.0, 1.0, 10000.0);
     96          }
     97          model.addColumns(buildObject2, true);
     98          printf("Time for 10000 addColumn using CoinBuild+-1 is %g\n", CoinCpuTime() - time1);
     99          model.dual();
     100          model = modelSave;
     101          // Now use build +-1
     102          model.deleteColumns(3, del);
     103          CoinModel modelObject2;
     104          time1 = CoinCpuTime();
     105          for ( k = 0; k < 10000; k++) {
     106               int column2Index[] = {0, 1};
     107               double column2Value[] = {1.0, 1.0, -1.0};
     108               int bias = k & 1;
     109               modelObject2.addColumn(2, column2Index, column2Value + bias,
     110                                      0.0, 1.0, 10000.0);
     111          }
     112          model.addColumns(modelObject2, true);
     113          printf("Time for 10000 addColumn using CoinModel+-1 is %g\n", CoinCpuTime() - time1);
     114          //model.writeMps("xx.mps");
     115          model.dual();
     116          model = modelSave;
     117          // Now use model
     118          CoinModel modelObject;
     119          time1 = CoinCpuTime();
     120          for ( k = 0; k < 100000; k++) {
     121               int column2Index[] = {0, 1};
     122               double column2Value[] = {1.0, -5.0};
     123               modelObject.addColumn(2, column2Index, column2Value,
     124                                     0.0, 1.0, 10000.0);
     125          }
     126          model.addColumns(modelObject);
     127          printf("Time for 100000 addColumn using CoinModel is %g\n", CoinCpuTime() - time1);
     128          model.dual();
     129          // Print column solution Just first 3 columns
     130          int numberColumns = model.numberColumns();
     131          numberColumns = CoinMin(3, numberColumns);
     132
     133          // Alternatively getColSolution()
     134          double * columnPrimal = model.primalColumnSolution();
     135          // Alternatively getReducedCost()
     136          double * columnDual = model.dualColumnSolution();
     137          // Alternatively getColLower()
     138          double * columnLower = model.columnLower();  // Alternatively getColUpper()
     139          double * columnUpper = model.columnUpper();
     140          // Alternatively getObjCoefficients()
     141          double * columnObjective = model.objective();
     142
     143          int iColumn;
     144
     145          std::cout << "               Primal          Dual         Lower         Upper          Cost"
     146                    << std::endl;
     147
     148          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     149               double value;
     150               std::cout << std::setw(6) << iColumn << " ";
     151               value = columnPrimal[iColumn];
     152               if (fabs(value) < 1.0e5)
     153                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     154               else
     155                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     156               value = columnDual[iColumn];
     157               if (fabs(value) < 1.0e5)
     158                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     159               else
     160                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     161               value = columnLower[iColumn];
     162               if (fabs(value) < 1.0e5)
     163                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     164               else
     165                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     166               value = columnUpper[iColumn];
     167               if (fabs(value) < 1.0e5)
     168                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     169               else
     170                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     171               value = columnObjective[iColumn];
     172               if (fabs(value) < 1.0e5)
     173                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     174               else
     175                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     176
     177               std::cout << std::endl;
     178          }
     179          std::cout << "--------------------------------------" << std::endl;
     180     }
     181     {
     182          // Now copy a model
     183          ClpSimplex  model;
     184          int status;
     185          if (argc < 2)
     186               status = model.readMps("../../Data/Sample/p0033.mps");
     187          else
     188               status = model.readMps(argv[1]);
     189          if (status) {
     190               printf("errors on input\n");
     191               exit(77);
     192          }
     193          model.initialSolve();
     194          int numberRows = model.numberRows();
     195          int numberColumns = model.numberColumns();
     196          const double * rowLower = model.rowLower();
     197          const double * rowUpper = model.rowUpper();
     198
     199          // Start off model2
     200          ClpSimplex model2;
     201          model2.addRows(numberRows, rowLower, rowUpper, NULL);
     202
     203          // Build object
     204          CoinBuild buildObject;
     205          // Add columns
     206          const double * columnLower = model.columnLower();
     207          const double * columnUpper = model.columnUpper();
     208          const double * objective = model.objective();
     209          CoinPackedMatrix * matrix = model.matrix();
     210          const int * row = matrix->getIndices();
     211          const int * columnLength = matrix->getVectorLengths();
     212          const CoinBigIndex * columnStart = matrix->getVectorStarts();
     213          const double * elementByColumn = matrix->getElements();
     214          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     215               CoinBigIndex start = columnStart[iColumn];
     216               buildObject.addColumn(columnLength[iColumn], row + start, elementByColumn + start,
     217                                     columnLower[iColumn], columnUpper[iColumn],
     218                                     objective[iColumn]);
     219          }
     220
     221          // add in
     222          model2.addColumns(buildObject);
     223          model2.initialSolve();
     224     }
     225     {
     226          // and again
     227          ClpSimplex  model;
     228          int status;
     229          if (argc < 2)
     230               status = model.readMps("../../Data/Sample/p0033.mps");
     231          else
     232               status = model.readMps(argv[1]);
     233          if (status) {
     234               printf("errors on input\n");
     235               exit(77);
     236          }
     237          model.initialSolve();
     238          int numberRows = model.numberRows();
     239          int numberColumns = model.numberColumns();
     240          const double * rowLower = model.rowLower();
     241          const double * rowUpper = model.rowUpper();
     242
     243          // Build object
     244          CoinModel buildObject;
     245          for (int iRow = 0; iRow < numberRows; iRow++)
     246               buildObject.setRowBounds(iRow, rowLower[iRow], rowUpper[iRow]);
     247
     248          // Add columns
     249          const double * columnLower = model.columnLower();
     250          const double * columnUpper = model.columnUpper();
     251          const double * objective = model.objective();
     252          CoinPackedMatrix * matrix = model.matrix();
     253          const int * row = matrix->getIndices();
     254          const int * columnLength = matrix->getVectorLengths();
     255          const CoinBigIndex * columnStart = matrix->getVectorStarts();
     256          const double * elementByColumn = matrix->getElements();
     257          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     258               CoinBigIndex start = columnStart[iColumn];
     259               buildObject.addColumn(columnLength[iColumn], row + start, elementByColumn + start,
     260                                     columnLower[iColumn], columnUpper[iColumn],
     261                                     objective[iColumn]);
     262          }
     263
     264          // add in
     265          ClpSimplex model2;
     266          model2.loadProblem(buildObject);
     267          model2.initialSolve();
     268     }
     269     return 0;
     270}

  • trunk/Clp/examples/addRows.cpp

    r1370 r1552  
    1515int main (int argc, const char *argv[])
    1616{
    17   try {
    18     // Empty model
    19     ClpSimplex  model;
    20    
    21     // Objective - just nonzeros
    22     int objIndex[]={0,2};
    23     double objValue[]={1.0,4.0};
    24     // Upper bounds - as dense vector
    25     double upper[]={2.0,COIN_DBL_MAX,4.0};
    26    
    27     // Create space for 3 columns
    28     model.resize(0,3);
    29     // Fill in
    30     int i;
    31     // Virtuous way
    32     // First objective
    33     for (i=0;i<2;i++)
    34       model.setObjectiveCoefficient(objIndex[i],objValue[i]);
    35     // Now bounds (lower will be zero by default but do again)
    36     for (i=0;i<3;i++) {
    37       model.setColumnLower(i,0.0);
    38       model.setColumnUpper(i,upper[i]);
    39     }
    40     /*
    41       We could also have done in non-virtuous way e.g.
    42       double * objective = model.objective();
    43       and then set directly
    44     */
    45     // Faster to add rows all at once - but this is easier to show
    46     // Now add row 1 as >= 2.0
    47     int row1Index[] = {0,2};
    48     double row1Value[]={1.0,1.0};
    49     model.addRow(2,row1Index,row1Value,
    50                  2.0,COIN_DBL_MAX);
    51     // Now add row 2 as == 1.0
    52     int row2Index[] = {0,1,2};
    53     double row2Value[]={1.0,-5.0,1.0};
    54     model.addRow(3,row2Index,row2Value,
    55                  1.0,1.0);
    56     // solve
    57     model.dual();
    58    
    59     /*
    60       Adding one row at a time has a significant overhead so let's
    61       try a more complicated but faster way
    62      
    63       First time adding in 10000 rows one by one
    64     */
    65     model.allSlackBasis();
    66     ClpSimplex modelSave=model;
    67     double time1 = CoinCpuTime();
    68     int k;
    69     for ( k=0;k<10000;k++) {
    70       int row2Index[] = {0,1,2};
    71       double row2Value[]={1.0,-5.0,1.0};
    72       model.addRow(3,row2Index,row2Value,
    73                    1.0,1.0);
    74     }
    75     printf("Time for 10000 addRow is %g\n",CoinCpuTime()-time1);
    76     model.dual();
    77     model=modelSave;
    78     // Now use build
    79     CoinBuild buildObject;
    80     time1 = CoinCpuTime();
    81     for ( k=0;k<10000;k++) {
    82       int row2Index[] = {0,1,2};
    83       double row2Value[]={1.0,-5.0,1.0};
    84       buildObject.addRow(3,row2Index,row2Value,
    85                          1.0,1.0);
    86     }
    87     model.addRows(buildObject);
    88     printf("Time for 10000 addRow using CoinBuild is %g\n",CoinCpuTime()-time1);
    89     model.dual();
    90     model=modelSave;
    91     int del[]={0,1,2};
    92     model.deleteRows(2,del);
    93     // Now use build +-1
    94     CoinBuild buildObject2;
    95     time1 = CoinCpuTime();
    96     for ( k=0;k<10000;k++) {
    97       int row2Index[] = {0,1,2};
    98       double row2Value[]={1.0,-1.0,1.0};
    99       buildObject2.addRow(3,row2Index,row2Value,
    100                           1.0,1.0);
    101     }
    102     model.addRows(buildObject2,true);
    103     printf("Time for 10000 addRow using CoinBuild+-1 is %g\n",CoinCpuTime()-time1);
    104     model.dual();
    105     model=modelSave;
    106     model.deleteRows(2,del);
    107     // Now use build +-1
    108     CoinModel modelObject2;
    109     time1 = CoinCpuTime();
    110     for ( k=0;k<10000;k++) {
    111       int row2Index[] = {0,1,2};
    112       double row2Value[]={1.0,-1.0,1.0};
    113       modelObject2.addRow(3,row2Index,row2Value,
    114                             1.0,1.0);
    115     }
    116     model.addRows(modelObject2,true);
    117     printf("Time for 10000 addRow using CoinModel+-1 is %g\n",CoinCpuTime()-time1);
    118     model.dual();
    119     model=ClpSimplex();
    120     // Now use build +-1
    121     CoinModel modelObject3;
    122     time1 = CoinCpuTime();
    123     for ( k=0;k<10000;k++) {
    124       int row2Index[] = {0,1,2};
    125       double row2Value[]={1.0,-1.0,1.0};
    126       modelObject3.addRow(3,row2Index,row2Value,
    127                             1.0,1.0);
    128     }
    129     model.loadProblem(modelObject3,true);
    130     printf("Time for 10000 addRow using CoinModel load +-1 is %g\n",CoinCpuTime()-time1);
    131     model.writeMps("xx.mps");
    132     model.dual();
    133     model=modelSave;
    134     // Now use model
    135     CoinModel modelObject;
    136     time1 = CoinCpuTime();
    137     for ( k=0;k<10000;k++) {
    138       int row2Index[] = {0,1,2};
    139       double row2Value[]={1.0,-5.0,1.0};
    140       modelObject.addRow(3,row2Index,row2Value,
    141                          1.0,1.0);
    142     }
    143     model.addRows(modelObject);
    144     printf("Time for 10000 addRow using CoinModel is %g\n",CoinCpuTime()-time1);
    145     model.dual();
    146     model.writeMps("b.mps");
    147     // Method using least memory - but most complicated
    148     time1 = CoinCpuTime();
    149     // Assumes we know exact size of model and matrix
    150     // Empty model
    151     ClpSimplex  model2;
    152     {
    153       // Create space for 3 columns and 10000 rows
    154       int numberRows=10000;
    155       int numberColumns=3;
    156       // This is fully dense - but would not normally be so
    157       int numberElements = numberRows*numberColumns;
    158       // Arrays will be set to default values
    159       model2.resize(numberRows,numberColumns);
    160       double * elements = new double [numberElements];
    161       CoinBigIndex * starts = new CoinBigIndex [numberColumns+1];
    162       int * rows = new int [numberElements];;
    163       int * lengths = new int[numberColumns];
    164       // Now fill in - totally unsafe but ....
    165       // no need as defaults to 0.0 double * columnLower = model2.columnLower();
    166       double * columnUpper = model2.columnUpper();
    167       double * objective = model2.objective();
    168       double * rowLower = model2.rowLower();
    169       double * rowUpper = model2.rowUpper();
    170       // Columns - objective was packed
    171       for (k=0;k<2;k++) {
    172         int iColumn = objIndex[k];
    173         objective[iColumn]=objValue[k];
    174       }
    175       for (k=0;k<numberColumns;k++)
    176         columnUpper[k]=upper[k];
    177       // Rows
    178       for (k=0;k<numberRows;k++) {
    179         rowLower[k]=1.0;
    180         rowUpper[k]=1.0;
    181       }
    182       // Now elements
    183       double row2Value[]={1.0,-5.0,1.0};
    184       CoinBigIndex put=0;
    185       for (k=0;k<numberColumns;k++) {
    186         starts[k]=put;
    187         lengths[k]=numberRows;
    188         double value=row2Value[k];
    189         for (int i=0;i<numberRows;i++) {
    190           rows[put]=i;
    191           elements[put]=value;
    192           put++;
    193         }
    194       }
    195       starts[numberColumns]=put;
    196       // assign to matrix
    197       CoinPackedMatrix * matrix = new CoinPackedMatrix(true,0.0,0.0);
    198       matrix->assignMatrix(true,numberRows,numberColumns,numberElements,
    199                            elements,rows,starts,lengths);
    200       ClpPackedMatrix * clpMatrix = new ClpPackedMatrix(matrix);
    201       model2.replaceMatrix(clpMatrix,true);
    202       printf("Time for 10000 addRow using hand written code is %g\n",CoinCpuTime()-time1);
    203       // If matrix is really big could switch off creation of row copy
    204       // model2.setSpecialOptions(256);
    205     }
    206     model2.dual();
    207     model2.writeMps("a.mps");
    208     // Print column solution
    209     int numberColumns = model.numberColumns();
    210    
    211     // Alternatively getColSolution()
    212     double * columnPrimal = model.primalColumnSolution();
    213     // Alternatively getReducedCost()
    214     double * columnDual = model.dualColumnSolution();
    215     // Alternatively getColLower()
    216     double * columnLower = model.columnLower();
    217     // Alternatively getColUpper()
    218     double * columnUpper = model.columnUpper();
    219     // Alternatively getObjCoefficients()
    220     double * columnObjective = model.objective();
    221    
    222     int iColumn;
    223    
    224     std::cout<<"               Primal          Dual         Lower         Upper          Cost"
    225              <<std::endl;
    226    
    227     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    228       double value;
    229       std::cout<<std::setw(6)<<iColumn<<" ";
    230       value = columnPrimal[iColumn];
    231       if (fabs(value)<1.0e5)
    232         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    233       else
    234         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    235       value = columnDual[iColumn];
    236       if (fabs(value)<1.0e5)
    237         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    238       else
    239         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    240       value = columnLower[iColumn];
    241       if (fabs(value)<1.0e5)
    242         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    243       else
    244         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    245       value = columnUpper[iColumn];
    246       if (fabs(value)<1.0e5)
    247         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    248       else
    249         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    250       value = columnObjective[iColumn];
    251       if (fabs(value)<1.0e5)
    252         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    253       else
    254         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    255      
    256       std::cout<<std::endl;
    257     }
    258     std::cout<<"--------------------------------------"<<std::endl;
    259     // Test CoinAssert
    260     std::cout<<"If Clp compiled with -g below should give assert, if with -O1 or COIN_ASSERT CoinError"<<std::endl;
    261     model=modelSave;
    262     model.deleteRows(2,del);
    263     // Deliberate error
    264     model.deleteColumns(1,del+2);
    265     // Now use build +-1
    266     CoinBuild buildObject3;
    267     time1 = CoinCpuTime();
    268     for ( k=0;k<10000;k++) {
    269       int row2Index[] = {0,1,2};
    270       double row2Value[]={1.0,-1.0,1.0};
    271       buildObject3.addRow(3,row2Index,row2Value,
    272                           1.0,1.0);
    273     }
    274     model.addRows(buildObject3,true);
    275   }
    276   catch (CoinError e) {
    277     e.print();
    278     if (e.lineNumber()>=0)
    279       std::cout<<"This was from a CoinAssert"<<std::endl;
    280   }
    281   return 0;
    282 }   
     17     try {
     18          // Empty model
     19          ClpSimplex  model;
     20
     21          // Objective - just nonzeros
     22          int objIndex[] = {0, 2};
     23          double objValue[] = {1.0, 4.0};
     24          // Upper bounds - as dense vector
     25          double upper[] = {2.0, COIN_DBL_MAX, 4.0};
     26
     27          // Create space for 3 columns
     28          model.resize(0, 3);
     29          // Fill in
     30          int i;
     31          // Virtuous way
     32          // First objective
     33          for (i = 0; i < 2; i++)
     34               model.setObjectiveCoefficient(objIndex[i], objValue[i]);
     35          // Now bounds (lower will be zero by default but do again)
     36          for (i = 0; i < 3; i++) {
     37               model.setColumnLower(i, 0.0);
     38               model.setColumnUpper(i, upper[i]);
     39          }
     40          /*
     41            We could also have done in non-virtuous way e.g.
     42            double * objective = model.objective();
     43            and then set directly
     44          */
     45          // Faster to add rows all at once - but this is easier to show
     46          // Now add row 1 as >= 2.0
     47          int row1Index[] = {0, 2};
     48          double row1Value[] = {1.0, 1.0};
     49          model.addRow(2, row1Index, row1Value,
     50                       2.0, COIN_DBL_MAX);
     51          // Now add row 2 as == 1.0
     52          int row2Index[] = {0, 1, 2};
     53          double row2Value[] = {1.0, -5.0, 1.0};
     54          model.addRow(3, row2Index, row2Value,
     55                       1.0, 1.0);
     56          // solve
     57          model.dual();
     58
     59          /*
     60            Adding one row at a time has a significant overhead so let's
     61            try a more complicated but faster way
     62
     63            First time adding in 10000 rows one by one
     64          */
     65          model.allSlackBasis();
     66          ClpSimplex modelSave = model;
     67          double time1 = CoinCpuTime();
     68          int k;
     69          for ( k = 0; k < 10000; k++) {
     70               int row2Index[] = {0, 1, 2};
     71               double row2Value[] = {1.0, -5.0, 1.0};
     72               model.addRow(3, row2Index, row2Value,
     73                            1.0, 1.0);
     74          }
     75          printf("Time for 10000 addRow is %g\n", CoinCpuTime() - time1);
     76          model.dual();
     77          model = modelSave;
     78          // Now use build
     79          CoinBuild buildObject;
     80          time1 = CoinCpuTime();
     81          for ( k = 0; k < 10000; k++) {
     82               int row2Index[] = {0, 1, 2};
     83               double row2Value[] = {1.0, -5.0, 1.0};
     84               buildObject.addRow(3, row2Index, row2Value,
     85                                  1.0, 1.0);
     86          }
     87          model.addRows(buildObject);
     88          printf("Time for 10000 addRow using CoinBuild is %g\n", CoinCpuTime() - time1);
     89          model.dual();
     90          model = modelSave;
     91          int del[] = {0, 1, 2};
     92          model.deleteRows(2, del);
     93          // Now use build +-1
     94          CoinBuild buildObject2;
     95          time1 = CoinCpuTime();
     96          for ( k = 0; k < 10000; k++) {
     97               int row2Index[] = {0, 1, 2};
     98               double row2Value[] = {1.0, -1.0, 1.0};
     99               buildObject2.addRow(3, row2Index, row2Value,
     100                                   1.0, 1.0);
     101          }
     102          model.addRows(buildObject2, true);
     103          printf("Time for 10000 addRow using CoinBuild+-1 is %g\n", CoinCpuTime() - time1);
     104          model.dual();
     105          model = modelSave;
     106          model.deleteRows(2, del);
     107          // Now use build +-1
     108          CoinModel modelObject2;
     109          time1 = CoinCpuTime();
     110          for ( k = 0; k < 10000; k++) {
     111               int row2Index[] = {0, 1, 2};
     112               double row2Value[] = {1.0, -1.0, 1.0};
     113               modelObject2.addRow(3, row2Index, row2Value,
     114                                   1.0, 1.0);
     115          }
     116          model.addRows(modelObject2, true);
     117          printf("Time for 10000 addRow using CoinModel+-1 is %g\n", CoinCpuTime() - time1);
     118          model.dual();
     119          model = ClpSimplex();
     120          // Now use build +-1
     121          CoinModel modelObject3;
     122          time1 = CoinCpuTime();
     123          for ( k = 0; k < 10000; k++) {
     124               int row2Index[] = {0, 1, 2};
     125               double row2Value[] = {1.0, -1.0, 1.0};
     126               modelObject3.addRow(3, row2Index, row2Value,
     127                                   1.0, 1.0);
     128          }
     129          model.loadProblem(modelObject3, true);
     130          printf("Time for 10000 addRow using CoinModel load +-1 is %g\n", CoinCpuTime() - time1);
     131          model.writeMps("xx.mps");
     132          model.dual();
     133          model = modelSave;
     134          // Now use model
     135          CoinModel modelObject;
     136          time1 = CoinCpuTime();
     137          for ( k = 0; k < 10000; k++) {
     138               int row2Index[] = {0, 1, 2};
     139               double row2Value[] = {1.0, -5.0, 1.0};
     140               modelObject.addRow(3, row2Index, row2Value,
     141                                  1.0, 1.0);
     142          }
     143          model.addRows(modelObject);
     144          printf("Time for 10000 addRow using CoinModel is %g\n", CoinCpuTime() - time1);
     145          model.dual();
     146          model.writeMps("b.mps");
     147          // Method using least memory - but most complicated
     148          time1 = CoinCpuTime();
     149          // Assumes we know exact size of model and matrix
     150          // Empty model
     151          ClpSimplex  model2;
     152          {
     153               // Create space for 3 columns and 10000 rows
     154               int numberRows = 10000;
     155               int numberColumns = 3;
     156               // This is fully dense - but would not normally be so
     157               int numberElements = numberRows * numberColumns;
     158               // Arrays will be set to default values
     159               model2.resize(numberRows, numberColumns);
     160               double * elements = new double [numberElements];
     161               CoinBigIndex * starts = new CoinBigIndex [numberColumns+1];
     162               int * rows = new int [numberElements];;
     163               int * lengths = new int[numberColumns];
     164               // Now fill in - totally unsafe but ....
     165               // no need as defaults to 0.0 double * columnLower = model2.columnLower();
     166               double * columnUpper = model2.columnUpper();
     167               double * objective = model2.objective();
     168               double * rowLower = model2.rowLower();
     169               double * rowUpper = model2.rowUpper();
     170               // Columns - objective was packed
     171               for (k = 0; k < 2; k++) {
     172                    int iColumn = objIndex[k];
     173                    objective[iColumn] = objValue[k];
     174               }
     175               for (k = 0; k < numberColumns; k++)
     176                    columnUpper[k] = upper[k];
     177               // Rows
     178               for (k = 0; k < numberRows; k++) {
     179                    rowLower[k] = 1.0;
     180                    rowUpper[k] = 1.0;
     181               }
     182               // Now elements
     183               double row2Value[] = {1.0, -5.0, 1.0};
     184               CoinBigIndex put = 0;
     185               for (k = 0; k < numberColumns; k++) {
     186                    starts[k] = put;
     187                    lengths[k] = numberRows;
     188                    double value = row2Value[k];
     189                    for (int i = 0; i < numberRows; i++) {
     190                         rows[put] = i;
     191                         elements[put] = value;
     192                         put++;
     193                    }
     194               }
     195               starts[numberColumns] = put;
     196               // assign to matrix
     197               CoinPackedMatrix * matrix = new CoinPackedMatrix(true, 0.0, 0.0);
     198               matrix->assignMatrix(true, numberRows, numberColumns, numberElements,
     199                                    elements, rows, starts, lengths);
     200               ClpPackedMatrix * clpMatrix = new ClpPackedMatrix(matrix);
     201               model2.replaceMatrix(clpMatrix, true);
     202               printf("Time for 10000 addRow using hand written code is %g\n", CoinCpuTime() - time1);
     203               // If matrix is really big could switch off creation of row copy
     204               // model2.setSpecialOptions(256);
     205          }
     206          model2.dual();
     207          model2.writeMps("a.mps");
     208          // Print column solution
     209          int numberColumns = model.numberColumns();
     210
     211          // Alternatively getColSolution()
     212          double * columnPrimal = model.primalColumnSolution();
     213          // Alternatively getReducedCost()
     214          double * columnDual = model.dualColumnSolution();
     215          // Alternatively getColLower()
     216          double * columnLower = model.columnLower();
     217          // Alternatively getColUpper()
     218          double * columnUpper = model.columnUpper();
     219          // Alternatively getObjCoefficients()
     220          double * columnObjective = model.objective();
     221
     222          int iColumn;
     223
     224          std::cout << "               Primal          Dual         Lower         Upper          Cost"
     225                    << std::endl;
     226
     227          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     228               double value;
     229               std::cout << std::setw(6) << iColumn << " ";
     230               value = columnPrimal[iColumn];
     231               if (fabs(value) < 1.0e5)
     232                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     233               else
     234                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     235               value = columnDual[iColumn];
     236               if (fabs(value) < 1.0e5)
     237                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     238               else
     239                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     240               value = columnLower[iColumn];
     241               if (fabs(value) < 1.0e5)
     242                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     243               else
     244                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     245               value = columnUpper[iColumn];
     246               if (fabs(value) < 1.0e5)
     247                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     248               else
     249                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     250               value = columnObjective[iColumn];
     251               if (fabs(value) < 1.0e5)
     252                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     253               else
     254                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     255
     256               std::cout << std::endl;
     257          }
     258          std::cout << "--------------------------------------" << std::endl;
     259          // Test CoinAssert
     260          std::cout << "If Clp compiled with -g below should give assert, if with -O1 or COIN_ASSERT CoinError" << std::endl;
     261          model = modelSave;
     262          model.deleteRows(2, del);
     263          // Deliberate error
     264          model.deleteColumns(1, del + 2);
     265          // Now use build +-1
     266          CoinBuild buildObject3;
     267          time1 = CoinCpuTime();
     268          for ( k = 0; k < 10000; k++) {
     269               int row2Index[] = {0, 1, 2};
     270               double row2Value[] = {1.0, -1.0, 1.0};
     271               buildObject3.addRow(3, row2Index, row2Value,
     272                                   1.0, 1.0);
     273          }
     274          model.addRows(buildObject3, true);
     275     } catch (CoinError e) {
     276          e.print();
     277          if (e.lineNumber() >= 0)
     278               std::cout << "This was from a CoinAssert" << std::endl;
     279     }
     280     return 0;
     281}

  • trunk/Clp/examples/decomp2.cpp

    r1370 r1552  
    99int main (int argc, const char *argv[])
    1010{
    11   /* Create a structured model by reading mps file and trying
    12      Dantzig-Wolfe decomposition (that's the 1 parameter)
    13   */
    14   // At present D-W rows are hard coded - will move stuff from OSL
    15   CoinStructuredModel model((argc<2) ? "../../Data/Netlib/czprob.mps"
    16                         : argv[1],1);
    17   if (!model.numberRows())
    18     exit(10);
    19   // Get default solver - could change stuff
    20   ClpSimplex solver;
    21   /*
    22     This driver does a simple Dantzig Wolfe decomposition
    23   */
    24   solver.solve(&model);
    25   // Double check
    26   solver.primal(1);
    27   return 0;
    28 }   
     11     /* Create a structured model by reading mps file and trying
     12        Dantzig-Wolfe decomposition (that's the 1 parameter)
     13     */
     14     // At present D-W rows are hard coded - will move stuff from OSL
     15     CoinStructuredModel model((argc < 2) ? "../../Data/Netlib/czprob.mps"
     16                               : argv[1], 1);
     17     if (!model.numberRows())
     18          exit(10);
     19     // Get default solver - could change stuff
     20     ClpSimplex solver;
     21     /*
     22       This driver does a simple Dantzig Wolfe decomposition
     23     */
     24     solver.solve(&model);
     25     // Double check
     26     solver.primal(1);
     27     return 0;
     28}

  • trunk/Clp/examples/decomp3.cpp

    r1370 r1552  
    1111int main (int argc, const char *argv[])
    1212{
    13   /* Create a structured model by reading mps file and trying
    14      Dantzig-Wolfe or Benders decomposition
    15   */
    16   // Get maximum number of blocks
    17   int maxBlocks=50;
    18   if (argc>2)
    19     maxBlocks=atoi(argv[2]);
    20   int decompose=1;
    21   if (maxBlocks<0) {
    22     maxBlocks=-maxBlocks;
    23     decompose=2;
    24   }
    25   if (maxBlocks<2) {
    26     printf("Second parameters is maximum number of blocks >=2)\n");
    27     exit(5);
    28   } else {
    29     printf("Allowing at most %d blocks\n",maxBlocks);
    30   }
    31   //#define PRESOLVE
     13     /* Create a structured model by reading mps file and trying
     14        Dantzig-Wolfe or Benders decomposition
     15     */
     16     // Get maximum number of blocks
     17     int maxBlocks = 50;
     18     if (argc > 2)
     19          maxBlocks = atoi(argv[2]);
     20     int decompose = 1;
     21     if (maxBlocks < 0) {
     22          maxBlocks = -maxBlocks;
     23          decompose = 2;
     24     }
     25     if (maxBlocks < 2) {
     26          printf("Second parameters is maximum number of blocks >=2)\n");
     27          exit(5);
     28     } else {
     29          printf("Allowing at most %d blocks\n", maxBlocks);
     30     }
     31     //#define PRESOLVE
    3232#ifndef PRESOLVE
    33   CoinStructuredModel model((argc<2) ? "../../Data/Netlib/czprob.mps"
    34                             : argv[1],decompose,maxBlocks);
    35   if (!model.numberRows())
    36     exit(10);
    37   // Get default solver - could change stuff
    38   ClpSimplex solver;
    39   // change factorization frequency from 200
    40   solver.setFactorizationFrequency(100+model.numberRows()/50);
    41   /*
    42     This driver does a simple Dantzig Wolfe decomposition
    43   */
    44   double time1 = CoinCpuTime() ;
    45   solver.solve(&model);
    46   std::cout<<"model took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    47   // Double check
    48   solver.primal(1);
     33     CoinStructuredModel model((argc < 2) ? "../../Data/Netlib/czprob.mps"
     34                               : argv[1], decompose, maxBlocks);
     35     if (!model.numberRows())
     36          exit(10);
     37     // Get default solver - could change stuff
     38     ClpSimplex solver;
     39     // change factorization frequency from 200
     40     solver.setFactorizationFrequency(100 + model.numberRows() / 50);
     41     /*
     42       This driver does a simple Dantzig Wolfe decomposition
     43     */
     44     double time1 = CoinCpuTime() ;
     45     solver.solve(&model);
     46     std::cout << "model took " << CoinCpuTime() - time1 << " seconds" << std::endl;
     47     // Double check
     48     solver.primal(1);
    4949#else
    50   ClpSimplex  model;
    51   int status = model.readMps((argc<2) ? "../../Data/Netlib/czprob.mps"
    52                              : argv[1]);
    53   if (status) {
    54     fprintf(stdout,"Bad readMps %s\n",argv[1]);
    55     exit(1);
    56   }
    57   ClpSimplex * model2=&model;
     50     ClpSimplex  model;
     51     int status = model.readMps((argc < 2) ? "../../Data/Netlib/czprob.mps"
     52                                : argv[1]);
     53     if (status) {
     54          fprintf(stdout, "Bad readMps %s\n", argv[1]);
     55          exit(1);
     56     }
     57     ClpSimplex * model2 = &model;
    5858
    59   CoinStructuredModel modelB;
    60   modelB.decompose(*model2->matrix(),
    61                    model2->rowLower(),model2->rowUpper(),
    62                    model2->columnLower(),model2->columnUpper(),
    63                    model2->objective(),1,maxBlocks,
    64                    model2->objectiveOffset());
    65   // change factorization frequency from 200
    66   model2->setFactorizationFrequency(100+modelB.numberRows()/50);
    67   /*
    68     This driver does a simple Dantzig Wolfe decomposition
    69   */
    70   double time1 = CoinCpuTime() ;
    71   model2->solve(&modelB);
    72   std::cout<<"model took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    73   // But resolve for safety
    74   model.primal(1);
     59     CoinStructuredModel modelB;
     60     modelB.decompose(*model2->matrix(),
     61                      model2->rowLower(), model2->rowUpper(),
     62                      model2->columnLower(), model2->columnUpper(),
     63                      model2->objective(), 1, maxBlocks,
     64                      model2->objectiveOffset());
     65     // change factorization frequency from 200
     66     model2->setFactorizationFrequency(100 + modelB.numberRows() / 50);
     67     /*
     68       This driver does a simple Dantzig Wolfe decomposition
     69     */
     70     double time1 = CoinCpuTime() ;
     71     model2->solve(&modelB);
     72     std::cout << "model took " << CoinCpuTime() - time1 << " seconds" << std::endl;
     73     // But resolve for safety
     74     model.primal(1);
    7575#endif
    76   return 0;
    77   ClpSimplex solver2;
    78   solver2.readMps((argc<2) ? "../../Data/Netlib/czprob.mps" : argv[1]);
    79   time1 = CoinCpuTime() ;
    80   solver2.dual();
    81   std::cout<<"second try took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
    82   return 0;
    83 }   
     76     return 0;
     77     ClpSimplex solver2;
     78     solver2.readMps((argc < 2) ? "../../Data/Netlib/czprob.mps" : argv[1]);
     79     time1 = CoinCpuTime() ;
     80     solver2.dual();
     81     std::cout << "second try took " << CoinCpuTime() - time1 << " seconds" << std::endl;
     82     return 0;
     83}

  • 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}

  • trunk/Clp/examples/defaults.cpp

    r1370 r1552  
    1010int main (int argc, const char *argv[])
    1111{
    12   ClpSimplex  model;
    13   int status;
    14   // Keep names
    15   if (argc<2)
    16     status=model.readMps("../../Data/Sample/p0033.mps",true);
    17   else
    18     status=model.readMps(argv[1],true);
    19   /*
    20     This driver is similar to minimum.cpp, but it
    21     sets all parameter values to their defaults.  The purpose of this
    22     is to give a list of most of the methods that the end user will need.
    23 
    24     There are also some more methods as for OsiSolverInterface e.g.
    25     some loadProblem methods and deleteRows and deleteColumns.
    26 
    27     Often two methods do the same thing, one having a name I like
    28     while the other adheres to OsiSolverInterface standards
    29   */
    30 
    31   // Use exact devex ( a variant of steepest edge)
    32   ClpPrimalColumnSteepest primalSteepest;
    33   model.setPrimalColumnPivotAlgorithm(primalSteepest);
    34 
    35 
    36   int integerValue;
    37   double value;
    38 
    39   // Infeasibility cost
    40   value = model.infeasibilityCost();
    41   std::cout<<"Default value of infeasibility cost is "<<value<<std::endl;
    42   model.setInfeasibilityCost(value);
    43 
    44   if (!status) {
    45     model.primal();
    46   }
    47   // Number of rows and columns - also getNumRows, getNumCols
    48   std::string modelName;
    49   model.getStrParam(ClpProbName,modelName);
    50   std::cout<<"Model "<<modelName<<" has "<<model.numberRows()<<" rows and "<<
    51     model.numberColumns()<<" columns"<<std::endl;
    52 
    53   /* Some parameters as in OsiSolverParameters.  ObjectiveLimits
    54      are not active yet.  dualTolerance, setDualTolerance,
    55      primalTolerance and setPrimalTolerance may be used as well */
    56 
    57   model.getDblParam(ClpDualObjectiveLimit,value);
    58   std::cout<<"Value of ClpDualObjectiveLimit is "<<value<<std::endl;
    59   model.getDblParam(ClpPrimalObjectiveLimit,value);
    60   std::cout<<"Value of ClpPrimalObjectiveLimit is "<<value<<std::endl;
    61   model.getDblParam(ClpDualTolerance,value);
    62   std::cout<<"Value of ClpDualTolerance is "<<value<<std::endl;
    63   model.getDblParam(ClpPrimalTolerance,value);
    64   std::cout<<"Value of ClpPrimalTolerance is "<<value<<std::endl;
    65   model.getDblParam(ClpObjOffset,value);
    66   std::cout<<"Value of ClpObjOffset is "<<value<<std::endl;
    67 
    68   // setDblParam(ClpPrimalTolerance) is same as this
    69   model.getDblParam(ClpPrimalTolerance,value);
    70   model.setPrimalTolerance(value);
    71 
    72   model.setDualTolerance( model.dualTolerance()) ;
    73 
    74   // Other Param stuff
    75 
    76   // Can also use maximumIterations
    77   model.getIntParam(ClpMaxNumIteration,integerValue);
    78   std::cout<<"Value of ClpMaxNumIteration is "<<integerValue<<std::endl;
    79   model.setMaximumIterations(integerValue);
    80 
    81   // Not sure this works yet
    82   model.getIntParam(ClpMaxNumIterationHotStart,integerValue);
    83   std::cout<<"Value of ClpMaxNumIterationHotStart is "
    84            <<integerValue<<std::endl;
    85 
    86   // Can also use getIterationCount and getObjValue
    87   /* Status of problem:
    88       0 - optimal
    89       1 - primal infeasible
    90       2 - dual infeasible
    91       3 - stopped on iterations etc
    92       4 - stopped due to errors
    93   */
    94   std::cout<<"Model status is "<<model.status()<<" after "
    95            <<model.numberIterations()<<" iterations - objective is "
    96            <<model.objectiveValue()<<std::endl;
    97 
    98   assert(!model.isAbandoned());
    99   assert(model.isProvenOptimal());
    100   assert(!model.isProvenPrimalInfeasible());
    101   assert(!model.isProvenDualInfeasible());
    102   assert(!model.isPrimalObjectiveLimitReached());
    103   assert(!model.isDualObjectiveLimitReached());
    104   assert(!model.isIterationLimitReached());
    105 
    106 
    107   // Things to help you determine if optimal
    108   assert(model.primalFeasible());
    109   assert(!model.numberPrimalInfeasibilities());
    110   assert(model.sumPrimalInfeasibilities()<1.0e-7);
    111   assert(model.dualFeasible());
    112   assert(!model.numberDualInfeasibilities());
    113   assert(model.sumDualInfeasibilities()<1.0e-7);
    114 
    115   // Save warm start and set to all slack
    116   unsigned char * basis1 = model.statusCopy();
    117   model.createStatus();
    118 
    119   // Now create another model and do hot start
    120   ClpSimplex model2=model;
    121   model2.copyinStatus(basis1);
    122   delete [] basis1;
    123 
    124   // Check model has not got basis (should iterate)
    125   model.dual();
    126  
    127   // Can use getObjSense
    128   model2.setOptimizationDirection(model.optimizationDirection());
    129 
    130   // Can use scalingFlag() to check if scaling on
    131   // But set up scaling
    132   model2.scaling();
    133 
    134   // Could play with sparse factorization on/off
    135   model2.setSparseFactorization(model.sparseFactorization());
    136 
    137   // Sets row pivot choice algorithm in dual
    138   ClpDualRowSteepest dualSteepest;
    139   model2.setDualRowPivotAlgorithm(dualSteepest);
    140 
    141   // Dual bound (i.e. dual infeasibility cost)
    142   value = model.dualBound();
    143   std::cout<<"Default value of dual bound is "<<value<<std::endl;
    144   model.setDualBound(value);
    145 
    146   // Do some deafult message handling
    147   // To see real use - see OsiOslSolverInterfaceTest.cpp
    148   CoinMessageHandler handler;
    149   model2.passInMessageHandler(& handler);
    150   model2.newLanguage(CoinMessages::us_en);
    151 
    152   //Increase level of detail
    153   model2.setLogLevel(4);
    154 
    155   // solve
    156   model2.dual();
    157   // flip direction twice and solve
    158   model2.setOptimizationDirection(-1);
    159   model2.dual();
    160   model2.setOptimizationDirection(1);
    161   //Decrease level of detail
    162   model2.setLogLevel(1);
    163   model2.dual();
    164 
    165   /*
    166     Now for getting at information.  This will not deal with:
    167 
    168     ClpMatrixBase * rowCopy() and ClpMatrixbase * clpMatrix()
    169 
    170     nor with
    171 
    172     double * infeasibilityRay() and double * unboundedRay()
    173     (NULL returned if none/wrong)
    174     Up to user to use delete [] on these arrays.
    175 
    176    */
    177 
    178 
    179   /*
    180     Now to print out solution.  The methods used return modifiable
    181     arrays while the alternative names return const pointers -
    182     which is of course much more virtuous
    183  
    184    */
    185 
    186   int numberRows = model2.numberRows();
    187 
    188   // Alternatively getRowActivity()
    189   double * rowPrimal = model2.primalRowSolution();
    190   // Alternatively getRowPrice()
    191   double * rowDual = model2.dualRowSolution();
    192   // Alternatively getRowLower()
    193   double * rowLower = model2.rowLower();
    194   // Alternatively getRowUpper()
    195   double * rowUpper = model2.rowUpper();
    196   // Alternatively getRowObjCoefficients()
    197   double * rowObjective = model2.rowObjective();
    198    
    199   // If we have not kept names (parameter to readMps) this will be 0
    200   assert(model2.lengthNames());
    201 
    202   // Row names
    203   const std::vector<std::string> * rowNames = model2.rowNames();
    204 
    205 
    206   int iRow;
    207 
    208   std::cout<<"                       Primal          Dual         Lower         Upper        (Cost)"
    209            <<std::endl;
    210 
    211   for (iRow=0;iRow<numberRows;iRow++) {
    212     double value;
    213     std::cout<<std::setw(6)<<iRow<<" "<<std::setw(8)<<(*rowNames)[iRow];
    214     value = rowPrimal[iRow];
    215     if (fabs(value)<1.0e5)
    216       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    217     else
    218       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    219     value = rowDual[iRow];
    220     if (fabs(value)<1.0e5)
    221       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    222     else
    223       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    224     value = rowLower[iRow];
    225     if (fabs(value)<1.0e5)
    226       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    227     else
    228       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    229     value = rowUpper[iRow];
    230     if (fabs(value)<1.0e5)
    231       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    232     else
    233       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    234     if (rowObjective) {
    235       value = rowObjective[iRow];
    236       if (fabs(value)<1.0e5)
    237         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    238       else
    239         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    240     }
    241     std::cout<<std::endl;
    242   }
    243   std::cout<<"--------------------------------------"<<std::endl;
    244 
    245   // Columns
    246 
    247   int numberColumns = model2.numberColumns();
    248 
    249   // Alternatively getColSolution()
    250   double * columnPrimal = model2.primalColumnSolution();
    251   // Alternatively getReducedCost()
    252   double * columnDual = model2.dualColumnSolution();
    253   // Alternatively getColLower()
    254   double * columnLower = model2.columnLower();
    255   // Alternatively getColUpper()
    256   double * columnUpper = model2.columnUpper();
    257   // Alternatively getObjCoefficients()
    258   double * columnObjective = model2.objective();
    259    
    260   // If we have not kept names (parameter to readMps) this will be 0
    261   assert(model2.lengthNames());
    262 
    263   // Column names
    264   const std::vector<std::string> * columnNames = model2.columnNames();
    265 
    266 
    267   int iColumn;
    268  
    269   std::cout<<"                       Primal          Dual         Lower         Upper          Cost"
    270            <<std::endl;
    271 
    272   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    273     double value;
    274     std::cout<<std::setw(6)<<iColumn<<" "<<std::setw(8)<<(*columnNames)[iColumn];
    275     value = columnPrimal[iColumn];
    276     if (fabs(value)<1.0e5)
    277       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    278     else
    279       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    280     value = columnDual[iColumn];
    281     if (fabs(value)<1.0e5)
    282       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    283     else
    284       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    285     value = columnLower[iColumn];
    286     if (fabs(value)<1.0e5)
    287       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    288     else
    289       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    290     value = columnUpper[iColumn];
    291     if (fabs(value)<1.0e5)
    292       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    293     else
    294       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    295     value = columnObjective[iColumn];
    296     if (fabs(value)<1.0e5)
    297       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    298     else
    299       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    300 
    301     std::cout<<std::endl;
    302   }
    303   std::cout<<"--------------------------------------"<<std::endl;
    304 
    305   std::cout<<resetiosflags(std::ios::fixed|std::ios::showpoint|std::ios::scientific);
    306 
    307   // Now matrix
    308   CoinPackedMatrix * matrix = model2.matrix();
    309 
    310   const double * element = matrix->getElements();
    311   const int * row = matrix->getIndices();
    312   const int * start = matrix->getVectorStarts();
    313   const int * length = matrix->getVectorLengths();
    314 
    315   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    316     std::cout<<"Column "<<iColumn;
    317     int j;
    318     for (j=start[iColumn];j<start[iColumn]+length[iColumn];j++)
    319       std::cout<<" ( "<<row[j]<<", "<<element[j]<<")";
    320     std::cout<<std::endl;
    321   }
    322   return 0;
    323 }   
     12     ClpSimplex  model;
     13     int status;
     14     // Keep names
     15     if (argc < 2)
     16          status = model.readMps("../../Data/Sample/p0033.mps", true);
     17     else
     18          status = model.readMps(argv[1], true);
     19     /*
     20       This driver is similar to minimum.cpp, but it
     21       sets all parameter values to their defaults.  The purpose of this
     22       is to give a list of most of the methods that the end user will need.
     23
     24       There are also some more methods as for OsiSolverInterface e.g.
     25       some loadProblem methods and deleteRows and deleteColumns.
     26
     27       Often two methods do the same thing, one having a name I like
     28       while the other adheres to OsiSolverInterface standards
     29     */
     30
     31     // Use exact devex ( a variant of steepest edge)
     32     ClpPrimalColumnSteepest primalSteepest;
     33     model.setPrimalColumnPivotAlgorithm(primalSteepest);
     34
     35
     36     int integerValue;
     37     double value;
     38
     39     // Infeasibility cost
     40     value = model.infeasibilityCost();
     41     std::cout << "Default value of infeasibility cost is " << value << std::endl;
     42     model.setInfeasibilityCost(value);
     43
     44     if (!status) {
     45          model.primal();
     46     }
     47     // Number of rows and columns - also getNumRows, getNumCols
     48     std::string modelName;
     49     model.getStrParam(ClpProbName, modelName);
     50     std::cout << "Model " << modelName << " has " << model.numberRows() << " rows and " <<
     51               model.numberColumns() << " columns" << std::endl;
     52
     53     /* Some parameters as in OsiSolverParameters.  ObjectiveLimits
     54        are not active yet.  dualTolerance, setDualTolerance,
     55        primalTolerance and setPrimalTolerance may be used as well */
     56
     57     model.getDblParam(ClpDualObjectiveLimit, value);
     58     std::cout << "Value of ClpDualObjectiveLimit is " << value << std::endl;
     59     model.getDblParam(ClpPrimalObjectiveLimit, value);
     60     std::cout << "Value of ClpPrimalObjectiveLimit is " << value << std::endl;
     61     model.getDblParam(ClpDualTolerance, value);
     62     std::cout << "Value of ClpDualTolerance is " << value << std::endl;
     63     model.getDblParam(ClpPrimalTolerance, value);
     64     std::cout << "Value of ClpPrimalTolerance is " << value << std::endl;
     65     model.getDblParam(ClpObjOffset, value);
     66     std::cout << "Value of ClpObjOffset is " << value << std::endl;
     67
     68     // setDblParam(ClpPrimalTolerance) is same as this
     69     model.getDblParam(ClpPrimalTolerance, value);
     70     model.setPrimalTolerance(value);
     71
     72     model.setDualTolerance( model.dualTolerance()) ;
     73
     74     // Other Param stuff
     75
     76     // Can also use maximumIterations
     77     model.getIntParam(ClpMaxNumIteration, integerValue);
     78     std::cout << "Value of ClpMaxNumIteration is " << integerValue << std::endl;
     79     model.setMaximumIterations(integerValue);
     80
     81     // Not sure this works yet
     82     model.getIntParam(ClpMaxNumIterationHotStart, integerValue);
     83     std::cout << "Value of ClpMaxNumIterationHotStart is "
     84               << integerValue << std::endl;
     85
     86     // Can also use getIterationCount and getObjValue
     87     /* Status of problem:
     88         0 - optimal
     89         1 - primal infeasible
     90         2 - dual infeasible
     91         3 - stopped on iterations etc
     92         4 - stopped due to errors
     93     */
     94     std::cout << "Model status is " << model.status() << " after "
     95               << model.numberIterations() << " iterations - objective is "
     96               << model.objectiveValue() << std::endl;
     97
     98     assert(!model.isAbandoned());
     99     assert(model.isProvenOptimal());
     100     assert(!model.isProvenPrimalInfeasible());
     101     assert(!model.isProvenDualInfeasible());
     102     assert(!model.isPrimalObjectiveLimitReached());
     103     assert(!model.isDualObjectiveLimitReached());
     104     assert(!model.isIterationLimitReached());
     105
     106
     107     // Things to help you determine if optimal
     108     assert(model.primalFeasible());
     109     assert(!model.numberPrimalInfeasibilities());
     110     assert(model.sumPrimalInfeasibilities() < 1.0e-7);
     111     assert(model.dualFeasible());
     112     assert(!model.numberDualInfeasibilities());
     113     assert(model.sumDualInfeasibilities() < 1.0e-7);
     114
     115     // Save warm start and set to all slack
     116     unsigned char * basis1 = model.statusCopy();
     117     model.createStatus();
     118
     119     // Now create another model and do hot start
     120     ClpSimplex model2 = model;
     121     model2.copyinStatus(basis1);
     122     delete [] basis1;
     123
     124     // Check model has not got basis (should iterate)
     125     model.dual();
     126
     127     // Can use getObjSense
     128     model2.setOptimizationDirection(model.optimizationDirection());
     129
     130     // Can use scalingFlag() to check if scaling on
     131     // But set up scaling
     132     model2.scaling();
     133
     134     // Could play with sparse factorization on/off
     135     model2.setSparseFactorization(model.sparseFactorization());
     136
     137     // Sets row pivot choice algorithm in dual
     138     ClpDualRowSteepest dualSteepest;
     139     model2.setDualRowPivotAlgorithm(dualSteepest);
     140
     141     // Dual bound (i.e. dual infeasibility cost)
     142     value = model.dualBound();
     143     std::cout << "Default value of dual bound is " << value << std::endl;
     144     model.setDualBound(value);
     145
     146     // Do some deafult message handling
     147     // To see real use - see OsiOslSolverInterfaceTest.cpp
     148     CoinMessageHandler handler;
     149     model2.passInMessageHandler(& handler);
     150     model2.newLanguage(CoinMessages::us_en);
     151
     152     //Increase level of detail
     153     model2.setLogLevel(4);
     154
     155     // solve
     156     model2.dual();
     157     // flip direction twice and solve
     158     model2.setOptimizationDirection(-1);
     159     model2.dual();
     160     model2.setOptimizationDirection(1);
     161     //Decrease level of detail
     162     model2.setLogLevel(1);
     163     model2.dual();
     164
     165     /*
     166       Now for getting at information.  This will not deal with:
     167
     168       ClpMatrixBase * rowCopy() and ClpMatrixbase * clpMatrix()
     169
     170       nor with
     171
     172       double * infeasibilityRay() and double * unboundedRay()
     173       (NULL returned if none/wrong)
     174       Up to user to use delete [] on these arrays.
     175
     176      */
     177
     178
     179     /*
     180       Now to print out solution.  The methods used return modifiable
     181       arrays while the alternative names return const pointers -
     182       which is of course much more virtuous
     183
     184      */
     185
     186     int numberRows = model2.numberRows();
     187
     188     // Alternatively getRowActivity()
     189     double * rowPrimal = model2.primalRowSolution();
     190     // Alternatively getRowPrice()
     191     double * rowDual = model2.dualRowSolution();
     192     // Alternatively getRowLower()
     193     double * rowLower = model2.rowLower();
     194     // Alternatively getRowUpper()
     195     double * rowUpper = model2.rowUpper();
     196     // Alternatively getRowObjCoefficients()
     197     double * rowObjective = model2.rowObjective();
     198
     199     // If we have not kept names (parameter to readMps) this will be 0
     200     assert(model2.lengthNames());
     201
     202     // Row names
     203     const std::vector<std::string> * rowNames = model2.rowNames();
     204
     205
     206     int iRow;
     207
     208     std::cout << "                       Primal          Dual         Lower         Upper        (Cost)"
     209               << std::endl;
     210
     211     for (iRow = 0; iRow < numberRows; iRow++) {
     212          double value;
     213          std::cout << std::setw(6) << iRow << " " << std::setw(8) << (*rowNames)[iRow];
     214          value = rowPrimal[iRow];
     215          if (fabs(value) < 1.0e5)
     216               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     217          else
     218               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     219          value = rowDual[iRow];
     220          if (fabs(value) < 1.0e5)
     221               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     222          else
     223               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     224          value = rowLower[iRow];
     225          if (fabs(value) < 1.0e5)
     226               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     227          else
     228               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     229          value = rowUpper[iRow];
     230          if (fabs(value) < 1.0e5)
     231               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     232          else
     233               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     234          if (rowObjective) {
     235               value = rowObjective[iRow];
     236               if (fabs(value) < 1.0e5)
     237                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     238               else
     239                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     240          }
     241          std::cout << std::endl;
     242     }
     243     std::cout << "--------------------------------------" << std::endl;
     244
     245     // Columns
     246
     247     int numberColumns = model2.numberColumns();
     248
     249     // Alternatively getColSolution()
     250     double * columnPrimal = model2.primalColumnSolution();
     251     // Alternatively getReducedCost()
     252     double * columnDual = model2.dualColumnSolution();
     253     // Alternatively getColLower()
     254     double * columnLower = model2.columnLower();
     255     // Alternatively getColUpper()
     256     double * columnUpper = model2.columnUpper();
     257     // Alternatively getObjCoefficients()
     258     double * columnObjective = model2.objective();
     259
     260     // If we have not kept names (parameter to readMps) this will be 0
     261     assert(model2.lengthNames());
     262
     263     // Column names
     264     const std::vector<std::string> * columnNames = model2.columnNames();
     265
     266
     267     int iColumn;
     268
     269     std::cout << "                       Primal          Dual         Lower         Upper          Cost"
     270               << std::endl;
     271
     272     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     273          double value;
     274          std::cout << std::setw(6) << iColumn << " " << std::setw(8) << (*columnNames)[iColumn];
     275          value = columnPrimal[iColumn];
     276          if (fabs(value) < 1.0e5)
     277               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     278          else
     279               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     280          value = columnDual[iColumn];
     281          if (fabs(value) < 1.0e5)
     282               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     283          else
     284               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     285          value = columnLower[iColumn];
     286          if (fabs(value) < 1.0e5)
     287               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     288          else
     289               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     290          value = columnUpper[iColumn];
     291          if (fabs(value) < 1.0e5)
     292               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     293          else
     294               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     295          value = columnObjective[iColumn];
     296          if (fabs(value) < 1.0e5)
     297               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     298          else
     299               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     300
     301          std::cout << std::endl;
     302     }
     303     std::cout << "--------------------------------------" << std::endl;
     304
     305     std::cout << resetiosflags(std::ios::fixed | std::ios::showpoint | std::ios::scientific);
     306
     307     // Now matrix
     308     CoinPackedMatrix * matrix = model2.matrix();
     309
     310     const double * element = matrix->getElements();
     311     const int * row = matrix->getIndices();
     312     const int * start = matrix->getVectorStarts();
     313     const int * length = matrix->getVectorLengths();
     314
     315     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     316          std::cout << "Column " << iColumn;
     317          int j;
     318          for (j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++)
     319               std::cout << " ( " << row[j] << ", " << element[j] << ")";
     320          std::cout << std::endl;
     321     }
     322     return 0;
     323}

  • trunk/Clp/examples/driver.cpp

    r1370 r1552  
    88int main (int argc, const char *argv[])
    99{
    10   ClpSimplex  model;
    11   int status;
    12   // Keep names when reading an mps file
    13   if (argc<2)
    14     status=model.readMps("../../Data/Sample/p0033.mps",true);
    15   else
    16     status=model.readMps(argv[1],true);
     10     ClpSimplex  model;
     11     int status;
     12     // Keep names when reading an mps file
     13     if (argc < 2)
     14          status = model.readMps("../../Data/Sample/p0033.mps", true);
     15     else
     16          status = model.readMps(argv[1], true);
    1717
    18   if (status) {
    19     fprintf(stderr,"Bad readMps %s\n",argv[1]);
    20     fprintf(stdout,"Bad readMps %s\n",argv[1]);
    21     exit(1);
    22   }
     18     if (status) {
     19          fprintf(stderr, "Bad readMps %s\n", argv[1]);
     20          fprintf(stdout, "Bad readMps %s\n", argv[1]);
     21          exit(1);
     22     }
    2323#ifdef STYLE1
    24   if (argc<3 ||!strstr(argv[2],"primal")) {
    25     // Use the dual algorithm unless user said "primal"
    26     model.initialDualSolve();
    27   } else {
    28     model.initialPrimalSolve();
    29   }
     24     if (argc < 3 || !strstr(argv[2], "primal")) {
     25          // Use the dual algorithm unless user said "primal"
     26          model.initialDualSolve();
     27     } else {
     28          model.initialPrimalSolve();
     29     }
    3030#else
    31   ClpSolve solvectl;
     31     ClpSolve solvectl;
    3232
    3333
    34   if (argc<3 ||(!strstr(argv[2],"primal")&&!strstr(argv[2],"barrier"))) {
    35     // Use the dual algorithm unless user said "primal" or "barrier"
    36     std::cout << std::endl << " Solve using Dual: " << std::endl;
    37     solvectl.setSolveType(ClpSolve::useDual);
    38     solvectl.setPresolveType(ClpSolve::presolveOn);
    39     model.initialSolve(solvectl);
    40   } else if (strstr(argv[2],"barrier")) {
    41     // Use the barrier algorithm if user said "barrier"
    42     std::cout << std::endl << " Solve using Barrier: " << std::endl;
    43     solvectl.setSolveType(ClpSolve::useBarrier);
    44     solvectl.setPresolveType(ClpSolve::presolveOn);
    45     model.initialSolve(solvectl);
    46   } else {
    47     std::cout << std::endl << " Solve using Primal: " << std::endl;
    48     solvectl.setSolveType(ClpSolve::usePrimal);
    49     solvectl.setPresolveType(ClpSolve::presolveOn);
    50     model.initialSolve(solvectl);
    51   }
     34     if (argc < 3 || (!strstr(argv[2], "primal") && !strstr(argv[2], "barrier"))) {
     35          // Use the dual algorithm unless user said "primal" or "barrier"
     36          std::cout << std::endl << " Solve using Dual: " << std::endl;
     37          solvectl.setSolveType(ClpSolve::useDual);
     38          solvectl.setPresolveType(ClpSolve::presolveOn);
     39          model.initialSolve(solvectl);
     40     } else if (strstr(argv[2], "barrier")) {
     41          // Use the barrier algorithm if user said "barrier"
     42          std::cout << std::endl << " Solve using Barrier: " << std::endl;
     43          solvectl.setSolveType(ClpSolve::useBarrier);
     44          solvectl.setPresolveType(ClpSolve::presolveOn);
     45          model.initialSolve(solvectl);
     46     } else {
     47          std::cout << std::endl << " Solve using Primal: " << std::endl;
     48          solvectl.setSolveType(ClpSolve::usePrimal);
     49          solvectl.setPresolveType(ClpSolve::presolveOn);
     50          model.initialSolve(solvectl);
     51     }
    5252#endif
    53   std::string modelName;
    54   model.getStrParam(ClpProbName,modelName);
    55   std::cout<<"Model "<<modelName<<" has "<<model.numberRows()<<" rows and "<<
    56     model.numberColumns()<<" columns"<<std::endl;
     53     std::string modelName;
     54     model.getStrParam(ClpProbName, modelName);
     55     std::cout << "Model " << modelName << " has " << model.numberRows() << " rows and " <<
     56               model.numberColumns() << " columns" << std::endl;
    5757
    58   // remove this to print solution
     58     // remove this to print solution
    5959
    60   exit(0);
     60     exit(0);
    6161
    62   /*
    63     Now to print out solution.  The methods used return modifiable
    64     arrays while the alternative names return const pointers -
    65     which is of course much more virtuous.
     62     /*
     63       Now to print out solution.  The methods used return modifiable
     64       arrays while the alternative names return const pointers -
     65       which is of course much more virtuous.
    6666
    67     This version just does non-zero columns
    68  
    69    */
     67       This version just does non-zero columns
     68
     69      */
    7070#if 0
    71   int numberRows = model.numberRows();
     71     int numberRows = model.numberRows();
    7272
    73   // Alternatively getRowActivity()
    74   double * rowPrimal = model.primalRowSolution();
    75   // Alternatively getRowPrice()
    76   double * rowDual = model.dualRowSolution();
    77   // Alternatively getRowLower()
    78   double * rowLower = model.rowLower();
    79   // Alternatively getRowUpper()
    80   double * rowUpper = model.rowUpper();
    81   // Alternatively getRowObjCoefficients()
    82   double * rowObjective = model.rowObjective();
    83    
    84   // If we have not kept names (parameter to readMps) this will be 0
    85   assert(model.lengthNames());
     73     // Alternatively getRowActivity()
     74     double * rowPrimal = model.primalRowSolution();
     75     // Alternatively getRowPrice()
     76     double * rowDual = model.dualRowSolution();
     77     // Alternatively getRowLower()
     78     double * rowLower = model.rowLower();
     79     // Alternatively getRowUpper()
     80     double * rowUpper = model.rowUpper();
     81     // Alternatively getRowObjCoefficients()
     82     double * rowObjective = model.rowObjective();
    8683
    87   // Row names
    88   const std::vector<std::string> * rowNames = model.rowNames();
     84     // If we have not kept names (parameter to readMps) this will be 0
     85     assert(model.lengthNames());
     86
     87     // Row names
     88     const std::vector<std::string> * rowNames = model.rowNames();
    8989
    9090
    91   int iRow;
     91     int iRow;
    9292
    93   std::cout<<"                       Primal          Dual         Lower         Upper        (Cost)"
    94            <<std::endl;
     93     std::cout << "                       Primal          Dual         Lower         Upper        (Cost)"
     94               << std::endl;
    9595
    96   for (iRow=0;iRow<numberRows;iRow++) {
    97     double value;
    98     std::cout<<std::setw(6)<<iRow<<" "<<std::setw(8)<<(*rowNames)[iRow];
    99     value = rowPrimal[iRow];
    100     if (fabs(value)<1.0e5)
    101       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    102     else
    103       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    104     value = rowDual[iRow];
    105     if (fabs(value)<1.0e5)
    106       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    107     else
    108       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    109     value = rowLower[iRow];
    110     if (fabs(value)<1.0e5)
    111       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    112     else
    113       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    114     value = rowUpper[iRow];
    115     if (fabs(value)<1.0e5)
    116       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    117     else
    118       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    119     if (rowObjective) {
    120       value = rowObjective[iRow];
    121       if (fabs(value)<1.0e5)
    122         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    123       else
    124         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    125     }
    126     std::cout<<std::endl;
    127   }
     96     for (iRow = 0; iRow < numberRows; iRow++) {
     97          double value;
     98          std::cout << std::setw(6) << iRow << " " << std::setw(8) << (*rowNames)[iRow];
     99          value = rowPrimal[iRow];
     100          if (fabs(value) < 1.0e5)
     101               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     102          else
     103               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     104          value = rowDual[iRow];
     105          if (fabs(value) < 1.0e5)
     106               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     107          else
     108               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     109          value = rowLower[iRow];
     110          if (fabs(value) < 1.0e5)
     111               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     112          else
     113               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     114          value = rowUpper[iRow];
     115          if (fabs(value) < 1.0e5)
     116               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     117          else
     118               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     119          if (rowObjective) {
     120               value = rowObjective[iRow];
     121               if (fabs(value) < 1.0e5)
     122                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     123               else
     124                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     125          }
     126          std::cout << std::endl;
     127     }
    128128#endif
    129   std::cout<<"--------------------------------------"<<std::endl;
     129     std::cout << "--------------------------------------" << std::endl;
    130130
    131   // Columns
     131     // Columns
    132132
    133   int numberColumns = model.numberColumns();
     133     int numberColumns = model.numberColumns();
    134134
    135   // Alternatively getColSolution()
    136   double * columnPrimal = model.primalColumnSolution();
    137   // Alternatively getReducedCost()
    138   double * columnDual = model.dualColumnSolution();
    139   // Alternatively getColLower()
    140   double * columnLower = model.columnLower();
    141   // Alternatively getColUpper()
    142   double * columnUpper = model.columnUpper();
    143   // Alternatively getObjCoefficients()
    144   double * columnObjective = model.objective();
    145    
    146   // If we have not kept names (parameter to readMps) this will be 0
    147   assert(model.lengthNames());
     135     // Alternatively getColSolution()
     136     double * columnPrimal = model.primalColumnSolution();
     137     // Alternatively getReducedCost()
     138     double * columnDual = model.dualColumnSolution();
     139     // Alternatively getColLower()
     140     double * columnLower = model.columnLower();
     141     // Alternatively getColUpper()
     142     double * columnUpper = model.columnUpper();
     143     // Alternatively getObjCoefficients()
     144     double * columnObjective = model.objective();
    148145
    149   // Column names
    150   const std::vector<std::string> * columnNames = model.columnNames();
     146     // If we have not kept names (parameter to readMps) this will be 0
     147     assert(model.lengthNames());
     148
     149     // Column names
     150     const std::vector<std::string> * columnNames = model.columnNames();
    151151
    152152
    153   int iColumn;
    154  
    155   std::cout<<"                       Primal          Dual         Lower         Upper          Cost"
    156            <<std::endl;
     153     int iColumn;
    157154
    158   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    159     double value;
    160     value = columnPrimal[iColumn];
    161     if (fabs(value)>1.0e-8) {
    162       std::cout<<std::setw(6)<<iColumn<<" "<<std::setw(8)<<(*columnNames)[iColumn];
    163       if (fabs(value)<1.0e5)
    164         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    165       else
    166         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    167       value = columnDual[iColumn];
    168       if (fabs(value)<1.0e5)
    169         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    170       else
    171         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    172       value = columnLower[iColumn];
    173       if (fabs(value)<1.0e5)
    174         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    175       else
    176         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    177       value = columnUpper[iColumn];
    178       if (fabs(value)<1.0e5)
    179         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    180       else
    181         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    182       value = columnObjective[iColumn];
    183       if (fabs(value)<1.0e5)
    184         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    185       else
    186         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    187      
    188       std::cout<<std::endl;
    189     }
    190   }
    191   std::cout<<"--------------------------------------"<<std::endl;
     155     std::cout << "                       Primal          Dual         Lower         Upper          Cost"
     156               << std::endl;
    192157
    193   return 0;
    194 }   
     158     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     159          double value;
     160          value = columnPrimal[iColumn];
     161          if (fabs(value) > 1.0e-8) {
     162               std::cout << std::setw(6) << iColumn << " " << std::setw(8) << (*columnNames)[iColumn];
     163               if (fabs(value) < 1.0e5)
     164                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     165               else
     166                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     167               value = columnDual[iColumn];
     168               if (fabs(value) < 1.0e5)
     169                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     170               else
     171                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     172               value = columnLower[iColumn];
     173               if (fabs(value) < 1.0e5)
     174                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     175               else
     176                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     177               value = columnUpper[iColumn];
     178               if (fabs(value) < 1.0e5)
     179                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     180               else
     181                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     182               value = columnObjective[iColumn];
     183               if (fabs(value) < 1.0e5)
     184                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     185               else
     186                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     187
     188               std::cout << std::endl;
     189          }
     190     }
     191     std::cout << "--------------------------------------" << std::endl;
     192
     193     return 0;
     194}

  • trunk/Clp/examples/driver2.cpp

    r1370 r1552  
    1515
    1616/** This just adds a model to CoinMessage and a void pointer so
    17     user can trap messages and do useful stuff. 
     17    user can trap messages and do useful stuff.
    1818    This is used in Clp/Test/unitTest.cpp
    1919
    … …  
    2424
    2525class MyMessageHandler : public CoinMessageHandler {
    26  
     26
    2727public:
    28   /**@name Overrides */
    29   //@{
    30   virtual int print();
    31   //@}
    32   /**@name set and get */
    33   //@{
    34   /// Model
    35   const ClpSimplex * model() const;
    36   void setModel(ClpSimplex * model);
    37   //@}
    38 
    39   /**@name Constructors, destructor */
    40   //@{
    41   /** Default constructor. */
    42   MyMessageHandler();
    43   /// Constructor with pointer to model
    44   MyMessageHandler(ClpSimplex * model,
    45                            FILE * userPointer=NULL);
    46   /** Destructor */
    47   virtual ~MyMessageHandler();
    48   //@}
    49 
    50   /**@name Copy method */
    51   //@{
    52   /** The copy constructor. */
    53   MyMessageHandler(const MyMessageHandler&);
    54   /** The copy constructor from an CoinSimplexMessageHandler. */
    55   MyMessageHandler(const CoinMessageHandler&);
    56  
    57   MyMessageHandler& operator=(const MyMessageHandler&);
    58   /// Clone
    59   virtual CoinMessageHandler * clone() const ;
    60   //@}
    61    
    62    
     28     /**@name Overrides */
     29     //@{
     30     virtual int print();
     31     //@}
     32     /**@name set and get */
     33     //@{
     34     /// Model
     35     const ClpSimplex * model() const;
     36     void setModel(ClpSimplex * model);
     37     //@}
     38
     39     /**@name Constructors, destructor */
     40     //@{
     41     /** Default constructor. */
     42     MyMessageHandler();
     43     /// Constructor with pointer to model
     44     MyMessageHandler(ClpSimplex * model,
     45                      FILE * userPointer = NULL);
     46     /** Destructor */
     47     virtual ~MyMessageHandler();
     48     //@}
     49
     50     /**@name Copy method */
     51     //@{
     52     /** The copy constructor. */
     53     MyMessageHandler(const MyMessageHandler&);
     54     /** The copy constructor from an CoinSimplexMessageHandler. */
     55     MyMessageHandler(const CoinMessageHandler&);
     56
     57     MyMessageHandler& operator=(const MyMessageHandler&);
     58     /// Clone
     59     virtual CoinMessageHandler * clone() const ;
     60     //@}
     61
     62
    6363protected:
    64   /**@name Data members
    65      The data members are protected to allow access for derived classes. */
    66   //@{
    67   /// Pointer back to model
    68   ClpSimplex * model_;
    69   //@}
     64     /**@name Data members
     65        The data members are protected to allow access for derived classes. */
     66     //@{
     67     /// Pointer back to model
     68     ClpSimplex * model_;
     69     //@}
    7070};
    7171
    … …  
    7676
    7777//-------------------------------------------------------------------
    78 // Default Constructor 
    79 //-------------------------------------------------------------------
    80 MyMessageHandler::MyMessageHandler () 
    81   : CoinMessageHandler(),
    82     model_(NULL)
    83 {
    84 }
    85 
    86 //-------------------------------------------------------------------
    87 // Copy constructor 
    88 //-------------------------------------------------------------------
    89 MyMessageHandler::MyMessageHandler (const MyMessageHandler & rhs) 
    90 : CoinMessageHandler(rhs),
    91     model_(rhs.model_)
    92 { 
    93 }
    94 
    95 MyMessageHandler::MyMessageHandler (const CoinMessageHandler & rhs) 
    96   : CoinMessageHandler(),
    97     model_(NULL)
    98 { 
     78// Default Constructor
     79//-------------------------------------------------------------------
     80MyMessageHandler::MyMessageHandler ()
     81     : CoinMessageHandler(),
     82       model_(NULL)
     83{
     84}
     85
     86//-------------------------------------------------------------------
     87// Copy constructor
     88//-------------------------------------------------------------------
     89MyMessageHandler::MyMessageHandler (const MyMessageHandler & rhs)
     90     : CoinMessageHandler(rhs),
     91       model_(rhs.model_)
     92{
     93}
     94
     95MyMessageHandler::MyMessageHandler (const CoinMessageHandler & rhs)
     96     : CoinMessageHandler(),
     97       model_(NULL)
     98{
    9999}
    100100
    101101// Constructor with pointer to model
    102102MyMessageHandler::MyMessageHandler(ClpSimplex * model,
    103                FILE * userPointer)
    104   : CoinMessageHandler(),
    105     model_(model)
    106 {
    107 }
    108 
    109 //-------------------------------------------------------------------
    110 // Destructor 
     103                                   FILE * userPointer)
     104     : CoinMessageHandler(),
     105       model_(model)
     106{
     107}
     108
     109//-------------------------------------------------------------------
     110// Destructor
    111111//-------------------------------------------------------------------
    112112MyMessageHandler::~MyMessageHandler ()
    … …  
    115115
    116116//----------------------------------------------------------------
    117 // Assignment operator 
     117// Assignment operator
    118118//-------------------------------------------------------------------
    119119MyMessageHandler &
    120120MyMessageHandler::operator=(const MyMessageHandler& rhs)
    121121{
    122   if (this != &rhs) {
    123     CoinMessageHandler::operator=(rhs);
    124     model_ = rhs.model_;
    125   }
    126   return *this;
     122     if (this != &rhs) {
     123          CoinMessageHandler::operator=(rhs);
     124          model_ = rhs.model_;
     125     }
     126     return *this;
    127127}
    128128//-------------------------------------------------------------------
    … …  
    131131CoinMessageHandler * MyMessageHandler::clone() const
    132132{
    133   return new MyMessageHandler(*this);
     133     return new MyMessageHandler(*this);
    134134}
    135135// Print out values from first 20 messages
    136 static int times=0;
    137 int 
     136static int times = 0;
     137int
    138138MyMessageHandler::print()
    139139{
    140   // You could have added a callback flag if you had wanted - see Clp_C_Interface.c
    141   times++;
    142   if (times<=20) {
    143     int messageNumber = currentMessage().externalNumber();
    144     if (currentSource()!="Clp")
    145       messageNumber += 1000000;
    146     int i;
    147     int nDouble=numberDoubleFields();
    148     printf("%d doubles - ",nDouble);
    149     for (i=0;i<nDouble;i++)
    150       printf("%g ",doubleValue(i));
    151     printf("\n");;
    152     int nInt=numberIntFields();
    153     printf("%d ints - ",nInt);
    154     for (i=0;i<nInt;i++)
    155       printf("%d ",intValue(i));
    156     printf("\n");;
    157     int nString=numberStringFields();
    158     printf("%d strings - ",nString);
    159     for (i=0;i<nString;i++)
    160       printf("%s ",stringValue(i).c_str());
    161     printf("\n");;
    162   }
    163   return CoinMessageHandler::print();
     140     // You could have added a callback flag if you had wanted - see Clp_C_Interface.c
     141     times++;
     142     if (times <= 20) {
     143          int messageNumber = currentMessage().externalNumber();
     144          if (currentSource() != "Clp")
     145               messageNumber += 1000000;
     146          int i;
     147          int nDouble = numberDoubleFields();
     148          printf("%d doubles - ", nDouble);
     149          for (i = 0; i < nDouble; i++)
     150               printf("%g ", doubleValue(i));
     151          printf("\n");;
     152          int nInt = numberIntFields();
     153          printf("%d ints - ", nInt);
     154          for (i = 0; i < nInt; i++)
     155               printf("%d ", intValue(i));
     156          printf("\n");;
     157          int nString = numberStringFields();
     158          printf("%d strings - ", nString);
     159          for (i = 0; i < nString; i++)
     160               printf("%s ", stringValue(i).c_str());
     161          printf("\n");;
     162     }
     163     return CoinMessageHandler::print();
    164164}
    165165const ClpSimplex *
    166166MyMessageHandler::model() const
    167167{
    168   return model_;
    169 }
    170 void 
     168     return model_;
     169}
     170void
    171171MyMessageHandler::setModel(ClpSimplex * model)
    172172{
    173   model_ = model;
     173     model_ = model;
    174174}
    175175
    176176int main (int argc, const char *argv[])
    177177{
    178   ClpSimplex  model;
    179   // Message handler
    180   MyMessageHandler messageHandler(&model);
    181   std::cout<<"Testing derived message handler"<<std::endl;
    182   model.passInMessageHandler(&messageHandler);
    183   int status;
    184   // Keep names when reading an mps file
    185   if (argc<2)
    186     status=model.readMps("../../Data/Sample/p0033.mps",true);
    187   else
    188     status=model.readMps(argv[1],true);
    189 
    190   if (status) {
    191     fprintf(stderr,"Bad readMps %s\n",argv[1]);
    192     fprintf(stdout,"Bad readMps %s\n",argv[1]);
    193     exit(1);
    194   }
    195 
    196   double time1 = CoinCpuTime();
    197   /*
    198     This driver shows how to do presolve.by hand (rather than with initialSolve)
    199   */
    200   ClpSimplex * model2;
    201   ClpPresolve pinfo;
    202   int numberPasses=5; // can change this
    203   /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
    204      not being integer, do "numberpasses" passes and throw away names
    205      in presolved model */
    206   model2 = pinfo.presolvedModel(model,1.0e-8,false,numberPasses,false);
    207   if (!model2) {
    208     fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    209             argv[1],1.0e-8);
    210     fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    211             argv[1],1.0e-8);
    212     // model was infeasible - maybe try again with looser tolerances
    213     model2 = pinfo.presolvedModel(model,1.0e-7,false,numberPasses,false);
    214     if (!model2) {
    215       fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    216               argv[1],1.0e-7);
    217       fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    218               argv[1],1.0e-7);
    219       exit(2);
    220     }
    221   }
    222   // change factorization frequency from 200
    223   model2->setFactorizationFrequency(100+model2->numberRows()/50);
    224   if (argc<3 ||!strstr(argv[2],"primal")) {
    225     // Use the dual algorithm unless user said "primal"
    226     /* faster if bounds tightened as then dual can flip variables
    227         to other bound to stay dual feasible.  We can trash the bounds as
    228         this model is going to be thrown away
    229     */
    230     int numberInfeasibilities = model2->tightenPrimalBounds();
    231     if (numberInfeasibilities)
    232       std::cout<<"** Analysis indicates model infeasible"
    233                <<std::endl;
    234     model2->crash(1000.0,2);
    235     ClpDualRowSteepest steep(1);
    236     model2->setDualRowPivotAlgorithm(steep);
    237     model2->dual();
    238   } else {
    239     ClpPrimalColumnSteepest steep(1);
    240     model2->setPrimalColumnPivotAlgorithm(steep);
    241     model2->primal();
    242   }
    243   pinfo.postsolve(true);
    244 
    245   int numberIterations=model2->numberIterations();;
    246   delete model2;
    247   /* After this postsolve model should be optimal.
    248      We can use checkSolution and test feasibility */
    249   model.checkSolution();
    250   if (model.numberDualInfeasibilities()||
    251       model.numberPrimalInfeasibilities())
    252     printf("%g dual %g(%d) Primal %g(%d)\n",
    253            model.objectiveValue(),
    254            model.sumDualInfeasibilities(),
    255            model.numberDualInfeasibilities(),
    256            model.sumPrimalInfeasibilities(),
    257            model.numberPrimalInfeasibilities());
    258   // But resolve for safety
    259   model.primal(1);
    260 
    261   numberIterations += model.numberIterations();;
    262   // for running timing tests
    263   std::cout<<argv[1]<<" Objective "<<model.objectiveValue()<<" took "<<
    264     numberIterations<<" iterations and "<<
    265     CoinCpuTime()-time1<<" seconds"<<std::endl;
    266 
    267   std::string modelName;
    268   model.getStrParam(ClpProbName,modelName);
    269   std::cout<<"Model "<<modelName<<" has "<<model.numberRows()<<" rows and "<<
    270     model.numberColumns()<<" columns"<<std::endl;
    271 
    272   // remove this to print solution
    273 
    274   exit(0);
    275 
    276   /*
    277     Now to print out solution.  The methods used return modifiable
    278     arrays while the alternative names return const pointers -
    279     which is of course much more virtuous.
    280 
    281     This version just does non-zero columns
    282  
    283    */
     178     ClpSimplex  model;
     179     // Message handler
     180     MyMessageHandler messageHandler(&model);
     181     std::cout << "Testing derived message handler" << std::endl;
     182     model.passInMessageHandler(&messageHandler);
     183     int status;
     184     // Keep names when reading an mps file
     185     if (argc < 2)
     186          status = model.readMps("../../Data/Sample/p0033.mps", true);
     187     else
     188          status = model.readMps(argv[1], true);
     189
     190     if (status) {
     191          fprintf(stderr, "Bad readMps %s\n", argv[1]);
     192          fprintf(stdout, "Bad readMps %s\n", argv[1]);
     193          exit(1);
     194     }
     195
     196     double time1 = CoinCpuTime();
     197     /*
     198       This driver shows how to do presolve.by hand (rather than with initialSolve)
     199     */
     200     ClpSimplex * model2;
     201     ClpPresolve pinfo;
     202     int numberPasses = 5; // can change this
     203     /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
     204        not being integer, do "numberpasses" passes and throw away names
     205        in presolved model */
     206     model2 = pinfo.presolvedModel(model, 1.0e-8, false, numberPasses, false);
     207     if (!model2) {
     208          fprintf(stderr, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     209                  argv[1], 1.0e-8);
     210          fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     211                  argv[1], 1.0e-8);
     212          // model was infeasible - maybe try again with looser tolerances
     213          model2 = pinfo.presolvedModel(model, 1.0e-7, false, numberPasses, false);
     214          if (!model2) {
     215               fprintf(stderr, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     216                       argv[1], 1.0e-7);
     217               fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     218                       argv[1], 1.0e-7);
     219               exit(2);
     220          }
     221     }
     222     // change factorization frequency from 200
     223     model2->setFactorizationFrequency(100 + model2->numberRows() / 50);
     224     if (argc < 3 || !strstr(argv[2], "primal")) {
     225          // Use the dual algorithm unless user said "primal"
     226          /* faster if bounds tightened as then dual can flip variables
     227          to other bound to stay dual feasible.  We can trash the bounds as
     228          this model is going to be thrown away
     229          */
     230          int numberInfeasibilities = model2->tightenPrimalBounds();
     231          if (numberInfeasibilities)
     232               std::cout << "** Analysis indicates model infeasible"
     233                         << std::endl;
     234          model2->crash(1000.0, 2);
     235          ClpDualRowSteepest steep(1);
     236          model2->setDualRowPivotAlgorithm(steep);
     237          model2->dual();
     238     } else {
     239          ClpPrimalColumnSteepest steep(1);
     240          model2->setPrimalColumnPivotAlgorithm(steep);
     241          model2->primal();
     242     }
     243     pinfo.postsolve(true);
     244
     245     int numberIterations = model2->numberIterations();;
     246     delete model2;
     247     /* After this postsolve model should be optimal.
     248        We can use checkSolution and test feasibility */
     249     model.checkSolution();
     250     if (model.numberDualInfeasibilities() ||
     251               model.numberPrimalInfeasibilities())
     252          printf("%g dual %g(%d) Primal %g(%d)\n",
     253                 model.objectiveValue(),
     254                 model.sumDualInfeasibilities(),
     255                 model.numberDualInfeasibilities(),
     256                 model.sumPrimalInfeasibilities(),
     257                 model.numberPrimalInfeasibilities());
     258     // But resolve for safety
     259     model.primal(1);
     260
     261     numberIterations += model.numberIterations();;
     262     // for running timing tests
     263     std::cout << argv[1] << " Objective " << model.objectiveValue() << " took " <<
     264               numberIterations << " iterations and " <<
     265               CoinCpuTime() - time1 << " seconds" << std::endl;
     266
     267     std::string modelName;
     268     model.getStrParam(ClpProbName, modelName);
     269     std::cout << "Model " << modelName << " has " << model.numberRows() << " rows and " <<
     270               model.numberColumns() << " columns" << std::endl;
     271
     272     // remove this to print solution
     273
     274     exit(0);
     275
     276     /*
     277       Now to print out solution.  The methods used return modifiable
     278       arrays while the alternative names return const pointers -
     279       which is of course much more virtuous.
     280
     281       This version just does non-zero columns
     282
     283      */
    284284#if 0
    285   int numberRows = model.numberRows();
    286 
    287   // Alternatively getRowActivity()
    288   double * rowPrimal = model.primalRowSolution();
    289   // Alternatively getRowPrice()
    290   double * rowDual = model.dualRowSolution();
    291   // Alternatively getRowLower()
    292   double * rowLower = model.rowLower();
    293   // Alternatively getRowUpper()
    294   double * rowUpper = model.rowUpper();
    295   // Alternatively getRowObjCoefficients()
    296   double * rowObjective = model.rowObjective();
    297    
    298   // If we have not kept names (parameter to readMps) this will be 0
    299   assert(model.lengthNames());
    300 
    301   // Row names
    302   const std::vector<std::string> * rowNames = model.rowNames();
    303 
    304 
    305   int iRow;
    306 
    307   std::cout<<"                       Primal          Dual         Lower         Upper        (Cost)"
    308            <<std::endl;
    309 
    310   for (iRow=0;iRow<numberRows;iRow++) {
    311     double value;
    312     std::cout<<std::setw(6)<<iRow<<" "<<std::setw(8)<<(*rowNames)[iRow];
    313     value = rowPrimal[iRow];
    314     if (fabs(value)<1.0e5)
    315       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    316     else
    317       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    318     value = rowDual[iRow];
    319     if (fabs(value)<1.0e5)
    320       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    321     else
    322       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    323     value = rowLower[iRow];
    324     if (fabs(value)<1.0e5)
    325       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    326     else
    327       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    328     value = rowUpper[iRow];
    329     if (fabs(value)<1.0e5)
    330       std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    331     else
    332       std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    333     if (rowObjective) {
    334       value = rowObjective[iRow];
    335       if (fabs(value)<1.0e5)
    336         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    337       else
    338         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    339     }
    340     std::cout<<std::endl;
    341   }
     285     int numberRows = model.numberRows();
     286
     287     // Alternatively getRowActivity()
     288     double * rowPrimal = model.primalRowSolution();
     289     // Alternatively getRowPrice()
     290     double * rowDual = model.dualRowSolution();
     291     // Alternatively getRowLower()
     292     double * rowLower = model.rowLower();
     293     // Alternatively getRowUpper()
     294     double * rowUpper = model.rowUpper();
     295     // Alternatively getRowObjCoefficients()
     296     double * rowObjective = model.rowObjective();
     297
     298     // If we have not kept names (parameter to readMps) this will be 0
     299     assert(model.lengthNames());
     300
     301     // Row names
     302     const std::vector<std::string> * rowNames = model.rowNames();
     303
     304
     305     int iRow;
     306
     307     std::cout << "                       Primal          Dual         Lower         Upper        (Cost)"
     308               << std::endl;
     309
     310     for (iRow = 0; iRow < numberRows; iRow++) {
     311          double value;
     312          std::cout << std::setw(6) << iRow << " " << std::setw(8) << (*rowNames)[iRow];
     313          value = rowPrimal[iRow];
     314          if (fabs(value) < 1.0e5)
     315               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     316          else
     317               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     318          value = rowDual[iRow];
     319          if (fabs(value) < 1.0e5)
     320               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     321          else
     322               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     323          value = rowLower[iRow];
     324          if (fabs(value) < 1.0e5)
     325               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     326          else
     327               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     328          value = rowUpper[iRow];
     329          if (fabs(value) < 1.0e5)
     330               std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     331          else
     332               std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     333          if (rowObjective) {
     334               value = rowObjective[iRow];
     335               if (fabs(value) < 1.0e5)
     336                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     337               else
     338                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     339          }
     340          std::cout << std::endl;
     341     }
    342342#endif
    343   std::cout<<"--------------------------------------"<<std::endl;
    344 
    345   // Columns
    346 
    347   int numberColumns = model.numberColumns();
    348 
    349   // Alternatively getColSolution()
    350   double * columnPrimal = model.primalColumnSolution();
    351   // Alternatively getReducedCost()
    352   double * columnDual = model.dualColumnSolution();
    353   // Alternatively getColLower()
    354   double * columnLower = model.columnLower();
    355   // Alternatively getColUpper()
    356   double * columnUpper = model.columnUpper();
    357   // Alternatively getObjCoefficients()
    358   double * columnObjective = model.objective();
    359    
    360   // If we have not kept names (parameter to readMps) this will be 0
    361   assert(model.lengthNames());
    362 
    363   // Column names
    364   const std::vector<std::string> * columnNames = model.columnNames();
    365 
    366 
    367   int iColumn;
    368  
    369   std::cout<<"                       Primal          Dual         Lower         Upper          Cost"
    370            <<std::endl;
    371 
    372   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    373     double value;
    374     value = columnPrimal[iColumn];
    375     if (fabs(value)>1.0e-8) {
    376       std::cout<<std::setw(6)<<iColumn<<" "<<std::setw(8)<<(*columnNames)[iColumn];
    377       if (fabs(value)<1.0e5)
    378         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    379       else
    380         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    381       value = columnDual[iColumn];
    382       if (fabs(value)<1.0e5)
    383         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    384       else
    385         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    386       value = columnLower[iColumn];
    387       if (fabs(value)<1.0e5)
    388         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    389       else
    390         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    391       value = columnUpper[iColumn];
    392       if (fabs(value)<1.0e5)
    393         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    394       else
    395         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    396       value = columnObjective[iColumn];
    397       if (fabs(value)<1.0e5)
    398         std::cout<<setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14)<<value;
    399       else
    400         std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;
    401      
    402       std::cout<<std::endl;
    403     }
    404   }
    405   std::cout<<"--------------------------------------"<<std::endl;
    406 
    407   return 0;
    408 }   
     343     std::cout << "--------------------------------------" << std::endl;
     344
     345     // Columns
     346
     347     int numberColumns = model.numberColumns();
     348
     349     // Alternatively getColSolution()
     350     double * columnPrimal = model.primalColumnSolution();
     351     // Alternatively getReducedCost()
     352     double * columnDual = model.dualColumnSolution();
     353     // Alternatively getColLower()
     354     double * columnLower = model.columnLower();
     355     // Alternatively getColUpper()
     356     double * columnUpper = model.columnUpper();
     357     // Alternatively getObjCoefficients()
     358     double * columnObjective = model.objective();
     359
     360     // If we have not kept names (parameter to readMps) this will be 0
     361     assert(model.lengthNames());
     362
     363     // Column names
     364     const std::vector<std::string> * columnNames = model.columnNames();
     365
     366
     367     int iColumn;
     368
     369     std::cout << "                       Primal          Dual         Lower         Upper          Cost"
     370               << std::endl;
     371
     372     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     373          double value;
     374          value = columnPrimal[iColumn];
     375          if (fabs(value) > 1.0e-8) {
     376               std::cout << std::setw(6) << iColumn << " " << std::setw(8) << (*columnNames)[iColumn];
     377               if (fabs(value) < 1.0e5)
     378                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     379               else
     380                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     381               value = columnDual[iColumn];
     382               if (fabs(value) < 1.0e5)
     383                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     384               else
     385                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     386               value = columnLower[iColumn];
     387               if (fabs(value) < 1.0e5)
     388                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     389               else
     390                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     391               value = columnUpper[iColumn];
     392               if (fabs(value) < 1.0e5)
     393                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     394               else
     395                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     396               value = columnObjective[iColumn];
     397               if (fabs(value) < 1.0e5)
     398                    std::cout << setiosflags(std::ios::fixed | std::ios::showpoint) << std::setw(14) << value;
     399               else
     400                    std::cout << setiosflags(std::ios::scientific) << std::setw(14) << value;
     401
     402               std::cout << std::endl;
     403          }
     404     }
     405     std::cout << "--------------------------------------" << std::endl;
     406
     407     return 0;
     408}

  • trunk/Clp/examples/driverC.c

    r1370 r1552  
    1313/* Call back function - just says whenever it gets Clp0005 or Coin0002 */
    1414static void callBack(Clp_Simplex * model, int messageNumber,
    15                      int nDouble, const double * vDouble,
    16                      int nInt, const int * vInt,
    17                      int nString, char ** vString)
     15                     int nDouble, const double * vDouble,
     16                     int nInt, const int * vInt,
     17                     int nString, char ** vString)
    1818{
    19   if (messageNumber==1000002) {
    20     /* Coin0002 */
    21     assert (nString==1&&nInt==3);
    22     printf("Name of problem is %s\n",vString[0]);
    23     printf("row %d col %d el %d\n",vInt[0],vInt[1],vInt[2]);
    24   } else if (messageNumber==5) {
    25     /* Clp0005 */
    26     int i;
    27     assert (nInt==4&&nDouble==3); /* they may not all print */
    28     for (i=0;i<3;i++)
    29       printf("%d %g\n",vInt[i],vDouble[i]);
    30   }
     19     if (messageNumber == 1000002) {
     20          /* Coin0002 */
     21          assert (nString == 1 && nInt == 3);
     22          printf("Name of problem is %s\n", vString[0]);
     23          printf("row %d col %d el %d\n", vInt[0], vInt[1], vInt[2]);
     24     } else if (messageNumber == 5) {
     25          /* Clp0005 */
     26          int i;
     27          assert (nInt == 4 && nDouble == 3); /* they may not all print */
     28          for (i = 0; i < 3; i++)
     29               printf("%d %g\n", vInt[i], vDouble[i]);
     30     }
    3131}
    32                      
     32
    3333
    3434
    3535int main (int argc, const char *argv[])
    3636{
    37   /* Get default model */
    38   Clp_Simplex  * model = Clp_newModel();
    39   int status;
    40   /* register callback */
    41   Clp_registerCallBack(model,callBack);
    42   /* Keep names when reading an mps file */
    43   if (argc<2)
    44     status=Clp_readMps(model,"../../Data/Sample/p0033.mps",1,0);
    45   else
    46     status=Clp_readMps(model,argv[1],1,0);
     37     /* Get default model */
     38     Clp_Simplex  * model = Clp_newModel();
     39     int status;
     40     /* register callback */
     41     Clp_registerCallBack(model, callBack);
     42     /* Keep names when reading an mps file */
     43     if (argc < 2)
     44          status = Clp_readMps(model, "../../Data/Sample/p0033.mps", 1, 0);
     45     else
     46          status = Clp_readMps(model, argv[1], 1, 0);
    4747
    48   if (status) {
    49     fprintf(stderr,"Bad readMps %s\n",argv[1]);
    50     fprintf(stdout,"Bad readMps %s\n",argv[1]);
    51     exit(1);
    52   }
     48     if (status) {
     49          fprintf(stderr, "Bad readMps %s\n", argv[1]);
     50          fprintf(stdout, "Bad readMps %s\n", argv[1]);
     51          exit(1);
     52     }
    5353
    54   if (argc<3 ||!strstr(argv[2],"primal")) {
    55     /* Use the dual algorithm unless user said "primal" */
    56     Clp_initialDualSolve(model);
    57   } else {
    58     Clp_initialPrimalSolve(model);
    59   }
     54     if (argc < 3 || !strstr(argv[2], "primal")) {
     55          /* Use the dual algorithm unless user said "primal" */
     56          Clp_initialDualSolve(model);
     57     } else {
     58          Clp_initialPrimalSolve(model);
     59     }
    6060
    61   {
    62     char modelName[80];
    63     Clp_problemName(model,80,modelName);
    64     printf("Model %s has %d rows and %d columns\n",
    65            modelName,Clp_numberRows(model),Clp_numberColumns(model));
    66   }
     61     {
     62          char modelName[80];
     63          Clp_problemName(model, 80, modelName);
     64          printf("Model %s has %d rows and %d columns\n",
     65                 modelName, Clp_numberRows(model), Clp_numberColumns(model));
     66     }
    6767
    68   /* remove this to print solution */
     68     /* remove this to print solution */
    6969
    70   /*exit(0); */
     70     /*exit(0); */
    7171
    72   {
    73     /*
    74       Now to print out solution.  The methods used return modifiable
    75       arrays while the alternative names return const pointers -
    76       which is of course much more virtuous.
    77      
    78       This version just does non-zero columns
    79      
    80     */
    81    
    82     /* Columns */
    83    
    84     int numberColumns = Clp_numberColumns(model);
    85     int iColumn;
    86    
    87    
    88     /* Alternatively getColSolution(model) */
    89     double * columnPrimal = Clp_primalColumnSolution(model);
    90     /* Alternatively getReducedCost(model) */
    91     double * columnDual = Clp_dualColumnSolution(model);
    92     /* Alternatively getColLower(model) */
    93     double * columnLower = Clp_columnLower(model);
    94     /* Alternatively getColUpper(model) */
    95     double * columnUpper = Clp_columnUpper(model);
    96     /* Alternatively getObjCoefficients(model) */
    97     double * columnObjective = Clp_objective(model);
    98    
    99     printf("--------------------------------------\n");
    100     /* If we have not kept names (parameter to readMps) this will be 0 */
    101     assert(Clp_lengthNames(model));
    102    
    103     printf("                       Primal          Dual         Lower         Upper          Cost\n");
    104    
    105     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    106       double value;
    107       value = columnPrimal[iColumn];
    108       if (value>1.0e-8||value<-1.0e-8) {
    109         char name[20];
    110         Clp_columnName(model,iColumn,name);
    111         printf("%6d %8s",iColumn,name);
    112         printf(" %13g",columnPrimal[iColumn]);
    113         printf(" %13g",columnDual[iColumn]);
    114         printf(" %13g",columnLower[iColumn]);
    115         printf(" %13g",columnUpper[iColumn]);
    116         printf(" %13g",columnObjective[iColumn]);
    117         printf("\n");
    118       }
    119     }
    120     printf("--------------------------------------\n");
    121   }
    122   return 0;
    123 }   
     72     {
     73          /*
     74            Now to print out solution.  The methods used return modifiable
     75            arrays while the alternative names return const pointers -
     76            which is of course much more virtuous.
     77
     78            This version just does non-zero columns
     79
     80          */
     81
     82          /* Columns */
     83
     84          int numberColumns = Clp_numberColumns(model);
     85          int iColumn;
     86
     87
     88          /* Alternatively getColSolution(model) */
     89          double * columnPrimal = Clp_primalColumnSolution(model);
     90          /* Alternatively getReducedCost(model) */
     91          double * columnDual = Clp_dualColumnSolution(model);
     92          /* Alternatively getColLower(model) */
     93          double * columnLower = Clp_columnLower(model);
     94          /* Alternatively getColUpper(model) */
     95          double * columnUpper = Clp_columnUpper(model);
     96          /* Alternatively getObjCoefficients(model) */
     97          double * columnObjective = Clp_objective(model);
     98
     99          printf("--------------------------------------\n");
     100          /* If we have not kept names (parameter to readMps) this will be 0 */
     101          assert(Clp_lengthNames(model));
     102
     103          printf("                       Primal          Dual         Lower         Upper          Cost\n");
     104
     105          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     106               double value;
     107               value = columnPrimal[iColumn];
     108               if (value > 1.0e-8 || value < -1.0e-8) {
     109                    char name[20];
     110                    Clp_columnName(model, iColumn, name);
     111                    printf("%6d %8s", iColumn, name);
     112                    printf(" %13g", columnPrimal[iColumn]);
     113                    printf(" %13g", columnDual[iColumn]);
     114                    printf(" %13g", columnLower[iColumn]);
     115                    printf(" %13g", columnUpper[iColumn]);
     116                    printf(" %13g", columnObjective[iColumn]);
     117                    printf("\n");
     118               }
     119          }
     120          printf("--------------------------------------\n");
     121     }
     122     return 0;
     123}

  • trunk/Clp/examples/dualCuts.cpp

    r1370 r1552  
    11/* $Id$ */
    2 /* Copyright (C) 2003, International Business Machines Corporation 
     2/* Copyright (C) 2003, International Business Machines Corporation
    33   and others.  All Rights Reserved.
    44
    5    This sample program is designed to illustrate programming 
     5   This sample program is designed to illustrate programming
    66   techniques using CoinLP, has not been thoroughly tested
    77   and comes without any warranty whatsoever.
    88
    9    You may copy, modify and distribute this sample program without 
     9   You may copy, modify and distribute this sample program without
    1010   any restrictions whatsoever and without any payment to anyone.
    1111*/
    … …  
    2222int main (int argc, const char *argv[])
    2323{
    24   ClpSimplex  model;
    25   int status;
    26   // Keep names
    27   if (argc<2) {
    28     status=model.readMps("small.mps",true);
    29   } else {
    30     status=model.readMps(argv[1],false);
    31   }
    32   if (status)
    33     exit(10);
    34   /*
    35     This driver implements a method of treating a problem as all cuts.
    36     So it adds in all E rows, solves and then adds in violated rows.
    37   */
    38 
    39   double time1 = CoinCpuTime();
    40   ClpSimplex * model2;
    41   ClpPresolve pinfo;
    42   int numberPasses=5; // can change this
    43   /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
    44      not being integer, do "numberpasses" passes and throw away names
    45      in presolved model */
    46   model2 = pinfo.presolvedModel(model,1.0e-8,false,numberPasses,false);
    47   if (!model2) {
    48     fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    49             argv[1],1.0e-8);
    50     fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    51             argv[1],1.0e-8);
    52     // model was infeasible - maybe try again with looser tolerances
    53     model2 = pinfo.presolvedModel(model,1.0e-7,false,numberPasses,false);
    54     if (!model2) {
    55       fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    56               argv[1],1.0e-7);
    57       fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
    58               argv[1],1.0e-7);
    59       exit(2);
    60     }
    61   }
    62   // change factorization frequency from 200
    63   model2->setFactorizationFrequency(100+model2->numberRows()/50);
    64  
    65   int numberColumns = model2->numberColumns();
    66   int originalNumberRows = model2->numberRows();
    67  
    68   // We will need arrays to choose rows to add
    69   double * weight = new double [originalNumberRows];
    70   int * sort = new int [originalNumberRows];
    71   int numberSort=0;
    72   char * take = new char [originalNumberRows];
    73  
    74   const double * rowLower = model2->rowLower();
    75   const double * rowUpper = model2->rowUpper();
    76   int iRow,iColumn;
    77   // Set up initial list
    78   numberSort=0;
    79   for (iRow=0;iRow<originalNumberRows;iRow++) {
    80     weight[iRow]=1.123e50;
    81     if (rowLower[iRow]==rowUpper[iRow]) {
    82       sort[numberSort++] = iRow;
    83       weight[iRow]=0.0;
    84     }
    85   }
    86   numberSort /= 2;
    87   // Just add this number of rows each time in small problem
    88   int smallNumberRows = 2*numberColumns;
    89   smallNumberRows=CoinMin(smallNumberRows,originalNumberRows/20);
    90   // and pad out with random rows
    91   double ratio = ((double)(smallNumberRows-numberSort))/((double) originalNumberRows);
    92   for (iRow=0;iRow<originalNumberRows;iRow++) {
    93     if (weight[iRow]==1.123e50&&CoinDrand48()<ratio)
    94       sort[numberSort++] = iRow;
    95   }
    96   /* This is optional.
    97      The best thing to do is to miss out random rows and do a set which makes dual feasible.
    98      If that is not possible then make sure variables have bounds.
    99 
    100      One way that normally works is to automatically tighten bounds.
    101   */
     24     ClpSimplex  model;
     25     int status;
     26     // Keep names
     27     if (argc < 2) {
     28          status = model.readMps("small.mps", true);
     29     } else {
     30          status = model.readMps(argv[1], false);
     31     }
     32     if (status)
     33          exit(10);
     34     /*
     35       This driver implements a method of treating a problem as all cuts.
     36       So it adds in all E rows, solves and then adds in violated rows.
     37     */
     38
     39     double time1 = CoinCpuTime();
     40     ClpSimplex * model2;
     41     ClpPresolve pinfo;
     42     int numberPasses = 5; // can change this
     43     /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
     44        not being integer, do "numberpasses" passes and throw away names
     45        in presolved model */
     46     model2 = pinfo.presolvedModel(model, 1.0e-8, false, numberPasses, false);
     47     if (!model2) {
     48          fprintf(stderr, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     49                  argv[1], 1.0e-8);
     50          fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     51                  argv[1], 1.0e-8);
     52          // model was infeasible - maybe try again with looser tolerances
     53          model2 = pinfo.presolvedModel(model, 1.0e-7, false, numberPasses, false);
     54          if (!model2) {
     55               fprintf(stderr, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     56                       argv[1], 1.0e-7);
     57               fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
     58                       argv[1], 1.0e-7);
     59               exit(2);
     60          }
     61     }
     62     // change factorization frequency from 200
     63     model2->setFactorizationFrequency(100 + model2->numberRows() / 50);
     64
     65     int numberColumns = model2->numberColumns();
     66     int originalNumberRows = model2->numberRows();
     67
     68     // We will need arrays to choose rows to add
     69     double * weight = new double [originalNumberRows];
     70     int * sort = new int [originalNumberRows];
     71     int numberSort = 0;
     72     char * take = new char [originalNumberRows];
     73
     74     const double * rowLower = model2->rowLower();
     75     const double * rowUpper = model2->rowUpper();
     76     int iRow, iColumn;
     77     // Set up initial list
     78     numberSort = 0;
     79     for (iRow = 0; iRow < originalNumberRows; iRow++) {
     80          weight[iRow] = 1.123e50;
     81          if (rowLower[iRow] == rowUpper[iRow]) {
     82               sort[numberSort++] = iRow;
     83               weight[iRow] = 0.0;
     84          }
     85     }
     86     numberSort /= 2;
     87     // Just add this number of rows each time in small problem
     88     int smallNumberRows = 2 * numberColumns;
     89     smallNumberRows = CoinMin(smallNumberRows, originalNumberRows / 20);
     90     // and pad out with random rows
     91     double ratio = ((double)(smallNumberRows - numberSort)) / ((double) originalNumberRows);
     92     for (iRow = 0; iRow < originalNumberRows; iRow++) {
     93          if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio)
     94               sort[numberSort++] = iRow;
     95     }
     96     /* This is optional.
     97        The best thing to do is to miss out random rows and do a set which makes dual feasible.
     98        If that is not possible then make sure variables have bounds.
     99
     100        One way that normally works is to automatically tighten bounds.
     101     */
    102102#if 0
    103   // However for some we need to do anyway
    104   double * columnLower = model2->columnLower();
    105   double * columnUpper = model2->columnUpper();
    106   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    107     columnLower[iColumn]=CoinMax(-1.0e6,columnLower[iColumn]);
    108     columnUpper[iColumn]=CoinMin(1.0e6,columnUpper[iColumn]);
    109   }
     103     // However for some we need to do anyway
     104     double * columnLower = model2->columnLower();
     105     double * columnUpper = model2->columnUpper();
     106     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     107          columnLower[iColumn] = CoinMax(-1.0e6, columnLower[iColumn]);
     108          columnUpper[iColumn] = CoinMin(1.0e6, columnUpper[iColumn]);
     109     }
    110110#endif
    111   model2->tightenPrimalBounds(-1.0e4,true);
    112   printf("%d rows in initial problem\n",numberSort);
    113   double * fullSolution = model2->primalRowSolution();
    114  
    115   // Just do this number of passes
    116   int maxPass=50;
    117   // And take out slack rows until this pass
    118   int takeOutPass=30;
    119   int iPass;
    120  
    121   const int * start = model2->clpMatrix()->getVectorStarts();
    122   const int * length = model2->clpMatrix()->getVectorLengths();
    123   const int * row = model2->clpMatrix()->getIndices();
    124   int * whichColumns = new int [numberColumns];
    125   for (iColumn=0;iColumn<numberColumns;iColumn++)
    126     whichColumns[iColumn]=iColumn;
    127   int numberSmallColumns=numberColumns;
    128   for (iPass=0;iPass<maxPass;iPass++) {
    129     printf("Start of pass %d\n",iPass);
    130     // Cleaner this way
    131     std::sort(sort,sort+numberSort);
    132     // Create small problem
    133     ClpSimplex small(model2,numberSort,sort,numberSmallColumns,whichColumns);
    134     small.setFactorizationFrequency(100+numberSort/200);
    135     //small.setPerturbation(50);
    136     //small.setLogLevel(63);
    137     // A variation is to just do N iterations
    138     //if (iPass)
    139     //small.setMaximumIterations(100);
    140     // Solve
    141     small.factorization()->messageLevel(8);
    142     if (iPass) {
    143       small.dual();
    144     } else {
    145       small.writeMps("continf.mps");
    146       ClpSolve solveOptions;
    147       solveOptions.setSolveType(ClpSolve::useDual);
    148       //solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    149       //solveOptions.setSpecialOption(1,2,200); // idiot
    150       small.initialSolve(solveOptions);
    151     }
    152     bool dualInfeasible = (small.status()==2);
    153     // move solution back
    154     double * solution = model2->primalColumnSolution();
    155     const double * smallSolution = small.primalColumnSolution();
    156     for (int j=0;j<numberSmallColumns;j++) {
    157       iColumn = whichColumns[j];
    158       solution[iColumn]=smallSolution[j];
    159       model2->setColumnStatus(iColumn,small.getColumnStatus(j));
    160     }
    161     for (iRow=0;iRow<numberSort;iRow++) {
    162       int kRow = sort[iRow];
    163       model2->setRowStatus(kRow,small.getRowStatus(iRow));
    164     }
    165     // compute full solution
    166     memset(fullSolution,0,originalNumberRows*sizeof(double));
    167     model2->times(1.0,model2->primalColumnSolution(),fullSolution);
    168     if (iPass!=maxPass-1) {
    169       // Mark row as not looked at
    170       for (iRow=0;iRow<originalNumberRows;iRow++)
    171         weight[iRow]=1.123e50;
    172       // Look at rows already in small problem
    173       int iSort;
    174       int numberDropped=0;
    175       int numberKept=0;
    176       int numberBinding=0;
    177       int numberInfeasibilities=0;
    178       double sumInfeasibilities=0.0;
    179       for (iSort=0;iSort<numberSort;iSort++) {
    180         iRow=sort[iSort];
    181         //printf("%d %g %g\n",iRow,fullSolution[iRow],small.primalRowSolution()[iSort]);
    182         if (model2->getRowStatus(iRow)==ClpSimplex::basic) {
    183           // Basic - we can get rid of if early on
    184           if (iPass<takeOutPass&&!dualInfeasible) {
    185             // may have hit max iterations so check
    186             double infeasibility = CoinMax(fullSolution[iRow]-rowUpper[iRow],
    187                                        rowLower[iRow]-fullSolution[iRow]);
    188             weight[iRow]=-infeasibility;
    189             if (infeasibility>1.0e-8) {
    190               numberInfeasibilities++;
    191               sumInfeasibilities += infeasibility;
    192             } else {
    193               weight[iRow]=1.0;
    194               numberDropped++;
    195             }
    196           } else {
    197             // keep
    198             weight[iRow]=-1.0e40;
    199             numberKept++;
    200           }
    201         } else {
    202           // keep
    203           weight[iRow]=-1.0e50;
    204           numberKept++;
    205           numberBinding++;
    206         }
    207       }
    208       // Now rest
    209       for (iRow=0;iRow<originalNumberRows;iRow++) {
    210         sort[iRow]=iRow;
    211         if (weight[iRow]==1.123e50) {
    212           // not looked at yet
    213           double infeasibility = CoinMax(fullSolution[iRow]-rowUpper[iRow],
    214                                      rowLower[iRow]-fullSolution[iRow]);
    215           weight[iRow]=-infeasibility;
    216           if (infeasibility>1.0e-8) {
    217             numberInfeasibilities++;
    218             sumInfeasibilities += infeasibility;
    219           }
    220         }
    221       }
    222       // sort
    223       CoinSort_2(weight,weight+originalNumberRows,sort);
    224       numberSort = CoinMin(originalNumberRows,smallNumberRows+numberKept);
    225       memset(take,0,originalNumberRows);
    226       for (iRow=0;iRow<numberSort;iRow++)
    227       take[sort[iRow]]=1;
    228       numberSmallColumns=0;
    229       for (iColumn=0;iColumn<numberColumns;iColumn++) {
    230         int n=0;
    231         for (int j=start[iColumn];j<start[iColumn]+length[iColumn];j++) {
    232           int iRow = row[j];
    233           if (take[iRow])
    234             n++;
    235         }
    236         if (n)
    237           whichColumns[numberSmallColumns++]=iColumn;
    238       }
    239       printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n",
    240              numberBinding,numberKept,numberDropped,numberSort,numberSmallColumns);
    241       printf("%d rows are infeasible - sum is %g\n",numberInfeasibilities,
    242                sumInfeasibilities);
    243       if (!numberInfeasibilities) {
    244         printf("Exiting as looks optimal\n");
    245         break;
    246       }
    247       numberInfeasibilities=0;
    248       sumInfeasibilities=0.0;
    249       for (iSort=0;iSort<numberSort;iSort++) {
    250         if (weight[iSort]>-1.0e30&&weight[iSort]<-1.0e-8) {
    251           numberInfeasibilities++;
    252           sumInfeasibilities += -weight[iSort];
    253         }
    254       }
    255       printf("in small model %d rows are infeasible - sum is %g\n",numberInfeasibilities,
    256                sumInfeasibilities);
    257     }
    258   }
    259   delete [] weight;
    260   delete [] sort;
    261   delete [] whichColumns;
    262   delete [] take;
    263   // If problem is big you may wish to skip this
    264   model2->dual();
    265   int numberBinding=0;
    266   for (iRow=0;iRow<originalNumberRows;iRow++) {
    267     if (model2->getRowStatus(iRow)!=ClpSimplex::basic)
    268       numberBinding++;
    269   }
    270   printf("%d binding rows at end\n",numberBinding);
    271   pinfo.postsolve(true);
    272 
    273   int numberIterations=model2->numberIterations();;
    274   delete model2;
    275   /* After this postsolve model should be optimal.
    276      We can use checkSolution and test feasibility */
    277   model.checkSolution();
    278   if (model.numberDualInfeasibilities()||
    279       model.numberPrimalInfeasibilities())
    280     printf("%g dual %g(%d) Primal %g(%d)\n",
    281            model.objectiveValue(),
    282            model.sumDualInfeasibilities(),
    283            model.numberDualInfeasibilities(),
    284            model.sumPrimalInfeasibilities(),
    285            model.numberPrimalInfeasibilities());
    286   // But resolve for safety
    287   model.primal(1);
    288 
    289   numberIterations += model.numberIterations();;
    290   printf("Solve took %g seconds\n",CoinCpuTime()-time1);
    291   return 0;
    292 }   
     111     model2->tightenPrimalBounds(-1.0e4, true);
     112     printf("%d rows in initial problem\n", numberSort);
     113     double * fullSolution = model2->primalRowSolution();
     114
     115     // Just do this number of passes
     116     int maxPass = 50;
     117     // And take out slack rows until this pass
     118     int takeOutPass = 30;
     119     int iPass;
     120
     121     const int * start = model2->clpMatrix()->getVectorStarts();
     122     const int * length = model2->clpMatrix()->getVectorLengths();
     123     const int * row = model2->clpMatrix()->getIndices();
     124     int * whichColumns = new int [numberColumns];
     125     for (iColumn = 0; iColumn < numberColumns; iColumn++)
     126          whichColumns[iColumn] = iColumn;
     127     int numberSmallColumns = numberColumns;
     128     for (iPass = 0; iPass < maxPass; iPass++) {
     129          printf("Start of pass %d\n", iPass);
     130          // Cleaner this way
     131          std::sort(sort, sort + numberSort);
     132          // Create small problem
     133          ClpSimplex small(model2, numberSort, sort, numberSmallColumns, whichColumns);
     134          small.setFactorizationFrequency(100 + numberSort / 200);
     135          //small.setPerturbation(50);
     136          //small.setLogLevel(63);
     137          // A variation is to just do N iterations
     138          //if (iPass)
     139          //small.setMaximumIterations(100);
     140          // Solve
     141          small.factorization()->messageLevel(8);
     142          if (iPass) {
     143               small.dual();
     144          } else {
     145               small.writeMps("continf.mps");
     146               ClpSolve solveOptions;
     147               solveOptions.setSolveType(ClpSolve::useDual);
     148               //solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
     149               //solveOptions.setSpecialOption(1,2,200); // idiot
     150               small.initialSolve(solveOptions);
     151          }
     152          bool dualInfeasible = (small.status() == 2);
     153          // move solution back
     154          double * solution = model2->primalColumnSolution();
     155          const double * smallSolution = small.primalColumnSolution();
     156          for (int j = 0; j < numberSmallColumns; j++) {
     157               iColumn = whichColumns[j];
     158               solution[iColumn] = smallSolution[j];
     159               model2->setColumnStatus(iColumn, small.getColumnStatus(j));
     160          }
     161          for (iRow = 0; iRow < numberSort; iRow++) {
     162               int kRow = sort[iRow];
     163               model2->setRowStatus(kRow, small.getRowStatus(iRow));
     164          }
     165          // compute full solution
     166          memset(fullSolution, 0, originalNumberRows * sizeof(double));
     167          model2->times(1.0, model2->primalColumnSolution(), fullSolution);
     168          if (iPass != maxPass - 1) {
     169               // Mark row as not looked at
     170               for (iRow = 0; iRow < originalNumberRows; iRow++)
     171                    weight[iRow] = 1.123e50;
     172               // Look at rows already in small problem
     173               int iSort;
     174               int numberDropped = 0;
     175               int numberKept = 0;
     176               int numberBinding = 0;
     177               int numberInfeasibilities = 0;
     178               double sumInfeasibilities = 0.0;
     179               for (iSort = 0; iSort < numberSort; iSort++) {
     180                    iRow = sort[iSort];
     181                    //printf("%d %g %g\n",iRow,fullSolution[iRow],small.primalRowSolution()[iSort]);
     182                    if (model2->getRowStatus(iRow) == ClpSimplex::basic) {
     183                         // Basic - we can get rid of if early on
     184                         if (iPass < takeOutPass && !dualInfeasible) {
     185                              // may have hit max iterations so check
     186                              double infeasibility = CoinMax(fullSolution[iRow] - rowUpper[iRow],
     187                                                             rowLower[iRow] - fullSolution[iRow]);
     188                              weight[iRow] = -infeasibility;
     189                              if (infeasibility > 1.0e-8) {
     190                                   numberInfeasibilities++;
     191                                   sumInfeasibilities += infeasibility;
     192                              } else {
     193                                   weight[iRow] = 1.0;
     194                                   numberDropped++;
     195                              }
     196                         } else {
     197                              // keep
     198                              weight[iRow] = -1.0e40;
     199                              numberKept++;
     200                         }
     201                    } else {
     202                         // keep
     203                         weight[iRow] = -1.0e50;
     204                         numberKept++;
     205                         numberBinding++;
     206                    }
     207               }
     208               // Now rest
     209               for (iRow = 0; iRow < originalNumberRows; iRow++) {
     210                    sort[iRow] = iRow;
     211                    if (weight[iRow] == 1.123e50) {
     212                         // not looked at yet
     213                         double infeasibility = CoinMax(fullSolution[iRow] - rowUpper[iRow],
     214                                                        rowLower[iRow] - fullSolution[iRow]);
     215                         weight[iRow] = -infeasibility;
     216                         if (infeasibility > 1.0e-8) {
     217                              numberInfeasibilities++;
     218                              sumInfeasibilities += infeasibility;
     219                         }
     220                    }
     221               }
     222               // sort
     223               CoinSort_2(weight, weight + originalNumberRows, sort);
     224               numberSort = CoinMin(originalNumberRows, smallNumberRows + numberKept);
     225               memset(take, 0, originalNumberRows);
     226               for (iRow = 0; iRow < numberSort; iRow++)
     227                    take[sort[iRow]] = 1;
     228               numberSmallColumns = 0;
     229               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     230                    int n = 0;
     231                    for (int j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
     232                         int iRow = row[j];
     233                         if (take[iRow])
     234                              n++;
     235                    }
     236                    if (n)
     237                         whichColumns[numberSmallColumns++] = iColumn;
     238               }
     239               printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n",
     240                      numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns);
     241               printf("%d rows are infeasible - sum is %g\n", numberInfeasibilities,
     242                      sumInfeasibilities);
     243               if (!numberInfeasibilities) {
     244                    printf("Exiting as looks optimal\n");
     245                    break;
     246               }
     247               numberInfeasibilities = 0;
     248               sumInfeasibilities = 0.0;
     249               for (iSort = 0; iSort < numberSort; iSort++) {
     250                    if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) {
     251                         numberInfeasibilities++;
     252                         sumInfeasibilities += -weight[iSort];
     253                    }
     254               }
     255               printf("in small model %d rows are infeasible - sum is %g\n", numberInfeasibilities,
     256                      sumInfeasibilities);
     257          }
     258     }
     259     delete [] weight;
     260     delete [] sort;
     261     delete [] whichColumns;
     262     delete [] take;
     263     // If problem is big you may wish to skip this
     264     model2->dual();
     265     int numberBinding = 0;
     266     for (iRow = 0; iRow < originalNumberRows; iRow++) {
     267          if (model2->getRowStatus(iRow) != ClpSimplex::basic)
     268               numberBinding++;
     269     }
     270     printf("%d binding rows at end\n", numberBinding);
     271     pinfo.postsolve(true);
     272
     273     int numberIterations = model2->numberIterations();;
     274     delete model2;
     275     /* After this postsolve model should be optimal.
     276        We can use checkSolution and test feasibility */
     277     model.checkSolution();
     278     if (model.numberDualInfeasibilities() ||
     279               model.numberPrimalInfeasibilities())
     280          printf("%g dual %g(%d) Primal %g(%d)\n",
     281                 model.objectiveValue(),
     282                 model.sumDualInfeasibilities(),
     283                 model.numberDualInfeasibilities(),
     284                 model.sumPrimalInfeasibilities(),
     285                 model.numberPrimalInfeasibilities());
     286     // But resolve for safety
     287     model.primal(1);
     288
     289     numberIterations += model.numberIterations();;
     290     printf("Solve took %g seconds\n", CoinCpuTime() - time1);
     291     return 0;
     292}

  • trunk/Clp/examples/ekk.cpp

    r1370 r1552  
    4343*/
    4444extern "C" int ekk_preSolveClp(EKKModel * model, bool keepIntegers,
    45                                int pass);
     45                               int pass);
    4646
    4747#include "ClpSimplex.hpp"
    … …  
    5151int main (int argc, const char *argv[])
    5252{
    53   const char * name;
    54   // problem is actually scaled for osl, dynamically for clp (slows clp)
    55   // default is primal, no presolve, minimise and use clp
    56   bool    primal=true,presolve=false;
    57   int useosl=0;
    58   bool freeFormat=false;
    59   bool exportIt=false;
     53     const char * name;
     54     // problem is actually scaled for osl, dynamically for clp (slows clp)
     55     // default is primal, no presolve, minimise and use clp
     56     bool    primal = true, presolve = false;
     57     int useosl = 0;
     58     bool freeFormat = false;
     59     bool exportIt = false;
    6060
    61   EKKModel * model;
    62   EKKContext * context;
     61     EKKModel * model;
     62     EKKContext * context;
    6363
    64   if ( argc>1 ) {
    65     name=argv[1];
    66   } else {
    67     name="../../Data/Sample/p0033.mps";
    68   }
     64     if ( argc > 1 ) {
     65          name = argv[1];
     66     } else {
     67          name = "../../Data/Sample/p0033.mps";
     68     }
    6969
    70   /* initialize OSL environment */
    71   context=ekk_initializeContext();
    72   model=ekk_newModel(context,"");
     70     /* initialize OSL environment */
     71     context = ekk_initializeContext();
     72     model = ekk_newModel(context, "");
    7373
    74   int i;
    75   printf("*** Options ");
    76   for (i=2;i<argc;i++) {
    77     printf("%s ",argv[i]);
    78   }
    79   printf("\n");
     74     int i;
     75     printf("*** Options ");
     76     for (i = 2; i < argc; i++) {
     77          printf("%s ", argv[i]);
     78     }
     79     printf("\n");
    8080
    81   // see if free format needed
     81     // see if free format needed
    8282
    83   for (i=2;i<argc;i++) {
    84     if (!strncmp(argv[i],"free",4)) {
    85       freeFormat=true;
    86     }
    87   }
     83     for (i = 2; i < argc; i++) {
     84          if (!strncmp(argv[i], "free", 4)) {
     85               freeFormat = true;
     86          }
     87     }
    8888
    89   // create model from MPS file
     89     // create model from MPS file
    9090
    91   if (!freeFormat) {
    92     ekk_importModel(model,name);
    93   } else {
    94     ekk_importModelFree(model,name);
    95   }
     91     if (!freeFormat) {
     92          ekk_importModel(model, name);
     93     } else {
     94          ekk_importModelFree(model, name);
     95     }
    9696
    97   // other options
    98   for (i=2;i<argc;i++) {
    99     if (!strncmp(argv[i],"max",3)) {
    100       if (!strncmp(argv[i],"max2",4)) {
    101         // This is for testing - it just reverses signs and maximizes
    102         int i,n=ekk_getInumcols(model);
    103         double * objective=ekk_getObjective(model);
    104         for (i=0;i<n;i++) {
    105           objective[i]=-objective[i];
    106         }
    107         ekk_setObjective(model,objective);
    108         ekk_setMaximize(model);
    109       } else {
    110         // maximize
    111         ekk_setMaximize(model);
    112       }
    113     }
    114     if (!strncmp(argv[i],"dual",4)) {
    115       primal=false;
    116     }
    117     if (!strncmp(argv[i],"presol",6)) {
    118       presolve=true;
    119     }
    120     if (!strncmp(argv[i],"osl",3)) {
    121       useosl=1;
    122     }
    123     if (!strncmp(argv[i],"both",4)) {
    124       useosl=2;
    125     }
    126     if (!strncmp(argv[i],"export",6)) {
    127       exportIt=true;
    128     }
    129   }
    130   if (useosl) {
    131     // OSL
    132     if (presolve)
    133       ekk_preSolve(model,3,NULL);
    134     ekk_scale(model);
    135    
    136     if (primal)
    137       ekk_primalSimplex(model,1);
    138     else
    139       ekk_dualSimplex(model);
    140     if (presolve) {
    141       ekk_postSolve(model,NULL);
    142       ekk_primalSimplex(model,3);
    143     }
    144     if (useosl==2)
    145       ekk_allSlackBasis(model); // otherwise it would be easy
    146   }
    147   if ((useosl&2)==0) {
    148     // CLP
    149     if (presolve)
    150       ekk_preSolveClp(model,true,5);
    151     /* 3 is because it is ignored if no presolve, and we
    152        are forcing Clp to re-optimize */
    153     if (primal)
    154       ekk_primalClp(model,1,3);
    155     else
    156       ekk_dualClp(model,3);
    157   }
    158   if (exportIt) {
    159     ClpSimplex * clp = new ClpSimplex();;
    160     int numberRows=ekk_getInumrows(model);
    161     int numberColumns= ekk_getInumcols(model);
    162     clp->loadProblem(numberColumns,numberRows,ekk_blockColumn(model,0),
    163                      ekk_blockRow(model,0),ekk_blockElement(model,0),
    164                      ekk_collower(model),ekk_colupper(model),
    165                      ekk_objective(model),
    166                      ekk_rowlower(model),ekk_rowupper(model));
    167     // Do integer stuff
    168     int * which =  ekk_listOfIntegers(model);
    169     int numberIntegers =  ekk_getInumints(model);
    170     for (int i=0;i<numberIntegers;i++)
    171       clp->setInteger(which[i]);
    172     ekk_free(which);
    173     clp->writeMps("try1.mps");
    174     delete clp;
    175   }
    176   ekk_deleteModel(model);
    177   ekk_endContext(context);
    178   return 0;
    179 }   
     97     // other options
     98     for (i = 2; i < argc; i++) {
     99          if (!strncmp(argv[i], "max", 3)) {
     100               if (!strncmp(argv[i], "max2", 4)) {
     101                    // This is for testing - it just reverses signs and maximizes
     102                    int i, n = ekk_getInumcols(model);
     103                    double * objective = ekk_getObjective(model);
     104                    for (i = 0; i < n; i++) {
     105                         objective[i] = -objective[i];
     106                    }
     107                    ekk_setObjective(model, objective);
     108                    ekk_setMaximize(model);
     109               } else {
     110                    // maximize
     111                    ekk_setMaximize(model);
     112               }
     113          }
     114          if (!strncmp(argv[i], "dual", 4)) {
     115               primal = false;
     116          }
     117          if (!strncmp(argv[i], "presol", 6)) {
     118               presolve = true;
     119          }
     120          if (!strncmp(argv[i], "osl", 3)) {
     121               useosl = 1;
     122          }
     123          if (!strncmp(argv[i], "both", 4)) {
     124               useosl = 2;
     125          }
     126          if (!strncmp(argv[i], "export", 6)) {
     127               exportIt = true;
     128          }
     129     }
     130     if (useosl) {
     131          // OSL
     132          if (presolve)
     133               ekk_preSolve(model, 3, NULL);
     134          ekk_scale(model);
     135
     136          if (primal)
     137               ekk_primalSimplex(model, 1);
     138          else
     139               ekk_dualSimplex(model);
     140          if (presolve) {
     141               ekk_postSolve(model, NULL);
     142               ekk_primalSimplex(model, 3);
     143          }
     144          if (useosl == 2)
     145               ekk_allSlackBasis(model); // otherwise it would be easy
     146     }
     147     if ((useosl & 2) == 0) {
     148          // CLP
     149          if (presolve)
     150               ekk_preSolveClp(model, true, 5);
     151          /* 3 is because it is ignored if no presolve, and we
     152             are forcing Clp to re-optimize */
     153          if (primal)
     154               ekk_primalClp(model, 1, 3);
     155          else
     156               ekk_dualClp(model, 3);
     157     }
     158     if (exportIt) {
     159          ClpSimplex * clp = new ClpSimplex();;
     160          int numberRows = ekk_getInumrows(model);
     161          int numberColumns = ekk_getInumcols(model);
     162          clp->loadProblem(numberColumns, numberRows, ekk_blockColumn(model, 0),
     163                           ekk_blockRow(model, 0), ekk_blockElement(model, 0),
     164                           ekk_collower(model), ekk_colupper(model),
     165                           ekk_objective(model),
     166                           ekk_rowlower(model), ekk_rowupper(model));
     167          // Do integer stuff
     168          int * which =  ekk_listOfIntegers(model);
     169          int numberIntegers =  ekk_getInumints(model);
     170          for (int i = 0; i < numberIntegers; i++)
     171               clp->setInteger(which[i]);
     172          ekk_free(which);
     173          clp->writeMps("try1.mps");
     174          delete clp;
     175     }
     176     ekk_deleteModel(model);
     177     ekk_endContext(context);
     178     return 0;
     179}

  • trunk/Clp/examples/ekk_interface.cpp

    r1370 r1552  
    77#include "ekk_c_api.h"
    88//#include "ekk_c_api_undoc.h"
    9  extern "C"{
    10  OSLLIBAPI void * OSLLINKAGE ekk_compressModel(EKKModel * model);
    11   OSLLIBAPI void OSLLINKAGE ekk_decompressModel(EKKModel * model,
    12                                                 void * compressInfo);
    13  }
    14 
    15 static ClpPresolve * presolveInfo=NULL;
    16 
    17 static ClpSimplex * clpmodel(EKKModel * model,int startup)
    18 {
    19   ClpSimplex * clp = new ClpSimplex();;
    20   int numberRows=ekk_getInumrows(model);
    21   int numberColumns= ekk_getInumcols(model);
    22   clp->loadProblem(numberColumns,numberRows,ekk_blockColumn(model,0),
    23                   ekk_blockRow(model,0),ekk_blockElement(model,0),
    24                   ekk_collower(model),ekk_colupper(model),
    25                   ekk_objective(model),
    26                   ekk_rowlower(model),ekk_rowupper(model));
    27   clp->setOptimizationDirection((int) ekk_getRmaxmin(model));
    28   clp->setPrimalTolerance(ekk_getRtolpinf(model));
    29   if (ekk_getRpweight(model)!=0.1)
    30     clp->setInfeasibilityCost(1.0/ekk_getRpweight(model));
    31   clp->setDualTolerance(ekk_getRtoldinf(model));
    32   if (ekk_getRdweight(model)!=0.1)
    33     clp->setDualBound(1.0/ekk_getRdweight(model));
    34   clp->setDblParam(ClpObjOffset,ekk_getRobjectiveOffset(model));
    35   const int * rowStatus =  ekk_rowstat(model);
    36   const double * rowSolution = ekk_rowacts(model);
    37   int i;
    38   clp->createStatus();
    39   double * clpSolution;
    40   clpSolution = clp->primalRowSolution();
    41   memcpy(clpSolution,rowSolution,numberRows*sizeof(double));
    42   const double * rowLower = ekk_rowlower(model);
    43   const double * rowUpper = ekk_rowupper(model);
    44   for (i=0;i<numberRows;i++) {
    45     ClpSimplex::Status status;
    46     if ((rowStatus[i]&0x80000000)!=0) {
    47       status = ClpSimplex::basic;
    48     } else {
    49       if (!startup) {
    50         // believe bits
    51         int ikey = rowStatus[i]&0x60000000;
    52         if (ikey==0x40000000) {
    53           // at ub
    54           status = ClpSimplex::atUpperBound;
    55           clpSolution[i]=rowUpper[i];
    56         } else if (ikey==0x20000000) {
    57           // at lb
    58           status = ClpSimplex::atLowerBound;
    59           clpSolution[i]=rowLower[i];
    60         } else if (ikey==0x60000000) {
    61           // free
    62           status = ClpSimplex::isFree;
    63           clpSolution[i]=0.0;
    64         } else {
    65           // fixed
    66           status = ClpSimplex::atLowerBound;
    67           clpSolution[i]=rowLower[i];
    68         }
    69       } else {
    70         status = ClpSimplex::superBasic;
    71       }
    72     }
    73     clp->setRowStatus(i,status);
    74   }
    75  
    76   const int * columnStatus = ekk_colstat(model);
    77   const double * columnSolution =  ekk_colsol(model);
    78   clpSolution = clp->primalColumnSolution();
    79   memcpy(clpSolution,columnSolution,numberColumns*sizeof(double));
    80   const double * columnLower = ekk_collower(model);
    81   const double * columnUpper = ekk_colupper(model);
    82   for (i=0;i<numberColumns;i++) {
    83     ClpSimplex::Status status;
    84     if ((columnStatus[i]&0x80000000)!=0) {
    85       status = ClpSimplex::basic;
    86     } else {
    87       if (!startup) {
    88         // believe bits
    89         int ikey = columnStatus[i]&0x60000000;
    90         if (ikey==0x40000000) {
    91           // at ub
    92           status = ClpSimplex::atUpperBound;
    93           clpSolution[i]=columnUpper[i];
    94         } else if (ikey==0x20000000) {
    95           // at lb
    96           status = ClpSimplex::atLowerBound;
    97           clpSolution[i]=columnLower[i];
    98         } else if (ikey==0x60000000) {
    99           // free
    100           status = ClpSimplex::isFree;
    101           clpSolution[i]=0.0;
    102         } else {
    103           // fixed
    104           status = ClpSimplex::atLowerBound;
    105           clpSolution[i]=columnLower[i];
    106         }
    107       } else {
    108         status = ClpSimplex::superBasic;
    109       }
    110     }
    111     clp->setColumnStatus(i,status);
    112   }
    113   return clp;
    114 }
    115 
    116 static int solve(EKKModel * model, int  startup,int algorithm,
    117                 int presolve)
    118 {
    119   // values pass or not
    120   if (startup)
    121     startup=1;
    122   // if scaled then be careful
    123   bool scaled = ekk_scaling(model)==1;
    124   if (scaled)
    125     ekk_scaleRim(model,1);
    126   void * compressInfo=NULL;
    127   ClpSimplex * clp;
    128   if (!presolve||!presolveInfo) {
    129     // no presolve or osl presolve - compact columns
    130     compressInfo=ekk_compressModel(model);
    131     clp = clpmodel(model,startup);;
    132   } else {
    133     // pick up clp model
    134     clp = presolveInfo->model();
    135   }
    136    
    137   // don't scale if alreday scaled
    138   if (scaled)
    139     clp->scaling(false);
    140   if (clp->numberRows()>10000)
    141     clp->factorization()->maximumPivots(100+clp->numberRows()/100);
    142   if (algorithm>0)
    143     clp->primal(startup);
    144   else
    145     clp->dual();
    146 
    147   int numberIterations = clp->numberIterations();
    148   if (presolve&&presolveInfo) {
    149     // very wasteful - create a clp copy of osl model
    150     ClpSimplex * clpOriginal = clpmodel(model,0);
    151     presolveInfo->setOriginalModel(clpOriginal);
    152     // do postsolve
    153     presolveInfo->postsolve(true);
    154     delete clp;
    155     delete presolveInfo;
    156     presolveInfo=NULL;
    157     clp = clpOriginal;
    158     if (presolve==3||(presolve==2&&clp->status())) {
    159       printf("Resolving from postsolved model\n");
    160       clp->primal(1);
    161       numberIterations += clp->numberIterations();
    162     }
    163   }
    164 
    165   // put back solution
    166 
    167   double * rowDual = (double *) ekk_rowduals(model);
    168   int numberRows=ekk_getInumrows(model);
    169   int numberColumns= ekk_getInumcols(model);
    170   int * rowStatus =  (int *) ekk_rowstat(model);
    171   double * rowSolution = (double *) ekk_rowacts(model);
    172   int i;
    173   int * columnStatus = (int *) ekk_colstat(model);
    174   double * columnSolution =  (double *) ekk_colsol(model);
    175   memcpy(rowSolution,clp->primalRowSolution(),numberRows*sizeof(double));
    176   memcpy(rowDual,clp->dualRowSolution(),numberRows*sizeof(double));
    177   for (i=0;i<numberRows;i++) {
    178     if (clp->getRowStatus(i)==ClpSimplex::basic)
    179       rowStatus[i]=0x80000000;
    180     else
    181       rowStatus[i]=0;
    182   }
    183  
    184   double * columnDual = (double *) ekk_colrcosts(model);
    185   memcpy(columnSolution,clp->primalColumnSolution(),
    186          numberColumns*sizeof(double));
    187   memcpy(columnDual,clp->dualColumnSolution(),numberColumns*sizeof(double));
    188   for (i=0;i<numberColumns;i++) {
    189     if (clp->getColumnStatus(i)==ClpSimplex::basic)
    190       columnStatus[i]=0x80000000;
    191     else
    192       columnStatus[i]=0;
    193   }
    194   ekk_setIprobstat(model,clp->status());
    195   ekk_setRobjvalue(model,clp->objectiveValue());
    196   ekk_setInumpinf(model,clp->numberPrimalInfeasibilities());
    197   ekk_setInumdinf(model,clp->numberDualInfeasibilities());
    198   ekk_setIiternum(model,numberIterations);
    199   ekk_setRsumpinf(model,clp->sumPrimalInfeasibilities());
    200   ekk_setRsumdinf(model,clp->sumDualInfeasibilities());
    201   delete clp;
    202   if (compressInfo )
    203     ekk_decompressModel(model,compressInfo);
    204 
    205   if (scaled)
    206       ekk_scaleRim(model,2);
    207   return 0;
     9extern "C" {
     10     OSLLIBAPI void * OSLLINKAGE ekk_compressModel(EKKModel * model);
     11     OSLLIBAPI void OSLLINKAGE ekk_decompressModel(EKKModel * model,
     12               void * compressInfo);
     13}
     14
     15static ClpPresolve * presolveInfo = NULL;
     16
     17static ClpSimplex * clpmodel(EKKModel * model, int startup)
     18{
     19     ClpSimplex * clp = new ClpSimplex();;
     20     int numberRows = ekk_getInumrows(model);
     21     int numberColumns = ekk_getInumcols(model);
     22     clp->loadProblem(numberColumns, numberRows, ekk_blockColumn(model, 0),
     23                      ekk_blockRow(model, 0), ekk_blockElement(model, 0),
     24                      ekk_collower(model), ekk_colupper(model),
     25                      ekk_objective(model),
     26                      ekk_rowlower(model), ekk_rowupper(model));
     27     clp->setOptimizationDirection((int) ekk_getRmaxmin(model));
     28     clp->setPrimalTolerance(ekk_getRtolpinf(model));
     29     if (ekk_getRpweight(model) != 0.1)
     30          clp->setInfeasibilityCost(1.0 / ekk_getRpweight(model));
     31     clp->setDualTolerance(ekk_getRtoldinf(model));
     32     if (ekk_getRdweight(model) != 0.1)
     33          clp->setDualBound(1.0 / ekk_getRdweight(model));
     34     clp->setDblParam(ClpObjOffset, ekk_getRobjectiveOffset(model));
     35     const int * rowStatus =  ekk_rowstat(model);
     36     const double * rowSolution = ekk_rowacts(model);
     37     int i;
     38     clp->createStatus();
     39     double * clpSolution;
     40     clpSolution = clp->primalRowSolution();
     41     memcpy(clpSolution, rowSolution, numberRows * sizeof(double));
     42     const double * rowLower = ekk_rowlower(model);
     43     const double * rowUpper = ekk_rowupper(model);
     44     for (i = 0; i < numberRows; i++) {
     45          ClpSimplex::Status status;
     46          if ((rowStatus[i] & 0x80000000) != 0) {
     47               status = ClpSimplex::basic;
     48          } else {
     49               if (!startup) {
     50                    // believe bits
     51                    int ikey = rowStatus[i] & 0x60000000;
     52                    if (ikey == 0x40000000) {
     53                         // at ub
     54                         status = ClpSimplex::atUpperBound;
     55                         clpSolution[i] = rowUpper[i];
     56                    } else if (ikey == 0x20000000) {
     57                         // at lb
     58                         status = ClpSimplex::atLowerBound;
     59                         clpSolution[i] = rowLower[i];
     60                    } else if (ikey == 0x60000000) {
     61                         // free
     62                         status = ClpSimplex::isFree;
     63                         clpSolution[i] = 0.0;
     64                    } else {
     65                         // fixed
     66                         status = ClpSimplex::atLowerBound;
     67                         clpSolution[i] = rowLower[i];
     68                    }
     69               } else {
     70                    status = ClpSimplex::superBasic;
     71               }
     72          }
     73          clp->setRowStatus(i, status);
     74     }
     75
     76     const int * columnStatus = ekk_colstat(model);
     77     const double * columnSolution =  ekk_colsol(model);
     78     clpSolution = clp->primalColumnSolution();
     79     memcpy(clpSolution, columnSolution, numberColumns * sizeof(double));
     80     const double * columnLower = ekk_collower(model);
     81     const double * columnUpper = ekk_colupper(model);
     82     for (i = 0; i < numberColumns; i++) {
     83          ClpSimplex::Status status;
     84          if ((columnStatus[i] & 0x80000000) != 0) {
     85               status = ClpSimplex::basic;
     86          } else {
     87               if (!startup) {
     88                    // believe bits
     89                    int ikey = columnStatus[i] & 0x60000000;
     90                    if (ikey == 0x40000000) {
     91                         // at ub
     92                         status = ClpSimplex::atUpperBound;
     93                         clpSolution[i] = columnUpper[i];
     94                    } else if (ikey == 0x20000000) {
     95                         // at lb
     96                         status = ClpSimplex::atLowerBound;
     97                         clpSolution[i] = columnLower[i];
     98                    } else if (ikey == 0x60000000) {
     99                         // free
     100                         status = ClpSimplex::isFree;
     101                         clpSolution[i] = 0.0;
     102                    } else {
     103                         // fixed
     104                         status = ClpSimplex::atLowerBound;
     105                         clpSolution[i] = columnLower[i];
     106                    }
     107               } else {
     108                    status = ClpSimplex::superBasic;
     109               }
     110          }
     111          clp->setColumnStatus(i, status);
     112     }
     113     return clp;
     114}
     115
     116static int solve(EKKModel * model, int  startup, int algorithm,
     117                int presolve)
     118{
     119     // values pass or not
     120     if (startup)
     121          startup = 1;
     122     // if scaled then be careful
     123     bool scaled = ekk_scaling(model) == 1;
     124     if (scaled)
     125          ekk_scaleRim(model, 1);
     126     void * compressInfo = NULL;
     127     ClpSimplex * clp;
     128     if (!presolve || !presolveInfo) {
     129          // no presolve or osl presolve - compact columns
     130          compressInfo = ekk_compressModel(model);
     131          clp = clpmodel(model, startup);;
     132     } else {
     133          // pick up clp model
     134          clp = presolveInfo->model();
     135     }
     136
     137     // don't scale if alreday scaled
     138     if (scaled)
     139          clp->scaling(false);
     140     if (clp->numberRows() > 10000)
     141          clp->factorization()->maximumPivots(100 + clp->numberRows() / 100);
     142     if (algorithm > 0)
     143          clp->primal(startup);
     144     else
     145          clp->dual();
     146
     147     int numberIterations = clp->numberIterations();
     148     if (presolve && presolveInfo) {
     149          // very wasteful - create a clp copy of osl model
     150          ClpSimplex * clpOriginal = clpmodel(model, 0);
     151          presolveInfo->setOriginalModel(clpOriginal);
     152          // do postsolve
     153          presolveInfo->postsolve(true);
     154          delete clp;
     155          delete presolveInfo;
     156          presolveInfo = NULL;
     157          clp = clpOriginal;
     158          if (presolve == 3 || (presolve == 2 && clp->status())) {
     159               printf("Resolving from postsolved model\n");
     160               clp->primal(1);
     161               numberIterations += clp->numberIterations();
     162          }
     163     }
     164
     165     // put back solution
     166
     167     double * rowDual = (double *) ekk_rowduals(model);
     168     int numberRows = ekk_getInumrows(model);
     169     int numberColumns = ekk_getInumcols(model);
     170     int * rowStatus =  (int *) ekk_rowstat(model);
     171     double * rowSolution = (double *) ekk_rowacts(model);
     172     int i;
     173     int * columnStatus = (int *) ekk_colstat(model);
     174     double * columnSolution =  (double *) ekk_colsol(model);
     175     memcpy(rowSolution, clp->primalRowSolution(), numberRows * sizeof(double));
     176     memcpy(rowDual, clp->dualRowSolution(), numberRows * sizeof(double));
     177     for (i = 0; i < numberRows; i++) {
     178          if (clp->getRowStatus(i) == ClpSimplex::basic)
     179               rowStatus[i] = 0x80000000;
     180          else
     181               rowStatus[i] = 0;
     182     }
     183
     184     double * columnDual = (double *) ekk_colrcosts(model);
     185     memcpy(columnSolution, clp->primalColumnSolution(),
     186            numberColumns * sizeof(double));
     187     memcpy(columnDual, clp->dualColumnSolution(), numberColumns * sizeof(double));
     188     for (i = 0; i < numberColumns; i++) {
     189          if (clp->getColumnStatus(i) == ClpSimplex::basic)
     190               columnStatus[i] = 0x80000000;
     191          else
     192               columnStatus[i] = 0;
     193     }
     194     ekk_setIprobstat(model, clp->status());
     195     ekk_setRobjvalue(model, clp->objectiveValue());
     196     ekk_setInumpinf(model, clp->numberPrimalInfeasibilities());
     197     ekk_setInumdinf(model, clp->numberDualInfeasibilities());
     198     ekk_setIiternum(model, numberIterations);
     199     ekk_setRsumpinf(model, clp->sumPrimalInfeasibilities());
     200     ekk_setRsumdinf(model, clp->sumDualInfeasibilities());
     201     delete clp;
     202     if (compressInfo )
     203          ekk_decompressModel(model, compressInfo);
     204
     205     if (scaled)
     206          ekk_scaleRim(model, 2);
     207     return 0;
    208208}
    209209/* As ekk_primalSimplex + postsolve instructions:
    … …