- Timestamp:
- May 24, 2010 9:03:59 PM (11 years ago)
- Location:
- trunk/Clp/examples
- Files:
-
- 31 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Clp/examples/addBits.cpp
r1370 r1552 34 34 int main (int argc, const char *argv[]) 35 35 { 36 // Empty model37 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); 41 41 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(); 181 67 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 15 15 int main (int argc, const char *argv[]) 16 16 { 17 {18 // Empty model19 ClpSimplex model;20 21 // Bounds on rows - as dense vector22 double lower[]={2.0,1.0};23 double upper[]={COIN_DBL_MAX,1.0};24 25 // Create space for 2 rows26 model.resize(2,0);27 // Fill in28 int i;29 // Now row bounds30 for (i=0;i<2;i++) {31 model.setRowLower(i,lower[i]);32 model.setRowUpper(i,upper[i]);33 }34 // Now add column 135 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 240 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 345 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 // solve50 model.dual();51 52 /*53 Adding one column at a time has a significant overhead so let's54 try a more complicated but faster way55 56 First time adding in 10000 columns one by one57 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 build73 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 +-186 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 +-1102 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 model118 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 columns130 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 else155 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 else160 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 else165 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 else170 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 else175 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 model183 ClpSimplex model;184 int status;185 if (argc<2)186 status=model.readMps("../../Data/Sample/p0033.mps");187 else188 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 model2200 ClpSimplex model2;201 model2.addRows(numberRows,rowLower,rowUpper,NULL);202 203 // Build object204 CoinBuild buildObject;205 // Add columns206 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 in222 model2.addColumns(buildObject);223 model2.initialSolve();224 }225 {226 // and again227 ClpSimplex model;228 int status;229 if (argc<2)230 status=model.readMps("../../Data/Sample/p0033.mps");231 else232 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 object244 CoinModel buildObject;245 for (int iRow=0;iRow<numberRows;iRow++)246 buildObject.setRowBounds(iRow,rowLower[iRow],rowUpper[iRow]);247 248 // Add columns249 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 in265 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 15 15 int main (int argc, const char *argv[]) 16 16 { 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 9 9 int main (int argc, const char *argv[]) 10 10 { 11 /* Create a structured model by reading mps file and trying12 Dantzig-Wolfe decomposition (that's the 1 parameter)13 */14 // At present D-W rows are hard coded - will move stuff from OSL15 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 stuff20 ClpSimplex solver;21 /*22 This driver does a simple Dantzig Wolfe decomposition23 */24 solver.solve(&model);25 // Double check26 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 11 11 int main (int argc, const char *argv[]) 12 12 { 13 /* Create a structured model by reading mps file and trying14 Dantzig-Wolfe or Benders decomposition15 */16 // Get maximum number of blocks17 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 PRESOLVE13 /* 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 32 32 #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 stuff38 ClpSimplex solver;39 // change factorization frequency from 20040 solver.setFactorizationFrequency(100+model.numberRows()/50);41 /*42 This driver does a simple Dantzig Wolfe decomposition43 */44 double time1 = CoinCpuTime() ;45 solver.solve(&model);46 std::cout<<"model took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;47 // Double check48 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); 49 49 #else 50 ClpSimplex model;51 int status = model.readMps((argc<2) ? "../../Data/Netlib/czprob.mps"52 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; 58 58 59 CoinStructuredModel modelB;60 modelB.decompose(*model2->matrix(),61 model2->rowLower(),model2->rowUpper(),62 model2->columnLower(),model2->columnUpper(),63 model2->objective(),1,maxBlocks,64 65 // change factorization frequency from 20066 model2->setFactorizationFrequency(100+modelB.numberRows()/50);67 /*68 This driver does a simple Dantzig Wolfe decomposition69 */70 double time1 = CoinCpuTime() ;71 model2->solve(&modelB);72 std::cout<<"model took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;73 // But resolve for safety74 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); 75 75 #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 9 9 int main (int argc, const char *argv[]) 10 10 { 11 ClpSimplex model;12 int status;13 // Keep names14 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 decomposition24 */25 // Get master rows in some mysterious way26 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 rows32 if (numberRows==105127) {33 // ken-1834 for (iRow=104976;iRow<numberRows;iRow++)35 rowBlock[iRow]=-1;36 } else if (numberRows==2426) {37 // ken-738 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-6048 for (iRow=10198;iRow<10280;iRow++)49 rowBlock[iRow]=-1;50 } else if (numberRows==1503) {51 // degen352 for (iRow=0;iRow<561;iRow++)53 rowBlock[iRow]=-1;54 } else if (numberRows==929) {55 // czprob56 for (iRow=0;iRow<39;iRow++)57 rowBlock[iRow]=-1;58 }59 CoinPackedMatrix * matrix = model.matrix();60 // get row copy61 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 at74 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 allocated84 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; // mark91 92 93 }94 if (nstack) {95 96 97 columnBlock[iColumn]=numberBlocks-1;98 99 100 101 for (k=rowStart[iRow];k<rowStart[iRow]+rowLength[iRow];k++) {102 103 104 105 if (columnBlock[iColumn]==-2) {106 columnBlock[iColumn]=numberBlocks-1; // mark107 108 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 117 assert (columnBlock[iColumn]==numberBlocks-1);118 119 120 121 } else {122 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 problems144 CoinPackedMatrix * top = new CoinPackedMatrix [numberBlocks];145 ClpSimplex * sub = new ClpSimplex [numberBlocks];146 ClpSimplex master;147 // Create all sub problems148 // Could do much faster - but do that later149 int * whichRow = new int [numberRows];150 int * whichColumn = new int [numberColumns];151 // get top matrix152 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); 171 171 #if 0 172 // temp173 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; 176 176 #endif 177 // and top matrix178 CoinPackedMatrix matrix = topMatrix;179 // and delete bits180 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 master188 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; 201 201 202 // Overkill in terms of space203 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 proposals212 int numberMasterColumns = master.numberColumns();213 master.resize(numberMasterRows+numberBlocks,numberMasterColumns);214 // Arrays to say which block and when created215 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 artificials220 {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 238 }239 // and resize matrix to double check clp will be happy240 //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); 242 242 243 for (iPass=0;iPass<maxPass;iPass++) {244 printf("Start of pass %d\n",iPass);245 // Solve master - may be infeasible246 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 } 250 250 251 master.primal();252 int problemStatus = master.status(); // do here as can change (delcols)253 if (master.numberIterations()==0&&iPass)254 break; // finished255 if (master.objectiveValue()>lastObjective-1.0e-7&&iPass>555)256 break; // finished257 lastObjective = master.objectiveValue();258 // mark basic ones and delete if necessary259 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 // delete267 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 273 when[numberKeep]=when[iColumn];274 whichBlock[numberKeep++]=whichBlock[iColumn];275 276 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 objective290 dual = master.infeasibilityRay();291 deleteDual = true;292 printf("The sum of infeasibilities is %g\n",293 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 solve303 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 objective311 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 319 }320 sub[iBlock].primal();321 memcpy(objective2,saveObj,numberColumns2*sizeof(double));322 // get proposal323 if (sub[iBlock].numberIterations()||!iPass) {324 double objValue =0.0;325 326 327 328 329 top[iBlock].times(solution,elementAdd+start);330 for (i=0;i<numberColumns2;i++)331 objValue += solution[i]*saveObj[i];332 333 int number=start;334 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 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 350 351 rowAdd[number]=numberMasterRows+iBlock;352 elementAdd[number++]=1.0;353 354 355 356 iBlock,smallest,largest,dj);357 if (dj<-1.0e-6||!iPass) {358 359 objective[numberProposals]=objValue;360 columnAdd[++numberProposals]=number;361 when[numberColumnsGenerated]=iPass;362 whichBlock[numberColumnsGenerated++]=iBlock;363 364 365 366 367 top[iBlock].times(solution,elementAdd+start);368 for (i=0;i<numberColumns2;i++)369 objValue += solution[i]*saveObj[i];370 371 int number=start;372 373 double smallest=1.0e100;374 double largest=0.0;375 for (i=0;i<numberMasterRows;i++) {376 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 386 387 388 iBlock,smallest,largest,dj);389 if (dj<-1.0e-6) {390 391 objective[numberProposals]=objValue;392 columnAdd[++numberProposals]=number;393 when[numberColumnsGenerated]=iPass;394 whichBlock[numberColumnsGenerated++]=iBlock;395 396 397 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 solution409 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 rows437 top[iBlock].reverseOrdering();438 // zero solution439 memset(sol,0,numberColumnsGenerated*sizeof(double));440 int i;441 for (i=numberMasterColumns;i<numberColumnsGenerated;i++)442 if (whichBlock[i-numberMasterColumns]==iBlock)443 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 450 else451 upper[iRow]=COIN_DBL_MAX;452 if (masterLower[iRow]>-1.0e20)453 454 else455 lower[iRow]=-COIN_DBL_MAX;456 }457 sub[iBlock].addRows(numberMasterRows,lower,upper,458 459 460 461 462 sub[iBlock].primal();463 const double * subSolution = sub[iBlock].primalColumnSolution();464 const double * subRowSolution = sub[iBlock].primalRowSolution();465 // move solution466 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 basic489 model.setStatus(iColumn,ClpSimplex::atUpperBound);490 } else if (fullSolution[iColumn]<=fullLower[iColumn]+1.0e-8) {491 // may help to make rest non basic492 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 10 10 int main (int argc, const char *argv[]) 11 11 { 12 ClpSimplex model;13 int status;14 // Keep names15 if (argc<2)16 status=model.readMps("../../Data/Sample/p0033.mps",true);17 else18 status=model.readMps(argv[1],true);19 /*20 This driver is similar to minimum.cpp, but it21 sets all parameter values to their defaults. The purpose of this22 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 like28 while the other adheres to OsiSolverInterface standards29 */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 cost40 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, getNumCols48 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. ObjectiveLimits54 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 this69 model.getDblParam(ClpPrimalTolerance,value);70 model.setPrimalTolerance(value);71 72 model.setDualTolerance( model.dualTolerance()) ;73 74 // Other Param stuff75 76 // Can also use maximumIterations77 model.getIntParam(ClpMaxNumIteration,integerValue);78 std::cout<<"Value of ClpMaxNumIteration is "<<integerValue<<std::endl;79 model.setMaximumIterations(integerValue);80 81 // Not sure this works yet82 model.getIntParam(ClpMaxNumIterationHotStart,integerValue);83 std::cout<<"Value of ClpMaxNumIterationHotStart is "84 <<integerValue<<std::endl;85 86 // Can also use getIterationCount and getObjValue87 /* Status of problem:88 0 - optimal89 1 - primal infeasible90 2 - dual infeasible91 3 - stopped on iterations etc92 4 - stopped due to errors93 */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 optimal108 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 slack116 unsigned char * basis1 = model.statusCopy();117 model.createStatus();118 119 // Now create another model and do hot start120 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 getObjSense128 model2.setOptimizationDirection(model.optimizationDirection());129 130 // Can use scalingFlag() to check if scaling on131 // But set up scaling132 model2.scaling();133 134 // Could play with sparse factorization on/off135 model2.setSparseFactorization(model.sparseFactorization());136 137 // Sets row pivot choice algorithm in dual138 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 handling147 // To see real use - see OsiOslSolverInterfaceTest.cpp148 CoinMessageHandler handler;149 model2.passInMessageHandler(& handler);150 model2.newLanguage(CoinMessages::us_en);151 152 //Increase level of detail153 model2.setLogLevel(4);154 155 // solve156 model2.dual();157 // flip direction twice and solve158 model2.setOptimizationDirection(-1);159 model2.dual();160 model2.setOptimizationDirection(1);161 //Decrease level of detail162 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 with171 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 modifiable181 arrays while the alternative names return const pointers -182 which is of course much more virtuous183 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 0200 assert(model2.lengthNames());201 202 // Row names203 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 else218 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 else223 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 else228 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 else233 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 else239 std::cout<<setiosflags(std::ios::scientific)<<std::setw(14)<<value;240 }241 std::cout<<std::endl;242 }243 std::cout<<"--------------------------------------"<<std::endl;244 245 // Columns246 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 0261 assert(model2.lengthNames());262 263 // Column names264 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 else279 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 else284 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 else289 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 else294 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 else299 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 matrix308 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 8 8 int main (int argc, const char *argv[]) 9 9 { 10 ClpSimplex model;11 int status;12 // Keep names when reading an mps file13 if (argc<2)14 status=model.readMps("../../Data/Sample/p0033.mps",true);15 else16 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); 17 17 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 } 23 23 #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 } 30 30 #else 31 ClpSolve solvectl;31 ClpSolve solvectl; 32 32 33 33 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 } 52 52 #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; 57 57 58 // remove this to print solution58 // remove this to print solution 59 59 60 exit(0);60 exit(0); 61 61 62 /*63 Now to print out solution. The methods used return modifiable64 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. 66 66 67 This version just does non-zero columns68 69 */67 This version just does non-zero columns 68 69 */ 70 70 #if 0 71 int numberRows = model.numberRows();71 int numberRows = model.numberRows(); 72 72 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(); 86 83 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(); 89 89 90 90 91 int iRow;91 int iRow; 92 92 93 std::cout<<" Primal Dual Lower Upper (Cost)"94 <<std::endl;93 std::cout << " Primal Dual Lower Upper (Cost)" 94 << std::endl; 95 95 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 else103 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 else108 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 else113 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 else118 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 else124 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 } 128 128 #endif 129 std::cout<<"--------------------------------------"<<std::endl;129 std::cout << "--------------------------------------" << std::endl; 130 130 131 // Columns131 // Columns 132 132 133 int numberColumns = model.numberColumns();133 int numberColumns = model.numberColumns(); 134 134 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(); 148 145 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(); 151 151 152 152 153 int iColumn; 154 155 std::cout<<" Primal Dual Lower Upper Cost" 156 <<std::endl; 153 int iColumn; 157 154 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; 192 157 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 15 15 16 16 /** 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. 18 18 This is used in Clp/Test/unitTest.cpp 19 19 … … 24 24 25 25 class MyMessageHandler : public CoinMessageHandler { 26 26 27 27 public: 28 /**@name Overrides */29 //@{30 virtual int print();31 //@}32 /**@name set and get */33 //@{34 /// Model35 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 model44 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 /// Clone59 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 63 63 protected: 64 /**@name Data members65 The data members are protected to allow access for derived classes. */66 //@{67 /// Pointer back to model68 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 //@} 70 70 }; 71 71 … … 76 76 77 77 //------------------------------------------------------------------- 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 //------------------------------------------------------------------- 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 { 99 99 } 100 100 101 101 // Constructor with pointer to model 102 102 MyMessageHandler::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 111 111 //------------------------------------------------------------------- 112 112 MyMessageHandler::~MyMessageHandler () … … 115 115 116 116 //---------------------------------------------------------------- 117 // Assignment operator 117 // Assignment operator 118 118 //------------------------------------------------------------------- 119 119 MyMessageHandler & 120 120 MyMessageHandler::operator=(const MyMessageHandler& rhs) 121 121 { 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; 127 127 } 128 128 //------------------------------------------------------------------- … … 131 131 CoinMessageHandler * MyMessageHandler::clone() const 132 132 { 133 return new MyMessageHandler(*this);133 return new MyMessageHandler(*this); 134 134 } 135 135 // Print out values from first 20 messages 136 static int times =0;137 int 136 static int times = 0; 137 int 138 138 MyMessageHandler::print() 139 139 { 140 // You could have added a callback flag if you had wanted - see Clp_C_Interface.c141 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(); 164 164 } 165 165 const ClpSimplex * 166 166 MyMessageHandler::model() const 167 167 { 168 return model_;169 } 170 void 168 return model_; 169 } 170 void 171 171 MyMessageHandler::setModel(ClpSimplex * model) 172 172 { 173 model_ = model;173 model_ = model; 174 174 } 175 175 176 176 int main (int argc, const char *argv[]) 177 177 { 178 ClpSimplex model;179 // Message handler180 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 file185 if (argc<2)186 status=model.readMps("../../Data/Sample/p0033.mps",true);187 else188 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 this203 /* Use a tolerance of 1.0e-8 for feasibility, treat problem as204 not being integer, do "numberpasses" passes and throw away names205 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 tolerances213 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 200223 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 variables227 to other bound to stay dual feasible. We can trash the bounds as 228 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 254 255 256 257 258 // But resolve for safety259 model.primal(1);260 261 numberIterations += model.numberIterations();;262 // for running timing tests263 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 solution273 274 exit(0);275 276 /*277 Now to print out solution. The methods used return modifiable278 arrays while the alternative names return const pointers -279 which is of course much more virtuous.280 281 This version just does non-zero columns282 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 */ 284 284 #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 0299 assert(model.lengthNames());300 301 // Row names302 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 else317 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 else322 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 else327 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 else332 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 else338 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 } 342 342 #endif 343 std::cout<<"--------------------------------------"<<std::endl;344 345 // Columns346 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 0361 assert(model.lengthNames());362 363 // Column names364 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 else380 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 else385 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 else390 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 else395 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 else400 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 13 13 /* Call back function - just says whenever it gets Clp0005 or Coin0002 */ 14 14 static void callBack(Clp_Simplex * model, int messageNumber, 15 16 17 int nString, char ** vString) 15 int nDouble, const double * vDouble, 16 int nInt, const int * vInt, 17 int nString, char ** vString) 18 18 { 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 } 31 31 } 32 32 33 33 34 34 35 35 int main (int argc, const char *argv[]) 36 36 { 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 else46 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); 47 47 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 } 53 53 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 } 60 60 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 } 67 67 68 /* remove this to print solution */68 /* remove this to print solution */ 69 69 70 /*exit(0); */70 /*exit(0); */ 71 71 72 {73 /*74 Now to print out solution. The methods used return modifiable75 arrays while the alternative names return const pointers -76 which is of course much more virtuous.77 78 This version just does non-zero columns79 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 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 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 1 1 /* $Id$ */ 2 /* Copyright (C) 2003, International Business Machines Corporation 2 /* Copyright (C) 2003, International Business Machines Corporation 3 3 and others. All Rights Reserved. 4 4 5 This sample program is designed to illustrate programming 5 This sample program is designed to illustrate programming 6 6 techniques using CoinLP, has not been thoroughly tested 7 7 and comes without any warranty whatsoever. 8 8 9 You may copy, modify and distribute this sample program without 9 You may copy, modify and distribute this sample program without 10 10 any restrictions whatsoever and without any payment to anyone. 11 11 */ … … 22 22 int main (int argc, const char *argv[]) 23 23 { 24 ClpSimplex model;25 int status;26 // Keep names27 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 this43 /* Use a tolerance of 1.0e-8 for feasibility, treat problem as44 not being integer, do "numberpasses" passes and throw away names45 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 tolerances53 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 20063 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 add69 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 list78 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 problem88 int smallNumberRows = 2*numberColumns;89 smallNumberRows=CoinMin(smallNumberRows,originalNumberRows/20);90 // and pad out with random rows91 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 */ 102 102 #if 0 103 // However for some we need to do anyway104 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 } 110 110 #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 passes116 int maxPass=50;117 // And take out slack rows until this pass118 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 way131 std::sort(sort,sort+numberSort);132 // Create small problem133 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 iterations138 //if (iPass)139 //small.setMaximumIterations(100);140 // Solve141 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); // idiot150 small.initialSolve(solveOptions);151 }152 bool dualInfeasible = (small.status()==2);153 // move solution back154 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 solution166 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 at170 for (iRow=0;iRow<originalNumberRows;iRow++)171 weight[iRow]=1.123e50;172 // Look at rows already in small problem173 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 182 if (model2->getRowStatus(iRow)==ClpSimplex::basic) {183 184 if (iPass<takeOutPass&&!dualInfeasible) {185 186 double infeasibility = CoinMax(fullSolution[iRow]-rowUpper[iRow],187 rowLower[iRow]-fullSolution[iRow]);188 weight[iRow]=-infeasibility;189 if (infeasibility>1.0e-8) {190 191 192 193 weight[iRow]=1.0;194 195 196 197 198 weight[iRow]=-1.0e40;199 200 201 202 203 weight[iRow]=-1.0e50;204 205 206 207 }208 // Now rest209 for (iRow=0;iRow<originalNumberRows;iRow++) {210 sort[iRow]=iRow;211 if (weight[iRow]==1.123e50) {212 213 double infeasibility = CoinMax(fullSolution[iRow]-rowUpper[iRow],214 rowLower[iRow]-fullSolution[iRow]);215 weight[iRow]=-infeasibility;216 if (infeasibility>1.0e-8) {217 218 219 220 221 }222 // sort223 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 233 234 235 236 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 243 if (!numberInfeasibilities) {244 245 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 252 253 254 }255 printf("in small model %d rows are infeasible - sum is %g\n",numberInfeasibilities,256 257 }258 }259 delete [] weight;260 delete [] sort;261 delete [] whichColumns;262 delete [] take;263 // If problem is big you may wish to skip this264 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 282 283 284 285 286 // But resolve for safety287 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 43 43 */ 44 44 extern "C" int ekk_preSolveClp(EKKModel * model, bool keepIntegers, 45 45 int pass); 46 46 47 47 #include "ClpSimplex.hpp" … … 51 51 int main (int argc, const char *argv[]) 52 52 { 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 clp56 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; 60 60 61 EKKModel * model;62 EKKContext * context;61 EKKModel * model; 62 EKKContext * context; 63 63 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 } 69 69 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, ""); 73 73 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"); 80 80 81 // see if free format needed81 // see if free format needed 82 82 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 } 88 88 89 // create model from MPS file89 // create model from MPS file 90 90 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 } 96 96 97 // other options98 for (i=2;i<argc;i++) {99 if (!strncmp(argv[i],"max",3)) {100 if (!strncmp(argv[i],"max2",4)) {101 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 109 } else {110 111 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 // OSL132 if (presolve)133 ekk_preSolve(model,3,NULL);134 ekk_scale(model);135 136 if (primal)137 ekk_primalSimplex(model,1);138 else139 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 easy146 }147 if ((useosl&2)==0) {148 // CLP149 if (presolve)150 ekk_preSolveClp(model,true,5);151 /* 3 is because it is ignored if no presolve, and we152 are forcing Clp to re-optimize */153 if (primal)154 ekk_primalClp(model,1,3);155 else156 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 stuff168 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 7 7 #include "ekk_c_api.h" 8 8 //#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 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 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 51 int ikey = rowStatus[i]&0x60000000;52 if (ikey==0x40000000) {53 54 55 clpSolution[i]=rowUpper[i];56 } else if (ikey==0x20000000) {57 58 59 clpSolution[i]=rowLower[i];60 } else if (ikey==0x60000000) {61 62 63 clpSolution[i]=0.0;64 65 66 67 clpSolution[i]=rowLower[i];68 69 } else {70 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 89 int ikey = columnStatus[i]&0x60000000;90 if (ikey==0x40000000) {91 92 93 clpSolution[i]=columnUpper[i];94 } else if (ikey==0x20000000) {95 96 97 clpSolution[i]=columnLower[i];98 } else if (ikey==0x60000000) {99 100 101 clpSolution[i]=0.0;102 103 104 105 clpSolution[i]=columnLower[i];106 107 } else {108 109 }110 }111 clp->setColumnStatus(i,status);112 }113 return clp;114 } 115 116 static int solve(EKKModel * model, int startup, int algorithm,117 118 { 119 // values pass or not120 if (startup)121 startup=1;122 // if scaled then be careful123 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 columns130 compressInfo=ekk_compressModel(model);131 clp = clpmodel(model,startup);;132 } else {133 // pick up clp model134 clp = presolveInfo->model();135 }136 137 // don't scale if alreday scaled138 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 else145 clp->dual();146 147 int numberIterations = clp->numberIterations();148 if (presolve&&presolveInfo) {149 // very wasteful - create a clp copy of osl model150 ClpSimplex * clpOriginal = clpmodel(model,0);151 presolveInfo->setOriginalModel(clpOriginal);152 // do postsolve153 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 solution166 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 else181 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 else192 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;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; 208 208 } 209 209 /* As ekk_primalSimplex + postsolve instructions: … …